text:
# Autogenerated wrapper script for NOMAD_jll for i686-linux-musl-cxx11
export libnomadAlgos, libnomadEval, libnomadInterface, libnomadUtils, libsgtelib, nomad

using CompilerSupportLibraries_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"
LIBPATH_default = ""

# Relative path to `libnomadAlgos`
const libnomadAlgos_splitpath = ["lib", "libnomadAlgos.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libnomadAlgos_path = ""

# libnomadAlgos-specific global declaration
# This will be filled out by __init__()
libnomadAlgos_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libnomadAlgos = "libnomadAlgos.so"

# Relative path to `libnomadEval`
const libnomadEval_splitpath = ["lib", "libnomadEval.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libnomadEval_path = ""

# libnomadEval-specific global declaration
# This will be filled out by __init__()
libnomadEval_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libnomadEval = "libnomadEval.so"

# Relative path to `libnomadInterface`
const libnomadInterface_splitpath = ["lib", "libnomadInterface.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libnomadInterface_path = ""

# libnomadInterface-specific global declaration
# This will be filled out by __init__()
libnomadInterface_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libnomadInterface = "libnomadInterface.so"

# Relative path to `libnomadUtils`
const libnomadUtils_splitpath = ["lib", "libnomadUtils.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libnomadUtils_path = ""

# libnomadUtils-specific global declaration
# This will be filled out by __init__()
libnomadUtils_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libnomadUtils = "libnomadUtils.so"

# Relative path to `libsgtelib`
const libsgtelib_splitpath = ["lib", "libsgtelib.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libsgtelib_path = ""

# libsgtelib-specific global declaration
# This will be filled out by __init__()
libsgtelib_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libsgtelib = "libsgtelib.so"

# Relative path to `nomad`
const nomad_splitpath = ["bin", "nomad"]

# This will be filled out by __init__() for all products, as it must be done at runtime
nomad_path = ""

# nomad-specific global declaration
function nomad(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(nomad_path)
    end
end

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"NOMAD")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (CompilerSupportLibraries_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (CompilerSupportLibraries_jll.LIBPATH_list,))

    global libnomadAlgos_path = normpath(joinpath(artifact_dir, libnomadAlgos_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libnomadAlgos_handle = dlopen(libnomadAlgos_path)
    push!(LIBPATH_list, dirname(libnomadAlgos_path))

    global libnomadEval_path = normpath(joinpath(artifact_dir, libnomadEval_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libnomadEval_handle = dlopen(libnomadEval_path)
    push!(LIBPATH_list, dirname(libnomadEval_path))

    global libnomadInterface_path = normpath(joinpath(artifact_dir, libnomadInterface_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libnomadInterface_handle = dlopen(libnomadInterface_path)
    push!(LIBPATH_list, dirname(libnomadInterface_path))

    global libnomadUtils_path = normpath(joinpath(artifact_dir, libnomadUtils_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libnomadUtils_handle = dlopen(libnomadUtils_path)
    push!(LIBPATH_list, dirname(libnomadUtils_path))

    global libsgtelib_path = normpath(joinpath(artifact_dir, libsgtelib_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libsgtelib_handle = dlopen(libsgtelib_path)
    push!(LIBPATH_list, dirname(libsgtelib_path))

    global nomad_path = normpath(joinpath(artifact_dir, nomad_splitpath...))

    push!(PATH_list, dirname(nomad_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(vcat(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ':')

end  # __init__()


text:
#####################################################################################################
# Plots corresponding to analysis region (or subregions) properties
#####################################################################################################
# Deps
using Plots
using Random
using StatsBase
using StatsPlots
using KernelDensity
using LinearAlgebra
using DelimitedFiles
default(titlefont=(14,"arial"),guidefont=(14,"arial"),tickfont=(10,"arial"),legendfont=(12,"arial"))
blnd_frnd_col = ["#999999","#E69F00","#56B4E9","#009E73","#F0E442","#0072B2","#D55E00","#CC79A7"]

#####################################################################################################
# Functions
#####################################################################################################
function read_in_data(data_dir::String,sub_max_sizes::Vector{Int64})

    # Ranges
    num_cpgs_rng = 0:1:25
    sub_size_rng = 0:20:600
    sts_sub_size = collect(sub_size_rng)
    sts_num_cpgs = collect(num_cpgs_rng)
    cent_sub_size = sts_sub_size .- step(sub_size_rng)/2.0
    cent_num_cpgs = sts_num_cpgs .- step(num_cpgs_rng)/2.0

    # Read in data
    size_histos = []
    cpgs_histos = []
    obs_sub_size = []
    for sub_size in sub_max_sizes
        println("$(sub_size)")
        # Read in data
        sub_data = readdlm("$(data_dir)/subregion_$(sub_size)_study/subregion_table.txt")
        # Get histograms
        sub_size_histo = fit(Histogram,sub_data[:,1],sub_size_rng)
        num_cpgs_histo = fit(Histogram,sub_data[:,2],num_cpgs_rng)
        # Normalize them
        sub_size_histo = normalize(sub_size_histo,mode=:probability)
        num_cpgs_histo = normalize(num_cpgs_histo,mode=:probability)
        # Push histos
        push!(size_histos,sub_size_histo)
        push!(cpgs_histos,num_cpgs_histo)
    end

    # Return
    return (cent_sub_size,size_histos),(cent_num_cpgs,cpgs_histos)

end
#####################################################################################################
# Produce histogram subregion size
#####################################################################################################

# Data dir
base_dir = "/Users/jordiabante/OneDrive - Johns Hopkins"
data_dir = "$(base_dir)/CpelNano/Data/Real-Data/Analysis-Region-Properties-hg38/"

# Produce histos
sub_max_sizes = [100,150,200,250,300,350,400,450,500]
(cent_sub_size,size_histos),(cent_num_cpgs,cpgs_histos) = read_in_data(data_dir,sub_max_sizes)

# Plot subregion sizes
p1s = []
for (i,size) in enumerate(sub_max_sizes)
    p = plot(cent_sub_size,size_histos[i].weights*100,seriestype=:bar,alpha=0.5,label="",title="$(size) bp",
        ylim=(0,100),xlabel="Analysis region size (bp)",ylabel="Percentage (%)")
    push!(p1s,p)
end

# Store plot
plot(p1s...,size=(1000,1000),layout=(3,3))
savefig("$(data_dir)/Analysis-Region-Sizes.pdf")

# Plot number of CpG sites
p2s = []
for (i,size) in enumerate(sub_max_sizes)
    p = plot(cent_num_cpgs,cpgs_histos[i].weights*100,seriestype=:bar,alpha=0.5,label="",xticks=0:2:25,
        ylim=(0,100),xlabel="Number of CpG sites",ylabel="Percentage (%)",title="$(size) bp",tickfont=(10,"arial"))
    push!(p2s,p)
end

# Store plot
plot(p2s...,size=(1000,1000),layout=(3,3))
savefig("$(data_dir)/Analysis-Region-CpG-Histogram.pdf")

ext: jl


text:
module Environments

using CSV
using DataFrames
using Random
using Setfield
using Distributions
using Plots

export HexEnv, visualize, setTiles, setTerminalStates, reset, matchRow, step, getNextAction, STEP_COST, visualizeEpsLabel, OUT_OF_BOUNDS_PENALTY, visualizeOptimalPolicy

STEP_COST = -1
OUT_OF_BOUNDS_PENALTY = -1
MOVE_PROB = 0.7
HEX_ACTIONS_NUM = 6
ANNOTATION_FONT_SIZE = 7
MARKER_SIZE = 30
FIGURE_DPI = 300
FIGURE_LOW_DPI = 100
ACTIONS = [-1 0; 0 1; 1 1; 1 0; 1 -1; 0 -1]

fig = nothing

struct HexEnv
    tiles::Matrix{Int64}
    currentState::Vector{Int64}
    vizTiles::Matrix{Float64}
    HexEnv() = new()
    HexEnv(tiles::Matrix{Int64}, currentState::Vector{Int64}, vizTiles::Matrix{Float64}) = new(tiles, currentState, vizTiles)
end;

function HexEnv(tiles::Matrix{Int64}, currentState::Vector{Int64})
    vizTiles = convert(Matrix{Float64}, tiles)
    for row in eachrow(vizTiles)
        if (row[2] % 2 == 1)
            row[1] = row[1] - 0.5
        end
        row[2] = row[2] + 0.5
    end
    return HexEnv(tiles, currentState, vizTiles)
end

function HexEnv(tiles::Matrix)
    return HexEnv(tiles, Int64[])
end

# visualize the environment with paths
function visualize(env::HexEnv, path::Any = nothing, directions::Any = nothing, directionAlpha::Any = nothing;
                    show::Bool = true, lowDpi::Bool = false, optimalPolicy::Bool = false, filename = nothing)

    width = findmax(env.vizTiles[:, 1])[1] + 0.5
    height = findmax(env.vizTiles[:, 2])[1] + 1

    if (lowDpi)
        fDpi = FIGURE_LOW_DPI
    else
        fDpi = FIGURE_DPI
    end

    negTmn = getNegRewardTerminals(env, path);
    posTmn = getPosRewardTerminals(env, path);
    normalTiles = getNormalTiles(env, path);

    global fig

    if (path !== nothing)
        if optimalPolicy
            fig = scatter(path[:, 1], path[:, 2], m = (0, :h, MARKER_SIZE), bg = :gray98, color = :dodgerblue, ticks=false,
                showaxis=false, xlims=(0, width), aspect_ratio = :equal, lab = nothing, dpi = fDpi, size = (width*60, height*60),
                markeralpha = directionAlpha);
            annotate!(fig, path[:, 1], path[:, 2], directions, :royalblue4);

            scatter!(fig, posTmn[:, 1], posTmn[:, 2], color = :darkseagreen1, m = (0, :h, MARKER_SIZE), lab = nothing);
            scatter!(fig, negTmn[:, 1], negTmn[:, 2], color = :lightsalmon, m = (0, :h, MARKER_SIZE), lab = nothing);

            annotate!(fig, posTmn[:, 1], posTmn[:, 2], [text(ann, :mediumseagreen, ANNOTATION_FONT_SIZE) for ann in posTmn[:, 3]]);
            annotate!(fig, negTmn[:, 1], negTmn[:, 2], [text(ann, :snow, ANNOTATION_FONT_SIZE) for ann in negTmn[:, 3]]);
        else
            fig = scatter(path[:, 1], path[:, 2], m = (0, :h, MARKER_SIZE), bg = :gray98, color = :wheat1, ticks=false,
                showaxis=false, xlims=(0, width), aspect_ratio = :equal, lab = nothing, dpi = fDpi, size = (width*60, height*60));
            annotate!(fig, path[:, 1], path[:, 2], directions, :wheat3);
            scatter!(fig, posTmn[:, 1], posTmn[:, 2], color = :darkseagreen1, m = (0, :h, MARKER_SIZE), lab = nothing);
            scatter!(fig, negTmn[:, 1], negTmn[:, 2], color = :lightsalmon, m = (0, :h, MARKER_SIZE), lab = nothing);
            scatter!(fig, normalTiles[:, 1], normalTiles[:, 2], color = :azure, m = (0, :h, MARKER_SIZE), lab = nothing);

            annotate!(fig, normalTiles[:, 1], normalTiles[:, 2], [text(ann, :turquoise3, ANNOTATION_FONT_SIZE) for ann in normalTiles[:, 3]]);
            annotate!(fig, posTmn[:, 1], posTmn[:, 2], [text(ann, :mediumseagreen, ANNOTATION_FONT_SIZE) for ann in posTmn[:, 3]]);
            annotate!(fig, negTmn[:, 1], negTmn[:, 2], [text(ann, :snow, ANNOTATION_FONT_SIZE) for ann in negTmn[:, 3]]);
        end

    else
        fig = scatter(negTmn[:, 1], negTmn[:, 2], m = (0, :h, MARKER_SIZE), color = :lightsalmon, bg = :gray98, ticks=false,
            showaxis=false, xlims=(0, width), aspect_ratio = :equal, lab = nothing, dpi = fDpi, size = (width*60, height*60));

        scatter!(fig, posTmn[:, 1], posTmn[:, 2], color = :darkseagreen1, m = (0, :h, MARKER_SIZE), lab = nothing);
        scatter!(fig, normalTiles[:, 1], normalTiles[:, 2], color = :azure, m = (0, :h, MARKER_SIZE), lab = nothing);
        annotate!(fig, normalTiles[:, 1], normalTiles[:, 2], [text(ann, :turquoise3, ANNOTATION_FONT_SIZE) for ann in normalTiles[:, 3]]);
        annotate!(fig, posTmn[:, 1], posTmn[:, 2], [text(ann, :mediumseagreen, ANNOTATION_FONT_SIZE) for ann in posTmn[:, 3]]);
        annotate!(fig, negTmn[:, 1], negTmn[:, 2], [text(ann, :snow, ANNOTATION_FONT_SIZE) for ann in negTmn[:, 3]]);
    end

    if show
        display(fig);
    end

    if (filename !== nothing)
        print("\nSaving current figure");
        savefig(fig, filename);
    end
end;

function getNormalTiles(env::HexEnv, path::Any = nothing)
    n = size(env.vizTiles)[1]

    if (path !== nothing)
        vizTiles = Matrix{Float64}(undef, 0, 4)
        for i in 1:n
            if (env.vizTiles[i, 4] == 0 && matchRow(path, env.vizTiles[i, 1:2]) == 0)
                vizTiles = vcat(vizTiles, env.vizTiles[i, :]')
            end
        end
    else
        vizTiles = env.vizTiles[env.vizTiles[:, 4] .== 0, :]
    end

    return vizTiles
end

function getNegRewardTerminals(env::HexEnv, path::Any = nothing)

    n = size(env.vizTiles)[1]

    if (path !== nothing)
        vizTiles = Matrix{Float64}(undef, 0, 4)
        for i in 1:n
            if (env.vizTiles[i, 4] == 1 && matchRow(path, env.vizTiles[i, 1:2]) == 0)
                vizTiles = vcat(vizTiles, env.vizTiles[i, :]')
            end
        end
    else
        vizTiles = env.vizTiles[env.vizTiles[:, 4] .== 1, :]
    end

    return vizTiles[vizTiles[:, 3] .< 0, :]
end

function getPosRewardTerminals(env::HexEnv, path::Any = nothing)
    n = size(env.vizTiles)[1]

    if (path !== nothing)
        vizTiles = Matrix{Float64}(undef, 0, 4)
        for i in 1:n
            if (env.vizTiles[i, 4] == 1 && matchRow(path, env.vizTiles[i, 1:2]) == 0)
                vizTiles = vcat(vizTiles, env.vizTiles[i, :]')
            end
        end
    else
        vizTiles = env.vizTiles[env.vizTiles[:, 4] .== 1, :]
    end

    return vizTiles[vizTiles[:, 3] .>= 0, :]
end

# randomly choose a tile from the list of tiles to start from
function reset(env::HexEnv)
    normalTiles = env.tiles[env.tiles[:, 4] .== 0, :]
    env = @set env.currentState = normalTiles[rand(1:size(normalTiles)[1]), 1:2]
end;

# check if state is in the list of tiles
function matchRow(matrix, row)
    idx = 0
    for i in 1:size(matrix)[1]
        if (matrix[i, :] == row)
            idx = i
            break
        end
    end
    return idx
end;

# transfer from one state to another
function step(env::HexEnv, action::Int64)
    nextState = env.currentState + ACTIONS[action, :]
    if (action != 1 && action != 4)
        nextState[1] -= abs((nextState[2] + 1) % 2)
    end

    # check if the next state is in the list of tiles (available in the environment)
    position = matchRow(env.tiles[:, 1:2], nextState)

    if (position != 0)

        reward = STEP_COST

        reward += env.tiles[position, 3]

        if (env.tiles[position, 4] == 1)
            return env, nextState, reward, true
        end

        return env, nextState, reward, false

    else
        return env, env.currentState, OUT_OF_BOUNDS_PENALTY, false
    end
end;

# get next action with stochastic effect
function getNextActionWithProb(action::Int64, prob::Float64)

    neighborsProb = (1 - prob)
    randNum = rand(Uniform(0., 1.))

    # |=====|=====|===================|
    #    A     B            C

    # if the random number in A
    if (randNum < neighborsProb / 2)
        action = (action - 1) % HEX_ACTIONS_NUM
        if (action == 0)
            action = HEX_ACTIONS_NUM
        end
        return action;
    # if the random number in B
    elseif (randNum < neighborsProb)
        action = (action + 1) % HEX_ACTIONS_NUM
        if (action == 0)
            action = HEX_ACTIONS_NUM
        end
        return action;
    # if the random number in C
    else
        return action;
    end

end;

# get next action
function getNextAction(currentState::Vector{Int64}, epsilon, qTable; prob::Bool=true)
    if rand(Uniform(0., 1.)) > epsilon
        if prob
            return getNextActionWithProb(findmax(qTable[currentState[1], currentState[2], :])[2], MOVE_PROB)
        else
            return findmax(qTable[currentState[1], currentState[2], :])[2]
        end
    else
        if prob
            return getNextActionWithProb(rand(1:HEX_ACTIONS_NUM), MOVE_PROB)
        else
            return rand(1:HEX_ACTIONS_NUM)
        end
    end
end;

# get arrow corresponding to action
function getArrow(action::Int64)
    if (action == 1)
        return "←"
    elseif (action == 2)
        return "↖"
    elseif (action == 3)
        return "↗"
    elseif (action == 4)
        return "→"
    elseif (action == 5)
        return "↘"
    elseif (action == 6)
        return "↙"
    else
        return "?"
    end
end;

# find the best path from start to terminal
function visualize(env::HexEnv, startState::Vector{Int64}, qTable::Array{Float64, 3}; filename::Any = nothing, show::Bool = true, prob::Bool = true, render::Bool = true)

    env = @set env.currentState = startState

    maxMove = 2 * size(env.tiles)[1];
    move = 0;

    totalReward = 0;

    terminal = false;
    stepTiles = Matrix{Float64}(undef, 0,2);
    stepDirections = String[];

    while (!terminal)

        row = convert(Vector{Float64}, env.currentState)

        if (row[2] % 2 == 1)
            row[1] = row[1] - 0.5
        end
        row[2] = row[2] + 0.5

        nextAction = getNextAction(env.currentState, 0., qTable, prob=prob)
        env, nextState, reward, terminal = step(env, nextAction)
        totalReward += reward

        position = matchRow(stepTiles, row)
        if (position == 0)
            stepTiles = vcat(stepTiles, row')
            push!(stepDirections, getArrow(nextAction))
        else
            stepDirections[position] = getArrow(nextAction)
        end

        env = @set env.currentState = nextState

        move += 1
        if (move > maxMove)
            print("\nCan't find a path")
            return
        end
    end

    if (!render)
        return totalReward;
    end

    visualize(env, stepTiles, stepDirections, show = show)

    if (terminal)
        print("\nReward: $totalReward \n")
    end

    if filename !== nothing
        print("\nSaving current figure");
        savefig(fig, filename);
    end
end;

function visualizeEpsLabel(episode::Int)
    annotate!(fig, 0, 0.5, text("Episode: $episode", ANNOTATION_FONT_SIZE, :left));
end

# visualize blank map
function visualize(tile::Vector{Int64}; show = true)
    row = convert(Vector{Float64}, tile)
    if (row[2] % 2 == 1)
        row[1] = row[1] - 0.5
    end
    row[2] = row[2] + 0.5
    scatter!([row[1]], [row[2]], color = :wheat1, m = (0, [:h], MARKER_SIZE), lab = nothing)
    if show
        display(fig);
    end
end;

# visualize the optimal policy for the given map
function visualizeOptimalPolicy(env::HexEnv, qTable::Array{Float64,3}; show = true, filename = nothing)
    normalTiles = env.tiles[env.tiles[:, 4] .== 0, :]

    tiles = Matrix{Float64}(undef, 0,2);
    tileDirections = String[];
    directionAlpha = Float64[];

    minQ = findmin(qTable)[1];
    qTable = if (minQ < 0) qTable .+ abs(minQ) else qTable end;

    maxQ = findmax(qTable)[1];

    for i in 1:size(normalTiles)[1]
        row = convert(Vector{Float64}, normalTiles[i, 1:2])
        if (row[2] % 2 == 1)
            row[1] = row[1] - 0.5
        end
        row[2] = row[2] + 0.5

        maxAction = findmax(qTable[normalTiles[i, 1], normalTiles[i, 2], :])

        push!(directionAlpha, maxAction[1] / maxQ)
        tiles = vcat(tiles, row')
        tileDirections = push!(tileDirections, getArrow(maxAction[2]))
    end

    visualize(env, tiles, tileDirections, directionAlpha, show = show, optimalPolicy = true)

    if filename !== nothing
        print("\nSaving current figure");
        savefig(fig, filename);
    end
end;

end;


text:
# abstract type PipelineType end

# abstract type Compute <: PipelineType end
# abstract type Graphics <: PipelineType end
# abstract type RayTracing <: PipelineType end

@enum PipelineType::Int8 begin
    PIPELINE_TYPE_GRAPHICS
    PIPELINE_TYPE_COMPUTE
    PIPELINE_TYPE_ASYNC_COMPUTE
end

function Vk.PipelineBindPoint(type::PipelineType)
    @match type begin
        &PIPELINE_TYPE_GRAPHICS => Vk.PIPELINE_BIND_POINT_GRAPHICS
        &(PIPELINE_TYPE_COMPUTE | PIPELINE_TYPE_ASYNC_COMPUTE) => Vk.PIPELINE_BIND_POINT_COMPUTE
    end
end

function PipelineType(bind_point::Vk.PipelineBindPoint)
    @match bind_point begin
        &Vk.PIPELINE_BIND_POINT_GRAPHICS => PIPELINE_TYPE_GRAPHICS
        &Vk.PIPELINE_BIND_POINT_COMPUTE => PIPELINE_TYPE_COMPUTE
    end
end

struct PipelineLayout <: LavaAbstraction
    handle::Vk.PipelineLayout
    descriptor_set_layouts::Vector{Vk.DescriptorSetLayout}
    push_constant_ranges::Vector{Vk.PushConstantRange}
end

struct Pipeline <: LavaAbstraction
    handle::Vk.Pipeline
    type::PipelineType
    layout::PipelineLayout
end


text:
# Autogenerated wrapper script for Zlib_jll for aarch64-apple-darwin
export libz

JLLWrappers.@generate_wrapper_header("Zlib")
JLLWrappers.@declare_library_product(libz, "@rpath/libz.1.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libz,
        "lib/libz.1.2.11.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
for func in [:clblasSrotmg, :clblasDrotmg]
    @eval @blasfun $func(SD1::cl.CL_mem, offSD1::Csize_t,
                              SD2::cl.CL_mem, offSD2::Csize_t,
                              SX1::cl.CL_mem, offSX1::Csize_t,
                              SY1::cl.CL_mem, offSY1::Csize_t,
                              SPARAM::cl.CL_mem, offSparam::Csize_t,
                              n_queues::cl.CL_uint,
                              queues::Ptr{cl.CL_command_queue},
                              n_events_in_wait_list::cl.CL_uint,
                              event_wait_list::Ptr{cl.CL_event},
                              events::Ptr{cl.CL_event})

    @eval @blasfun2 $func(SD1::cl.CL_mem, offSD1::Csize_t,
                               SD2::cl.CL_mem, offSD2::Csize_t,
                               SX1::cl.CL_mem, offSX1::Csize_t,
                               SY1::cl.CL_mem, offSY1::Csize_t,
                               SPARAM::cl.CL_mem, offSparam::Csize_t)
end


text:
@testset "rcc_test" begin
    A = [0.1, 0.5, 3., 15.]
    S = collect(0:0.1:0.9)
    C = [0.01, 0.1, 0.5, 1., 10.]
    N = [1, 5, 10, 100]
    for a in A
        for s in S
            for c in C
                for n in N
                    @test typeof(rcc(a,s,c; N=n)) <: MatrixNetwork
                end
            end
        end
    end
end


text: Dict{String,Any}("nl_end" => Dict{String,Any}("value" => "value\n","type" => "string"),"lit_nl_mid" => Dict{String,Any}("value" => "val\\nue","type" => "string"),"nl_mid" => Dict{String,Any}("value" => "val\nue","type" => "string"),"lit_nl_uni" => Dict{String,Any}("value" => "val\\ue","type" => "string"),"lit_nl_end" => Dict{String,Any}("value" => "value\\n","type" => "string"))

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
module DataAndPlotting

using HDF5
using Plots

# Wave parameters
const ki = 3
const kj = 3

# Grid size
const ni = 300
const nj = 200

# i ranges from 1...ni
# i-1 ranges from 0...ni-1
# x ranges from 0...1
function index2coord(i, ni)
    @assert ni > 1
    @assert 1 <= i <= ni
    (i - 1) / (ni - 1)
end

# Set up standing wave
export initialize
function initialize()
    arr = zeros(Float64, ni, nj)
    for j in 1:nj, i in 1:ni
        x = index2coord(i, ni)
        y = index2coord(j, nj)
        arr[i,j] = sin(2π*ki*x) * sin(2π*kj*y)
    end
    arr
end

# Write data to file
export output
function output(arr::Array{Float64, 2}, filename::String)
    h5write(filename, "wave", arr)
end

# Read data from file
export input
function input(filename::String)
    arr = h5read(filename, "wave")
end

# Create a beautiful plot
export makeplot
function makeplot(arr::Array{Float64, 2}, filename::String)
    heatmap(arr, clim=(-1.0, +1.0), color=:viridis)
    savefig(filename)
end

end


text:
import MathOptInterface
import MathOptInterface: Utilities

const MOI = MathOptInterface
const MOIU = MOI.Utilities

# indices
const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex

# function aliases
const SV = MOI.SingleVariable
const SAF = MOI.ScalarAffineFunction{Float64}
const SQF = MOI.ScalarQuadraticFunction{Float64}

# set aliases
const Bounds{T} = Union{
    MOI.EqualTo{T},
    MOI.GreaterThan{T},
    MOI.LessThan{T},
    MOI.Interval{T}
}

mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::BaronModel
    nlp_block_data::Union{Nothing, MOI.NLPBlockData}
    options
end

MOIU.@model(Model, # modelname
    (MOI.ZeroOne, MOI.Integer), # scalarsets
    (MOI.Interval, MOI.LessThan, MOI.GreaterThan, MOI.EqualTo), # typedscalarsets
    (), # vectorsets
    (), # typedvectorsets
    (MOI.SingleVariable,), # scalarfunctions
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction), # typedscalarfunctions
    (), # vectorfunctions
    () # typedvectorfunctions
)

Optimizer(; options...) = Optimizer(BaronModel(; options...), nothing, options)

# empty
function MOI.is_empty(model::Optimizer)
    BARON.is_empty(model.inner) && model.nlp_block_data === nothing
end

function MOI.empty!(model::Optimizer)
    model.inner = BaronModel(; model.options...)
    model.nlp_block_data = nothing
end

# copy
function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; kws...)
    return MOIU.automatic_copy_to(dest, src; kws...)
end

# allocate-load interface
MOIU.supports_allocate_load(model::Optimizer, copy_names::Bool) = true

MOIU.allocate(model::Optimizer, args...) = nothing

function MOIU.allocate_constraint(model::Optimizer, f::MOI.AbstractFunction, s::MOI.AbstractSet)
    constraint_info = model.inner.constraint_info
    push!(constraint_info, ConstraintInfo())
    return CI{typeof(f), typeof(s)}(length(constraint_info))
end

function MOIU.allocate_constraint(model::Optimizer, f::SV, s::MOI.AbstractSet)
    # use negative indices for variable bounds
    CI{typeof(f), typeof(s)}(-f.variable.value)
end

# optimize
function MOI.optimize!(model::Optimizer)
    write_bar_file(model.inner)
    run(`$baron_exec $(model.inner.problem_file_name)`)
    read_results(model.inner)
end

include(joinpath("moi", "util.jl"))
include(joinpath("moi", "variables.jl"))
include(joinpath("moi", "constraints.jl"))
include(joinpath("moi", "objective.jl"))
include(joinpath("moi", "results.jl"))


text:
mutable struct SelectOp <: AbstractSelectOp
    p::libgb.GxB_SelectOp
    function SelectOp(p::libgb.GxB_SelectOp)
        d = new(p)
        function f(selectop)
            libgb.GxB_SelectOp_free(Ref(selectop.p))
        end
        return finalizer(f, d)
    end
end

const SelectUnion = Union{AbstractSelectOp, libgb.GxB_SelectOp}

function SelectOp()
    return SelectOp(libgb.GxB_SelectOp_new())
end

Base.unsafe_convert(::Type{libgb.GxB_SelectOp}, selectop::SelectOp) = selectop.p

baremodule SelectOps
import ..SSGrB: load_global, SelectOp
import ..Types
import ..libgb: GB_SelectOp_opaque
end

function loadselectops()
    builtins = ["GxB_TRIL",
    "GxB_TRIU",
    "GxB_DIAG",
    "GxB_OFFDIAG",
    "GxB_NONZERO",
    "GxB_EQ_ZERO",
    "GxB_GT_ZERO",
    "GxB_GE_ZERO",
    "GxB_LT_ZERO",
    "GxB_LE_ZERO",
    "GxB_NE_THUNK",
    "GxB_EQ_THUNK",
    "GxB_GT_THUNK",
    "GxB_GE_THUNK",
    "GxB_LT_THUNK",
    "GxB_LE_THUNK"]
    for name ∈ builtins
        simple = Symbol(string(name[5:end]))
        constquote = quote
            const $simple = SelectOp(load_global($name, GB_SelectOp_opaque))
        end
        @eval(SelectOps, $constquote)
    end
end
isloaded(d::AbstractSelectOp) = d.p !== C_NULL
function Base.getindex(o::AbstractSelectOp)
    o.p
end
function validtypes(::AbstractSelectOp)
    return nothing
end

function printtest(io::IO, M::SelectOp)
    str = mktemp() do _, f
        cf = Libc.FILE(f)
        libgb.GxB_SelectOp_fprint(M, "", libgb.GxB_SHORT, cf)
        flush(f)
        seekstart(f)
        x = read(f, String)
        close(cf)
        x
    end
    print(io, str)
end

function Base.show(io::IO, d::SelectOp)
    printtest(io, d)
end


text:
#InspectDR: Operation on datasets
#-------------------------------------------------------------------------------

#==Types declarations
===============================================================================#

#Input dataset data:
#TODO: make x an immutable (read-only) type so data cannot be changed once identified.
type IDataset{DATAF1} #DATAF1::Bool: Data is function of 1 argument (optimization possible)
	x::Vector
	y::Vector
end

#==Accessors
===============================================================================#

Point(ds::IDataset, i::Int) = Point2D(ds.x[i], ds.y[i])
Point(ds::Vector{Point2D}, i::Int) = ds[i]

#TODO: Specialize for sorted datasets?
function getextents(ds::IDataset)
	(xmin, xmax) = extrema_nan(ds.x)
	(ymin, ymax) = extrema_nan(ds.y)
	return PExtents2D(xmin, xmax, ymin, ymax)
end

#==Pre-processing display data
===============================================================================#

#_reduce: Obtain reduced dataset by limiting to the extents & max resolution:
#-------------------------------------------------------------------------------
#Generic algorithm... Just transfer all data for now
#    xres_max: Max number of x-points in data window.
#returns: Vector{Point2D}
#TODO: clip data beyond extents.
#WARNING: not clipping might cause display issues when applying the transform
function _reduce(input::IDataset, ext::PExtents2D, xres_max::Integer)
	x = input.x; y = input.y
	n_ds = length(x) #numer of points of input dataset
	result = Array(Point2D, n_ds)
	for i in 1:n_ds
		result[i] = Point2D(x[i], y[i])
	end
	return result
end

#Optimized for functions of 1 argument:
#    xres_max: Max number of x-points in data window.
#returns: Vector{Point2D}
function _reduce(input::IDataset{true}, ext::PExtents2D, xres_max::Integer)
	xres = (ext.xmax - ext.xmin)/ xres_max
	const min_lookahead = 2
	const thresh_xres = (min_lookahead+1)*xres
	x = input.x; y = input.y
	n_ds = length(x) #numer of points of input dataset
	sz = min(n_ds, xres_max)+4 #Add 4 pts, just in case (TODO: fix)
	result = Array(Point2D, sz)
	n = 0 #Number of points in reduced dataset
	i = 1 #Index into input dataset

	if length(x) != length(y)
		error("x & y - vector length mismatch.")
	elseif length(x) < 1
		resize!(result, 0)
		return result
	end

	#Discard data before visible extents:
	while i < n_ds
		if x[i] > ext.xmin; break; end
		i+=1
	end

	i = max(i-1, 1)
	lastx = x[i]
	lasty = y[i]
	n+=1
	result[n] = Point2D(lastx, lasty)
	i+=1

	while i <= n_ds
		if lastx >= ext.xmax; break; end
		xthresh = min(lastx+thresh_xres, ext.xmax)
		ilahead_start = i+min_lookahead
		ilahead = min(ilahead_start, n_ds)
		while ilahead < n_ds
			if x[ilahead] >= xthresh; break; end
			ilahead += 1
		end

		#@assert(ilahead<=nds)

		if ilahead > ilahead_start #More data than xres allows
			#TODO: make this a function??
			#"Internal limits":
			(ymin_i, ymax_i) = extrema(y[i:(ilahead-1)])
			p1 = Point2D(x[ilahead-1], y[ilahead-1])
			p2 = Point2D(x[ilahead], y[ilahead])
			nexty = interpolate(p1, p2, xthresh)

			#"External limits:
			(ymin, ymax) = minmax(lasty, nexty)

			#Add points representing "internal limits"
			#(if they exceed external):
			yint = [ymin_i, ymax_i]
			ysel = Bool[false, false]
			ysel[1] = ymin_i < ymin
			ysel[2] = ymax_i > ymax
#ysel = Bool[true, true] #Debug: add points no matter what
			curx = lastx + 1.5*xres #TODO: resize to have steps of 1xres?
#@show ilahead-i
			offset = lasty < nexty ?0 :1 #Add min or max first??
			for j in (offset+(1:2))
				idx = 1+(j&0x1)
				if !ysel[idx]; continue; end #Only add points if desired
				n+=1;
				result[n] = Point2D(curx,yint[idx])
			end
			#Done adding points

			n+=1;
			lastx = xthresh
			lasty = nexty
			result[n] = Point2D(lastx, lasty)
			i = ilahead
		else #Plot actual data points
			while i <= ilahead
				n += 1
				lastx = x[i]
				lasty = y[i]
				result[n] = Point2D(lastx, lasty)
				i += 1
			end
		end
	end

	resize!(result, n)
	return result
end

#Last line


text:
# Hello world example, similar to the Boost.Python hello world

using CxxWrap
using Base.Test

# Wrap the functions defined in C++
wrap_modules("build/cxxwrap/libhello")

# Output:
Hello.hello()


text:
using LinearAlgebra
using Polynomials

include("../tools/plotting.jl")
"""
Use Taylor's Theorem with Remainder to check
validity of computed Jacobian of f: R^n -> R^m.

More specifically, check that error of linear approximation decays like O(h^2)

last: smallest power of exp^(-x) we wish to compute in scheme

Requirements: derivative/jacobian must be in the right 2D shape, and not be a 1D list

Hessian is only for functions from R^n -> R^1
"""
function checkDerivative(f, df, x0, hessian = nothing, first = 3, last = 12, num = 10)
    if size(x0, 2) != 1
        x0 = reshape(x0, length(x0), 1) #reshape x0 into 2D column vector or 1D Array
    end
    f0 = f(x0)
    df0 = df(x0)
    n = length(x0)
    m = length(f0)
    @assert size(f0, 2)==1 #no reason for f0 to be 2D
    #println(size(df0))
    #println("m: ", m)
    #println("n: ", n)
    #@assert size(df0) == (m, n) ||  size(df0) == (m, ) #keep 1D list as possibility
    #println("size df0, 1: ", size(df0, 1))
    #println("size df0, 2: ", size(df0, 2))
    try
        @assert size(df0, 1) == m && size(df0, 2) == n
    catch e
        #@info "df0", df0
        #@info "m", m
        #@info "n", n
    end

    try
        df0 = reshape(df0, m, n) #reshapes to m x n if non-scalar
    catch e
    end

    if hessian!=nothing
        @assert m == 1 #we can only compute the Hessian for functions to R^1
        d2f0 = hessian(x0)
        #println(size(d2f0))
        @assert size(d2f0, 1) == n && size(d2f0, 1) ==n
    end
    #what if x0 or f0 are constants?
    dx = rand(size(x0, 1), size(x0, 2))
    if length(f0) >1
        f0 = reshape(f0, length(f0), 1) #reshape into col vector
    end
    h = collect(first:(last-first)/num:last)
    for i=1:length(h)
        h[i] = exp(-h[i])
    end
    A = zeros(length(h))
    for i = 1:length(h)
        fi = f(x0 .+ h[i]*dx) #dx has the same shape x0, so we can safely do .+, which allows us to add scalar and array
                               #if dx did not have same shape as x0, then .+ could have unintended broadcasting effects
        @assert size(fi, 1) == size(f0, 1) && size(fi, 2) == size(f0, 2) #this lets us do fi .- f0 safely
        if false#debug
            println("x0: ", x0)
            println("dx: ", dx)
            println("fi: ", fi)
            println("f0: ", f0)
            println("df0: ", df0)
            println("increment", h[i]*dx)
            println("h[i]: ", h[i])
        end
            if hessian!=nothing
                inc = h[i] .* dx
                expr1 = fi .- f0
                expr2 = df0 * inc .+ 0.5* inc' * d2f0 * inc
                #println("size expr1: ", size(expr1))
                #println("size expr2: ", size(expr2))
                @assert size(expr1, 1) == size(expr2, 1) && size(expr1, 2) == size(expr2, 2)
                A[i] = norm(expr1 .- expr2)
                #A[i] = norm((fi - f0) - df0' * inc .- 0.5* inc' * d2f0 * inc)
            else
                expr1 = fi .- f0
                expr2 = df0 * (h[i] .* dx)
                #println("size expr1: ", size(expr1))
                #println("size expr2: ", size(expr2))
                @assert size(expr1, 1) == size(expr2, 1) && size(expr1, 2)==size(expr2, 2)
                A[i] = norm(expr1 .- expr2)
                #catch DimensionMismatch #we catch the case when f: R^1 -> R^n, in which case  df0'*dx will yield an error
                #    println("Warning: finitedifference dimension mismatch")
                #    @assert size(expr1) == size(expr2)
                #    A[i] = norm(expr1 - expr2)
                #end
        end
    end
    r1 = log.(h)
    r2 = log.(A)
    plt = Plots.plot(r1, r2, title = "Finite Difference Derivative Checker", xlabel = "log of h", ylabel = "log of error")#, fontfamily = font(48, "Courier") , reuse = false)
    return (r1, r2, plt, fit(r1, r2, 1))
end

"""
df is an in-place function (has no return value)
"""
function checkDerivative_in_place(f, df, x0, store, hessian = nothing, first = 3, last = 12, num = 10)
    if size(x0, 2) != 1
        x0 = reshape(x0, length(x0), 1) #reshape x0 into 2D column vector or 1D Array
    end
    f0 = f(x0)
    #df0 = df(x0)
    df(store, x0) #result in store
    n = length(x0)
    m = length(f0)
    @assert size(f0, 2)==1 #no reason for f0 to be 2D
    #println(size(df0))
    #println("m: ", m)
    #println("n: ", n)
    #@assert size(df0) == (m, n) ||  size(df0) == (m, ) #keep 1D list as possibility
    #println("size df0, 1: ", size(df0, 1))
    #println("size df0, 2: ", size(df0, 2))
    try
        @assert size(store, 1) == m && size(store, 2) == n
    catch e
        #@info "df0", store
        #@info "m", m
        #@info "n", n
        try
            store = reshape(store, m, n) # force reshape to m x n if non-scalar
        catch e
        end
    end
    if hessian!=nothing
        @assert m == 1 #we can only compute the Hessian for functions to R^1
        d2f0 = hessian(x0)
        #println(size(d2f0))
        @assert size(d2f0, 1) == n && size(d2f0, 1) ==n
    end
    #what if x0 or f0 are constants?
    dx = rand(size(x0, 1), size(x0, 2))
    if length(f0) >1
        f0 = reshape(f0, length(f0), 1) #reshape into col vector
    end
    h = collect(first:(last-first)/num:last)
    for i=1:length(h)
        h[i] = exp(-h[i])
    end
    A = zeros(length(h))
    for i = 1:length(h)
        fi = f(x0 .+ h[i]*dx) #dx has the same shape x0, so we can safely do .+, which allows us to add scalar and array
                               #if dx did not have same shape as x0, then .+ could have unintended broadcasting effects
        @assert size(fi, 1) == size(f0, 1) && size(fi, 2) == size(f0, 2) #this lets us do fi .- f0 safely
        if false#debug
            println("x0: ", x0)
            println("dx: ", dx)
            println("fi: ", fi)
            println("f0: ", f0)
            println("df0: ", store)
            println("increment", h[i]*dx)
            println("h[i]: ", h[i])
        end
            if hessian!=nothing
                inc = h[i] .* dx
                expr1 = fi .- f0
                expr2 = store * inc .+ 0.5* inc' * d2f0 * inc
                #println("size expr1: ", size(expr1))
                #println("size expr2: ", size(expr2))
                @assert size(expr1, 1) == size(expr2, 1) && size(expr1, 2) == size(expr2, 2)
                A[i] = norm(expr1 .- expr2)
                #A[i] = norm((fi - f0) - df0' * inc .- 0.5* inc' * d2f0 * inc)
            else
                expr1 = fi .- f0
                expr2 = store * (h[i] .* dx)
                #println("size expr1: ", size(expr1))
                #println("size expr2: ", size(expr2))
                @assert size(expr1, 1) == size(expr2, 1) && size(expr1, 2)==size(expr2, 2)
                A[i] = norm(expr1 .- expr2)
                #catch DimensionMismatch #we catch the case when f: R^1 -> R^n, in which case  df0'*dx will yield an error
                #    println("Warning: finitedifference dimension mismatch")
                #    @assert size(expr1) == size(expr2)
                #    A[i] = norm(expr1 - expr2)
                #end
        end
    end
    r1 = log.(h)
    r2 = log.(A)
    plt = Plots.plot(r1, r2, title = "Finite Difference Derivative Checker", xlabel = "log of h", ylabel = "log of error")#, fontfamily = font(48, "Courier") , reuse = false)
    return (r1, r2, plt, fit(r1, r2, 1))
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# Autogenerated wrapper script for Doxygen_jll for x86_64-unknown-freebsd
export doxygen

using Libiconv_jll
JLLWrappers.@generate_wrapper_header("Doxygen")
JLLWrappers.@declare_executable_product(doxygen)
function __init__()
    JLLWrappers.@generate_init_header(Libiconv_jll)
    JLLWrappers.@init_executable_product(
        doxygen,
        "bin/doxygen",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module TestJoin
    using Base.Test
    using DataFrames

    name = DataFrame(ID = [1, 2, 3], Name = ["John Doe", "Jane Doe", "Joe Blogs"])
    job = DataFrame(ID = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])

    # Join on symbols or vectors of symbols
    join(name, job, on = :ID)
    join(name, job, on = [:ID])

    # Soon we won't allow natural joins
    #@test_throws join(name, job)

    # Test output of various join types
    outer = DataFrame(ID = [1, 2, 3, 4],
                      Name = @data(["John Doe", "Jane Doe", "Joe Blogs", NA]),
                      Job = @data(["Lawyer", "Doctor", NA, "Farmer"]))

    # (Tests use current column ordering but don't promote it)
    right = outer[!isna(outer[:Job]), [:Name, :ID, :Job]]
    left = outer[!isna(outer[:Name]), :]
    inner = left[!isna(left[:Job]), :]
    semi = inner[:, [:ID, :Name]]
    anti = left[isna(left[:Job]), [:ID, :Name]]

    @test isequal(join(name, job, on = :ID), inner)
    @test isequal(join(name, job, on = :ID, kind = :inner), inner)
    @test isequal(join(name, job, on = :ID, kind = :outer), outer)
    @test isequal(join(name, job, on = :ID, kind = :left), left)
    @test isequal(join(name, job, on = :ID, kind = :right), right)
    @test isequal(join(name, job, on = :ID, kind = :semi), semi)
    @test isequal(join(name, job, on = :ID, kind = :anti), anti)

    # Join on multiple keys
    df1 = DataFrame(A = 1, B = 2, C = 3)
    df2 = DataFrame(A = 1, B = 2, D = 4)

    join(df1, df2, on = [:A, :B])

    # Test output of cross joins
    df1 = DataFrame(A = 1:2, B = 'a':'b')
    df2 = DataFrame(A = 1:3, C = 3:5)

    cross = DataFrame(A = [1, 1, 1, 2, 2, 2],
                      B = ['a', 'a', 'a', 'b', 'b', 'b'],
                      C = [3, 4, 5, 3, 4, 5])

    @test join(df1, df2[[:C]], kind = :cross) == cross

    # Cross joins handle naming collisions
    @test size(join(df1, df1, kind = :cross)) == (4, 4)

    # Cross joins don't take keys
    @test_throws ArgumentError join(df1, df2, on = :A, kind = :cross)
end


text:
module ParticlePlots

using Plots

# functions from partViz
export plotWalls!,
       plotParticles!

using StokesParticles

function plotWalls!(p,wallList::Vector)
    N = length(wallList)

    for i=1:N
      plotWall!(p,wallList[i])
    end

    return p
  end

function plotParticles!(plotObj,pList)
    colorInd = [:black,:red,:green,:yellow,:purple,:orange,:blue]
    N = length(pList)

    for i=1:N
        if mod(i,20) == 0
            tc = mod(i,length(colorInd))+1
            plotParticle!(plotObj,pList[i],colorInd[tc])
        else
            plotParticle!(plotObj,pList[i],:white)
        end
    end

    nothing
end

function plotParticle!(plotObj::Plots.Plot,particle,colorSym::Symbol)
    θ = range(0,stop=2*pi,length=40)
    xs = particle.pos.x .+ particle.radius*cos.(θ)
    ys = particle.pos.y .+ particle.radius*sin.(θ)

    plot!(plotObj,xs,ys,label="",
        seriestype=[:shape,],color=colorSym,
        linewidth=1,linecolor=:black)
end

function plotWall!(plotObj::Plots.Plot,wall::StokesParticles.LineWall)
    radius = wall.thickness/2

    x1 = wall.nodes[1].x; x2 = wall.nodes[2].x
    y1 = wall.nodes[1].y; y2 = wall.nodes[2].y
    θWALL = atan(y2-y1,x2-x1)

    Nθ = 15
    θ1 = range(-pi/2 + θWALL,stop=pi/2 + θWALL,length=Nθ)
    θ2 = range( pi/2 + θWALL,stop = 3*pi/2 + θWALL,length=Nθ)

    nx = wall.n[1]
    ny = wall.n[2]

    # 1 arc
    x = [x2 + radius*cos(θ1[i]) for i=1:Nθ]
    y = [y2 + radius*sin(θ1[i]) for i=1:Nθ]
    plot!(plotObj,x,y,label="",color=:black)

    # 2 line
    x = [x2 + radius*nx,x1 + radius*nx]
    y = [y2 + radius*ny,y1 + radius*ny]
    plot!(plotObj,x,y,label="",color=:black)

    # 3 arc
    x = [x1 + radius*cos(θ2[i]) for i=1:Nθ]
    y = [y1 + radius*sin(θ2[i]) for i=1:Nθ]
    plot!(plotObj,x,y,label="",color=:black)

    # 4 line
    x = [x1 - radius*nx,x2 - radius*nx]
    y = [y1 - radius*ny,y2 - radius*ny]
    plot!(plotObj,x,y,label="",color=:black)
end
function plotWall!(plotObj::Plots.Plot,wall::StokesParticles.CircleWall)
    radius = wall.thickness/2
    Nplot = 50; Δθ = 2*pi/(Nplot-1)
    curvex = [wall.center.x + (wall.radius + radius)*cos((i-1)*Δθ) for i=1:Nplot]
    curvey = [wall.center.y + (wall.radius + radius)*sin((i-1)*Δθ) for i=1:Nplot]
    plot!(plotObj,curvex,curvey,label="",color=:black)
    curvex = [wall.center.x + (wall.radius - radius)*cos((i-1)*Δθ) for i=1:Nplot]
    curvey = [wall.center.y + (wall.radius - radius)*sin((i-1)*Δθ) for i=1:Nplot]
    plot!(plotObj,curvex,curvey,label="",color=:black)
end
function plotWall!(plotObj::Plots.Plot,wall::StokesParticles.ArcWall)
    radius = wall.thickness/2

    pa = [wall.nodes[1].x,wall.nodes[1].y]
    c  = [wall.nodes[2].x,wall.nodes[2].y]
    pb = [wall.nodes[3].x,wall.nodes[3].y]
    r = sqrt((pa[1]-c[1])^2 + (pa[2]-c[2])^2)
    θ1 = atan(pa[2]-c[2],pa[1]-c[1])
    θ2 = atan(pb[2]-c[2],pb[1]-c[1])

    if θ1 > θ2
        θ2 += 2*pi
    end

    Nplot = 50; Δθ = (θ2 - θ1)/(Nplot-1)

    # 1 big arc -- done
    curvex = [c[1] + (r+radius)*cos((i-1)*Δθ+θ1) for i=1:Nplot]
    curvey = [c[2] + (r+radius)*sin((i-1)*Δθ+θ1) for i=1:Nplot]
    plot!(plotObj,curvex,curvey,label="",color=:black)

    # 2 far cap
    θWALL = θ2 + pi/2
    θ2cap = range(-pi/2 + θWALL, stop=pi/2 + θWALL,length=Nplot)
    x = [pb[1] + radius*cos(θ2cap[i]) for i=1:Nplot]
    y = [pb[2] + radius*sin(θ2cap[i]) for i=1:Nplot]
    plot!(plotObj,x,y,label="",color=:black)

    # 3 small arc
    curvex = [c[1] + (r-radius)*cos((i-1)*Δθ+θ1) for i=1:Nplot]
    curvey = [c[2] + (r-radius)*sin((i-1)*Δθ+θ1) for i=1:Nplot]
    plot!(plotObj,curvex,curvey,label="",color=:black)

    # 4 close cap
    θWALL = θ1 + pi/2
    θ1cap = range( pi/2 + θWALL,stop = 3*pi/2 + θWALL,length=Nplot)
    x = [pa[1] + radius*cos(θ1cap[i]) for i=1:Nplot]
    y = [pa[2] + radius*sin(θ1cap[i]) for i=1:Nplot]
    plot!(plotObj,x,y,label="",color=:black)
end

end # ParticlePlots


text:
using Waterfall
import Plots

# This file generates a histogram for each step in the cascade. It shows how a step HEIGHT
# depends on plotting order.
include(joinpath(WATERFALL_DIR,"src","includes.jl"))
include(joinpath(WATERFALL_DIR,"bin","io.jl"))

nperm = 1000
x = define_permute(df[[1:11;15],:], nperm; kwargs...)

data = collect_data(x)
N = length(data)

v = get_value(data)
lab = get_label.(data)
sublab = uppercase.(get_sublabel.(data))
title = ["$(sublab[ii]) ($(lab[ii]))" for ii in 1:N]

ii=6
for ii in 2:11
    Plots.histogram(v[ii,:];
        title=title[ii],
        xlabel="Efficiency (%)",
        color=COLORCYCLE[ii-1],
        legend=false,
        ylim=(0,1000),
        # xlim=(-2.,-1.,),
    )
    Plots.png(joinpath("fig","hist",name(ii-1)))
end

for ii in [1,N]
    Plots.histogram(v[ii,:];
            title=title[ii],
            xlabel="Efficiency (%)",
            color="black",
            legend=false,
            ylim=(0,1000),
            # xlim=(-2.,-1.,),
        )
    Plots.png(joinpath("fig","hist",name(ii-1)))
end


text:
n_features = (28, 28)
n_targets = 1

@testset "trainset" begin
    d = MNIST()

    @test d.split == :train
    @test extrema(d.features) == (0, 1)
    @test convert2image(d, 1) isa AbstractMatrix{<:Gray}
    @test convert2image(d, 1:2) isa AbstractArray{<:Gray, 3}

    test_supervised_array_dataset(d;
        n_features, n_targets, n_obs=60000,
        Tx=Float32, Ty=Int,
        conv2img=true)
end

@testset "testset" begin
    d = MNIST(split=:test, Tx=UInt8)

    @test d.split == :test
    @test extrema(d.features) == (0, 255)
    @test convert2image(d, 1) isa AbstractMatrix{<:Gray}

    test_supervised_array_dataset(d;
        n_features, n_targets, n_obs=10000,
        Tx=UInt8, Ty=Int,
        conv2img=true)

end


text:
using OpticalPropagation
using Documenter

makedocs(;
    modules=[OpticalPropagation],
    authors="Yong-an Lu <miroox@outlook.com>",
    repo="https://github.com/HIT-UOI-SR/OpticalPropagation.jl/blob/{commit}{path}#L{line}",
    sitename="OpticalPropagation.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://HIT-UOI-SR.github.io/OpticalPropagation.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "2D Light Field" => "lightfield2d.md",
        "Optical Propagation Methods" => "methods.md",
    ],
)

deploydocs(;
    repo="github.com/HIT-UOI-SR/OpticalPropagation.jl",
)


text:
# implement julia_commander for showing docstrings

function parse_doc(doc::AbstractString)
    doc = replace(doc, r"\n\n\n+" => "\n\n")
    capture_ranges = [
        (m.offset):(m.offset + ncodeunits(m.match) - 1) for
        m in eachmatch.(r"(```.+?```)"s, doc)
    ]
    for m in eachmatch(r"(^|\n)(#+ |!!! )(.*)\n", doc)
        if all(rg -> m.offset ∉ rg, capture_ranges)
            if m.captures[2] == "# "
                doc = replace(
                    doc,
                    m.match =>
                        m.captures[1] *
                        "**" *
                        m.captures[3] *
                        "**\n" *
                        "≡"^(length(m.captures[3]));
                    count=1,
                )
            elseif m.captures[2] == "!!! "
                doc = replace(
                    doc, m.match => m.captures[1] * "__" * m.captures[3] * "__\n"; count=1
                )
            else
                doc = replace(
                    doc,
                    m.match =>
                        m.captures[1] *
                        "*" *
                        m.captures[3] *
                        "*\n" *
                        "-"^(length(m.captures[3]));
                    count=1,
                )
            end
        end
    end
    doc = replace(doc, r"(```.+)\n" => "```julia\n")
    doc = replace(doc, "```\n\n" => "```\n")
    for m in eachmatch(r"\[([^ ]*)\]\(@ref\)", doc)
        doc = replace(doc, m.match => "`" * m.captures[1] * "`"; count=1)
    end
    return doc
end

function commander(c::Client, m::Message, ::Val{:julia_doc})
    # @info "julia_commander called"
    # @info "Message content" m.content m.author.username m.author.discriminator
    startswith(m.content, COMMAND_PREFIX * "j ") ||
        startswith(m.content, COMMAND_PREFIX * "j? ") ||
        m.content == COMMAND_PREFIX * "j" ||
        return nothing
    regex = Regex(
        COMMAND_PREFIX *
        raw"j(\?| help| doc| packages| package| stats| top| bottom)? *(.*)$",
    )
    matches = match(regex, m.content)
    if matches === nothing || matches.captures[1] in (" help", nothing)
        help_commander(c, m, Val(:julia_doc))
    elseif matches.captures[1] in ("?", " doc")
        handle_julia_help_commander(c, m, matches.captures[2])
    elseif matches.captures[1] == " packages"
        handle_julia_package_list(c, m)
    elseif matches.captures[1] == " package"
        handle_julia_names_in_package(c, m, matches.captures[2])
    elseif matches.captures[1] in (" stats", " top", " bottom")
        handle_doc_stats(c, m, matches.captures[1], matches.captures[2])
    else
        reply(c, m, "Sorry, I don't understand the request; use `j help` for help")
    end
    return nothing
end

function help_commander(c::Client, m::Message, ::Val{:julia_doc})
    # @info "Sending help for message" m.id m.author
    reply(
        c,
        m,
        """
        The `j` command shows the docstring of keywords, names in Base and names in other selected packages, as well as statistics on the names that have been queried. Only names with meaningful help information have been recorded from each package.

        To see which packages have been added, use `j packages`. If you would like a particular package to be added, just let us know.

        Here is the list of all available `j` commands and their use:
        ```
        j help
        j? <name>
        j doc <name>
        j packages
        j package <PkgName>
        j stats [<name>]
        j top [<number>]
        j bottom [<number>]
        ```
        * `j help` returns this help.
        * `j packages` shows which packages are available with names.
        * `j package <PkgName>` shows which names have bee recorded from to package `<PkgName>`.
        * `j? <name>` and `j doc <name>` return the documentation for `<name>`
        * `j stats <name>` return how many times the docstring for `name` has been queried.
        * `j top [<number>]` return the top 10 names that have been queried (or top `number`, if given).
        * `j bottom [<number>]` return the bottom 10 names that have been queried (or bottom `number`, if given).
        """,
    )
    return nothing
end

function handle_julia_help_commander(c::Client, m::Message, name::AbstractString)
    # @info "julia_help_commander called"
    try
        name = strip(name)
        name = replace(name, r"\s{2,}" => " ")
        channel = @discord get_channel(c, m.channel_id)
        all_names = load_names(
            joinpath(@__DIR__, "..", "..", "data", "docs", "all_names.json")
        )
        pkg_in_name, name =
            occursin('.', name) ? u -> (first(u), last(u))(split(name, '.')) : ("", name)
        all_docs = load_docs(
            joinpath(@__DIR__, "..", "..", "data", "docs", "all_docs.json")
        )
        if name in keys(all_names)
            doc_in_pkgs = Dict{String,String}()
            for pkg in keys(all_names[name])
                if pkg_in_name in ("", pkg)
                    push!(doc_in_pkgs, pkg => all_docs[pkg][name][2])
                end
            end

            if length(doc_in_pkgs) > 0
                for (k, v) in doc_in_pkgs
                    if k in ("Base", "Keywords")
                        doc = v
                    else
                        doc = "*In Package* `$k`:\n" * v
                    end
                    doc = parse_doc(doc)
                    docs = split_message(
                        doc; extrastyles=[r"\n.*\n≡.+\n", r"\n.*\n-+\n", r"[^\s]+"]
                    )
                    for doc_chunk in docs
                        # @info doc_chunk
                        reply(c, m, doc_chunk)
                    end
                end

                count = update_names_count("namescount", name, m.channel_id, channel.name)
                reply(c, m, "*(Count number for `$name`: $count)*")
            else
                reply(c, m, "No documentation for `$name` found in package `$pkg`")
            end
        else
            reply(c, m, "No documentation found.")
        end
    catch ex
        @info ex
        reply(c, m, "Sorry, it didn't work.")
    end
    return nothing
end

function handle_doc_stats(
    c::Client, m::Message, captured1::AbstractString, captured2::AbstractString
)
    # @info "handle_doc_stats called"
    captured2 = strip(captured2)
    captured2 = replace(captured2, r"\s{2,}" => " ")
    try
        if captured1 == " stats"
            if captured2 == ""
                statsmgs = stats_namescount("namescount")
            else
                statsmgs = stats_namescount("namescount"; name=captured2)
            end
        elseif captured1 in (" top", " bottom")
            if captured2 != "" && all(isdigit, captured2)
                statsmgs = stats_namescount(
                    "namescount";
                    place=strip(captured1),
                    number=max(1, parse(Int, captured2)),
                )
            else
                statsmgs = stats_namescount("namescount"; place=strip(captured1))
            end
        end
        reply(c, m, statsmgs)
    catch ex
        @info ex
        reply(c, m, "Sorry, it didn't work.")
    end
    return nothing
end

function handle_julia_package_list(c::Client, m::Message)
    all_docs_filename = joinpath(@__DIR__, "..", "..", "data", "docs", "all_docs.json")
    if isfile(all_docs_filename)
        all_docs = load_docs(all_docs_filename)
        msg =
            "Besides the **keywords** and **Base**, " *
            "there are $(length(all_docs)-2) packages " *
            "available with recorded names:\n\n" *
            join(sort(collect(keys(all_docs))), ", ") *
            "."
    else
        msg = "No packages available."
    end
    return reply(c, m, msg)
end

function handle_julia_names_in_package(c::Client, m::Message, pkg::AbstractString)
    pkg = strip(pkg)
    pkg = replace(pkg, r"\s{2,}" => " ")
    # @info pkg
    all_names_filename = joinpath(@__DIR__, "..", "..", "data", "docs", "all_names.json")
    all_docs_filename = joinpath(@__DIR__, "..", "..", "data", "docs", "all_docs.json")
    if isfile(all_names_filename) && isfile(all_docs_filename)
        all_names = load_names(all_names_filename)
        pkg_list = [name for (name, d) in all_names if pkg in keys(d)]
        pkg_list_exported = [
            name for name in pkg_list if all_names[name][pkg] == "exported"
        ]
        pkg_list_nonexported = [
            name for name in pkg_list if all_names[name][pkg] == "nonexported"
        ]
        # @info pkg_list
        if length(pkg_list) > 0
            msg = ""
            if length(pkg_list_exported) == 0
                msg *= "**There are no exported names on record from package `$pkg`:**\n\n"
            elseif length(pkg_list_exported) == 1
                msg *= "**There is 1 exported name on record from package `$pkg`:**\n\n"
            else
                msg *=
                    "**There are $(length(pkg_list_exported)) exported names " *
                    "on record from package `$pkg`:**\n\n"
            end
            if length(pkg_list_exported) > 0
                msg *= join(sort(pkg_list_exported), ", ") * ".\n\n"
            end
            if length(pkg_list_nonexported) == 0
                msg *= "**There aren't any nonexported names:**\n\n"
            elseif length(pkg_list_nonexported) == 1
                msg *= "**There is 1 nonexported name:**\n\n"
            else
                msg *=
                    "**There are $(length(pkg_list_nonexported)) " *
                    "non-exported names:**\n\n"
            end
            if length(pkg_list_nonexported) > 1
                msg *= join(sort(pkg_list_nonexported), ", ") * "."
            end
        else
            msg = "No names found in `$pkg`."
        end
    else
        msg = "No packages available."
    end

    msg = split_message(msg)
    for msg_chunk in msg
        reply(c, m, msg_chunk)
    end
end


text:
xmin(g::GridTopology) = minimum(g.xs)
ymin(g::GridTopology) = minimum(g.ys)
xcellsize(g::GridTopology) = step(g.xs)
ycellsize(g::GridTopology) = step(g.ys)
xcells(g::GridTopology) = length(g.xs)
ycells(g::GridTopology) = length(g.ys)
indices(g::SEGrid) = g.indices
indices(g::SEGrid, idx) = g.indices[:, idx]
boundingbox(g::GridTopology) = Bbox(xmin(g), xmax(g), ymin(g), ymax(g))

@forward_func GridData.grid xmin, ymin, xcellsize, ycellsize, xcells, ycells, boundingbox
@forward_func Locations{GridData}.coords xmin, ymin, xcellsize, ycellsize, xcells, ycells, boundingbox, indices
@forward_func SubGridData.grid xmin, ymin, xcellsize, ycellsize, xcells, ycells, boundingbox
@forward_func SubLocations{SubGridData}.coords xmin, ymin, xcellsize, ycellsize, xcells, ycells, boundingbox, indices
@forward_func SEAssemblage.site xmin, ymin, xcellsize, ycellsize, xcells, ycells, boundingbox, indices

show(io::IO, b::Bbox) = println(io, "xmin:\t$(b.xmin)\nxmax:\t$(b.xmax)\nymin:\t$(b.ymin)\nymax:\t$(b.ymax)")
show(io::IO, g::GridData) = (println(io, "Spatial Grid\n", gridline(g.grid)); println(io, size(g.indices,1), " sites");)
show(io::IO, g::GridTopology) = println(io, "Spatial GridTopology\n", gridline(g))
gridline(g::GridTopology) =
    """
       lower left : $(xmin(g)), $(ymin(g))
       cellsizes  : $(xcellsize(g)), $(ycellsize(g))
       size       : $(xcells(g)), $(ycells(g))
    """


text:
module js

import GR

pxwidth = 640
pxheight = 480

@static if VERSION < v"0.7.0-DEV.4762"
    macro cfunction(f, rt, tup)
        :(Base.cfunction($(esc(f)), $(esc(rt)), Tuple{$(esc(tup))...}))
    end
end

id_count = 0
js_running = false
dispid = 0

mutable struct JSTermWidget
    identifier::Int
    width::Int
    height::Int
    display::Int
    visible::Bool
end

function inject_js()
  global comm
  comm = nothing
  _js_fallback = "https://gr-framework.org/downloads/gr-latest.js"
  _gr_js = if isfile(joinpath(ENV["GRDIR"], "lib", "gr.js"))
    _gr_js = try
      _gr_js = open(joinpath(ENV["GRDIR"], "lib", "gr.js")) do f
        _gr_js = read(f, String)
        _gr_js
      end
    catch e
      nothing
    end
    _gr_js
  end
  if _gr_js === nothing
      error(string("Unable to open '", joinpath(ENV["GRDIR"], "lib", "gr.js"), "'."))
  else
      display(HTML(string("""

      """)))
  end
end

function JSTermWidget(id::Int, width::Int, height::Int, display::Int)
  global id_count, js_running
  if GR.isijulia()
    id_count += 1
    if !js_running
      inject_js()
      js_running = true
    end
    JSTermWidget(id, width, height, display, false)
  else
    error("JSTermWidget is only available in IJulia environments")
  end
end

comm = nothing

evthandler = Dict()
global_evthandler = nothing

function register_evthandler(f::Function, device, port)
  global evthandler
  if GR.isijulia()
    send_command(Dict("command" => "enable_events"), "cmd", string(device, port))
    evthandler[string(device, port)] = f
  else
    error("register_evthandler is only available in IJulia environments")
  end
end

function unregister_evthandler(device, port)
  global evthandler
  if GR.isijulia()
    if global_evthandler === nothing
      send_command(Dict("command" => "disable_events"), "cmd", string(device, port))
    end
    evthandler[string(device, port)] = nothing
  else
    error("unregister_evthandler is only available in IJulia environments")
  end
end

function register_evthandler(f::Function)
  global global_evthandler
  if GR.isijulia()
    send_command(Dict("command" => "enable_events"), "cmd", nothing)
    global_evthandler = f
  else
    error("register_evthandler is only available in IJulia environments")
  end
end

function unregister_evthandler()
  global global_evthandler, evthandler
  if GR.isijulia()
    send_command(Dict("command" => "disable_events"), "cmd", nothing)
    for key in keys(evthandler)
      if evthandler[key] !== nothing
        send_command(Dict("command" => "enable_events"), "cmd", key)
      end
    end
    global_evthandler = nothing
  else
    error("unregister_evthandler is only available in IJulia environments")
  end
end

function send_command(msg, msgtype, id=nothing)
  global comm
  if GR.isijulia()
    if comm === nothing
      error("JSTerm comm not initialized.")
    else
      if id !== nothing
        Main.IJulia.send_comm(comm, merge(msg, Dict("type" => msgtype, "id" => id)))
      else
        Main.IJulia.send_comm(comm, merge(msg, Dict("type" => msgtype)))
      end
    end
  else
    error("send_command is only available in IJulia environments")
  end
end

function send_evt(msg, device, port)
  if GR.isijulia()
    send_command(msg, "evt", string(device, port))
  else
    error("send_evt is only available in IJulia environments")
  end
end

function send_evt(msg, identifier)
  if GR.isijulia()
    send_command(msg, "evt", identifier)
  else
    error("send_evt is only available in IJulia environments")
  end
end

function disable_jseventhandling(device, port)
  if GR.isijulia()
    send_command(Dict("command" => "disable_jseventhandling"), "cmd", string(device, port))
  else
    error("disable_jseventhandling is only available in IJulia environments")
  end
end

function enable_jseventhandling(device, port)
  if GR.isijulia()
    send_command(Dict("command" => "enable_jseventhandling"), "cmd", string(device, port))
  else
    error("enable_jseventhandling is only available in IJulia environments")
  end
end

function disable_jseventhandling()
  if GR.isijulia()
    send_command(Dict("command" => "disable_jseventhandling"), "cmd", nothing)
  else
    error("disable_jseventhandling is only available in IJulia environments")
  end
end

function enable_jseventhandling()
  if GR.isijulia()
    send_command(Dict("command" => "enable_jseventhandling"), "cmd", nothing)
  else
    Main.IJulia.send_comm(comm, Dict("json" => f, "type" => "evt"))
    error("enable_jseventhandling is only available in IJulia environments")
  end
end

function comm_msg_callback(msg)
  data = msg.content["data"]
  if haskey(data, "type")
    if data["type"] == "save"
      display(HTML(string("
", data["content"]["data"], "
")))
    elseif data["type"] == "evt"
      global_evthandler(data["content"])
      if haskey(evthandler, data["id"]) && evthandler[data["id"]] !== nothing
        evthandler[data["id"]](data["content"])
      end
    end
  end
end

function jsterm_send(data::String)
  global js_running, draw_end_condition, comm, dispid
  if GR.isijulia()
    if comm === nothing
      comm = Main.IJulia.Comm("jsterm_comm")
      comm.on_close = function comm_close_callback(msg)
        global js_running
        js_running = false
      end
      comm.on_msg = comm_msg_callback
    end
    display(HTML(string("

Processing Input
")))
    Main.IJulia.send_comm(comm, Dict("json" => data, "type"=>"draw", "display"=>dispid))
    dispid += 1
  else
    error("jsterm_send is only available in IJulia environments")
  end
end

function recv(name::Cstring, id::Int32, msg::Cstring)
    # receives string from C and sends it to JS via Comm
    global js_running, draw_end_condition
    if !js_running
      inject_js()
      js_running = true
    end
    jsterm_send(unsafe_string(msg))
    return convert(Int32, 1)
end

function send(name::Cstring, id::Int32)
    # Dummy function, not in use
    return convert(Cstring, "String")
end

jswidgets = nothing
send_c = nothing
recv_c = nothing

function initjs()
    global jswidgets, send_c, recv_c

    jswidgets = Dict{Int32, JSTermWidget}()
    send_c = @cfunction(send, Cstring, (Cstring, Int32))
    recv_c = @cfunction(recv, Int32, (Cstring, Int32, Cstring))
    send_c, recv_c
end

end # module


text:
function maketable4(panel::DataFrame; trialrunonly::Bool=trialrunonly,
    runpanela::Bool=true,
    )::Nothing
  #load the data we will work with
  #panel::DataFrame = constructportfolios(refreshportfolios=false, trialrunonly=trialrunonly)
  sort!(panel, [:permno, :date])

  ################common data
  cqports::PortfolioConstructions = leverageconstructions()
  T4_FNAMES::Vector{Symbol} = CQUART_FNAMES_FD
  T4_FVALS::Vector{Symbol} = CQUART_FVALS_D
  T4_FGROUPS::Vector{Symbol} = CQUART_FGROUPS_FD

  #########################PANEL A#########################
  T4_PANEL_A_ROWKEYS::Vector{Symbol} = [:constant, :plusxmkt, :plussmbhml, :plusrmwcma, :plusumd]

  T4_PANEL_A_ROWLABEL::Dict = Dict(
    :constant=>"Constant",
    :plusxmkt=>"+XMKT",
    :plussmbhml=>"+SMB+HML",
    :plusrmwcma=>"+RMW+CMA",
    :plusumd=>"+UMD")

  T4_PANEL_A_CONTROLS::Dict = Dict(
    :constant=>CONTROL0,
    :plusxmkt=>[:mkt],
    :plussmbhml=>[:mkt, :P_ff3m_smb,	:P_ff3m_hml],
    :plusrmwcma=>[:mkt, :P_ff5_smb,	:P_ff5_hml,	:P_ff5_rmw,	:P_ff5_cma],
    :plusumd=>[:mkt, :P_ff5_smb,	:P_ff5_hml,	:P_ff5_rmw,	:P_ff5_cma, :P_ff3m_umd])

  T4_PANEL_A_FOCALS::Dict = Dict(
    :constant=>T4_FNAMES,
    :plusxmkt=>T4_FNAMES,
    :plussmbhml=>T4_FNAMES,
    :plusrmwcma=>T4_FNAMES,
    :plusumd=>T4_FNAMES)

  T4_PANEL_A_NAME::String = "t4panela-$(REPLICATION_TYPE[])-$(OUT_SUFFIX[])"

  T4_PANEL_A_MIN_ROWS::Int = 24

  runpanela && tablexpanelC(panel,
    fvals=T4_FVALS,
    fgroups=T4_FGROUPS,
    Fnames = T4_FNAMES,
    rowkeys=T4_PANEL_A_ROWKEYS,
    rowlabel=T4_PANEL_A_ROWLABEL,
    controls=T4_PANEL_A_CONTROLS,
    focals=T4_PANEL_A_FOCALS,
    tablename=T4_PANEL_A_NAME,
    minrowspergroup = T4_PANEL_A_MIN_ROWS
  )
  return nothing
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
julia> @elapsed bigstr = string(BigInt(5)^4^3^2)
0.017507363

julia> length(bigstr)
183231

julia> bigstr[1:20]
"62060698786608744707"

julia> bigstr[end-20:end]
"892256259918212890625"

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export QPEPSConfig, mbit, vbit, nspins, get_circuit, QMPSConfig, nbath
export Bag, setcontent!

abstract type AbstractQPEPSConfig end
struct QPEPSConfig <: AbstractQPEPSConfig
    nm::Int
    nv::Int
    nrepeat::Int
    depth::Int
end

struct QMPSConfig <: AbstractQPEPSConfig
    nv::Int
    nrepeat::Int
    depth::Int
end

nm(c::QMPSConfig) = 1
nm(c::QPEPSConfig) = c.nm
# index physical qubits
mbit(c::AbstractQPEPSConfig) = 1:nm(c)
mbit(c::AbstractQPEPSConfig, ibit::Int) = mbit(c)[ibit]

# index virtual qubits
vbit(c::AbstractQPEPSConfig, ibit::Int, iv::Int) = vbit(c, ibit)[iv]
vbit(c::AbstractQPEPSConfig, ibit::Int) = (nm(c)+c.nv*(ibit-1)+1):(nm(c)+c.nv*ibit)
Yao.nqubits(c::AbstractQPEPSConfig) = (c.nv+1)*nm(c)

# total system spins
nspins(c::QPEPSConfig) = c.nm*c.nrepeat + c.nv*c.nm
nspins(c::QMPSConfig) = c.nrepeat + c.nv-1  # one of them is an ancilla

# number of bath qubits
nbath(c::QPEPSConfig) = c.nv
nbath(c::QMPSConfig) = c.nv-1

"""get a SU(2) symmetric PEPS ansatz."""
function get_circuit(c::QPEPSConfig)
    nbit = nqubits(c)
    circ = chain(nbit)
    add!(block) = push!(circ, block)
    for i=1:c.nrepeat
        # for physical qubits, form a singlet
        chain(nbit, singlet_block(nbit, mbit(c,ibit), mbit(c,ibit+1)) for ibit=1:2:c.nm) |> add!
        if i==1
            # create singlets if it is the initial repeatition.
            for iv=1:c.nv
                chain(nbit, singlet_block(nbit, vbit(c,ibit,iv), vbit(c,ibit+1,iv)) for ibit=1:2:c.nm) |> add!
            end
        end
        for ir=1:c.depth
            for ibit=1:c.nm
                # for each physical qubit, interact with its own bath.
                chain(nbit, pswap(nbit,mbit(c,ibit),vbit(c,ibit,iv)) for iv=1:c.nv) |> add!
                # interact neghboring physical qubits
                pswap(nbit,mbit(c,ibit),mbit(c,ibit%c.nm+1)) |> add!
                # interact neghboring virtual qubits
                chain(nbit, pswap(nbit,vbit(c,ibit,iv),vbit(c,ibit%c.nm+1,iv)) for iv=1:c.nv) |> add!
            end
        end
        # measure physical qubits
        locs = (mbit(c)...,)
        Bag(basis_rotor(Z, nbit, locs)) |> add!
        Measure(nbit; locs=locs, resetto=0) |> add!
        if i==c.nrepeat
            for iv=1:c.nv
                locs = map(ibit->c.nm+c.nv*(ibit-1)+iv, (1:c.nm...,))
                Bag(basis_rotor(Z, nbit, locs)) |> add!
                Measure(nbit; locs=locs, resetto = 0) |> add!
            end
        end
    end
    return circ
end

"""get a SU(2) symmetric MPS ansatz."""
function get_circuit(c::QMPSConfig)
    @assert c.nrepeat%2 == 0
    @assert c.nv%2 == 1
    nbit = nqubits(c)
    circ = chain(nbit)
    add!(block) = push!(circ, block)
    for i=1:c.nrepeat
        # for physical qubits, form a singlet
        if i%2==1
            singlet_block(nbit, 1, nbit) |> add!
        else
            swap(nbit, 1, nbit) |> add!
        end
        if i==1
            # create singlets if it is the initial repeatition.
            chain(nbit, singlet_block(nbit, vbit(c,1,iv), vbit(c,1,iv+1)) for iv=1:2:nbath(c)) |> add!
        end
        for ir=1:c.depth
            chain(nbit, pswap(nbit,iv,iv%c.nv+1) for iv=1:c.nv) |> add!
        end
        # measure physical qubits
        locs = (1,)
        Bag(basis_rotor(Z, nbit, locs)) |> add!
        Measure(nbit; locs=locs, resetto=0) |> add!
        if i==c.nrepeat
            for iv=1:nbath(c)
                Bag(basis_rotor(Z, nbit, (1+iv,))) |> add!
                Measure(nbit; locs=(1+iv,), resetto = 0) |> add!
            end
        end
    end
    return circ
end


text:
##
using IFGF
using LinearAlgebra
using StaticArrays
using Random

include(joinpath(IFGF.PROJECT_ROOT,"test","testutils.jl"))

# parameters
p        = (3,5,5)
ppw      = 10 # point per wavelength
T        = ComplexF64
r        = 1 # sphere radius
area     = 4π*r^2
nmax     = 100
Δs       = (0.5,π/2,π/2)

# loop over wavenumbers
kvec = [2π,2π,4π,8π]
# kvec = [2π,2π,4π,8π,16π,32π,64π]

for k in kvec
    n        = ceil(Int,ppw^2*k^2*area/(4π^2))
    K        = HelmholtzKernel(k)
    Xpts = Ypts = sphere_uniform(n,r)
    I   = collect(1:1000)
    B   = ones(T,n)
    C   = zeros(T,n)
    exa = [sum(K( Xpts[i],Ypts[j])*B[j] for j in 1:n) for i in I]
    ta  = @elapsed A = assemble_ifgf(K,Xpts,Ypts;nmax,order=p.-1,Δs,threads=true)
    tp = @elapsed mul!(C,A,B,1,0;threads=true)
    er = norm(C[I]-exa,2) / norm(exa,2)
    println("="^80)
    println(A)
    println("lambda=$(2*pi/k),N=$n,checksum=$(exa[1]),tassemble = $ta, tprod=$tp,
    error=$er")
end

max_stars_repo_name: WaveProp/IFGF

max_issues_repo_name: WaveProp/IFGF

max_forks_repo_name: WaveProp/IFGF


text:
begin
    local b = startboard()
    local list = MoveList(100)
    local list2 = MoveList(100)
    local positions_tested = 0

    while positions_tested < 100000
        positions_tested += 1
        recycle!(list)
        moves(b, list)
        if list.count == 0 || squarecount(occupiedsquares(b)) < 10
            b = startboard()
            continue
        else
            if !ischeck(b)
                recycle!(list2)
                pseudocaptures(b, list2)
                for m ∈ list
                    if moveiscapture(b, m) || ispromotion(m)
                        @test m ∈ list2
                    else
                        @test m ∉ list2
                    end
                end
            end
            domove!(b, rand(list))
        end
    end
end

begin
    local b = startboard()
    local list = MoveList(100)
    local list2 = MoveList(100)
    local positions_tested = 0

    while positions_tested < 100000
        positions_tested += 1
        recycle!(list)
        moves(b, list)
        if list.count == 0 || squarecount(occupiedsquares(b)) < 10
            b = startboard()
            continue
        else
            if !ischeck(b)
                recycle!(list2)
                pseudoquiets(b, list2)
                for m ∈ list
                    if moveiscapture(b, m) || ispromotion(m)
                        @test m ∉ list2
                    else
                        @test m ∈ list2
                    end
                end
            end
            domove!(b, rand(list))
        end
    end
end

begin
    local b = startboard()
    local list = MoveList(100)
    local list2 = MoveList(100)
    local positions_tested = 0

    while positions_tested < 100000
        positions_tested += 1
        recycle!(list)
        moves(b, list)
        if list.count == 0 || squarecount(occupiedsquares(b)) < 10
            b = startboard()
            continue
        else
            if !ischeck(b)
                recycle!(list2)
                pseudochecks(b, list2)
                for m ∈ list
                    if moveiscapture(b, m) || ispromotion(m) || moveiscastle(b, m)
                        @test m ∉ list2
                    elseif moveischeck(b, m)
                        @test m ∈ list2
                    else
                        @test m ∉ list2
                    end
                end
            end
            domove!(b, rand(list))
        end
    end
end

begin
    local b = startboard()
    local list = MoveList(100)
    local list2 = MoveList(100)
    local positions_tested = 0

    while positions_tested < 100000
        positions_tested += 1
        recycle!(list)
        moves(b, list)
        if list.count == 0 || squarecount(occupiedsquares(b)) < 10
            b = startboard()
            continue
        else
            if !ischeck(b)
                recycle!(list2)
                pseudomoves(b, list2)
                for m ∈ list2
                    @test moveispseudo(b, m)
                    @test pseudoislegal(b, m) || m ∉ list
                end
            end
            domove!(b, rand(list))
        end
    end
end


text:
#
# Version 1.0 (SDM-16)
#
#     30.10.2016: - Tested with Octave/Matlab version.
#                   Outputs are the same (very similar at least)
#                   for the artificially generated and Human Activity signals.
#
#                 - Added function scale01() to scale PCCF values to
#                   the range [0, 1].
#
# include("./probs.jl")
module Pccf
# SAC 2016 - Improved Matrix based version
using Distributions
export pccf, write_pccf, scale01

normalize = vec -> vec / sum(vec)

function pccf(T, mu, sigma)
    # exact analytical solution
    t = 1:T
    chp_probs = zeros(T)
    for dist_i in 1:T
        dist = Normal(dist_i * mu, sqrt(dist_i*sigma^2))
        chp_probs += pdf(dist, t)
    end
    # normalize(chp_probs)
    chp_probs
end

function weights_matrix(T)
  n = T
  M = zeros(n,n)
  M[1,:] = 1:n
  for i = 2:n
      for j = i:n
          M[i,j] = M[i-1, j-1]
      end
  end
  M'
end

function gaussPdf(x, mu, sigma)
    k = 1/(sigma * sqrt(2*pi))
    p = k * exp(-((x-mu)^2/(2*sigma^2)))
end

function pccf_mp(T, mu, sigma)
    # message passing version
    gaussDist = Normal(mu, sigma)
    P         = zeros(T, T)
    P[:,1]    = pdf(gaussDist, 1:T)
    W         = pdf(gaussDist, weights_matrix(T))
    for i = 1:T-1
        P[i+1:T, i+1] = W[1:T-i, 1:T-i] * P[i:end-1, i]
    end
    #normalize( sum(P,2) ) # change probabilities
    sum(P,2) # change probabilities
end

function pccfWithConfirmation(probs,delayInMu,changes)
    # PCCF with confirmation
    T=length(probs)
    confirmed=changes[1:delayInMu:length(changes)]
    out=Float64[]
    nchps=length(confirmed)
    for i=1:nchps
        if i==1
            out=vcat(out,probs[1:(confirmed[i]-1)])
        else
            out=vcat(out,probs[1:(confirmed[i]-confirmed[i-1])])
        end
        if i==nchps
            out=vcat(out,probs[1:(T-confirmed[i]+1)])
        end
    end
    @assert length(out)==T
    out
end

function write_pccf(filepath, T = 100, mu = 25.0, sigma = 3.2)
    # r = pccf(T, mu, sigma)
    r = pccf_exact(T, mu, sigma)
    f = open(filepath, "w")
    writedlm(f, r)
    close(f)
    println("... PCCF function is written into $filepath")
end

"Scale values of the vector to the range [0, 1]"
function scale01(v)
    minv = minimum(v)
    maxv = maximum(v)
    (v - minv) / (maxv - minv)
end

end
#------------------------- End module Pccf -----------------------------------#

# module PccfIcdm
# # using Probs
# export pccf, write_pccf
# # Old version based on the init window of size w
#
# function gaussPdf(x, mu, sigma)
#     k = 1/(sigma * sqrt(2*pi))
#     p = k * exp(-((x-mu)^2/(2*sigma^2)))
# end
#
# function initProbs(w, mu, sigma)
# 	f(x) = gaussPdf(x, mu, sigma)
# 	init = map(f, collect(1:w))
# 	init / sum(init)
# end
#
# function possiblePrevPositions(prev_len, delta, i)
# 	if i >= 2 && i <= prev_len + delta
# 		collect(max(1, i-delta) : 1 : min(i-1, prev_len))
# 	else
# 		nothing
# 	end
# end
#
# function forwardProbs(prevProbs, w, mu, sigma)
# 	delta = w
# 	n_prev = length( prevProbs )
# 	nextProbs = zeros(n_prev + w - 1)
# 	n_next = length(nextProbs)
# 	for i = 2 : (n_next+1)
# 		prevPositions = possiblePrevPositions( n_prev, delta, i)
# 		if prevPositions != nothing
# 			prevProbs1 = prevProbs[prevPositions]
# 			deltas = map(x -> i - x, prevPositions)
# 			prob_deltas = map(x -> gaussPdf(x, mu, sigma), deltas)
# 			nextProbs[i-1] = sum(prob_deltas .* prevProbs1)
# 		end
# 	end
# 	nextProbs = nextProbs / sum(nextProbs)
# end
#
# function pccf(num_events, w, mu, sigma)
# 	tic()
# 	allM = zeros(num_events * w, num_events)
# 	allM[1:w, 1] = initProbs(w, mu, sigma)
# 	for i = 2:num_events
# 		prev_vect = allM[(i-1):((i-1)*w), i-1]
# 		next_vec_of_probs = forwardProbs(prev_vect, w, mu, sigma)
# 		allM[(i:(i*w)), i] = next_vec_of_probs
# 	end
# 	nn, mm = size(allM)
# 	cc = zeros(nn)
# 	for i = 1:nn
# 		s = 0.
# 		for j = 1:mm
# 			s += allM[i,j]
# 		end
# 		cc[i] = s
# 	end
# 	toc()
# 	cc = cc / sum(cc)
# end
#
# function write_pccf(filepath; num_of_events = 30, w = 150, mu = 25.0, sigma = 3.2)
# 	r = pccf(num_of_events, w, mu, sigma)
# 	f = open(filepath, "w")
# 	writedlm(f, r)
# 	close(f)
# 	println("... PCCF function is written into $filepath")
# end
# end
# #----------------------- END MODULE PccfIcdm ---------------------------------#

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
using StaticArrays
using Base: @_inline_meta
export reconstruct, DelayEmbedding, embed, τrange
export AbstractEmbedding

#####################################################################################
# Univariate Delay Coordinates
#####################################################################################
"Super-type of embedding methods."
abstract type AbstractEmbedding end

"""
    DelayEmbedding(γ, τ, h = nothing) → `embedding`
Return a delay coordinates embedding structure to be used as a function-like-object,
given a timeseries and some index. Calling
```julia
embedding(s, n)
```
will create the `n`-th delay vector of the embedded space, which has `γ`
temporal neighbors with delay(s) `τ`.
`γ` is the embedding dimension minus 1, `τ` is the delay time(s) while `h` are
extra weights, as in [`embed`](@ref) for more.

**Be very careful when choosing `n`, because `@inbounds` is used internally.**
Use [`τrange`](@ref)!
"""
struct DelayEmbedding{γ, H} <: AbstractEmbedding
    delays::SVector{γ, Int}
    h::SVector{γ, H}
end
@inline DelayEmbedding(γ, τ, h = nothing) = DelayEmbedding(Val{γ}(), τ, h)
@inline function DelayEmbedding(g::Val{γ}, τ::Int, h::H) where {γ, H}
    idxs = [k*τ for k in 1:γ]
    hs = hweights(g, h)
    htype = H <: Union{Nothing, Real} ? H : eltype(H)
    return DelayEmbedding{γ, htype}(SVector{γ, Int}(idxs...), hs)
end
@inline function DelayEmbedding(g::Val{γ}, τ::AbstractVector, h::H) where {γ, H}
    γ != length(τ) && throw(ArgumentError(
        "Delay time vector length must equal the embedding dimension minus 1."
    ))
    hs = hweights(g, h)
    htype = H <: Union{Nothing, Real} ? H : eltype(H)
    return DelayEmbedding{γ, htype}(SVector{γ, Int}(τ...), hs)
end
hweights(::Val{γ}, h::Nothing) where {γ} = SVector{γ, Nothing}(fill(nothing, γ))
hweights(::Val{γ}, h::T) where {γ, T<:Real} = SVector{γ, T}([h^b for b in 1:γ]...)
hweights(::Val{γ}, h::AbstractVector) where {γ} = SVector{γ}(h)

@generated function (r::DelayEmbedding{γ, Nothing})(s::AbstractArray{T}, i) where {γ, T}
    gens = [:(s[i + r.delays[$k]]) for k=1:γ]
    quote
        @_inline_meta
        @inbounds return SVector{$γ+1,T}(s[i], $(gens...))
    end
end

# This is the version with weights
@generated function (r::DelayEmbedding{γ})(s::AbstractArray{T}, i) where {γ, T, R<:Real}
    gens = [:(r.h[$k]*(s[i + r.delays[$k]])) for k=1:γ]
    quote
        @_inline_meta
        @inbounds return SVector{$γ+1,T}(s[i], $(gens...))
    end
end

"""
    embed(s, d, τ [, h])
Embed `s` using delay coordinates with embedding dimension `d` and delay time `τ`
and return the result as a [`Dataset`](@ref). Optionally use weight `h`, see below.

Here `τ > 0`, use [`genembed`](@ref) for a generalized version.

## Description
If `τ` is an integer, then the ``n``-th entry of the embedded space is
```math
(s(n), s(n+\\tau), s(n+2\\tau), \\dots, s(n+(d-1)\\tau))
```
If instead `τ` is a vector of integers, so that `length(τ) == d-1`,
then the ``n``-th entry is
```math
(s(n), s(n+\\tau[1]), s(n+\\tau[2]), \\dots, s(n+\\tau[d-1]))
```

The resulting set can have same
invariant quantities (like e.g. lyapunov exponents) with the original system
that the timeseries were recorded from, for proper `d` and `τ`.
This is known as the Takens embedding theorem [^Takens1981] [^Sauer1991].
The case of different delay times allows embedding systems with many time scales,
see[^Judd1998].

If provided, `h` can be weights to multiply the entries of the embedded space.
If `h isa Real` then the embedding is
```math
(s(n), h \\cdot s(n+\\tau), h^2 \\cdot s(n+2\\tau), \\dots,h^{d-1} \\cdot s(n+γ\\tau))
```
Otherwise `h` can be a vector of length `d-1`, which the decides the weights of each
entry directly.

## References
[^Takens1981] : F. Takens, *Detecting Strange Attractors in Turbulence — Dynamical
Systems and Turbulence*, Lecture Notes in Mathematics **366**, Springer (1981)

[^Sauer1991] : T. Sauer *et al.*, J. Stat. Phys. **65**, pp 579 (1991)

[^Judd1998]: K. Judd & A. Mees, [Physica D **120**, pp 273 (1998)](https://www.sciencedirect.com/science/article/pii/S0167278997001188)

[^Farmer1988]: Farmer & Sidorowich, [Exploiting Chaos to Predict the Future and Reduce Noise"](http://www.nzdl.org/gsdlmod?e=d-00000-00---off-0cltbibZz-e--00-1----0-10-0---0---0direct-10---4-------0-1l--11-en-50---20-home---00-3-1-00-0--4--0--0-0-11-10-0utfZz-8-00&a=d&cl=CL3.16&d=HASH013b29ffe107dba1e52f1a0c_1245)

"""
function embed(s::AbstractVector{T}, d, τ, h::H = nothing) where {T, H}
    if d == 1
        return Dataset(s)
    end
    htype = H <: Union{Nothing, Real} ? H : eltype(H)
    de::DelayEmbedding{d-1, htype} = DelayEmbedding(d-1, τ, h)
    return embed(s, de)
end

@inline function embed(s::AbstractVector{T}, de::DelayEmbedding{γ}) where {T, γ}
    r = τrange(s, de)
    data = Vector{SVector{γ+1, T}}(undef, length(r))
    @inbounds for i in r
        data[i] = de(s, i)
    end
    return Dataset{γ+1, T}(data)
end

"""
    τrange(s, de::AbstractEmbedding)
Return the range `r` of valid indices `n` to create delay vectors out of `s` using `de`.
"""
τrange(s, de::AbstractEmbedding) = 1:(length(s) - maximum(de.delays))

#####################################################################################
# Generalized embedding (arbitrary combination of timeseries and delays)
#####################################################################################
export GeneralizedEmbedding, genembed

"""
    GeneralizedEmbedding(τs, js = ones(length(τs)), ws = nothing) -> `embedding`
Return a delay coordinates embedding structure to be used as a function.
Given a timeseries *or* trajectory (i.e. `Dataset`) `s` and calling
```julia
embedding(s, n)
```
will create the delay vector of the `n`-th point of `s` in the embedded space using
generalized embedding (see [`genembed`](@ref)).

`js` is ignored for timeseries input `s` (since all entries of `js` must be `1` in
this case) and in addition `js` defaults to `(1, ..., 1)` for all `τ`.

**Be very careful when choosing `n`, because `@inbounds` is used internally.**
Use [`τrange`](@ref)!
"""
struct GeneralizedEmbedding{D, W} <: AbstractEmbedding
    τs::NTuple{D, Int}
    js::NTuple{D, Int}
    ws::NTuple{D, W}
end
function GeneralizedEmbedding(τs::NTuple{D, Int}) where {D} # type stable version
    GeneralizedEmbedding{D, Nothing}(
        τs, NTuple{D, Int}(ones(D)), NTuple{D, Nothing}(fill(nothing, D))
    )
end

function GeneralizedEmbedding(τs, js = ones(length(τs)), ws = nothing)
    D = length(τs)
    a = NTuple{D, Int}(τs)
    b = NTuple{D, Int}(js)
    W = isnothing(ws) ? Nothing : typeof(promote(ws...)[1])
    c = isnothing(ws) ? NTuple{D, W}(fill(nothing, D)) : NTuple{D, W}(promote(ws...))
    return GeneralizedEmbedding{D, W}(a, b, c)
end

function Base.show(io::IO, g::GeneralizedEmbedding{D, W}) where {D, W}
    print(io, "$D-dimensional generalized embedding\n")
    print(io, "  τs: $(g.τs)\n")
    print(io, "  js: $(g.js)\n")
    wp = W == Nothing ? "nothing" : g.ws
    print(io, "  ws: $(wp)")
end

const Data{T} = Union{Dataset{D, T}, AbstractVector{T}} where {D}

# dataset version
@generated function (g::GeneralizedEmbedding{D, W})(s::AbstractDataset{Z, T}, i::Int) where {D,W,Z,T}
    gens = if W == Nothing
        [:(s[i + g.τs[$k], g.js[$k]]) for k=1:D]
    else
        [:(g.ws[$k]*s[i + g.τs[$k], g.js[$k]]) for k=1:D]
    end
    X = W == Nothing ? T : promote_type(T, W)
    quote
        @_inline_meta
        @inbounds return SVector{$D,$X}($(gens...))
    end
end

# vector version
@generated function (g::GeneralizedEmbedding{D, W})(s::AbstractVector{T}, i::Int) where {D,W,T}
    gens = if W == Nothing
        [:(s[i + g.τs[$k]]) for k=1:D]
    else
        [:(g.ws[$k]*s[i + g.τs[$k]]) for k=1:D]
    end
    X = W == Nothing ? T : promote_type(T, W)
    quote
        @_inline_meta
        @inbounds return SVector{$D,$X}($(gens...))
    end
end

τrange(s, ge::GeneralizedEmbedding) =
max(1, (-minimum(ge.τs) + 1)):min(length(s), length(s) - maximum(ge.τs))

"""
    genembed(s, τs, js = ones(...); ws = nothing) → dataset
Create a generalized embedding of `s` which can be a timeseries or arbitrary `Dataset`,
and return the result as a new `Dataset`.

The generalized embedding works as follows:
- `τs` denotes what delay times will be used for each of the entries
  of the delay vector. It is recommended that `τs[1] = 0`.
  `τs` is allowed to have *negative entries* as well.
- `js` denotes which of the timeseries contained in `s`
  will be used for the entries of the delay vector. `js` can contain duplicate indices.
- `ws` are optional weights that weight each embedded entry (the i-th entry of the
    delay vector is weighted by `ws[i]`). If provided, it is recommended that `ws[1] == 1`.

`τs, js, ws` are tuples (or vectors) of length `D`, which also coincides with the embedding
dimension. For example, imagine input trajectory ``s = [x, y, z]`` where ``x, y, z`` are
timeseries (the columns of the `Dataset`).
If `js = (1, 3, 2)` and `τs = (0, 2, -7)` the created delay vector at
each step ``n`` will be
```math
(x(n), z(n+2), y(n-7))
```
Using `ws = (1, 0.5, 0.25)` as well would create
```math
(x(n), \\frac{1}{2} z(n+2), \\frac{1}{4} y(n-7))
```

`js` can be skipped, defaulting to index 1 (first timeseries) for all delay entries, while
it has no effect if `s` is a timeseries instead of a `Dataset`.

See also [`embed`](@ref). Internally uses [`GeneralizedEmbedding`](@ref).
"""
function genembed(s, τs, js = ones(length(τs)); ws = nothing)
    D = length(τs)
    W = isnothing(ws) ? Nothing : typeof(promote(ws...)[1])
    ge::GeneralizedEmbedding{D, W} = GeneralizedEmbedding(τs, js, ws)
    return genembed(s, ge)
end

function genembed(s, ge::GeneralizedEmbedding{D, W}) where {D, W}
    r = τrange(s, ge)
    T = eltype(s)
    X = W == Nothing ? T : promote_type(T, W)
    data = Vector{SVector{D, X}}(undef, length(r))
    @inbounds for (i, n) in enumerate(r)
        data[i] = ge(s, n)
    end
    return Dataset{D, X}(data)
end