text:
abstract AbstractVegaType

function Base.call{T<:AbstractVegaType}(::Type{T}; kwargs...)
    out = T()
    for (sym, val) in kwargs
        setfield!(out, sym, val)
    end
    out
end

# TODO: Fill me in
type VegaTransform
	x::Dict{Any, Any}
end

tojs(x::VegaTransform) = x.x

Base.copy(x::VegaTransform) = VegaTransform()

type VegaStream <: AbstractVegaType
    _type::Union{AbstractString,Void}
    expr::Union{AbstractString,Void}
    scale::Union{AbstractString,Void}
end

type VegaDataRef <: AbstractVegaType
    data::Union{AbstractString,Void}
    field::Union{AbstractString,AbstractVector,Dict{Any,Any},Void}
    fields::Union{AbstractString,AbstractVector,Dict{Any,Any},Void}
    sort::Union{Bool,Dict{Any,Any},Void}
end

type VegaValueRef <: AbstractVegaType
    value::Any
    field::Union{AbstractString,Dict{Any,Any},Void}
    group::Union{AbstractString,Bool,Void}
    scale::Union{AbstractString,Dict{Any,Any},Void}
    mult::Union{Number,Void}
    offset::Union{Number,Void}
    band::Union{Bool,Void}
    signal::Union{AbstractString,Void}
end

type VegaPadding <: AbstractVegaType
    top::Union{Number,Void}
    left::Union{Number,Void}
    bottom::Union{Number,Void}
    right::Union{Number,Void}
end

type VegaMarkPropertySet <: AbstractVegaType
    x::Union{VegaValueRef,Void}
    x2::Union{VegaValueRef,Void}
    xc::Union{VegaValueRef,Void}
    width::Union{VegaValueRef,Void}
    y::Union{VegaValueRef,Void}
    y2::Union{VegaValueRef,Void}
    yc::Union{VegaValueRef,Void}
    height::Union{VegaValueRef,Void}
    opacity::Union{VegaValueRef,Void}
    fill::Union{VegaValueRef,Void}
    fillOpacity::Union{VegaValueRef,Void}
    stroke::Union{VegaValueRef,Void}
    strokeWidth::Union{VegaValueRef,Void}
    strokeOpacity::Union{VegaValueRef,Void}
    strokeDash::Union{VegaValueRef,Void}
    strokeDashOffset::Union{VegaValueRef,Void}
    cursor::Union{VegaValueRef,Void}
    clip::Union{VegaValueRef,Void}
    size::Union{VegaValueRef,Void}
    shape::Union{VegaValueRef,Void}
    path::Union{VegaValueRef,Void}
    innerRadius::Union{VegaValueRef,Void}
    outerRadius::Union{VegaValueRef,Void}
    startAngle::Union{VegaValueRef,Void}
    endAngle::Union{VegaValueRef,Void}
    interpolate::Union{VegaValueRef,Void}
    tension::Union{VegaValueRef,Void}
    orient::Union{VegaValueRef,Void}
    url::Union{VegaValueRef,Void}
    align::Union{VegaValueRef,Void}
    baseline::Union{VegaValueRef,Void}
    text::Union{VegaValueRef,Void}
    dx::Union{VegaValueRef,Void}
    dy::Union{VegaValueRef,Void}
    radius::Union{VegaValueRef,Void}
    theta::Union{VegaValueRef,Void}
    angle::Union{VegaValueRef,Void}
    font::Union{VegaValueRef,Void}
    fontSize::Union{VegaValueRef,Void}
    fontWeight::Union{VegaValueRef,Void}
    fontStyle::Union{VegaValueRef,Void}
end

type VegaAxis <: AbstractVegaType
    _type::Union{AbstractString,Void}
    scale::Union{AbstractString,Void}
    orient::Union{AbstractString,Void}
    title::Union{AbstractString,Void}
    titleOffset::Union{Number,Void}
    format::Union{AbstractString,Void}
    formatType::Union{AbstractString,Void}
    ticks::Union{Number,Void}
    values::Union{AbstractVector,Void}
    subdivide::Union{Number,Void}
    tickPadding::Union{Number,Void}
    tickSize::Union{Number,Void}
    tickSizeMajor::Union{Number,Void}
    tickSizeMinor::Union{Number,Void}
    tickSizeEnd::Union{Number,Void}
    offset::Union{Dict{Any,Any},Number,Void}
    layer::Union{AbstractString,Void}
    grid::Union{Bool,Void}
    properties::Union{Dict{Any,Any},VegaMarkPropertySet,Void}
end

type VegaFormat <: AbstractVegaType
    _type::Union{AbstractString,Void}
    parse::Union{Dict{Any,Any},Void}
    property::Union{AbstractString,Void}
    feature::Union{AbstractString,Void}
    mesh::Union{AbstractString,Void}
    children::Union{Dict{Any,Any},Void}
end

type VegaData <: AbstractVegaType
    name::Union{AbstractString,Void}
    format::Union{VegaFormat,Void}
    values::Union{AbstractVector,Void}
    source::Union{AbstractString,Void}
    url::Union{AbstractString,Void}
    transform::Union{Array{VegaTransform,1},Void}
    modify::Union{Array{VegaStream,1},Void}
end

type VegaMarkProperties <: AbstractVegaType
    enter::Union{VegaMarkPropertySet,Void}
    exit::Union{VegaMarkPropertySet,Void}
    update::Union{VegaMarkPropertySet,Void}
    hover::Union{VegaMarkPropertySet,Void}
end

type VegaMarkFrom <: AbstractVegaType
    data::Union{AbstractString,Void}
    mark::Union{AbstractString,Void}
    transform::Union{Array{VegaTransform,1},Void}
    value::Union{AbstractString,Void}
end

type VegaSignal <: AbstractVegaType
    name::Union{AbstractString,Void}
    init::Any
    verbose::Union{Bool,Void}
    expr::Union{AbstractString,Void}
    scale::Any
    streams::Any
end

type VegaScale <: AbstractVegaType
    name::Union{AbstractString,Void}
    _type::Union{AbstractString,Void}
    domain::Union{AbstractVector,VegaDataRef,Void}
    domainMin::Union{Number,VegaDataRef,VegaSignal,Void}
    domainMax::Union{Number,VegaDataRef,VegaSignal,Void}
    range::Union{AbstractString,AbstractVector,VegaDataRef,Void}
    rangeMin::Union{Number,VegaDataRef,VegaSignal,Void}
    rangeMax::Union{Number,VegaDataRef,VegaSignal,Void}
    reverse::Union{Bool,VegaDataRef,Void}
    round::Union{Bool,Void}
    points::Union{Bool,Void}
    padding::Union{Number,Void}
    outerPadding::Union{Number,Void}
    bandWidth::Union{Number,VegaSignal,Void}
    sort::Union{Bool,Dict{Any,Any},Void}
    clamp::Union{Bool,Void}
    nice::Union{AbstractString,Bool,Void}
    exponent::Union{Number,Void}
    zero::Union{Bool,Void}
end

type VegaMark <: AbstractVegaType
    name::Union{AbstractString,Void}
    key::Union{AbstractString,Void}
    _type::Union{AbstractString,Void}
    from::Union{VegaMarkFrom,Void}
    properties::Union{VegaMarkProperties,Void}
    delay::Union{VegaValueRef,Void}
    ease::Union{AbstractString,Void}
    marks::Union{AbstractVector,Void}
    scales::Union{Array{VegaScale,1},Void}
    axes::Union{Array{VegaAxis,1},Void}
    interactive::Union{Bool,Void}
end

type VegaLegend <: AbstractVegaType
    size::Union{AbstractString,Void}
    shape::Union{AbstractString,Void}
    fill::Union{AbstractString,Void}
    stroke::Union{AbstractString,Void}
    orient::Union{AbstractString,Void}
    offset::Union{Number,Void}
    title::Union{AbstractString,Void}
    format::Union{AbstractString,Void}
    formatType::Union{AbstractString,Void}
    values::Union{AbstractVector,Void}
    properties::Union{Dict{Any,Any},VegaMarkPropertySet,Void}
end

type VegaScene <: AbstractVegaType
    fill::Any
    fillOpacity::Any
    stroke::Any
    strokeOpacity::Any
    strokeWidth::Any
    strokeDash::Any
    strokeDashOffset::Any
end

type VegaContainer <: AbstractVegaType
    scene::Union{Array{VegaScene,1},Void}
    scales::Union{Array{VegaScale,1},Void}
    axes::Union{Array{VegaAxis,1},Void}
    legends::Union{Array{VegaLegend,1},Void}
    marks::Union{Array{VegaMark,1},Void}
end

type VegaPredicate <: AbstractVegaType
    name::Union{AbstractString,Void}
    _type::Union{AbstractString,Void}
    operands::Union{AbstractVector,Void}
    range::Union{AbstractVector,Void}
    scale::Union{AbstractVector,Void}
    item::Any
end

type VegaVisualization <: AbstractVegaType
    name::Union{AbstractString,Void}
    width::Union{Int64,Void}
    height::Union{Int64,Void}
    padding::Union{AbstractString,Number,VegaPadding,Void}
    viewport::Union{Array{Int64,1},Void}
    background::Union{AbstractString,Void}
    data::Union{Array{VegaData,1},Void}
    scales::Union{Array{VegaScale,1},Void}
    axes::Union{Array{VegaAxis,1},Void}
    marks::Union{Array{VegaMark,1},Void}
    legends::Union{Array{VegaLegend,1},Void}
    signals::Union{Array{Dict{Any,Any},1},Void}
    predicates::Union{Array{Dict{Any,Any},1},Void}
    modify::Union{Array{Dict{Any,Any},1},Void}
end

max_stars_repo_name: jpthor/vega


text:
# ==============================================================================
# MARK: - Packages
# ==============================================================================
# The required packages to run the backend.

import JSON # Package requirement: add JSON
import Dates # Already included in standard library
using Base64 # Already included in standard library
using WAV # Package requirement: add WAV
using Dash # Package requirement: add https://github.com/plotly/Dash.jl.git
using Plotly # Package requirement: add https://github.com/plotly/Plotly.jl

# ==============================================================================
# MARK: - Global Outlets
# ==============================================================================
# The current values of the UI displayed on the screen.

plotlyBackend_viewDropDown = nothing # Holds the string value of the view drop down.
plotlyBackend_projectionChecklist = nothing # Holds an array of the checks that are selected.
plotlyBackend_lineWidthSlider = nothing # Holds a number that the slider is currently on.
plotlyBackend_XSlider = nothing # Holds an array of size 2 representing the [min, max] that the slider is currently on.
plotlyBackend_YSlider = nothing # Holds an array of size 2 representing the [min, max] that the slider is currently on.
plotlyBackend_ZSlider = nothing # Holds an array of size 2 representing the [min, max] that the slider is currently on.
plotlyBackend_XSliderPushable = nothing # Holds [1] if checked
plotlyBackend_YSliderPushable = nothing # Holds [1] if checked
plotlyBackend_ZSliderPushable = nothing # Holds [1] if checked
plotlyBackend_fileUploader = Dict("fileName" => nothing, "dateModified" => nothing) # Holds a dictionary with the keys: fileName and dateModified

# ==============================================================================
# MARK: - Global Variables
# ==============================================================================
# Variables saved in memory so calculations do not need to calculated everytime.

plotlyBackend_dict = nothing # A dictionary of the raw data.
plotlyBackend_data = nothing # An array of Plotly line objects representing the data of the graph.
plotlyBackend_layout = nothing # A dictionary representing the layout of the graph.
plotlyBackend_freqUnits = nothing # A string value representing the freqency units.
plotlyBackend_XExtremities = nothing # An array of size 2 representing the [min, max] x extremities of the data.
plotlyBackend_YExtremities = nothing # An array of size 2 representing the [min, max] y extremities of the data.
plotlyBackend_ZExtremities = nothing # An array of size 2 representing the [min, max] z extremities of the data.
plotlyBackend_ZExtremitiesWithXMinProjection = nothing # An array of size 2 representing the [min, max] z extremities while showing FFT projection on X MIN.
plotlyBackend_ZExtremitiesWithYMinProjection = nothing # An array of size 2 representing the [min, max] z extremities while showing real projection on Y MIN.
plotlyBackend_soundData = nothing # A dictionary containing the key: "y" and "samplingFq" to create WAV files.

# ==============================================================================
# MARK: - Global Constants
# ==============================================================================
# Constants that will remain the same during the lifecycle of the app.

const plotlyBackend_3DGraphTitle = "ISA Plot" # The graph title when 3D mode is toggled.
const plotlyBackend_XYGraphTitle = "Orthographic XY" # The graph title when XY mode is toggled.
const plotlyBackend_XZGraphTitle = "Orthographic XZ" # The graph title when XZ mode is toggled.
const plotlyBackend_YZGraphTitle = "Orthographic YZ" # The graph title when YZ mode is toggled.

# ==============================================================================
# MARK: - Public Functions
# ==============================================================================
# The public functions that outside files should access. The rest of the file
# contains private functions that outside files should not access.

# A function used to run a dash server with default data.
function runPlotlyBackend()
    defaultData = Dict("freqUnits" => "rad/s", "components" => [Dict("x" => [1,2,3], "y" => [4,5,6], "z" => [7,8,9], "colorMap" => ["rgb(255,255,0)", "rgb(0,255,0)", "rgb(255,255,0)"])])
    runPlotlyBackend(defaultData)
end

# A function used to run a dash server with specified data.
function runPlotlyBackend(dict)
    println("Loading Dash Server...")
    plotlyBackend_assignAllGlobalVariablesSafely(dict)
    plotlyBackend_hostDashServer()
end

# ==============================================================================
# MARK: - Dash Server
# ==============================================================================
# Deploys Dash server and matches interactions with the UI to the proper function.

function plotlyBackend_hostDashServer()

    app = dash("ISA app")
    app.layout = plotlyBackend_createHTML()

    callback!(app, callid"viewDropDown.value, lineWidthSlider.value, projectionChecklist.value, xSlider.value, ySlider.value, zSlider.value, uploader.contents, uploader.filename, uploader.last_modified => isaGraph.figure") do viewType, lineWidthSlider, projectionChecklist, xSlider, ySlider, zSlider, fileContents, fileName, fileDateModified
        if (plotlyBackend_viewDropDown != viewType) plotlyBackend_viewDropDownInteracted(viewType) end
        if (plotlyBackend_lineWidthSlider != lineWidthSlider) plotlyBackend_lineWidthSliderInteracted(lineWidthSlider) end
        if (plotlyBackend_projectionChecklist != projectionChecklist) plotlyBackend_projectionCheckListInteracted(projectionChecklist) end
        if (plotlyBackend_XSlider != xSlider) plotlyBackend_xSliderInteracted(xSlider) end
        if (plotlyBackend_YSlider != ySlider) plotlyBackend_ySliderInteracted(ySlider) end
        if (plotlyBackend_ZSlider != zSlider) plotlyBackend_zSliderInteracted(zSlider) end
        if (fileContents != nothing && (plotlyBackend_fileUploader["fileName"] != fileName || plotlyBackend_fileUploader["dateModified"] != fileDateModified)) plotlyBackend_uploaderInteracted(fileName, fileDateModified, fileContents) end
        return (data = plotlyBackend_data, layout = plotlyBackend_layout)
    end

    callback!(app, callid"zSliderPushable.value, projectionChecklist.value, uploader.contents, uploader.filename, uploader.last_modified => zSliderDiv.children") do zSliderPushable, projectionChecklist, fileContents, fileName, fileDateModified
        zPushData = nothing
        if (plotlyBackend_ZSliderPushable != zSliderPushable)
            plotlyBackend_ZSliderPushableInteracted(zSliderPushable)
            zPushData = (plotlyBackend_ZSliderPushable == [1])
        end
        if plotlyBackend_projectionChecklist != projectionChecklist plotlyBackend_projectionCheckListInteracted(projectionChecklist) end
        if (fileContents != nothing && (plotlyBackend_fileUploader["fileName"] != fileName || plotlyBackend_fileUploader["dateModified"] != fileDateModified)) plotlyBackend_uploaderInteracted(fileName, fileDateModified, fileContents) end
        return plotlyBackend_determineZSliderDashComponent(zPushData)
    end

    callback!(app, callid"downloadSoundButton.n_clicks => downloadSoundButtonOutletDiv.children") do numberOfClicks
        if (numberOfClicks != -1) plotlyBackend_downloadSoundButtonInteracted() end
    end

    callback!(app, callid"saveButton.n_clicks => downloadStatusDiv.children") do numberOfClicks
        if numberOfClicks != -1 return plotlyBackend_saveButtonInteracted() end
    end

    callback!(app, callid"xSliderPushable.value, ySliderPushable.value, uploader.contents, uploader.filename, uploader.last_modified => uploadFileDiv.children, xSliderDiv.children, ySliderDiv.children") do xSliderPushable, ySliderPushable, fileContents, fileName, fileDateModified
        if (plotlyBackend_XSliderPushable != xSliderPushable) plotlyBackend_XSliderPushableInteracted(xSliderPushable) end
        if (plotlyBackend_YSliderPushable != ySliderPushable) plotlyBackend_YSliderPushableInteracted(ySliderPushable) end
        xPushData = (plotlyBackend_XSliderPushable == [1])
        yPushData = (plotlyBackend_YSliderPushable == [1])
        if (fileContents != nothing && (plotlyBackend_fileUploader["fileName"] != fileName || plotlyBackend_fileUploader["dateModified"] != fileDateModified))
            plotlyBackend_uploaderInteracted(fileName, fileDateModified, fileContents)
            sleep(1) # Purposely make it 1 second slower to always show loading icon
            xPushData = nothing
            yPushData = nothing
        end
        return ["", plotlyBackend_determineXSliderDashComponent(xPushData), plotlyBackend_determineYSliderDashComponent(yPushData)]
    end

    println("Deployed. CMD-click this link: http://127.0.0.1:8080")
    run_server(app, "0.0.0.0", 8080)

end

# ==============================================================================
# MARK: - Interactions
# ==============================================================================
# These functions trigger when the corresponding UI element is interacted with.
# In general, these functions update the global outlets with the newest set of
# values, then updates anything else that needs updating.

function plotlyBackend_xSliderInteracted(xSlider)
    global plotlyBackend_XSlider = xSlider
    plotlyBackend_updateXMinProjectionLocation()
    plotlyBackend_updateLayout()
end

function plotlyBackend_ySliderInteracted(ySlider)
    global plotlyBackend_YSlider = ySlider
    plotlyBackend_updateYMinProjectionLocation()
    plotlyBackend_updateLayout()
end

function plotlyBackend_zSliderInteracted(zSlider)
    global plotlyBackend_ZSlider = zSlider
    plotlyBackend_updateZMinProjectionLocation()
    plotlyBackend_updateLayout()
end

function plotlyBackend_XSliderPushableInteracted(xSliderPushable)
    global plotlyBackend_XSliderPushable = xSliderPushable
end

function plotlyBackend_YSliderPushableInteracted(ySliderPushable)
    global plotlyBackend_YSliderPushable = ySliderPushable
end

function plotlyBackend_ZSliderPushableInteracted(zSliderPushable)
    global plotlyBackend_ZSliderPushable = zSliderPushable
end

function plotlyBackend_lineWidthSliderInteracted(lineWidthSlider)
    global plotlyBackend_lineWidthSlider = lineWidthSlider
    plotlyBackend_updateGlobalDataLineWidth()
end

function plotlyBackend_viewDropDownInteracted(viewType)
    global plotlyBackend_viewDropDown = viewType
    plotlyBackend_updateLayout()
end

function plotlyBackend_projectionCheckListInteracted(projectionChecklist)
    global plotlyBackend_projectionChecklist = projectionChecklist
    plotlyBackend_updateXMinProjectionVisibility() # FFT Projection
    plotlyBackend_updateYMinProjectionVisibility() # Real Projection
    plotlyBackend_updateZMinProjectionVisibility() # STFT Projection
    plotlyBackend_updateXMinProjectionLocation()
    plotlyBackend_updateYMinProjectionLocation()
    plotlyBackend_updateZMinProjectionLocation()
    plotlyBackend_updateLayout()
end

function plotlyBackend_downloadSoundButtonInteracted()
    println("WAV Button Pressed")
    values = plotlyBackend_data[end - 2]["z"] # Get the projection on the X MIN axis's z values
    maxAbsoluteValue = max(abs(plotlyBackend_ZExtremitiesWithXMinProjection[1]), abs(plotlyBackend_ZExtremitiesWithXMinProjection[2]))
    y = map(value -> (value / maxAbsoluteValue), values) # Formula: y = y divided by max(abs(Y))
    samplingFq = trunc(Int, float(size(plotlyBackend_data[end - 2]["z"])[1]) / (plotlyBackend_XExtremities[2] - plotlyBackend_XExtremities[1]))
    fileName = "ISA Sound " * Dates.format(Dates.now(), "mm-dd-yyyy HH.MM.SS")
    filePath = "./downloads/$(fileName).wav"
    wavwrite(y, filePath, Fs=samplingFq)
end

function plotlyBackend_saveButtonInteracted()
    sleep(1) # Purposely make it 1 second slower to show loading icon
    fileName = "ISA Graph " * Dates.format(Dates.now(), "mm-dd-yyyy HH.MM.SS")
    filePath = "./downloads/$(fileName).json"
    if isfile(filePath)
        println("file already exists")
        return "Save failed."
    end
    open(filePath, "w") do f
        write(f, JSON.json(plotlyBackend_dict))
    end
    return "Saved in ISA/downloads"
end

function plotlyBackend_uploaderInteracted(fileName, fileDateModified, fileContents)
    global plotlyBackend_fileUploader = Dict("fileName" => fileName, "dateModified" => fileDateModified)
    description = split(fileContents, ",")[1]
    data = split(fileContents, ",")[2]
    if !occursin("json", description)
        println("Must be a JSON file")
        return
    end
    println("Uploading file...")
    bytes = base64decode(data)
    jsonString = String(bytes)
    try
        dict = JSON.parse(jsonString)
        plotlyBackend_assignAllGlobalVariablesSafely(dict)
    catch
        println("JSON file not formatted correctly.")
        return
    end
end

# ==============================================================================
# MARK: - Updates
# ==============================================================================
# Partially mutates an already assigned global variable. Make sure all global
# variables are already assigned when using.

function plotlyBackend_updateXMinProjectionLocation()
    if plotlyBackend_data[end - 2]["visible"] && plotlyBackend_data[end - 2]["x"][1] != plotlyBackend_XSlider[1]
        plotlyBackend_data[end - 2]["x"] .= plotlyBackend_XSlider[1]
    end
end

function plotlyBackend_updateYMinProjectionLocation()
    if plotlyBackend_data[end - 1]["visible"] && plotlyBackend_data[end - 1]["y"][1] != plotlyBackend_YSlider[1]
        plotlyBackend_data[end - 1]["y"] .= plotlyBackend_YSlider[1]
    end
end

function plotlyBackend_updateZMinProjectionLocation()
    if plotlyBackend_data[end]["visible"] && plotlyBackend_data[end]["z"][1][1] != plotlyBackend_ZSlider[1]
        onePercentBuffer = float(plotlyBackend_ZSlider[2] + plotlyBackend_ZSlider[1]) / float(100) # Sometimes if we put the plane exactly on Z-MIN, it gets hidden
        for dimension in plotlyBackend_data[end]["z"]
            dimension .= (plotlyBackend_ZSlider[1] + abs(onePercentBuffer))
        end
    end
end

function plotlyBackend_updateXMinProjectionVisibility()
    plotlyBackend_data[end - 2]["visible"] = plotlyBackend_determineChecklistStatus()["FFT"]
end

function plotlyBackend_updateYMinProjectionVisibility()
    plotlyBackend_data[end - 1]["visible"] = plotlyBackend_determineChecklistStatus()["real"]
end

function plotlyBackend_updateZMinProjectionVisibility()
    plotlyBackend_data[end]["visible"] = plotlyBackend_determineChecklistStatus()["STFT"]
end

function plotlyBackend_updateLayout()
    scene = nothing
    titleText = nothing
    if plotlyBackend_viewDropDown == "XY"
        scene = plotlyBackend_determineXYScene()
        titleText = plotlyBackend_XYGraphTitle
    elseif plotlyBackend_viewDropDown == "XZ"
        scene = plotlyBackend_determineXZScene()
        titleText = plotlyBackend_XZGraphTitle
    elseif plotlyBackend_viewDropDown == "YZ"
        scene = plotlyBackend_determineYZScene()
        titleText = plotlyBackend_YZGraphTitle
    else
        scene = plotlyBackend_determine3DScene()
        titleText = plotlyBackend_3DGraphTitle
    end
    plotlyBackend_layout["scene"] = scene
    plotlyBackend_layout["title"] = titleText
end

function plotlyBackend_updateGlobalDataLineWidth()
    for component in plotlyBackend_data[1:end-1] # Make sure not to mess with Z-projection plane
        component["line"]["width"] = plotlyBackend_lineWidthSlider
    end
end

# ==============================================================================
# MARK: - Assign Globals
# ==============================================================================
# Functions that completely assigns a global variable with a new value. Note that
# force assign functions have dependencies, so please be careful when using them.

function plotlyBackend_assignAllGlobalVariablesSafely(dict)
    # Must be done in order
    plotlyBackend_assignGlobalDict(dict)
    plotlyBackend_assignGlobalLineWidth(5)
    plotlyBackend_assignGlobalData(dict, plotlyBackend_lineWidthSlider)
    plotlyBackend_assignGlobalXYZExtremities(plotlyBackend_data)
    plotlyBackend_assignGlobalFreqUnits(dict)
    plotlyBackend_assignGlobalViewDropDown("3D")
    plotlyBackend_forceAssignGlobalXMinProjectionStuff(plotlyBackend_lineWidthSlider) # depends on global data & XYZ Extremities being set
    plotlyBackend_forceAssignGlobalYMinProjectionStuff(plotlyBackend_lineWidthSlider) # depends on global data & XYZ Extremities being set & X Min projection
    plotlyBackend_forceAssignGlobalSliderValues() # depends on global XYZ & XY min projection being set
    plotlyBackend_forceAssignGlobalZMinProjectionStuff(plotlyBackend_lineWidthSlider) # depends on global data & XYZ Extremities & Z Slider being set & XY Min projection
    plotlyBackend_forceAssignSoundData() # depends on global data & XYZ Extremities being set & XYZ Min projection
    plotlyBackend_forceAssignGlobalLayout() # depends on global sliders & freqUnits & global view drop down being set
end

# REQUIRED DEPENDENCY: Global data & XYZ MIN Projection's extremities must be set
function plotlyBackend_forceAssignSoundData()
    values = plotlyBackend_data[end - 2]["z"] # Get the projection on the X MIN axis's z values
    maxAbsoluteValue = max(abs(plotlyBackend_ZExtremitiesWithXMinProjection[1]), abs(plotlyBackend_ZExtremitiesWithXMinProjection[2]))
    y = map(value -> (value / maxAbsoluteValue), values) # Formula: y = y divided by max(abs(Y))
    samplingFq = trunc(Int, float(size(plotlyBackend_data[end - 2]["z"])[1]) / (plotlyBackend_XExtremities[2] - plotlyBackend_XExtremities[1]))
    global plotlyBackend_soundData = Dict("y" => y, "samplingFq" => samplingFq)
end

# REQUIRED DEPENDENCY: Global XYZ & XYZ min projection must be assigned
function plotlyBackend_forceAssignGlobalSliderValues()
    global plotlyBackend_XSlider = plotlyBackend_XExtremities
    global plotlyBackend_YSlider = plotlyBackend_YExtremities
    global plotlyBackend_ZSlider = plotlyBackend_determineZSliderExtremities()
end

# REQUIRED DEPENDENCY: Global sliders & freqUnits & global view drop down must be assigned
function plotlyBackend_forceAssignGlobalLayout()
    scene = nothing
    titleText = nothing
    if plotlyBackend_viewDropDown == "XY"
        scene = plotlyBackend_determineXYScene()
        titleText = plotlyBackend_XYGraphTitle
    elseif plotlyBackend_viewDropDown == "XZ"
        scene = plotlyBackend_determineXZScene()
        titleText = plotlyBackend_XZGraphTitle
    elseif plotlyBackend_viewDropDown == "YZ"
        scene = plotlyBackend_determineYZScene()
        titleText = plotlyBackend_YZGraphTitle
    else
        scene = plotlyBackend_determine3DScene()
        titleText = plotlyBackend_3DGraphTitle
    end
    global plotlyBackend_layout = Dict(
        "scene" => scene,
        "title" => titleText
    )
end

# REQUIRED DEPENDENCY: Global data & XYZ Extremities must be assigned
function plotlyBackend_forceAssignGlobalXMinProjectionStuff(lineWidth)
    numberOfPointsInProjection = size(plotlyBackend_data[1]["x"])[1]
    x = repeat([plotlyBackend_XExtremities[1]], numberOfPointsInProjection)
    y = LinRange(plotlyBackend_YExtremities[1], plotlyBackend_YExtremities[2], numberOfPointsInProjection)
    z = repeat([0.0], size(x)[1])
    zMin = plotlyBackend_data[1]["z"][1]
    zMax = plotlyBackend_data[1]["z"][1]
    for component in plotlyBackend_data
        for (index, zValue) in enumerate(component["z"])
            z[index] += zValue
            z[index] += zValue
            zMin = min(zMin, min(z[index], zValue))
            zMax = max(zMax, max(z[index], zValue))
        end
    end
    projectionLine = Dict(
        "line" => Dict(
            "color" => "rgb(255,0,0)",
            "width" => lineWidth
        ),
        "mode" => "lines",
        "type" => "scatter3d",
        "opacity" => 1.0,
        "x" => x,
        "y" => y,
        "z" => z,
        "visible" => plotlyBackend_determineChecklistStatus()["FFT"],
        "showlegend" => true
    )
    push!(plotlyBackend_data, projectionLine)
    global plotlyBackend_ZExtremitiesWithXMinProjection = [zMin, zMax]
end

# REQUIRED DEPENDENCY: Global data & XYZ Extremities & X Min Projection must be assigned
function plotlyBackend_forceAssignGlobalYMinProjectionStuff(lineWidth)
    x = plotlyBackend_data[1]["x"]
    y = repeat([plotlyBackend_YExtremities[1]], size(x)[1])
    z = repeat([0.0], size(x)[1])
    zMin = plotlyBackend_data[1]["z"][1]
    zMax = plotlyBackend_data[1]["z"][1]
    for component in plotlyBackend_data[1:end-1] # We can cut of the last element since its just a projection
        for (index, zValue) in enumerate(component["z"])
            z[index] += zValue
            zMin = min(zMin, min(z[index], zValue))
            zMax = max(zMax, max(z[index], zValue))
        end
    end
    projectionLine = Dict(
        "line" => Dict(
            "color" => "rgb(192,192,192)",
            "width" => lineWidth
        ),
        "mode" => "lines",
        "type" => "scatter3d",
        "opacity" => 1.0,
        "x" => x,
        "y" => y,
        "z" => z,
        "visible" => plotlyBackend_determineChecklistStatus()["real"],
        "showlegend" => true
    )
    push!(plotlyBackend_data, projectionLine)
    global plotlyBackend_ZExtremitiesWithYMinProjection = [zMin, zMax]
end

# REQUIRED DEPENDENCY: Global data & XYZ Extremities & XY MIN Projection & Z Slider must be assigned
# If you need to loop through plotlyBackend_data make sure you cut the last 2 components
# because they are just projections on the X MIN and Y MIN axis. This projection
# is has O(N^2) data points so make sure it is reasonable in size or the UI will freeze
function plotlyBackend_forceAssignGlobalZMinProjectionStuff(lineWidth)
    sizeOfProjection = 64
    x = LinRange(plotlyBackend_XExtremities[1], plotlyBackend_XExtremities[2], sizeOfProjection) # must be length of z[1]
    y = LinRange(plotlyBackend_YExtremities[1], plotlyBackend_YExtremities[2], sizeOfProjection) # must be length of z
    z = [repeat([plotlyBackend_ZExtremities[1]], size(x)[1]) for _ in 1:size(y)[1]] # A 2D array of size: length(y) by length(x)
    minIntensity = 1.0
    maxIntensity = 0.0
    intensities = [LinRange(0.0, 1.0, size(x)[1]) for _ in 1:size(z)[1]] # A 2D array that correlatates to each color a square is
    projectionPlane = Dict(
        "type" => "surface",
        "opacity" => 1.0,
        "cmax" => maxIntensity,
        "cmin" => minIntensity,
        "showlegend" => false,
        "showscale" => false,
        "surfacecolor" => intensities,
        "colorscale" => "Viridis",
        "x" => x,
        "y" => y,
        "z" => z,
        "visible" => plotlyBackend_determineChecklistStatus()["STFT"],
        "showlegend" => false
    )
    push!(plotlyBackend_data, projectionPlane)
end

function plotlyBackend_assignGlobalDict(dict)
    global plotlyBackend_dict = dict
end

function plotlyBackend_assignGlobalLineWidth(value)
    global plotlyBackend_lineWidthSlider = value
end

function plotlyBackend_assignGlobalData(dict, lineWidth)
    data = []
    for component in dict["components"]
        myLine = Dict(
            "line" => Dict(
                "color" => component["colorMap"],
                "width" => lineWidth
            ),
            "mode" => "lines",
            "type" => "scatter3d",
            "opacity" => 1.0,
            "visible" => true,
            "x" => Float64.(component["x"]),
            "y" => Float64.(component["y"]),
            "z" => Float64.(component["z"]),
        )
        push!(data, myLine)
    end
    global plotlyBackend_data = data
end

function plotlyBackend_assignGlobalXYZExtremities(data)
    xMin = data[1]["x"][1]
    yMin = data[1]["y"][1]
    zMin = data[1]["z"][1]
    xMax = data[1]["x"][1]
    yMax = data[1]["y"][1]
    zMax = data[1]["z"][1]
    for line in data
        xMin = min(xMin, minimum(line["x"]))
        yMin = min(yMin, minimum(line["y"]))
        zMin = min(zMin, minimum(line["z"]))
        xMax = max(xMax, maximum(line["x"]))
        yMax = max(yMax, maximum(line["y"]))
        zMax = max(zMax, maximum(line["z"]))
    end
    global plotlyBackend_XExtremities = [xMin, xMax]
    global plotlyBackend_YExtremities = [yMin, yMax]
    global plotlyBackend_ZExtremities = [zMin, zMax]
end

function plotlyBackend_assignGlobalFreqUnits(dict)
    global plotlyBackend_freqUnits = haskey(dict, "freqUnits") ? dict["freqUnits"] : "rad/s"
end

function plotlyBackend_assignGlobalViewDropDown(value)
    global plotlyBackend_viewDropDown = value
end

# ==============================================================================
# MARK: - Determines
# ==============================================================================
# Returns specific values based on the currently set global variables. Make sure
# the global variables are assigned before using.

function plotlyBackend_determineIfSignalIsAudiable()
    # https://discourse.julialang.org/t/how-to-play-an-audio-in-julia/31881
    lengthInSeconds = float(size(plotlyBackend_soundData["y"])[1]) / float(plotlyBackend_soundData["samplingFq"])
    if lengthInSeconds < 0.5 return false end
    if plotlyBackend_soundData["samplingFq"] < 4000 || plotlyBackend_soundData["samplingFq"] > 50000 return false end
    return true
end

function plotlyBackend_determineXSliderDashComponent(pushData = nothing)
    dcc_rangeslider(
        id = "xSlider",
        min = plotlyBackend_XExtremities[1],
        max = plotlyBackend_XExtremities[2],
        tooltip = (
            alwaysvisible = true,
            placement = "bottom",
        ),
        pushable = (pushData == true) ? (plotlyBackend_XSlider[2] - plotlyBackend_XSlider[1]) : false,
        step = (plotlyBackend_XExtremities[2] - plotlyBackend_XExtremities[1]) / 100,
        value = (pushData == nothing) ? [plotlyBackend_XExtremities[1], plotlyBackend_XExtremities[2]] : [plotlyBackend_XSlider[1], plotlyBackend_XSlider[2]]
    )
end

function plotlyBackend_determineYSliderDashComponent(pushData = nothing)
    dcc_rangeslider(
        id = "ySlider",
        min = plotlyBackend_YExtremities[1],
        max = plotlyBackend_YExtremities[2],
        tooltip = (
            alwaysvisible = true,
            placement = "bottom",
        ),
        pushable = (pushData == true) ? (plotlyBackend_YSlider[2] - plotlyBackend_YSlider[1]) : false,
        step = (plotlyBackend_YExtremities[2] - plotlyBackend_YExtremities[1]) / 100,
        value = (pushData == nothing) ? [plotlyBackend_YExtremities[1], plotlyBackend_YExtremities[2]] : [plotlyBackend_YSlider[1], plotlyBackend_YSlider[2]]
    )
end

function plotlyBackend_determineZSliderDashComponent(pushData = nothing)
    zSliderExtremities = plotlyBackend_determineZSliderExtremities()
    return dcc_rangeslider(
        id = "zSlider",
        min = zSliderExtremities[1],
        max = zSliderExtremities[2],
        tooltip = (
            alwaysvisible = true,
            placement = "bottom",
        ),
        pushable = (pushData == true) ? (plotlyBackend_ZSlider[2] - plotlyBackend_ZSlider[1]) : false,
        step = (zSliderExtremities[2] - zSliderExtremities[1]) / 100,
        value = (pushData == nothing) ? zSliderExtremities : [plotlyBackend_ZSlider[1], plotlyBackend_ZSlider[2]]
    )
end

# Returns a [min, max] extremity that should be assigned to the z slider
function plotlyBackend_determineZSliderExtremities()
    result = [plotlyBackend_ZExtremities[1], plotlyBackend_ZExtremities[2]]
    extraConsideration = []
    checklistStatus = plotlyBackend_determineChecklistStatus()
    if checklistStatus["FFT"] push!(extraConsideration, plotlyBackend_ZExtremitiesWithXMinProjection) end
    if checklistStatus["real"] push!(extraConsideration, plotlyBackend_ZExtremitiesWithYMinProjection) end
    for extremity in extraConsideration
        result[1] = min(result[1], extremity[1])
        result[2] = max(result[2], extremity[2])
    end
    return result
end

# Returns a dictionary that determines whether or not a checkbox in the checklist is checked
function plotlyBackend_determineChecklistStatus()
    checklistStatus = Dict()
    checklistStatus["FFT"] = (plotlyBackend_projectionChecklist != nothing && 1 in plotlyBackend_projectionChecklist)
    checklistStatus["real"] = (plotlyBackend_projectionChecklist != nothing && 2 in plotlyBackend_projectionChecklist)
    checklistStatus["STFT"] = (plotlyBackend_projectionChecklist != nothing && 3 in plotlyBackend_projectionChecklist)
    return checklistStatus
end

function plotlyBackend_determine3DScene()
    return (
        aspectmode = "manual",
        aspectratio = (x = 1, y = 1, z = 0.25),
        bgcolor = "rgb(255, 255, 255)",
        camera = (
            center = (x = 0, y = 0, z = 0),
            eye = (x = -0.5, y = -2.25, z = 1.25),
            projection = (type = "perspective",),
            up = (x = 0, y = 0, z = 1),
        ),
        dragmode = true,
        xaxis = (
            backgroundcolor = "rgba(0, 0, 128, 0.35)",
            range = plotlyBackend_XSlider,
            showbackground = true,
            showgrid = false,
            showticklabels = true,
            tickangle = "auto",
            title = (text = "Time", font = (color = 0, size = 10))
        ),
        yaxis = (
            backgroundcolor = "rgba(0, 0, 0, 0.2)",
            range = plotlyBackend_YSlider,
            showbackground = true,
            showgrid = false,
            showticklabels = true,
            tickangle = "auto",
            title = (text = string("Frequency in ", plotlyBackend_freqUnits), font = (color = 0, size = 10))
        ),
        zaxis = (
            backgroundcolor = "rgba(140, 60, 250, 1)",
            range = plotlyBackend_ZSlider,
            showbackground = true,
            showgrid = false,
            showticklabels = true,
            tickangle = "auto",
            title = (text = "Real", font = (color = 0, size = 10))
        )
    )
end

function plotlyBackend_determineXYScene()
    return (
        aspectmode = "manual",
        aspectratio = (x = 1, y = 1, z = 0.25),
        bgcolor = "rgb(255, 255, 255)",
        camera = (
            center = (x = 0, y = 0, z = 0),
            eye = (x = 0, y = 0, z = 1),
            projection = (type = "orthographic",),
            up = (x = 0, y = 1, z = 0),
        ),
        dragmode = false,
        xaxis = (
            backgroundcolor = "rgba(0, 0, 128, 0.35)",
            range = plotlyBackend_XSlider,
            showbackground = true,
            showgrid = false,
            showticklabels = true,
            tickangle = 0,
            title = (text = "Time", font = (color = 0, size = 10))
        ),
        yaxis = (
            backgroundcolor = "rgba(0, 0, 0, 0.2)",
            range = plotlyBackend_YSlider,
            showbackground = true,
            showgrid = false,
            showticklabels = true,
            tickangle = "auto",
            title = (text = string("Frequency in ", plotlyBackend_freqUnits), font = (color = 0, size = 10))
        ),
        zaxis = (
            backgroundcolor = "rgba(140, 60, 250, 1)",
            range = plotlyBackend_ZSlider,
            showbackground = true,
            showgrid = false,
            showticklabels = false,
            tickangle = "auto",
            title = (text = "",)
        )
    )
end

function plotlyBackend_determineXZScene()
    return (
        aspectmode = "manual",
        aspectratio = (x = 1, y = 1, z = 0.25),
        bgcolor = "rgb(255, 255, 255)",
        camera = (
            center = (x = 0, y = 0, z = 0),
            eye = (x = 0, y = -1, z = 0),
            projection = (type = "orthographic",),
            up = (x = 0, y = 0, z = 1),
        ),
        dragmode = false,
        xaxis = (
            backgroundcolor = "rgba(0, 0, 128, 0.35)",
            range = plotlyBackend_XSlider,
            showbackground = true,
            showgrid = false,
            showticklabels = true,
            tickangle = 0,
            title = (text = "Time", font = (color = 0, size = 10))
        ),
        yaxis = (
            backgroundcolor = "rgba(0, 0, 0, 0.2)",
            range = plotlyBackend_YSlider,
            showbackground = true,
            showgrid = false,
            showticklabels = false,
            tickangle = "auto",
            title = (text = "")
        ),
        zaxis = (
            backgroundcolor = "rgba(140, 60, 250, 1)",
            range = plotlyBackend_ZSlider,
            showbackground = true,
            showgrid = false,
            showticklabels = true,
            tickangle = "auto",
            title = (text = "Real", font = (color = 0, size = 10))
        )
    )
end

function plotlyBackend_determineYZScene()
    return (
        aspectmode = "manual",
        aspectratio = (x = 1, y = 1, z = 0.25),
        bgcolor = "rgb(255, 255, 255)",
        camera = (
            center = (x = 0, y = 0, z = 0),
            eye = (x = 1, y = 0, z = 0),
            projection = (type = "orthographic",),
            up = (x = 0, y = 0, z = 1),
        ),
        dragmode = false,
        xaxis = (
            backgroundcolor = "rgba(0, 0, 128, 0.35)",
            range = plotlyBackend_XSlider,
            showbackground = true,
            showgrid = false,
            showticklabels = false,
            tickangle = "auto",
            title = (text = "")
        ),
        yaxis = (
            backgroundcolor = "rgba(0, 0, 0, 0.2)",
            range = plotlyBackend_YSlider,
            showbackground = true,
            showgrid = false,
            showticklabels = true,
            tickangle = 0,
            title = (text = string("Frequency in ", plotlyBackend_freqUnits), font = (color = 0, size = 10))
        ),
        zaxis = (
            backgroundcolor = "rgba(140, 60, 250, 1)",
            range = plotlyBackend_ZSlider,
            showbackground = true,
            showgrid = false,
            showticklabels = true,
            tickangle = "auto",
            title = (text = "Real", font = (color = 0, size = 10))
        )
    )
end

# ==============================================================================
# MARK: - Dash HMTL
# ==============================================================================
# The HTML the dash server will present.

function plotlyBackend_createHTML()
    return html_div() do
        html_div() do
            html_div(style = (width = "275px", display = "inline-block", marginRight = "20px"),) do
                dcc_dropdown(
                    id = "viewDropDown",
                    options = [
                        (label = "3D", value = "3D"),
                        (label = "2D X & Y (Time and Frequency)", value = "XY"),
                        (label = "2D X & Z (Time and Real)", value = "XZ"),
                        (label = "2D Y & Z (Frequency and Real)", value = "YZ")
                    ],
                    value = "3D",
                    searchable = false,
                    clearable = false
                )
            end,
            html_div(style = (width = "250px", display = "inline-block", transform = "translateY(-20%)", marginRight = "20px")) do
                html_div() do
                    html_div(style = (display = "inline-block", marginRight = "20px", transform = "translateY(-5%)"),) do
                        html_p("Line Width")
                    end,
                    html_div(style = (width = "60%", display = "inline-block"),) do
                        dcc_slider(
                            id = "lineWidthSlider",
                            min = 1,
                            max = 10,
                            step = 1,
                            tooltip = (
                                alwaysvisible = true,
                                placement = "bottom",
                            ),
                            value = plotlyBackend_lineWidthSlider
                        )
                    end
                end
            end,
            html_div(style = (display = "inline-block", transform = "translateY(-75%)"),) do
                dcc_checklist(
                    id = "projectionChecklist",
                    options = [
                        (label = "FFT Projection", value = 1),
                        (label = "Real Projection", value = 2),
                        (label = "STFT Projection", value = 3)
                    ],
                    labelStyle = (display = "inline-block", marginRight = "20px")
                )
            end,
            html_div(style = (display = "inline-block", transform = "translateY(-75%)"),) do
                html_button(
                    id = "downloadSoundButton",
                    n_clicks = -1,
                    children = ["Download Sound"],
                    disabled = !plotlyBackend_determineIfSignalIsAudiable()
                ),
                html_div(id = "downloadSoundButtonOutletDiv")
            end
        end,
        dcc_graph(
            id = "isaGraph",
            figure = (
                data = plotlyBackend_data,
                layout = plotlyBackend_layout
            )
        ),
        html_div(style = (textAlign = "center",)) do
            html_div() do
                html_div(style = (display = "inline-block", marginRight = "20px", transform = "translateY(-5%)"),) do
                    html_p("Time Axis Traversal")
                end,
                html_div(id = "xSliderDiv", style = (width = "50%", display = "inline-block", marginRight = "20px"),) do
                    plotlyBackend_determineXSliderDashComponent()
                end,
                html_div(id = "xSliderPushableDiv", style = (width = "100px", display = "inline-block", transform = "translateY(-10%)"),) do
                    dcc_checklist(
                        id = "xSliderPushable",
                        options = [(label = "Push Mode", value = 1)],
                    )
                end
            end,
            html_div() do
                html_div(style = (display = "inline-block", marginRight = "20px", transform = "translateY(-5%)"),) do
                    html_p("Frequency Axis Traversal")
                end,
                html_div(id = "ySliderDiv", style = (width = "50%", marginRight = "20px", display = "inline-block"),) do
                    plotlyBackend_determineYSliderDashComponent()
                end,
                html_div(id = "ySliderPushableDiv", style = (width = "100px", display = "inline-block", transform = "translateY(-10%)"),) do
                    dcc_checklist(
                        id = "ySliderPushable",
                        options = [(label = "Push Mode", value = 1)],
                    )
                end
            end,
            html_div() do
                html_div(style = (display = "inline-block", marginRight = "20px", transform = "translateY(-5%)"),) do
                    html_p("Real Axis Traversal")
                end,
                html_div(id = "zSliderOutterDiv", style = (width = "50%", marginRight = "20px", display = "inline-block"),) do
                    html_div(id = "zSliderDiv") do
                        plotlyBackend_determineZSliderDashComponent()
                    end
                end,
                html_div(id = "zSliderPushableDiv", style = (width = "100px", display = "inline-block", transform = "translateY(-10%)"),) do
                    dcc_checklist(
                        id = "zSliderPushable",
                        options = [(label = "Push Mode", value = 1)],
                    )
                end
            end
        end,
        html_br(),
        html_br(),
        html_div(style = (display = "inline",),) do
            html_div(style = (width = "19%", display = "inline-block", float = "left", textAlign = "center", transform = "translateY(50%)", top = "50%", position = "relative")) do
                html_h3("Save Current Signal"),
                dcc_loading(
                    type = "default",
                    children = [html_div(id = "downloadStatusDiv",)]
                ),
                html_button(
                    id = "saveButton",
                    n_clicks = -1,
                    children = ["Save!"]
                )
            end,
            html_div(style = (width = "80%", display = "inline-block")) do
                html_h3("Visualize a Different Signal"),
                html_p("Instructions: This can take several seconds. Try not to interact with the UI while uploading a file. Seldomly, the uploader will have troubles with layout. If this happens, just toggle the projection buttons and the layout should correct itself. If there are further issues, please reupload until correct. Make sure to upload a correctly formatted JSON file with valid data. Note: the key \'freqUnits\' is optional and is defaulted to \'rad/s\'. See valid file examples in ISA/downloads."),
                dcc_loading(
                    type = "default",
                    children = [html_div(id = "uploadFileDiv",)]
                ),
                dcc_upload(
                    id = "uploader",
                    children = [
                        "Drag a JSON file or click to select"
                    ],
                    style = (
                        width = "100%",
                        height = "60px",
                        lineHeight = "60px",
                        borderWidth = "1px",
                        borderStyle = "dashed",
                        borderRadius = "5px",
                        textAlign = "center",
                    ),
                    # Do not multiple files to be uploaded
                    multiple = false
                )
            end
        end
    end
end

ext: jl

max_stars_repo_name: kricheso/ISA

max_issues_repo_name: kricheso/ISA

max_forks_repo_name: kricheso/ISA


text:
module MonteCarloFun

"
arxiv doc :     1008.4686
                1205.4446

plan:
        1) Fitting a model to data
                Fitting a straight line
                magic of Gaussian

"

#2018_10_05
#fitting a straight line

using Plots
using LinearAlgebra

x = [1., 2., 3., 4.]
y = [14, 26, 37, 43]
sigma = [1, 1, 2, 2]

plot(x, y, seriestype=:scatter, yerr=sigma, legend=false)

Y=y
N = length(Y)
A = zeros(N,2)

A[:,1] .= 1.;
A[:,2] .= x;

W = Diagonal(1. ./ sigma.^2)

MB1 = A \ Y
MB2 = (sqrt(W)*A) \ (sqrt(W)*Y)

plot(x, y, seriestype=:scatter, yerr=sigma, legend=:topleft)
plot!(x, A * MB1, label="Unweighted")
plot!(x, A * MB2, label="Weighted")

C = inv(A' * W * A)

dMB = randn((100,2)) * sqrt(C)

Angle = collect(range(0,stop=2π,length=100))
xa = [sin.(Angle) cos.(Angle)]
cMB = xa * sqrt(C)

plot([MB1[1]], [MB1[2]], seriestype=:scatter, label="Unweighted",xlabel="b", ylabel="m")
plot!([MB2[1]], [MB2[2]], seriestype=:scatter, label="Weighted",xlabel="b", ylabel="m")
plot!(MB2[1] .+ dMB[:,1], MB2[2] .+ dMB[:,2], seriestype=:scatter, label="Weighted")
plot!(MB2[1] .+ cMB[:,1], MB2[2] .+ cMB[:,2], seriestype=:scatter, label="Weighted")

plot(x, y, seriestype=:scatter, yerr=sigma, legend=:topleft)
plot!(x, A * MB2, line=(:blue, 1, 3), label="Weighted")
for i in 1:50
    dmb = MB2 + dMB[i,:]
    plot!(x, A * dmb, label= "" , line=(:black, 0.2, 2))
end
plot!()

function line_log_likehood(x,y,sigma,m,b)
    return -0.5 * sum(@. (y -(m*x+b))^2 / sigma^2)
end

bbest,mbest = MB2
line_log_likehood(x,y,sigma,mbest,bbest)

function line_lnl(b,m)
    return line_log_likehood(x, y, sigma, m, b)
end

function make_line_lnl_function(x,y,sigma)
    function f(m,b)
        return line_log_likehood(x,y,sigma,m,b)
    end
    return f
end

LineLnl = make_line_lnl_function(x,y,sigma)
LineLnl(mbest,bbest)

end # module


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

makedocs(sitename="EmpiricalModeDecomposition.jl Documentation",
        pages = [
                 "Library" => [
                    "Reference" => "reference.md",
                    "Extras" => "extras.md",
                    "Internals" => "internals.md"
                    ],
                "Methods" => "methods.md"
              ]
       )
deploydocs(
           repo = "github.com/atmnpatel/EmpiricalModeDecomposition.jl.git"
)


text:
function Base.merge!(
    combine::Function,
    d::MuttsDict,
    others::Union{MuttsDict,AbstractDict}...,
)
    for other in others
        for (k,v) in other
            if haskey(d, k)
                d[k] = combine(d[k], v)
            else
                d[k] = v
            end
        end
    end
    return d
end

function Base.merge!(d::MuttsDict, others::Union{MuttsDict,AbstractDict}...)
    for other in others
        for (k,v) in other
            d[k] = v
        end
    end
    return d
end

function Base.merge(d::MuttsDict, others::Union{MuttsDict,AbstractDict}...)
    result = branch!(d)
    Base.merge!(result, others...)
    mark_immutable!(result)
    return result
end


text:
using Revise
using JuMP,Gurobi,BilevelJuMP
includet("./bilevelJuMPExtension.jl")
# min -4sum(x) -3sum(y)
# s.t.
# y = argmin_y sum(y)
#      2x + y <= 4
#       x +2y <= 4
#       x     >= 0
#           y >= 0
#
# sol: x = 2, y = 0

# opt is the orginal problem
opt =  quote
    function some_model()
        blm = BilevelModel()
        sz = 10
        um = Upper(blm)
        lm = Lower(blm)
        @variable(um, x[1:sz])
        @variable(lm, y[1:sz])

        @objective(um, Min, -4sum(x) -3sum(y))

        @objective(lm, Min, sum(y))

        for i in 1:sz
            @constraint(lm, 2x[i]+y[i] <= 4)
            @constraint(lm, x[i]+2y[i] <= 4)
            @constraint(lm, x[i] >= 0)
            @constraint(lm, y[i] >= 0)
        end
        @constraint(lm,con[i in 1:sz],y[i]==1) #sign
        return blm
    end
end

# test.1 simple removal and addition of constarint named "con"
printstyled("****REMOVAL OF CONSTRAINT CON****\n";color=:blue)
opt = @delete_cons(opt,con)
mdl = eval(opt)()
optimize!(mdl,Gurobi.Optimizer())
@assert value(mdl[:y][1]) ≈ 0
printstyled("****ADDITION OF CONSTRAINT CON****\n";color=:blue)
opt = @add_lower_cons(opt,con,opt[:y][1] == 1)
mdl = eval(opt)()
optimize!(mdl,Gurobi.Optimizer())
@assert value(mdl[:y][1]) ≈ 1
printstyled("****SIMPLE TEST PASSED****\n";color=:green)

#test.2 complex removal and addition of constarint named "con" within loops
opt = [opt,opt]
n = 10
@time begin
for ii = 1:2
    for i = 1:2
        for j = 1:2
            printstyled("****REMOVAL OF CONSTRAINT CON****\n";color=:blue)
            opt[ii] = @delete_cons(opt[ii],con)
            f = eval(opt[ii])
            optimize!(f(),Gurobi.Optimizer(OutputFlag=0))
            sth = [0.5,0.3]
            idx = 2
            printstyled("****ADDITION OF CONSTRAINT CON****\n";color=:blue)
            opt[ii] = @add_lower_cons(opt[ii],con,sum(opt[ii][:y][i] + sth[idx]^2 for i in 1:n) == 10 + 0.1*(i+j))
            # println(opt[ii])
            f = eval(opt[ii])
            optimize!(f(),Gurobi.Optimizer(OutputFlag=0))
        end
    end
end
end
printstyled("****COMPLEX ADDITION TEST WITH LOOPS PASSED****\n";color=:green)

#test.3 batch addtion of constraints
@time begin
    for t in 1:2
        printstyled("****REMOVAL OF CONSTRAINT CON****\n";color=:blue)
        opt[t] = @delete_cons(opt[t],con)
        f = eval(opt[t])
        optimize!(f(),Gurobi.Optimizer(OutputFlag=0))
        sth = [1 for i in 1:n]
        printstyled("****ADDITION OF CONSTRAINT CON****\n";color=:blue)
        opt[t] = @add_lower_cons(opt[t],con[i = 1:n],opt[t][:y][i] == sth[i])
        # println(opt[t])
        f = eval(opt[t])
        optimize!(f(),Gurobi.Optimizer(OutputFlag=0))
    end
end
printstyled("****BATCH CONSTRAINTS ADDITION PASSED****\n";color=:green)


text:
import .BandedMatrices: BandedMatrix, isbanded, bandwidth

const BandedDerivativeOperator= Union{DerivativeOperator,DissipationOperator,VarCoefDerivativeOperator}

@inline function bandwidth(D::BandedDerivativeOperator, k::Int)
    if k == 1
        lower_bandwidth(D)
     else
        upper_bandwidth(D)
    end
end

isbanded(D::BandedDerivativeOperator) = true

function BandedMatrix(D::BandedDerivativeOperator)
    T = eltype(D)
    l = lower_bandwidth(D)
    u = upper_bandwidth(D)
    n,m = size(D)

    # create uninitialised matrix
    B = BandedMatrix(BandedMatrices.Zeros{T}(n, m), (l, u))

    # B.data[i,:] is the ith band, starting with the upper ones
    # B.data[:,j] is the jth column, entries out of range are undefined
    e = fill(zero(T), length(grid(D)))
    dest = similar(e)
    # left boundary
    for j in 1:l
        e[j] = 1
        mul!(dest, D, e)
        e[j] = 0
        B.data[end-j-u+1:end,j] = dest[1:j+u]
    end
    # inner part
    for j in l+1:length(e)-u
        e[j] = 1
        mul!(dest, D, e)
        e[j] = 0
        B.data[:,j] = dest[j-l:j+u]
    end
    # right boundary
    for j in length(e)-u+1:length(e)
        e[j] = 1
        mul!(dest, D, e)
        e[j] = 0
        B.data[1:length(e)+1-j+l,j] = dest[j-l:end]
     end

    B
end


text:
using Test
using BayesianTools
using StatsFuns
using Random
const T = typeof(1.0)

using BayesianTools.ProductDistributions

d1 = ProductDistribution(Normal(0,1), Uniform(0,1), InverseGamma(0.01, 0.01))

@test length(d1) == 3
@test insupport(d1, [0., 1., 1.])
@test insupport(d1, [0., .1, 21.])
@test !insupport(d1, [0., 2., .1])
@test !insupport(d1, [0., -2., .1])
@test !insupport(d1, [0., .1, -2])

@test logpdf(d1, [0., 0.5, 0.5]) .== logpdf(d1.marginals[1], 0.) + logpdf(d1.marginals[2], 0.5) + logpdf(d1.marginals[3], 0.5)
@test pdf(d1, [0., 0.5, 0.5]) ≈ pdf(d1.marginals[1], 0.)*pdf(d1.marginals[2], 0.5)*pdf(d1.marginals[3], 0.5)

d1 = Improper(0, +Inf)
d2 = Improper(-Inf, +Inf)
d3 = Improper(-Inf, 5)
d4 = Improper(5, +Inf)
d5 = Improper()

@test d2===d5
@test_throws AssertionError Improper(0,1)
@test_throws AssertionError Improper(+Inf,-Inf)

@test !insupport(d1, -1)
@test insupport(d1, 1)

@test minimum(d1) == 0
@test minimum(d2) == -Inf
@test minimum(d3) == -Inf
@test minimum(d4) == 5

@test maximum(d1) == +Inf
@test maximum(d2) == +Inf
@test maximum(d3) == 5.0
@test maximum(d4) == +Inf

@test pdf(d1, 100) == 1.0
@test logpdf(d1, 100) == 0.0
@test pdf(d1, -100) == 0.0
@test logpdf(d1, -100) == Inf

@test pdf(d2, 100) == 1.0
@test logpdf(d2, 100) == 0.0
@test pdf(d2, -100) == 1.0
@test logpdf(d2, -100) == 0.0

@test pdf(d4, 7) == 1.0
@test logpdf(d4, 7) == 0.0
@test pdf(d4, 4) == 0.0
@test logpdf(d4, 4) == Inf

using BayesianTools.Links

for d in [Uniform(0,1), Beta(2,2)]
    @test link(d, 0.1) == StatsFuns.logit(0.1)
    @test invlink(d, link(d, 0.1)) ≈ 0.1
end

for d in [Normal(0,1), TDist(2)]
    @test link(d, 0.1) == 0.1
    @test invlink(d, link(d, 0.1)) ≈ 0.1
end
d = ProductDistribution(Normal(0,1), Uniform(0,1), InverseGamma(0.01, 0.01))
@test link(d, invlink(d, [0.0, -1., 2.])) ≈ [0.0, -1., 2.] atol=1e-15

d = ProductDistribution(Normal(0,1), Uniform(0,1), InverseGamma(0.01, 0.01), Improper(0,+Inf))
@test link(d, invlink(d, [0.0, -1., 2., -2.])) ≈ [0.0, -1., 2., -2.] atol=1e-15

d = ProductDistribution(Normal(), Uniform())
Random.seed!(1)
x2 = hcat([rand( d ) for i in 1:1000]...)
@test maximum(x2[2,:])<=1.0
@test maximum(x2[1,:])>=0.0

@test(params(d) == (0.0, 1.0, 0.0, 1.0))

max_stars_repo_name: gragusa/ProductDistributions

max_issues_repo_name: gragusa/ProductDistributions

max_forks_repo_name: gragusa/ProductDistributions


text:
using FactCheck, TimeSeriesTools

include("data_quality.jl")
include("visualizations.jl")

exitstatus()


text:
module ProgressiveHedging

import DataFrames
import JuMP
import MathOptInterface
const MOI = MathOptInterface

using Distributed
using Printf
using TimerOutputs

#### Exports ####

# High-level Functions
export solve, solve_extensive

# Callbacks
export Callback
export cb, mean_deviation, variable_fixing
export apply_to_subproblem
export SubproblemCallback, spcb

# Exceptions
export UnimplementedError

# ID types and functions
export Index, NodeID, ScenarioID, StageID, VariableID, XhatID
export index, scenario, stage, value
export scid, stid, index
export convert_to_variable_ids, convert_to_xhat_id, stage_id

# Penalty Parameter Interface
export AbstractPenaltyParameter
export get_penalty_value
export is_initial_value_dependent, is_subproblem_dependent, is_variable_dependent
export penalty_map
export process_penalty_initial_value
export process_penalty_subproblem

# Penalty Parameter (Concrete)
export ProportionalPenaltyParameter, ScalarPenaltyParameter, SEPPenaltyParameter

# PHData interaction function
export PHData
export consensus_variables, probability, scenario_tree, scenarios
export get_callback, get_callback_ext

# Result retrieval functions
export lower_bounds
export print_timing
export residuals
export retrieve_aug_obj_value, retrieve_obj_value
export retrieve_soln, retrieve_no_hats, retrieve_w
export retrieve_xhat_history, retrieve_no_hat_history, retrieve_w_history

# Scenario tree types and functions
export ScenarioNode, ScenarioTree
export add_node, add_leaf, root, two_stage_tree

# Subproblem Interface
export AbstractSubproblem, VariableInfo
export add_ph_objective_terms, objective_value, report_values, report_variable_info
export solve_subproblem, update_ph_terms
export warm_start
export ef_copy_model
export ef_node_dict_constructor
export report_penalty_info
export add_lagrange_terms, update_lagrange_terms

# Subproblems Types (Concrete)
export JuMPSubproblem

# (Consensus) Variable interaction functions
export HatVariable
export is_leaf, name, value, branch_value, leaf_value, w_value, xhat_value
export is_integer, scenario_bundle, value, variables

#### Includes ####

include("id_types.jl")
include("penalty_parameter_types.jl")
include("scenario_tree.jl")

include("subproblem.jl")
include("jumpsubproblem.jl")

include("subproblem_callback.jl")
include("message.jl")
include("worker.jl")
include("worker_management.jl")

include("structs.jl")
include("penalty_parameter_functions.jl")

include("utils.jl")

include("algorithm.jl")
include("setup.jl")

include("callbacks.jl")

#### Functions ####

"""
    solve(tree::ScenarioTree,
          subproblem_constructor::Function,
          r<:Real,
          other_args...;
          max_iter::Int=1000,
          atol::Float64=1e-6,
          rtol::Float64=1e-6,
          gap_tol::Float64=-1.0,
          lower_bound::Int=0,
          report::Int=0,
          save_iterates::Int=0,
          save_residuals::Bool=false,
          timing::Bool=true,
          warm_start::Bool=false,
          callbacks::Vector{Callback}=Vector{Callback}(),
          worker_assignments::Dict{Int,Set{ScenarioID}}=Dict{Int,Set{ScenarioID}}(),
          args::Tuple=(),
          kwargs...)

Solve the stochastic programming problem described by `tree` and the models created by `subproblem_constructor` using Progressive Hedging.

**Arguments**

* `tree::ScenararioTree` : Scenario tree describing the structure of the problem to be solved.
* `subproblem_constructor::Function` : User created function to construct a subproblem. Should accept a `ScenarioID` (a unique identifier for each scenario subproblem) as an argument and returns a subtype of `AbstractSubproblem`.
* `r<:AbstractPenaltyParameter` : PH penalty parameter
* `other_args` : Other arguments that should be passed to `subproblem_constructor`. See also keyword arguments `args` and `kwargs`.

**Keyword Arguments**

* `max_iter::Int` : Maximum number of iterations to perform before returning. Defaults to 1000.
* `atol::Float64` : Absolute error tolerance. Defaults to 1e-6.
* `rtol::Float64` : Relative error tolerance. Defaults to 1e-6.
* `gap_tol::Float64` : Relative gap tolerance. Terminate when the relative gap between the lower bound and objective are smaller than `gap_tol`. Any value < 0.0 disables this termination condition. Defaults to -1.0. See also the `lower_bound` keyword argument.
* `lower_bound::Int` : Compute and save a lower-bound using (Gade, et. al. 2016) every `lower_bound` iterations. Any value <= 0 disables lower-bound computation. Defaults to 0.
* `report::Int` : Print progress to screen every `report` iterations. Any value <= 0 disables printing. Defaults to 0.
* `save_iterates::Int` : Save PH iterates every `save_iterates` steps. Any value <= 0 disables saving iterates. Defaults to 0.
* `save_residuals::Int` : Save PH residuals every `save_residuals` steps. Any value <= 0 disables saving residuals. Defaults to 0.
* `timing::Bool` : Print timing info after solving if true. Defaults to true.
* `warm_start::Bool` : Flag indicating that solver should be "warm started" by using the previous solution as the starting point (not compatible with all solvers)
* `callbacks::Vector{Callback}` : Collection of `Callback` structs to call after each PH iteration. Callbacks will be executed in the order they appear. See `Callback` struct for more info. Defaults to empty vector.
* `subproblem_callbacks::Vector{SubproblemCallback}` : Collection of `SubproblemCallback` structs to call before solving each subproblem. Each callback is called on each subproblem but does not affect other subproblems. See `SubproblemCallback` struct for more info. Defaults to empty vector.
* `worker_assignments::Dict{Int,Set{ScenarioID}}` : Dictionary specifying which scenario subproblems a worker will create and solve. The key values are worker ids as given by Distributed (see `Distributed.workers()`). The user is responsible for ensuring the specified workers exist and that every scenario is assigned to a worker. If no dictionary is given, scenarios are assigned to workers in round robin fashion. Defaults to empty dictionary.
* `args::Tuple` : Tuple of arguments to pass to `model_cosntructor`. Defaults to (). See also `other_args` and `kwargs`.
* `kwargs` : Any keyword arguments not specified here that need to be passed to `subproblem_constructor`.  See also `other_args` and `args`.

**Return Arguments**

* `niter` : Number of iterations performed.
* `abs_res` : Absolute residual. The sum of the L2 norm of movement between consensus variables from one iteration to the next and the L2 norm of the lack of consensus divided by the number of consensus variables.
* `rel_res` : Relative residual. The absolute residual divided by the maximum consensus variable value.
* `obj` : Objective value with progressive hedging terms removed.
* `soln_df` : `DataFrame` containing variables, their values and their stage and scenario IDs.
* `phd` : A data structure used for the progressive hedging algorithm. See `PHData` for more information.
"""
function solve(tree::ScenarioTree,
               subproblem_constructor::Function,
               r::R,
               other_args...;
               max_iter::Int=1000,
               atol::Float64=1e-6,
               rtol::Float64=1e-6,
               gap_tol::Float64=-1.0,
               lower_bound::Int=0,
               report::Int=0,
               save_iterates::Int=0,
               save_residuals::Int=0,
               timing::Bool=false,
               warm_start::Bool=false,
               callbacks::Vector{Callback}=Vector{Callback}(),
               subproblem_callbacks::Vector{SubproblemCallback}=Vector{SubproblemCallback}(),
               worker_assignments::Dict{Int,Set{ScenarioID}}=Dict{Int,Set{ScenarioID}}(),
               args::Tuple=(),
               kwargs...
               ) where {R <: AbstractPenaltyParameter}

    timo = TimerOutputs.TimerOutput()

    reset_penalty_map(r)

    if length(scenarios(tree)) == 1
        @warn("Given scenario tree indicates a deterministic problem (only one scenario).")
    elseif length(scenarios(tree)) <= 0
        error("Given scenario tree has no scenarios specified. Make sure 'add_leaf' is being called on leaves of the scenario tree.")
    end

    psum = sum(values(tree.prob_map))
    if !isapprox(psum, 1.0, atol=1e-8)
        error("Total probability of scenarios in given scenario tree is $psum.")
    end

    # Initialization
    if report > 0
        println("Initializing...")
    end

    (phd, winf) = @timeit(timo,
                          "Intialization",
                          initialize(tree,
                                     subproblem_constructor,
                                     r,
                                     worker_assignments,
                                     warm_start,
                                     timo,
                                     report,
                                     lower_bound,
                                     subproblem_callbacks,
                                     (other_args...,args...);
                                     kwargs...)
                          )

    for cb in callbacks
        add_callback(phd, cb)
    end

    # Solution
    if report > 0
        println("Solving...")
    end
    (niter, abs_res, rel_res) = @timeit(timo,
                                        "Solution",
                                        hedge(phd,
                                              winf,
                                              max_iter,
                                              atol,
                                              rtol,
                                              gap_tol,
                                              report,
                                              save_iterates,
                                              save_residuals,
                                              lower_bound)
                                        )

    # Post Processing
    if report > 0
        println("Done.")
    end

    soln_df = retrieve_soln(phd)
    obj = retrieve_obj_value(phd)

    if timing
        print_timing(phd)
    end

    return (niter, abs_res, rel_res, obj, soln_df, phd)
end

"""
    solve_extensive(tree::ScenarioTree,
          subproblem_constructor::Function,
          optimizer::Function,
          other_args...;
          opt_args::NamedTuple=NamedTuple(),
          subproblem_type::Type{S}=JuMPSubproblem,
          args::Tuple=(),
          kwargs...)

Solve given problem using Progressive Hedging.

**Arguments**

* `tree::ScenararioTree` : Scenario tree describing the structure of the problem to be solved.
* `subproblem_constructor::Function` : User created function to construct a subproblem. Should accept a `ScenarioID` (a unique identifier for each scenario subproblem) as an argument and returns a subtype of `AbstractSubproblem` specified by `subproblem_type`.
* `optimizer::Function` : Function which works with `JuMP.set_optimizer`
* `other_args` : Other arguments that should be passed to `subproblem_constructor`. See also keyword arguments `args` and `kwargs`

**Keyword Arguments**

* `subproblem_type<:JuMP.AbstractModel` : Type of model to create or created by `subproblem_constructor` to represent the subproblems. Defaults to JuMPSubproblem
* `opt_args::NamedTuple` : arguments passed to function given by `optimizer`
* `args::Tuple` : Tuple of arguments to pass to `model_cosntructor`. Defaults to (). See also `other_args` and `kwargs`.
* `kwargs` : Any keyword arguments not specified here that need to be passed to `subproblem_constructor`.  See also `other_args` and `args`.
"""
function solve_extensive(tree::ScenarioTree,
                         subproblem_constructor::Function,
                         optimizer::Function,
                         other_args...;
                         opt_args::NamedTuple=NamedTuple(),
                         subproblem_type::Type{S}=JuMPSubproblem,
                         args::Tuple=(),
                         kwargs...
                         ) where {S <: AbstractSubproblem}

    if length(scenarios(tree)) == 1
        @warn("Given scenario tree indicates a deterministic problem (only one scenario).")
    elseif length(scenarios(tree)) <= 0
        error("Given scenario tree has no scenarios specified. Make sure 'add_leaf' is being called on leaves of the scenario tree.")
    end

    psum = sum(values(tree.prob_map))
    if !isapprox(psum, 1.0, atol=1e-8)
        error("Total probability of scenarios in given scenario tree is $psum.")
    end

    model = build_extensive_form(tree,
                                 subproblem_constructor,
                                 Tuple([other_args...,args...]),
                                 subproblem_type;
                                 kwargs...)

    JuMP.set_optimizer(model, optimizer)
    for (key, value) in pairs(opt_args)
        JuMP.set_optimizer_attribute(model, string(key), value)
    end
    JuMP.optimize!(model)

    return model
end

end # module

#### TODOs ####

# 1. Add flag to solve call to turn off parallelism
# 2. Add handling of nonlinear objective functions
# 3. Add tests for nonlinear constraints (these should work currently but testing is needed)


text:
import Base: length, first, last, getindex

abstract IndexRepr

immutable LineRepr <: IndexRepr
    id::NTuple{2, Int}
end
interior_indices(::LineRepr) = ()

immutable CircleArcRepr <: IndexRepr
    id::NTuple{3, Int}
end
interior_indices(::CircleArcRepr) = (2,)

immutable EllipseArcRepr <: IndexRepr
    id::NTuple{3, Int}
end
interior_indices(::EllipseArcRepr) = (2,)

for (T, size) in zip((:LineRepr, :CircleArcRepr, :EllipseArcRepr), (2, 3, 3))
    eval(quote
            first(r::$T) = first(r.id)
            last(r::$T) = last(r.id)
            getindex(r::$T, i::Int) = r.id[i]
         end)
end

# Convert to primitives  -----------------------------------------------------------------

function primitives(l::Line)
    points = [l.p0, l.p1]
    indices = LineRepr((1, 2))
    (points, indices)
end

function primitives(c::CircleArc)
    points = [c.p0, c.center, c.p1]
    indices = CircleArcRepr((1, 2, 3))
    (points, indices)
end

function primitives(c::EllipseArc)
    points = [c.p0, c.center, c.p1]
    indices = EllipseArcRepr((1, 2, 3))
    (points, indices)
end

function primitives(shape::BoundingBox)
    dx = shape.ur.x - shape.ll.x
    dy = shape.ur.y - shape.ll.y

    points = [shape.ll, shape.ll + Point(dx, 0), shape.ur, shape.ll + Point(0, dy)]
    primitives(ClosedPolygon(points))
end

function primitives(shape::Square)
    points = [shape.ll,
              shape.ll + Point(shape.size, 0),
              shape.ll + Point(shape.size, shape.size),
              shape.ll + Point(0, shape.size)]

    primitives(ClosedPolygon(points))
end

function primitives(shape::Rectangle)
    points = [shape.ll,
              shape.ll + Point(shape.size_x, 0),
              shape.ll + Point(shape.size_x, shape.size_y),
              shape.ll + Point(0, shape.size_y)]

    primitives(ClosedPolygon(points))
end

# Physical points, curves defined in terms of LOCAL point numbering, orientation
function primitives(shape::ClosedPolygon)
    # Point part
    points = shape.points

    npoints = length(points)
    pindices = 1:npoints  # Local
    curves = map(LineRepr, zip([pindices[2:end]..., first(pindices)],
                               [pindices[1:end-1]..., last(pindices)]))

    (points, [(c, 1) for c in curves])
end

function primitives(shape::Circle)
    points = [shape.center,
              shape.center + Point(shape.radius, 0),
              shape.center + Point(0, shape.radius),
              shape.center + Point(-shape.radius, 0),
              shape.center + Point(0, -shape.radius)]

    curves = map(CircleArcRepr, [(2, 1, 3), (3, 1, 4), (4, 1, 5), (5, 1, 2)])

    (points, [(c, 1) for c in curves])
end

function primitives(shape::Ellipse)
    if shape.size_x > shape.size_y
        points = [shape.center,
                  shape.center + Point(shape.size_x, 0),
                  shape.center + Point(0, shape.size_y),
                  shape.center + Point(-shape.size_x, 0),
                  shape.center + Point(0, -shape.size_y)]
    else
        points = [shape.center,
                  shape.center + Point(0, shape.size_y),
                  shape.center + Point(shape.size_x, 0),
                  shape.center + Point(0, -shape.size_y),
                  shape.center + Point(-shape.size_x, 0)]
    end
    curves = map(EllipseArcRepr, [(2, 1, 3), (4, 1, 3), (4, 1, 5), (2, 1, 5)])
    orientation = [1, -1, 1, -1]

    (points, collect(zip(curves, orientation)))
end

function primitives(shape::Loop)
    # NOTE: loop as it comes here has curves defined in terms of points which can be not can
    # coincide. Here we try to keep only the uniqe points - the check for duplicates,
    # hoever, involves only the interior points of curves.

    curves = shape.curves
    orientation = shape.orientation
    @assert first(orientation) == 1

    points, indices = primitives(first(curves))
    # First curve for sure can add all its points
    all_points = Vector{Point}(points)
    # And also its indices
    all_indices = Vector{IndexRepr}([indices])
    # The connection over last point

    link = last(points)
    link_index = last(indices)
    for (curve, sign) in zip(curves[2:end-1], orientation[2:end-1])
        counter = length(all_points)
        # Connection over first
        points, indices = primitives(curve)
        # First, interior, Last
        if sign == 1
            # First, is known
            global_indices = [link_index]
            # Need to check interior points
            for interior_id in interior_indices(indices)
                point = points[interior_id]
                index = findfirst(q -> q == point, all_points)
                # Not seen
                if index == 0
                    push!(all_points, point)

                    counter += 1
                    index = counter
                end
                # Seen with index
                push!(global_indices, index)
            end
            # The last gut is new
            link = last(points)
            push!(all_points, link)
            counter += 1

            link_index = counter
            push!(global_indices, link_index)

            repr = typeof(indices)(tuple(global_indices...))
        else
            # Last, is known so remember to add it later
            # First is new
            link = first(points)
            push!(all_points, link)
            counter +=1

            global_indices = [counter]
            # Need to check interior points
            for interior_id in interior_indices(indices)
                point = points[interior_id]
                index = findfirst(q -> q == point, all_points)
                # Not seen
                if index == 0
                    push!(all_points, point)

                    counter += 1
                    index = counter
                end
                # Seen with index
                push!(global_indices, index)
            end
            # The last point is there already as link_index
            push!(global_indices, link_index)

            repr = typeof(indices)(tuple(global_indices...))
            link_index = first(repr)
        end
        push!(all_indices, repr)
    end

    # The last curve closes the loop so if 1/-1 oriented we have seen its last/first point
    # That point was seen as 1, last is given by link_index
    counter = length(all_points)
    curve, sign = last(curves), last(orientation)
    points, indices = primitives(curve)

    global_indices = []
    fi, li = (sign == 1) ? (link_index, 1) : (1, link_index)

    push!(global_indices, fi)
    # Interior guys
    for interior_id in interior_indices(indices)
        point = points[interior_id]
        index = findfirst(q -> q == point, all_points)
        # Not seen
        if index == 0
            push!(all_points, point)

            counter += 1
            index = counter
        end
        # Seen with index
        push!(global_indices, index)
    end
    push!(global_indices, li)

    repr = typeof(indices)(tuple(global_indices...))
    push!(all_indices, repr)

    (all_points, collect(zip(all_indices, orientation)))
end

# Convert to gmsh ------------------------------------------------------------------------

function mesh(canvas::Canvas, size::Vector, file::AbstractString)
    base, ext = splitext(file)
    for (index, size_i) in enumerate(size)
        file_i = base * "$(index)" * ext
        mesh(canvas, size_i, file_i)
    end
end

function mesh(canvas::Canvas, size::Real, file::AbstractString)
    base, ext = splitext(file)

    geo_file = "$(base).geo"
    if ext == ".msh"
        @assert gmsh_script(canvas, "$(geo_file)", size)
        gmsh = success(`gmsh -2 $(geo_file) -o $(file)`)
        return gmsh && success(`rm $(geo_file)`)
    end
    # XML mesh, mesh itself & physical & facet
    msh_file = "$(base).msh"
    if ext == ".xml"
        @assert mesh(canvas, size, msh_file)
        convert = success(`dolfin-convert $(msh_file) $(file)`)
        return convert && success(`rm $(msh_file)`)
    end
    # HDF5
    xml_file = "$(base).xml"
    if ext == ".h5"
        @assert mesh(canvas, size, xml_file)

const python = "
from dolfin import Mesh, MeshFunction, mpi_comm_world, HDF5File;
import os;

xml_mesh = '$(xml_file)'
base, ext = os.path.splitext(xml_mesh);
assert os.path.exists(xml_mesh);

xml_volumes = '$(base)_physical_region.xml';
assert os.path.exists(xml_volumes);

xml_facets = '$(base)_facet_region.xml';
has_facet_regions = os.path.exists(xml_facets);

mesh = Mesh(xml_mesh);
cell_f = MeshFunction('size_t', mesh, xml_volumes);

h5_file = '$(base).h5';
out = HDF5File(mesh.mpi_comm(), h5_file, 'w');
out.write(mesh, '/mesh');
out.write(cell_f, '/cell_markers');

if has_facet_regions: out.write(MeshFunction('size_t', mesh, xml_facets), '/facet_markers')
#"
        convert = success(`python -c "$(python)"`)
        xmls = [xml_file, "$(base)_physical_region.xml", "$(base)_facet_region.xml"]
        return convert && all(success(`rm $(file)`) for file in filter(isfile, xmls))
    end
end

gmsh_script(canvas::Canvas, size::Real) = write(STDOUT, gmsh_code(canvas, size))

function gmsh_script(canvas::Canvas, file::AbstractString, size::Real)
    @assert last(splitext(file)) == ".geo"
    f = open(file, "w")
    count = write(f, gmsh_code(canvas, size))
    close(f)
    count > 0
end

function gmsh_code(canvas::Canvas, size::Real)
    counter = 1  # Of any GMSH entities

    code = Vector{AbstractString}(["SIZE = $(size);\n"])
    loops = Vector{Int}()  # Each shape is a closed loop, their collection together with the
                           # canvas loop defines the outer domain

    # We first write shapes
    for (tag, shape) in enumerate(canvas.shapes)
        points, curves = primitives(shape)
        # Will need to convert entities from local to global
        local_to_global = collect(counter:counter+length(points)-1)
        # Code for points
        points_gmsh, counter = gmsh_code(points, local_to_global, counter)
        # In case of lines we want a loop as well; for loops and to define physical ...
        curves_gmsh, shape_loop, counter = gmsh_code(curves, local_to_global, counter)
        # Line Loop is another entity
        # NOTE, all stuff below is what gmsh_code of a curve should do .... Anyways, there
        # will be rewrite of this because of 3d
        loop_gmsh = "Line Loop($(counter)) = {$(join(map(x -> "$(x)", shape_loop), ", "))};"
        counter += 1
        # Plane surface in terms of line loop
        surf_gmsh = "Plane Surface($(counter)) = {$(counter-1)};"
        # Physical, not entities
        physurf_gmsh = "Physical Surface($(tag)) = {$(counter)};"
        physline_gmsh = "Physical Line($(tag)) = {$(join(map(x-> "$(x)", shape_loop), ", "))};"

        loops = vcat(loops, counter-1)  # Add to 'inner' boundary
        # Prepare for next round
        counter += 1  # Counter so that it is ready to use
        shape_code = join([points_gmsh,
                           curves_gmsh,
                           loop_gmsh,
                           surf_gmsh,
                           physurf_gmsh,
                           physline_gmsh,
                           "//--------------------\n"], "\n")  # For readability
        push!(code, shape_code)
    end
    # Now the canvas
    tag = length(canvas.shapes) + 1

    points, curves = primitives(canvas.bbox)
    local_to_global = collect(counter:counter+length(points)-1)

    points_gmsh, counter = gmsh_code(points, local_to_global, counter)

    curves_gmsh, shape_loop, counter = gmsh_code(curves, local_to_global, counter)
    loop_gmsh = "Line Loop($(counter)) = {$(join(map(x -> "$(x)", shape_loop), ", "))};"
    counter += 1

    surf_gmsh = "Plane Surface($(counter)) = {$(counter-1), $(join(map(x -> "$(x)", loops), ", "))};"
    physurf_gmsh = "Physical Surface($(tag)) = {$(counter)};"
    physline_gmsh = "Physical Line($(tag)) = {$(join(map(x-> "$(x)", shape_loop), ", "))};"

    canvas_code = join([points_gmsh,
                        curves_gmsh,
                        loop_gmsh,
                        surf_gmsh,
                        physurf_gmsh,
                        physline_gmsh,
                        "//--------------------\n"], "\n")
    push!(code, canvas_code)

    # Done
    join(code, "\n")
end

function gmsh_code(points::Vector{Point}, local_to_global::Vector{Int}, counter::Int)
    # Code for points
    points_gmsh = Vector{AbstractString}()

    for (local_index, point) in enumerate(points)
        str = "Point($(local_to_global[local_index])) = {$(point.x), $(point.y), 0., SIZE};"
        push!(points_gmsh, str)
        counter +=1
    end
    points_gmsh = join(points_gmsh, "\n")

    (points_gmsh, counter)
end

function gmsh_code(line::Tuple{LineRepr, Int}, local_to_global::Vector{Int}, counter::Int)
    line, orientation = line
    v0, v1 = local_to_global[first(line)], local_to_global[last(line)]
    line_gmsh = "Line($(counter)) = {$(v0), $(v1)};"
    loop_piece = orientation*counter

    counter += 1
    (line_gmsh, loop_piece, counter)
end

function gmsh_code(curve::Tuple{CircleArcRepr, Int}, local_to_global::Vector{Int}, counter::Int)
    curve, orientation = curve
    v0, v1 = local_to_global[first(curve)], local_to_global[last(curve)]
    center = local_to_global[curve[2]]
    curve_gmsh = "Circle($(counter)) = {$(v0), $(center), $(v1)};"
    loop_piece = orientation*counter

    counter += 1
    (curve_gmsh, loop_piece, counter)
end

function gmsh_code(curve::Tuple{EllipseArcRepr, Int}, local_to_global::Vector{Int}, counter::Int)
    curve, orientation = curve
    v0, center, v1 = local_to_global[curve[1]], local_to_global[curve[2]], local_to_global[curve[3]]
    curve_gmsh = "Ellipse($(counter)) = {$(v0), $(center), $(v0), $(v1)};"
    loop_piece = orientation*counter

    counter += 1
    (curve_gmsh, loop_piece, counter)
end

function gmsh_code{T<:IndexRepr}(curves::Vector{Tuple{T, Int}}, local_to_global::Vector{Int}, counter::Int)
    curves_gmsh = Vector{AbstractString}()
    shape_loop = Vector{Int}()

    for curve in curves
        curve_gmsh, loop_piece, counter = gmsh_code(curve, local_to_global, counter)
        push!(curves_gmsh, curve_gmsh)
        push!(shape_loop, loop_piece)
    end
    curves_gmsh = join(curves_gmsh, "\n")

    (curves_gmsh, shape_loop, counter)
end

max_stars_repo_name: MiroK/emi-cylinders

max_issues_repo_name: MiroK/emi-cylinders

max_forks_repo_name: MiroK/emi-cylinders


text:
# Faster B-spline manifold

"""
B-spline manifold for lower polynomial degree
TODO: make the field `bsplinespaces` to be conposite type, not abstract type, for performance
"""
struct FastBSplineManifold{T} <: AbstractBSplineManifold
    bsplinespaces::Vector{<:FastBSplineSpace}
    controlpoints::Array{T} where T<:Point
    function FastBSplineManifold(Ps::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{T}) where T
        Ps = FastBSplineSpace.(Ps)
        if collect(size(a)) ≠ dim.(Ps)
            throw(DimensionMismatch())
        else
            P = convert(Vector{FastBSplineSpace}, Ps)
            new{T}(P, float(a))
        end
    end
end

struct BSplineCurve{p1,T} <: AbstractBSplineManifold
    bsplinespace1::FastBSplineSpace{p1}
    controlpoints::Array{T,1} where T<:Point
    function BSplineCurve(P1::AbstractBSplineSpace, a::AbstractArray{T,1}) where T
        p1 = degree(P1)
        if size(a) ≠ (dim(P1),)
            throw(DimensionMismatch())
        else
            new{p1,T}(P1, float(a))
        end
    end
end
function BSplineCurve(P::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{<:Point,1})
    if length(P) ≠ 1
        throw(DimensionMismatch())
    else
        BSplineCurve(P[1], a)
    end
end

struct BSplineSurface{p1,p2,T} <: AbstractBSplineManifold
    bsplinespace1::FastBSplineSpace{p1}
    bsplinespace2::FastBSplineSpace{p2}
    controlpoints::Array{T,2} where T<:Point
    function BSplineSurface(P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, a::AbstractArray{T,2}) where T
        p1, p2 = degree(P1), degree(P2)
        if size(a) ≠ (dim(P1), dim(P2))
            throw(DimensionMismatch())
        else
            new{p1,p2,T}(P1, P2, float(a))
        end
    end
end
function BSplineSurface(P::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{<:Point,2})
    if length(P) ≠ 2
        throw(DimensionMismatch())
    else
        BSplineSurface(P[1], P[2], a)
    end
end

struct BSplineSolid{p1,p2,p3,T} <: AbstractBSplineManifold
    bsplinespace1::FastBSplineSpace{p1}
    bsplinespace2::FastBSplineSpace{p2}
    bsplinespace3::FastBSplineSpace{p3}
    controlpoints::Array{T,3} where T<:Point
    function BSplineSolid(P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, P3::AbstractBSplineSpace, a::AbstractArray{T,3}) where T
        p1, p2, p3 = degree(P1), degree(P2), degree(P3)
        if size(a) ≠ (dim(P1), dim(P2), dim(P3))
            throw(DimensionMismatch())
        else
            new{p1,p2,p3,T}(P1, P2, P3, float(a))
        end
    end
end
function BSplineSolid(P::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{<:Point,3})
    if length(P) ≠ 3
        throw(DimensionMismatch())
    else
        BSplineSolid(P[1], P[2], P[3], a)
    end
end

function BSplineManifold(M::BSplineCurve{p1}) where {p1}
    P1 = M.bsplinespace1
    a = M.controlpoints
    return BSplineManifold([P1], a)
end

function BSplineManifold(M::BSplineSurface{p1,p2}) where {p1} where {p2}
    k1 = M.bsplinespace1
    P2 = M.bsplinespace2
    a = M.controlpoints
    return BSplineManifold([P1, P2], a)
end

function BSplineManifold(M::BSplineSolid{p1,p2,p3}) where {p1} where {p2} where {p3}
    P1 = M.bsplinespace1
    P2 = M.bsplinespace2
    P3 = M.bsplinespace3
    a = M.controlpoints
    return BSplineManifold([P1, P2, P3], a)
end

"""
convert AbstractBSplineManifold to FastBSplineManifold
"""
function FastBSplineManifold(M::AbstractBSplineManifold)
    return FastBSplineManifold(BSplineSpace.(M.bsplinespaces), M.controlpoints)
end

function (M::FastBSplineManifold)(t::AbstractVector{<:Real})
    # TODO: faster
    return BSplineManifold(M)(t)
end

function (M::BSplineCurve{p1})(t1::Real) where {p1}
    P1, = bsplinespaces(M)
    a = controlpoints(M)
    n1 = dim(P1)
    k1 = P1.knots
    j1 = _knotindex(P1, t1)
    b1 = _bsplinebasis(P1,t1,j1)

    f1 = j1-p1
    l1 = j1
    S1 = b1[1]*a[f1]
    for i1 in f1+1:l1
        S1 += b1[i1-f1+1]*a[i1]
    end

    return S1
end

function (M::BSplineSurface{p1,p2})(t1::Real,t2::Real) where {p1} where {p2}
    P1, P2 = bsplinespaces(M)
    a = controlpoints(M)
    n1, n2 = dim(P1), dim(P2)
    k1, k2 = P1.knots, P2.knots
    j1, j2 = _knotindex(P1, t1), _knotindex(P2, t2)
    b1 = _bsplinebasis(P1,t1,j1)
    b2 = _bsplinebasis(P2,t2,j2)

    f1, f2 = j1-p1, j2-p2
    l1, l2 = j1, j2
    S2 = b2[1]*a[f1,f2]
    for i2 in f2+1:l2
        S2 += b2[i2-f2+1]*a[f1,i2]
    end
    S1 = b1[1]*S2
    for i1 in f1+1:l1
        S2 = b2[1]*a[i1,f2]
        for i2 in f2+1:l2
            S2 += b2[i2-f2+1]*a[i1,i2]
        end
        S1 += b1[i1-f1+1]*S2
    end

    return S1
end

function (M::BSplineSolid{p1,p2,p3})(t1::Real,t2::Real,t3::Real) where {p1} where {p2} where {p3}
    P1, P2, P3 = M.bsplinespace1, M.bsplinespace2, M.bsplinespace3
    a = controlpoints(M)
    n1, n2, n3 = dim(P1), dim(P2), dim(P3)
    k1, k2, k3 = P1.knots, P2.knots, P3.knots
    j1, j2, j3 = _knotindex(P1, t1), _knotindex(P2, t2), _knotindex(P3, t3)
    b1 = _bsplinebasis(P1,t1,j1)
    b2 = _bsplinebasis(P2,t2,j2)
    b3 = _bsplinebasis(P3,t3,j3)

    f1, f2, f3 = j1-p1, j2-p2, j3-p3
    l1, l2, l3 = j1, j2, j3
    S3 = b3[1]*a[f1,f2,f3]
    for i3 in f3+1:l3
        S3 += b3[i3-f3+1]*a[f1,f2,i3]
    end
    S2 = b1[1]*b2[1]*S3
    for i2 in f2+1:l2
        S3 = b3[1]*a[f1,i2,f3]
        for i3 in f3+1:l3
            S3 += b3[i3-f3+1]*a[f1,i2,i3]
        end
        S2 += b1[1]*b2[i2-f2+1]*S3
    end
    S1 = S2
    for i1 in f1+1:l1
        S3 = b3[1]*a[i1,f2,f3]
        for i3 in f3+1:l3
            S3 += b3[i3-f3+1]*a[i1,f2,i3]
        end
        S2 = b1[i1-f1+1]*b2[1]*S3
        for i2 in f2+1:l2
            S3 = b3[1]*a[i1,i2,f3]
            for i3 in f3+1:l3
                S3 += b3[i3-f3+1]*a[i1,i2,i3]
            end
            S2 += b1[i1-f1+1]*b2[i2-f2+1]*S3
        end
        S1 += S2
    end

    return S1
end

# function (M::Union{BSplineCurve,BSplineSurface,BSplineSolid})(t::AbstractVector{<:Real})
#     return M(t...)
# end

function (M::BSplineCurve)(t::AbstractVector{<:Real})
    return M(t...)
end

function (M::BSplineSurface)(t::AbstractVector{<:Real})
    return M(t...)
end

function (M::BSplineSolid)(t::AbstractVector{<:Real})
    return M(t...)
end

dim(M::FastBSplineManifold) = length(M.bsplinespaces)
dim(M::BSplineCurve) = 1
dim(M::BSplineSurface) = 2
dim(M::BSplineSolid) = 3

function bsplinespaces(M::FastBSplineManifold)
    return M.bsplinespaces
end

function bsplinespaces(M::BSplineCurve{p1}) where {p1}
    return (M.bsplinespace1, )
end

function bsplinespaces(M::BSplineSurface{p1,p2}) where {p1} where {p2}
    return (M.bsplinespace1, M.bsplinespace2)
end

function bsplinespaces(M::BSplineSolid{p1,p2,p3}) where {p1} where {p2} where {p3}
    return (M.bsplinespace1, M.bsplinespace2, M.bsplinespace3)
end

function controlpoints(M::FastBSplineManifold)
    return M.controlpoints
end

function controlpoints(M::BSplineCurve{p1}) where {p1}
    return M.controlpoints
end

function controlpoints(M::BSplineSurface{p1,p2}) where {p1} where {p2}
    return M.controlpoints
end

function controlpoints(M::BSplineSolid{p1,p2,p3}) where {p1} where {p2} where {p3}
    return M.controlpoints
end


text:
using Pkg

Pkg.develop(PackageSpec(url="https://github.com/carlobaldassi/GaussDCA.jl"))

msa = map(Residue, rand(1:21,100,20))
dca = gaussdca(msa, min_separation=2)

@test  isnan(dca[1,1])
@test  isnan(dca[1,2])
@test !isnan(dca[1,3])


text:
export NeuralPBC, unstack, controller
export SetDistanceLoss, gradient
export CustomDagger

struct NeuralPBC{N,HD}
    Hd::HD
    function NeuralPBC(N::Int, Hd)
        new{N,typeof(Hd)}(Hd)
    end
end

(P::NeuralPBC)(x, ps) = controller(P, x, ps)

function DiffEqFlux.initial_params(P::NeuralPBC{N,HD}) where {N,HD<:FastChain}
    gains = ones(Float32,N)
    vcat(DiffEqFlux.initial_params(P.Hd), gains)
end

function unstack(P::NeuralPBC{N}, ps::AbstractVector) where {N}
    return (ps[1:end-N], ps[end-N+1:end])
end

function controller(P::NeuralPBC{N}, x, ps) where {N}
    Hdθ, gains = unstack(P,ps)
    dot(gains, jacobian(P.Hd, x, Hdθ))
end

abstract type NeuralPBCLoss end

function gradient(l::NeuralPBCLoss, prob::ODEProblem, x0::VX, θ; dt=0.1) where
    {T<:Real, VX<:Union{Vector{T}, SubArray{T}}}
    loss(ps) = l(trajectory(prob, x0, ps; saveat=dt, sensealg=DiffEqFlux.ReverseDiffAdjoint()))
    val, back = Zygote.pullback(loss, θ)
    return first(back(1)), val
end
function gradient(l::NeuralPBCLoss, prob::ODEProblem, batch::VB, θ; dt=0.1) where
    {T<:Real, VX<:Union{Vector{T}, SubArray{T}}, VB<:Vector{VX}}
    batchsize = length(batch)
    gs = Vector{typeof(θ)}(undef, batchsize)
    ls = Vector{eltype(batch[1])}(undef, batchsize)
    Threads.@threads for id in 1:batchsize
        thisgs, thisls = gradient(l, prob, batch[id], θ, dt=dt)
        if !isnothing(thisgs)
            gs[id] = thisgs
        else
            gs[id] = 0θ
        end
        ls[id] = thisls
    end
    (sum(gs)/batchsize, sum(ls)/batchsize)
end
function gradient(ls::Tuple{Vararg{L}}, prob::ODEProblem, x0, θ; dt=0.1) where {L<:NeuralPBCLoss}
    loss(ps) = begin
        x = trajectory(prob, x0, ps; saveat=dt, sensealg=DiffEqFlux.ReverseDiffAdjoint())
        sum(l(x) for l in ls)
    end
    val, back = Zygote.pullback(loss, θ)
    return first(back(1)), val
end

struct SetDistanceLoss{T,F} <: NeuralPBCLoss
    xstar::Vector{T}
    radius::T
    f::F
end
function SetDistanceLoss(f::Function,xstar,r)
    SetDistanceLoss{eltype(xstar),typeof(f)}(xstar,r,f)
end
function (l::SetDistanceLoss)(x::AbstractVector)
    delta = l.f(x)
    return delta < l.radius ? zero(eltype(x)) : delta - l.radius
end
function (l::SetDistanceLoss)(x::AbstractMatrix)
    delta = mapreduce(l.f, min, eachcol(x))
    return delta < l.radius ? zero(eltype(x)) : delta - l.radius
end