text:
module OBJFileIO

export Shapeobject, loadOBJFile, createCylinderVertices

############### Begin Type declaration for Shape Object ################
type Shapeobject
    #A new type declaration. This is a shape object that is composed of
    #An array of vertices and an array of faces
    vertices::Array{Float64, 2}
    faces::Array{Int64, 2}

    #Constructor that creates an object from an obj file
    function Shapeobject(filename::String)
        #Create instance of the object
        shapeOBJ = new()
        #Populate the object fields using the loadOBJFile method
        shapeOBJ.vertices, shapeOBJ.faces = loadOBJFile(filename)
        #Return object
        return shapeOBJ
    end #end constructor

    #Constructor that creates an object from an obj file
    function Shapeobject(cylDiamter::Float64, cylLength::Float64,
                         numVertices::Int64, filename::String)
        #Create instance of the object
        shapeOBJ = new()
        #Populate the object.faces field using the loadOBJFile method
        x, shapeOBJ.faces = loadOBJFile(filename)
        #Populate the object.vertices field using the createCylinderVertices method
        shapeOBJ.vertices = createCylinderVertices(cylDiamter, cylLength, numVertices)
        #Return object
        return shapeOBJ
    end #end constructor

end # end of type declaration

############### Begin Functions ###############

#FUNCTION
#SUMMARY:
#This function loads an obj file and returns an array of vertices and an array of faces
#
#INPUTS:
#    objFile::String -
#    Path to the obj file containing the object that you want
#
#OUTPUTS:
#    Tuple of arrays
#
#    objectVertices::Array{Float64,2} -
#    array containing the vertices of the object
#
#    objectFaces::Array{Float64,2} -
#    array containing the faces of the object
function loadOBJFile(objFile::String)
    #Initiate variables for the object vertices and faces
    objectVertices = zeros(Float64,1,3)
    objectFaces = zeros(Int64,1,3)

    #Initialise variables to temporarily store the data for each line.
    lineVertex = zeros(Float64,1,3)
    lineFace = zeros(Int64,1,3)

    #Initialise counter for the vertex and face line number
    vertexLineCounter = 1
    faceLineCounter = 1

    #Open OBJ file for reading
    objFileContent = open(objFile)

    #Loop through each line of the file and look for vertex lines and face lines
    #Those lines begin with a 'v' or and 'f' respectively.
    #Foor those lines store the values in the repsective objectvertices or object faces
    #array.
    for line in eachline(objFileContent)
        if line[1] == 'v'
            vertexLine = split(line)
            #Initialise counter for the vertex line number
            for i = 1:length(vertexLine)
                if i != 1
                    if vertexLineCounter == 1
                        objectVertices[vertexLineCounter,i-1] = parsefloat(vertexLine[i])
                    else
                        lineVertex[i-1] = parsefloat(vertexLine[i])
                    end
                end
            end
            if vertexLineCounter != 1
                #Concatenate the vertex from the parsed line to the object vertices array
                objectVertices = vcat(objectVertices,lineVertex)
            end
            #Increment counter
            vertexLineCounter += 1

        elseif line[1] == 'f'
            faceLine = split(line)
            for i = 1:length(faceLine)
                if i != 1
                    if faceLineCounter == 1
                        objectFaces[faceLineCounter,i-1] = parseint(faceLine[i])
                    else
                        lineFace[i-1] = parsefloat(faceLine[i])
                    end
                end
            end
            if faceLineCounter != 1
                #Concatenate the face from the parsed line to the object vertices array
                objectFaces = vcat(objectFaces,lineFace)
            end
            #Increment counter
            faceLineCounter += 1
        end
    end

    #Close the file now we are done with it
    close(objFileContent)

    #return Object vertices and faces
    return objectVertices, objectFaces
end #end loadOBJFile function

#FUNCTION
#SUMMARY:
#Takes in parameters and returns a cylinder object from those parameters
#
#INPUTS:
#    cylDiameter::Float64 -
#    Diameter of cylinder
#
#    cylLength::Float64 -
#    Axial length of the cylinder
#
#    numVertices::Int64 -
#    number of vertices around each end of cylinder i.e. this is half of the total number
#    of vertices in the cylinder object
#
#OUTPUTS:
#    vertices::Array{Float64,2}
#    An 2*numVertices x 3 array containing the coordinates of each vertex of the cylinder
function createCylinderVertices(cylDiamter::Float64, cylLength::Float64, numVertices::Int64)
    #Calculate radius
    cylRadius = cylDiamter/2

    #Create xCoordinates
    cylBase = cylLength/2
    xCoordBase = -cylBase
    xCoordTop = cylBase

    #Calculate angular step around circle
    angleToVertex = -2π/numVertices

    #Preallocate vertices array
    vertices = zeros(Float64, 2 * numVertices, 3)

    #Loop through each vertex
    for vertex = 1:numVertices
        #Calculate points around the circle
        yCoord = cylRadius * cos((vertex - 1) * angleToVertex)
        ZCoord = cylRadius * sin((vertex - 1) * angleToVertex)

        #Add points to the vertices array for the Base
        vertices[2 * vertex - 1, 1] = xCoordBase
        vertices[2 * vertex - 1, 2] = yCoord
        vertices[2 * vertex - 1, 3] = ZCoord

        #Add points to the vertices array for the Top
        vertices[2 * vertex, 1] = xCoordTop
        vertices[2 * vertex, 2] = yCoord
        vertices[2 * vertex, 3] = ZCoord
    end

    #return vertices
    return vertices
end #end createCylinderVertices method

end #end module


text:
#=
Figures are supposed to fill the gap that Scenes in combination with Layoutables leave.
A scene is supposed to be a generic canvas on which plot objects can be placed and drawn.
Layoutables always require one specific type of scene, with a PixelCamera, in order to draw
their visual components there.
Figures also have layouts, which scenes do not have.
This is because every figure needs a layout, while not every scene does.
Figures keep track of the Layoutables that are created inside them, which scenes don't.

The idea is there are three types of plotting commands.
They can return either:
    - a FigureAxisPlot (figure, axis and plot)
    - an AxisPlot (axis and plot)
    - or an AbstractPlot (only plot)

This is needed so that users can without much boilerplate create the necessary structures and access them accordingly.

A normal plotting command creates a FigureAxisPlot, which can be splatted into figure, axis and plot. It displays like a figure, so that simple plots show up immediately.

A non-mutating plotting command that references a FigurePosition or a FigureSubposition creates an axis at that position and returns an AxisPlot, which can be splatted into axis and plot.

All mutating plotting commands return AbstractPlots.
They can either reference a FigurePosition or a FigureSubposition, in which case it is looked up
if an axis is placed at that position (if not it errors) or it can reference an axis directly to plot into.
=#

get_scene(fig::Figure) = fig.scene

const _current_figure = Ref{Union{Nothing, Figure}}(nothing)
current_figure() = _current_figure[]
current_figure!(fig) = (_current_figure[] = fig)

current_axis() = current_axis(current_figure())
current_axis(fig::Figure) = fig.current_axis[]
function current_axis!(fig::Figure, ax)
    if ax.parent !== fig
        error("This axis' parent is not the given figure")
    end
    fig.current_axis[] = ax
    ax
end
function current_axis!(fig::Figure, ::Nothing)
    fig.current_axis[] = nothing
end
function current_axis!(ax)
    fig = ax.parent
    if !(fig isa Figure)
        error("Axis parent is not a figure but a $(typeof(ax.parent)). Only axes in figures can have current_axis! called on them.")
    end

    current_axis!(fig, ax)
    # if the current axis is in a different figure, we switch to that as well
    # so that current_axis and current_figure are not out of sync
    current_figure!(fig)
    ax
end

function Figure(; kwargs...)
    scene, layout = layoutscene(; kwargs...)
    Figure(
        scene,
        layout,
        [],
        Attributes(),
        Ref{Any}(nothing)
    )
end

export Figure, current_axis, current_figure, current_axis!, current_figure!

# the FigurePosition is used to plot into a specific part of a figure and create
# an axis there, like `scatter(fig[1, 2], ...)`
struct FigurePosition
    fig::Figure
    gp::GridLayoutBase.GridPosition
end

function Base.getindex(fig::Figure, rows, cols, side = GridLayoutBase.Inner())
    FigurePosition(fig, fig.layout[rows, cols, side])
end

function Base.setindex!(fig::Figure, obj, rows, cols, side = GridLayoutBase.Inner())
    fig.layout[rows, cols, side] = obj
    obj
end

function Base.setindex!(fig::Figure, obj::AbstractArray, rows, cols)
    fig.layout[rows, cols] = obj
    obj
end

Base.lastindex(f::Figure, i) = lastindex(f.layout, i)
Base.lastindex(f::FigurePosition, i) = lastindex(f.fig, i)

# for now just redirect figure display/show to the internal scene
Base.show(io::IO, fig::Figure) = show(io, fig.scene)
Base.show(io::IO, ::MIME"text/plain", fig::Figure) = print(io, "Figure()")
# Base.show(io::IO, ::MIME"image/svg+xml", fig::Figure) = show(io, MIME"image/svg+xml"(), fig.scene)

# a FigureSubposition is just a deeper nested position in a figure's layout, and it doesn't
# necessarily have to refer to an existing layout either, because those can be created
# when it becomes necessary
struct FigureSubposition{T}
    parent::T
    rows
    cols
    side::GridLayoutBase.Side
end

# fp[1, 2] creates a FigureSubposition, a nested version of FigurePosition
# currently, because at the time of creation there is not necessarily a gridlayout
# at all nested locations, the end+1 syntax doesn't work, because it's not known how many
# rows/cols the nested grid has if it doesn't exist yet
function Base.getindex(parent::Union{FigurePosition,FigureSubposition},
        rows, cols, side = GridLayoutBase.Inner())
    FigureSubposition(parent, rows, cols, side)
end

function Base.setindex!(parent::Union{FigurePosition,FigureSubposition}, obj,
    rows, cols, side = GridLayoutBase.Inner())
    layout = get_layout_at!(parent, createmissing = true)
    isnothing(layout) && error("No single GridLayout could be found or created at this FigureSubposition. This means that there are at least two GridLayouts at this position already and it's unclear which one is meant.")
    layout[rows, cols, side] = obj
    obj
end

function Base.setindex!(parent::FigurePosition, obj)
    parent.gp[] = obj
    obj
end

function Base.setindex!(parent::FigureSubposition, obj)
    layout = get_layout_at!(parent.parent, createmissing = true)
    isnothing(layout) && error("No single GridLayout could be found or created at this FigureSubposition. This means that there are at least two GridLayouts at this position already and it's unclear which one is meant.")
    layout[parent.rows, parent.cols, parent.side] = obj
    obj
end

# to power a simple syntax for plotting into nested grids like
# `scatter(fig[1, 1][2, 3], ...)` we need to either find the only gridlayout that
# sits at that position, or we create all the ones that are missing along the way
# as a convenience, so that users don't have to manually create gridlayouts all that often
function get_layout_at!(fp::FigurePosition; createmissing = false)
    c = contents(fp.gp, exact = true)
    layouts = filter(x -> x isa GridLayoutBase.GridLayout, c)
    if isempty(layouts) && createmissing
        fp.gp[] = GridLayoutBase.GridLayout()
    elseif length(layouts) == 1
        first(layouts)::GridLayoutBase.GridLayout
    else
        nothing
    end
end

function get_layout_at!(fsp::FigureSubposition; createmissing = false)
    layout = get_layout_at!(fsp.parent; createmissing = createmissing)

    if isnothing(layout)
        return nothing
    end

    gp = layout[fsp.rows, fsp.cols, fsp.side]

    c = contents(gp, exact = true)
    layouts = filter(x -> x isa GridLayoutBase.GridLayout, c)
    if isempty(layouts) && createmissing
        gp[] = GridLayoutBase.GridLayout()
    elseif length(layouts) == 1
        first(layouts)::GridLayoutBase.GridLayout
    else
        nothing
    end
end

get_figure(fsp::FigureSubposition) = get_figure(fsp.parent)
get_figure(fp::FigurePosition) = fp.fig

function GridLayoutBase.contents(f::FigurePosition; exact = false)
    GridLayoutBase.contents(f.gp, exact = exact)
end

function GridLayoutBase.contents(f::FigureSubposition; exact = false)
    layout = get_layout_at!(f.parent, createmissing = false)
    isnothing(layout) && return []
    GridLayoutBase.contents(layout[f.rows, f.cols, f.side], exact = exact)
end

function GridLayoutBase.content(f::FigurePosition)
    content(f.gp)
end

function GridLayoutBase.content(f::FigureSubposition)
    cs = contents(f, exact = true)
    if length(cs) == 1
        return cs[1]
    else
        error("There is not exactly one object at the given FigureSubposition")
    end
end

ext: jl


text:
using StochasticProcesses

bm = BrownianMotion()
t = linspace(0, 1, 10000)
cumsim(bm, t, 1000)

Profile.clear()

println("cumsim(BrownianMotion) 1000x1000:")
for i in 1:5
  @time cumsim(bm, t, 1000)
end


text:
using Distributions
using LinearAlgebra
using Plots

function gen_uniform_randmat(N,type)
    if type == :complex
        Re = rand(Uniform(-1, 1), N, N)
        Im = rand(Uniform(-1, 1), N, N)im
        A = Re + Im
    elseif type==:real
        A = rand(Uniform(-1, 1), N, N)
    else
        error("type should be :real or :complex")
    end
end

function gen_randnmat(N, type)
    if type == :complex
        randn(Complex{Float64}, N, N)
    elseif type==:real
        randn(N, N)
    else
        error("type should be :real or :complex")
    end
end

function main(N)
    A = gen_uniform_randmat(N, :complex)
    λs = eigvals(A)
    scatter(
        real(λs),
        imag(λs),
        aspect_ratio=:equal
    )
end

main(500)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
function solve()
    n = parse(Int, readline())
    d = Dict{String, Int}()
    for i=1:n
        s = readline()
        if s[1] == '!'
            s = s[2:end]
            v = -1
        else
            v = 1
        end
        if s in keys(d) && d[s] != v
            return s
        else
            d[s] = v
        end
    end
    "satisfiable"
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
function match_begin(line)
    reg = REGEXS[:TEX_BEGIN_REGEX]
    return match(reg, line)
end

function match_end(line)
    reg = REGEXS[:TEX_END_REGEX]
    return match(reg, line)
end


text:
using Dates
using JSON

DFG_VERSION = "0.18.1";

_variableTypeConvert = Dict{Symbol, String}(
    :Point2 => "RoME.Point2",
    :Pose2 => "RoME.Pose2",
    :ContinuousScalar => "IncrementalInference.ContinuousScalar",
    # TBD - https://github.com/JuliaRobotics/Caesar.jl/discussions/807
    :Position1 => "IncrementalInference.ContinuousScalar",
    #TODO deprecate
    :Pose1 => "IncrementalInference.ContinuousScalar"
)

struct Variable
    label::String
    dataEntry::String
    nstime::String
    variableType::String
    dataEntryType::String
    ppeDict::String
    solverDataDict::String
    smallData::String
    solvable::Int
    tags::String
    timestamp:: String
    _version::String
end

struct SolverDataDict
    vecval::Vector{Float64}
    dimval::Int
    vecbw::Vector{Float64}
    dimbw::Int
    BayesNetOutVertIDs::Vector{Int}
    dimIDs::Vector{Int}
    dims::Int
    eliminated::Bool
    BayesNetVertID::String
    separator::Vector{Int}
    variableType::String
    initialized::Bool
    infoPerCoord::Vector{Float64}
    ismargin::Bool
    dontmargin::Bool
    solveInProgress::Int
    solvedCount::Int
    solveKey::String
end
function SolverDataDict(variableType::String, solveKey::String, dims::Int)
    return SolverDataDict(
        zeros(dims*100),
        dims,
        zeros(dims),
        dims,
        [],
        collect(range(0,dims-1, step=1)),
        dims,
        false,
        "_null",
        [],
        variableType,
        false,
        zeros(dims),
        false,
        false,
        0,
        0,
        solveKey)
end

"""
Internal utility function to create the correct solver data (variable data)
given a variable type.
"""
function _getSolverDataDict(variableType::String, solveKey::String)::SolverDataDict
    if variableType == "RoME.Point2"
        return SolverDataDict(variableType, solveKey, 2)
    end
    if variableType == "RoME.Pose2"
        return SolverDataDict(variableType, solveKey, 3)
    end
    if variableType == "IncrementalInference.ContinuousScalar"
        return SolverDataDict(variableType, solveKey, 1)
    end
    throw(error("Variable type '$(variableType)' not supported."))
end

function Variable(label::String, type::Union{String, Symbol}, tags::Vector{String} = ["VARIABLE"], timestamp::String = string(now(Dates.UTC))*"Z")::Variable
    variableType = type isa Symbol ? get(_variableTypeConvert, type, Nothing) : type
    type == Nothing && error("Variable type '$(type) is not supported")

    solverDataDict = Dict("default" => _getSolverDataDict(variableType, "default"))
    result = Variable(
        label,
        "{}",
        "0",
        variableType,
        "{}",
        "{}",
        json(solverDataDict),
        "{}",
        1,
        json(tags),
        timestamp,
        DFG_VERSION
    )
    return result
end


text:
# Auxiliary functions for location processing
function get_ENU(node::Int, model::Simulation)
	"Get ENU coordinates of a specific node in the Simulation"
    model.m.nodes[m.n[node]]
end

function get_coordinates(agent::Commuter, model::Simulation)
	"Get XY coordinates of a specific node in the Simulation"
    pos = get_ENU(agent.pos, model)
    (getX(pos),
     getY(pos))
end

# Auxiliary functions for colors
agentcolor(agent) = agent.infected ? :red : :black
agentsize(agent) = agent.infected ? 6 : 5

# Plot simulation state
function plot_simstate(model::Simulation, background=Plots.current())
	"""
	Description:
		Plots simulation state on specific moment in time (after run).
	Example:
		plotmap(s.m, width=500, height=500)
		plot_simstate(s)
	"""
	# Plots black/red(if infected) agents routes
    for agent in model.agents
        path = [agent.pos]
        while path[end] in keys(agent.path)
            push!(path, agent.path[path[end]])
        end
        length(path) < 2 && continue
        a_x, a_y = get_coordinates(agent, model)
        enu_coords = get_ENU.(path[2:end], Ref(model))

        a_x2, a_y2 = (getX(enu_coords[1]), getY(enu_coords[1]))

        !(a_x == a_x2 && a_y == a_y2) &&  plot!(
                background, [a_x, a_x2], [a_y, a_y2],
                color=agentcolor(agent), arrow =arrow(:open))

        length(enu_coords) > 2 &&  plot!(
                background, getX.(enu_coords),  getY.(enu_coords),
                color=agentcolor(agent), arrow =arrow(:closed))
    end

	# Plots black wagons routes
    for wagon in model.wagons
        path = [wagon.pos]
        while path[end] in keys(wagon.path)
            push!(path, wagon.path[path[end]])
        end
        length(path) < 2 && continue
        a_x, a_y = get_coordinates(wagon, model)
        enu_coords = get_ENU.(path[2:end], Ref(model))

        a_x2, a_y2 = (getX(enu_coords[1]), getY(enu_coords[1]))

        !(a_x == a_x2 && a_y == a_y2) &&  plot!(
                background, [a_x, a_x2], [a_y, a_y2],
                color=:black, arrow =arrow(:open), linewidth = 3)

        length(enu_coords) > 2 &&  plot!(
                background, getX.(enu_coords),  getY.(enu_coords),
                color=:black, arrow =arrow(:closed), linewidth = 3)
    end

	# Define simulation objects appearance
    colors = vcat([agentcolor(agent) for agent in s.agents], [:yellow for wagon in s.wagons])
    sizes = vcat([agentsize(agent) for agent in model.agents], [10 for wagon in model.wagons])
    markers = vcat([:circle for i in 1:length(model.agents)], [:square for i in 1:length(model.wagons)])
    pos_agents = [get_coordinates(agent, model) for agent in model.agents]
    pos_wagons = [get_coordinates(wagon, model) for wagon in model.wagons]
    pos = vcat(pos_agents, pos_wagons)

	# Plot scatter
    scatter!(background,
        pos;
        markercolor = colors,
        markersize = sizes,
        markershapes = markers,
        label = "",
        markerstrokewidth = 0.5,
        markerstrokecolor = :black,
        markeralpha = 0.9
    )

    # Annotate how many agents are inside a TTC car.
    pos_wagons_x = [w[1] for w in pos_wagons]
    pos_wagons_y = [w[2] for w in pos_wagons]
    annot_text = [text(string(length(wagon.passengers)),8,:black) for wagon in model.wagons]
    annotate!(background, pos_wagons_x, pos_wagons_y, annot_text)

end

max_stars_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_issues_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_forks_repo_name: NykPol/EpidemicInUrbanNetworkToronto


text:
# diebold_mariano.jl
# Diebold-Mariano
#
# Copyright (C) 2020   Guilherme Bodin
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

export DieboldMarianoTest

struct DieboldMarianoTest <: TTest
    n::Int         # number of observations
    xbar::Real     # estimated mean
    df::Int        # degrees of freedom
    t::Real        # test statistic
    stderr::Real   # empirical standard error
    μ0::Real       # mean under h_0
end

"""
    DieboldMarianoTest(e1::AbstractVector{<:Real}, e2::AbstractVector{<:Real}; loss=abs2, lookahead=1)

Perform the modified Diebold-Mariano test proposed by Harvey, Leybourne and Newbold of the null
hypothesis that the two methods have the same forecast accuracy. `loss` is the loss function described
in Diebold, F.X. and Mariano, R.S. (1995) Comparing predictive accuracy. Journal of Business and
Economic Statistics, 13, 253-263. and `lookahead` is the number of steps ahead of the forecast.

# References

  * Diebold, F.X. and Mariano, R.S. (1995) Comparing predictive accuracy.
    Journal of Business and Economic Statistics, 13, 253-263.

  * Harvey, D., Leybourne, S., & Newbold, P. (1997). Testing the equality of prediction
    mean squared errors. International Journal of forecasting, 13(2), 281-291.

"""
function DieboldMarianoTest(e1::AbstractVector{<:Real}, e2::AbstractVector{<:Real};
                            loss::Function=abs2, lookahead::Integer=1)
    length(e1) == length(e2) || throw(DimensionMismatch("inputs must have the same length"))
    n = length(e1)
    # Calculate the loss diferential series based on the loss function g
    d = loss.(e1) .- loss.(e2)
    dm_cov = autocov(d, collect(0:lookahead-1))
    dm_var = (dm_cov[1] + 2 * sum(dm_cov[2:end]))/n
    # Statistic from the original Diebold-Mariano test
    xbar_dm = mean(d)
    stderr_dm = sqrt(dm_var)
    statistic_dm = xbar_dm/stderr_dm
    k = sqrt((1 + (1 - 2*lookahead + (lookahead/n)*(lookahead - 1))/n))
    # Statistic from the modified Diebold-Mariano test proposed by Harvey, Leybourne and Newbold
    statistic_hln = statistic_dm * k
    xbar = xbar_dm * k
    stderr = stderr_dm
    return DieboldMarianoTest(n, xbar, n - 1, statistic_hln, stderr, 0.0)
end

testname(::DieboldMarianoTest) = "Diebold-Mariano test"
population_param_of_interest(x::DieboldMarianoTest) = ("Mean", x.μ0, x.xbar)
default_tail(test::DieboldMarianoTest) = :both

function show_params(io::IO, x::DieboldMarianoTest, ident)
    println(io, ident, "number of observations: $(x.n)")
    println(io, ident, "DM statistic:           $(x.t)")
    println(io, ident, "degrees of freedom:     $(x.df)")
end


text:
using Convex
using SDPA_GMP
using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges
using SemidefiniteOptInterface
const SDOI = SemidefiniteOptInterface

using SparseArrays
# using SDPA

using LinearAlgebra
using GenericLinearAlgebra

eye(n) = Matrix(big"1.0"*I, n, n)

using SCS

# # x = [1, 2, 3]
# # P = Variable(3, 3)
# # p1 = Problem{BigFloat}(:minimize, matrixfrac(x, P), [P <= big(2)*eye(3), P >= big(0.5) * eye(3)])

# x = Semidefinite(3)
# p = Problem{BigFloat}(:minimize, sumlargesteigs(x, 2), [x >= big(1)])
# #
# # mock1 = SDPA_GMP.sdpa_gmp_binary_solve(p1)
# #
# Y = Variable(5,5)
# X = -1*rand(BigFloat, 5, 5)
# p2 = Problem{BigFloat}(:minimize, tr(Y), [ diag(Y)[2:5] == diag(X)[2:5], Y[1,1] == big(0.0) ])
# #

# x = Variable(1)
# y = Variable(1)
# p2 = Problem{BigFloat}(:maximize, x + y, [ x <= big"1.0" , y <= big"2.0"])

# mock = Convex.solve!(p2, SDPA_GMP.SDPAGMPoptimizer(BigFloat, verbose = true));

# # mock2 = SDPA_GMP.sdpa_gmp_binary_solve(p2)

E12, E21 = ComplexVariable(2, 2), ComplexVariable(2, 2)
s1, s2 = [big"0.25" big"-0.25"*im; big"0.25"*im big"0.25"], [big"0.5" big"0.0"; big"0.0" big"0.0"]
p = Problem{BigFloat}(:minimize, real(tr(E12 * (s1 + 2 * s2) + E21 * (s2 + 2 * s1))), [E12 ⪰ 0, E21 ⪰ 0, E12 + E21 == Diagonal(ones(2)) ])
solve!(p, SDPA_GMP.Optimizer{BigFloat}());

# #
# # mock3 = SDPA_GMP.sdpa_gmp_binary_solve(p3)
x = Variable(Positive())
y = Semidefinite(3)
p = Problem{BigFloat}(:minimize, nuclearnorm(y), y[2,1]<=4, y[2,2]>=3, y[3,3]<=2)
# solve!(p, SDPA_GMP.SDPAGMPoptimizer(BigFloat, verbose = true));

# x = Variable(4)
# p = Problem{BigFloat}(:minimize, sum(Diagonal(x)), x == [1, 2, 3, 4])
# @test vexity(p) == AffineVexity()
# solve!(p, solver)
# @test p.optval ≈ 10 atol=TOL
# @test all(abs.(evaluate(Diagonal(x)) - Diagonal([1, 2, 3, 4])) .<= TOL)

# x = Variable(4,1)
# p = Problem{BigFloat}(:minimize, dotsort(x, [1,2,3,4]), sum(x) >= 7, x >=0, x<=2, x[4]<=1)

x = Variable(3)
p = Problem{BigFloat}(:minimize, norm_1(x), [-2 <= x, x <= 1])

# x = Variable(4, 4)
# p = Problem{BigFloat}(:minimize, sumlargest(x, 2), sumsmallest(x, 4) >= 1)

# A = Semidefinite(2)
# B = [1 0; 0 0]
# ρ = kron(B, A)
# constraints = [partialtrace(ρ, 1, [2; 2]) == BigFloat[0.09942819 0.29923607; 0.29923607 0.90057181], ρ in :SDP]
# p = Problem{BigFloat}(:minimize, Constant(0), constraints)

a = big"2.214123515"+big"0.8143999"*im
x = ComplexVariable()
objective = norm2(a-x)
c1 = real(x)>=0
p = Problem{BigFloat}(:minimize, objective,c1)

mock = solve!(p, SDPA_GMP.Optimizer{BigFloat}(presolve = false));
solve!(p, ProxSDP.Optimizer(log_verbose = true));
# d = solve!(p, e); # g = SDOI.mockSDoptimizer(Float64);
# f = SCS.Optimizer();

# g = [1.2; 4.3];
# u = Variable(2);
# objective = norm2(u-g);
# c2 = [];
# r = minimize(objective);

n = 5
M = rand(n,n)+ im*rand(n,n)
M = M + M'
 # now M is hermitian
x = ComplexVariable(n,n)
objective = norm2(M - x)
c1 = x in :SDP
p = Problem{Float64}(:minimize, objective, c1)
mmo = solve!(p, SDPA_GMP.Optimizer{Float64}());
c, A, b, cones, var_to_ranges, vartypes, conic_constraints = Convex.conic_problem(p3);


text:
## aggregate!{K,N<:Number}(aggr::Aggregator, key::K, values::Vector{N}) =
##     error("The method is not implemented for this type of aggregator")

## aggregate!{N<:Number}(aggr::Aggregator, values::Vector{N}) =
##     error("The method is not implemented for this type of aggregator")


text:
# Sigma can be used to approximate the volume of an n-dimensional ball
using Sigma
using Base.Test

## Ground Truth
# unit n-sphere volume (n-ball surface area)
S(n) = n == 0 ? 2 : 2π*V(n-1)

# unit n-ball volume
V(n) = n == 0 ? 1 : S(n-1) / n

# ratio of n_ball volume to enclosing unit_sphere
ratio_V(n) = V(n) / 2^n

## Returns the probability that any x,y,z,.. within a unit
## n dim box will be within the enclosed sphere.
## Approximates volume of nball ground-truth
function unit_n_ball(num_dims::Integer)
  ncube = mvuniform(-1.0, 1.0, num_dims)
  sphere = sqrt(sum([x^2 for x in ncube])) < 1.0
  ncube, sphere
end

## Prob that point within n dim box will be element of smaller cube
## of side length 1 10th
function unit_n_box(lenratio::Float64, num_dims::Integer)
  @assert 0 < lenratio <= 1
  bigcube = mvuniform(-1.0, 1.0, num_dims)
  smallcube  = [(x > -lenratio) & (x < lenratio) for x in bigcube]
  smallcube_cond = apply(&, smallcube)
  bigcube, smallcube_cond
end

model, condition = unit_n_ball(10)
rand(model, condition)

model, condition = unit_n_box(0.5, 10)
rand(model, condition)


text:
using OmegaCore, OmegaGrad, Distributions
# using OmegaTestModels

# import ForwardDiff
# using Flux

# const ΩTs = [LinearΩ{Vector{Int}, UnitRange{Int64}, Vector{Float64}},
#              Omega.SimpleΩ{Vector{Int}, Float64},
#              LinearΩ{Vector{Int}, UnitRange{Int64}, TrackedArray{Float64, 1, Array{Float64,1}}}]

# const gradalgs = [Omega.TrackerGrad,
#                   Omega.ForwardDiffGrad,
#                   # Omega.ZygoteGrad
#                   ]

# function testgrad(OT = Omega.SimpleΩ{Vector{Int}, Float64})
#   μ = uniform(0.0, 1.0)
#   x = normal(μ, 1.0) + normal(μ, 1.0)
#   y = logerr(x ==ₛ 1.0)
#   ω = OT()
#   y(ω)
#   Omega.lineargradient(y, ω, Omega.ForwardDiffGrad)
# end

# testgrad(LinearΩ{Vector{Int}, UnitRange{Int64}, Vector{Float64}})

# # Test them all
# istracked(::Type{LinearΩ{I, SEG, T}}) where {I, SEG, T <: TrackedArray} = true
# istracked(ΩT) = false

# function testgrad2(; modelcond, ΩT, gradalg)
#   @show ΩT
#   println()
#   Omega.lineargradient(logerr(modelcond), ΩT(), gradalg)
# end

# # Filter to differentiable models
# models = filter(hascond ∧ isdiff, allmodels)

# # Dont use TrackerGrad on wrong Omegas, etc
# f1(m, ΩT, gradalg) = gradalg == Omega.TrackerGrad ?  istracked(ΩT) : true
# f2(m, ΩT, gradalg) = gradalg == Omega.ForwardDiffGrad ? !istracked(ΩT) : true

# function runtests()
#   for model in models, gradalg in gradalgs, ΩT in ΩTs
#     if f1(model, ΩT, gradalg) && f2(model, ΩT, gradalg)
#       @show model.name
#       @show ΩT
#       @show gradalg
#       testgrad2(; modelcond = model.cond, ΩT = ΩT, gradalg = gradalg)
#     end
#   end
# end

# runtests()

include("zygotetests.jl")

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
include("FunField/DegreeLocalization.jl")
include("FunField/Round2.jl")
include("FunField/Factor.jl")
include("FunField/HessQR.jl")
include("FunField/IntClsZx.jl")


text:
#
# Ask for molecular formula
#

function askmform()
    println("Enter molecular formula:")
    return readline()
end


text:
module TestSkLearn

using Random
using Test
using Lale
using Statistics
using DataFrames: DataFrame

const IRIS = getiris()
const X = IRIS[:,1:3] |> DataFrame
const XC = IRIS[:,1:4] |> DataFrame
const YC = IRIS[:,5] |> Vector
const Y = IRIS[:,4] |> Vector

const classifiers = [
    "LinearSVC","QuadraticDiscriminantAnalysis","MLPClassifier",
    "RandomForestClassifier",
    "SVC","LinearSVC",
    "MLPClassifier",
    "SGDClassifier","KNeighborsClassifier",
    "DecisionTreeClassifier",
    "PassiveAggressiveClassifier","RidgeClassifier",
    "ExtraTreesClassifier","GradientBoostingClassifier",
    "BaggingClassifier","AdaBoostClassifier","GaussianNB","MultinomialNB",
 ]

const regressors = [
    "SVR",
    "Ridge",
    "SGDRegressor",
    "KNeighborsRegressor",
    "DecisionTreeRegressor",
    "RandomForestRegressor",
    "ExtraTreesRegressor",
    "GradientBoostingRegressor",
    "AdaBoostRegressor"
]
function fit_test(learner::String,in::DataFrame,out::Vector)
   _learner=LaleOp(learner)
   fit!(_learner,in,out)
   @test _learner.model != Dict()
   return _learner
end

@testset "lale classifiers" begin
   Random.seed!(123)
   for cl in classifiers
      println(cl)
      fit_test(cl,XC,YC)
   end
end

function fit_transform_reg(model::LaleOp,in::DataFrame,out::Vector)
   @test sum((transform(model,in) .- out).^2)/length(out) < 2.0
end
@testset "lale regressors" begin
   Random.seed!(123)
   for rg in regressors
      println(rg)
      model=fit_test(rg,X,Y)
      fit_transform_reg(model,X,Y)
   end
end

end


text:
function raw_ref(model::BilevelModel, idx::Int)
    if haskey(model.ctr_upper, idx)
        return model.ctr_upper[idx]
    elseif haskey(model.ctr_lower, idx)
        return model.ctr_lower[idx]
    else
        error("Index $(idx) does not belong to BilevelModel")
    end
end
raw_ref(cref::BilevelConstraintRef) = raw_ref(cref.model, cref.index)

function BilevelConstraintRef(model, idx)
    raw = raw_ref(model, idx)
    JuMP.ConstraintRef(model, idx, raw.shape)
end
JuMP.constraint_type(::AbstractBilevelModel) = BilevelConstraintRef
level(cref::BilevelConstraintRef) = cref.model.ctr_level[cref.index]
function JuMP.add_constraint(::BilevelModel, ::JuMP.AbstractConstraint, ::String="")
    error(
        "Can't add constraint directly to the bilevel model `m`, "*
        "attach the constraint to the upper or lower model "*
        "with @constraint(Upper(m), ...) or @constraint(Lower(m), ...)")
end
function JuMP.constraint_object(con_ref::ConstraintRef{BilevelModel, Int})
    raw = raw_ref(con_ref)
    return JuMP.constraint_object(raw)
end
# JuMP.add_constraint(m::UpperModel, c::JuMP.VectorConstraint, name::String="") =
#     error("no vec ctr")
function JuMP.add_constraint(m::InnerBilevelModel, c::Union{JuMP.ScalarConstraint{F,S},JuMP.VectorConstraint{F,S}}, name::String="") where {F,S}
    blm = bilevel_model(m)
    blm.nextconidx += 1
    cref = JuMP.ConstraintRef(blm, blm.nextconidx, JuMP.shape(c))
    func = JuMP.jump_function(c)
    level_func = replace_variables(func, bilevel_model(m), mylevel_model(m), mylevel_var_list(m), level(m))
    level_c = JuMP.build_constraint(error, level_func, c.set)
    level_cref = JuMP.add_constraint(mylevel_model(m), level_c, name)
    blm.ctr_level[cref.index] = level(m)
    mylevel_ctr_list(m)[cref.index] = level_cref
    blm.constraints[cref.index] = c
    blm.ctr_info[cref.index] = empty_info(c)
    JuMP.set_name(cref, name)
    cref
end
function JuMP.delete(m::AbstractBilevelModel, cref::BilevelConstraintRef)
    error("can't delete")
    # m.need_rebuild = true
    # delete!(m.constraints, cref.index)
    # delete!(m.connames, cref.index)
end
JuMP.is_valid(m::BilevelModel, cref::BilevelConstraintRef) = cref.index in keys(m.constraints)
JuMP.is_valid(m::InnerBilevelModel, cref::BilevelConstraintRef) =
    JuMP.is_valid(bilevel_model(m), cref) && level(cref) == level(m)
function JuMP.constraint_object(cref::BilevelConstraintRef, F::Type, S::Type)
    c = cref.model.constraints[cref.index]
    # `TypeError` should be thrown is `F` and `S` are not correct
    # This is needed for the tests in `constraints.jl`
    c.func::F
    c.set::S
    c
end
function empty_info(c::JuMP.ScalarConstraint{F,S}) where {F,S}
    return ConstraintInfo{Float64}()
end
function empty_info(c::JuMP.VectorConstraint{F,S}) where {F,S}
    return ConstraintInfo{Vector{Float64}}(MOI.dimension(c.set))
end

function set_dual_start(cref::BilevelConstraintRef, value::T) where T<:Number
    cref.model.ctr_info[cref.index].start = value
end
function set_dual_start(cref::BilevelConstraintRef, value::T) where T<:Vector{S} where S
    array = cref.model.ctr_info[cref.index].start
    @assert length(array) == length(value)
    copyto!(array, value)
end
function get_dual_start(cref::BilevelConstraintRef)
    cref.model.ctr_info[cref.index].start
end
function set_dual_upper_bound(cref::BilevelConstraintRef, value::T) where T<:Number
    cref.model.ctr_info[cref.index].upper = value
end
function set_dual_upper_bound(cref::BilevelConstraintRef, value::T) where T<:Vector{S} where S
    array = cref.model.ctr_info[cref.index].upper
    @assert length(array) == length(value)
    copyto!(array, value)
end
function get_dual_upper_bound(cref::BilevelConstraintRef)
    cref.model.ctr_info[cref.index].upper
end
function set_dual_lower_bound(cref::BilevelConstraintRef, value::T) where T<:Number
    cref.model.ctr_info[cref.index].lower = value
end
function set_dual_lower_bound(cref::BilevelConstraintRef, value::T) where T<:Vector{S} where S
    array = cref.model.ctr_info[cref.index].lower
    @assert length(array) == length(value)
    copyto!(array, value)
end
function get_dual_lower_bound(cref::BilevelConstraintRef)
    cref.model.ctr_info[cref.index].lower
end

function set_primal_upper_bound_hint(vref::BilevelVariableRef, value::T) where T<:Number
    vref.model.var_info[vref.idx].upper = value
end
function get_primal_upper_bound_hint(vref::BilevelVariableRef)
    vref.model.var_info[vref.idx].upper
end
function set_primal_lower_bound_hint(vref::BilevelVariableRef, value::T) where T<:Number
    vref.model.var_info[vref.idx].lower = value
end
function get_primal_lower_bound_hint(vref::BilevelVariableRef)
    vref.model.var_info[vref.idx].lower
end

function JuMP.value(cref::BilevelConstraintRef; result::Int = 1)
    if level(cref) == BilevelJuMP.LOWER_ONLY
        # Constraint index on the lower model
        con_lower_idx = cref.model.ctr_lower[cref.index].index
        # Single bilevel model constraint associated with the lower level constraint
        con_sblm_idx = cref.model.lower_to_sblm[con_lower_idx]
    else
        # Constraint index on the lower model
        con_upper_idx = cref.model.ctr_upper[cref.index].index
        # Single bilevel model constraint associated with the lower level constraint
        con_sblm_idx = cref.model.upper_to_sblm[con_upper_idx]
    end
    # Solver constraint associated with the single bilevel model constraint
    con_solver_idx = cref.model.sblm_to_solver[con_sblm_idx]
    return MOI.get(cref.model.solver, MOI.ConstraintPrimal(result), con_solver_idx)
end
# variables again (duals)
# code for using dual variables associated with lower level constraints
# in the upper level

function JuMP.num_constraints(model::BilevelModel)
    return length(model.ctr_level)
end
function JuMP.num_constraints(model::LowerModel)
    return length(model.m.ctr_lower)
end
function JuMP.num_constraints(model::UpperModel)
    return length(model.m.ctr_upper)
end

struct DualOf
    ci::BilevelConstraintRef
end
struct DualVariableInfo
    info::JuMP.VariableInfo
    ci::BilevelConstraintRef
end
function JuMP.build_variable(
    _error::Function,
    info::JuMP.VariableInfo,
    dual_of::DualOf;
    extra_kw_args...,
)

    if level(dual_of.ci) != LOWER_ONLY
        error("Variables can only be tied to LOWER level constraints, got $(dual_of.ci.level) level")
    end
    for (kwarg, _) in extra_kw_args
        _error("Unrecognized keyword argument $kwarg")
    end

    if info.has_lb
        set_dual_lower_bound(dual_of.ci, info.lower_bound)
        # info.has_lb = false
        # info.lower_bound = NaN
    end
    if info.has_ub
        set_dual_upper_bound(dual_of.ci, info.upper_bound)
        # info.has_ub = false
        # info.upper_bound = NaN
    end
    info.has_fix   && _error("Dual variable does not support fixing")
    if info.has_start
        set_dual_start(dual_of.ci, info.start)
        # info.has_start = false
        # info.start = NaN
    end
    info.binary    && _error("Dual variable cannot be binary")
    info.integer   && _error("Dual variable cannot be integer")

    info = JuMP.VariableInfo(false, NaN, false, NaN,
                             false, NaN, false, NaN,
                             false, false)

    return DualVariableInfo(
        info,
        dual_of.ci
    )
end
function JuMP.add_variable(inner::UpperModel, dual_info::DualVariableInfo, name::String="")
    # TODO vector version
    m = bilevel_model(inner)
    m.nextvaridx += 1
    vref = BilevelVariableRef(m, m.nextvaridx, DUAL_OF_LOWER)
    v_upper = JuMP.add_variable(m.upper, JuMP.ScalarVariable(dual_info.info), name)
    m.var_upper[vref.idx] = v_upper
    m.var_level[vref.idx] = DUAL_OF_LOWER
    m.upper_var_to_lower_ctr_link[v_upper] = m.ctr_lower[dual_info.ci.index] # TODO improve this
    m.variables[vref.idx] = JuMP.ScalarVariable(dual_info.info)
    JuMP.set_name(vref, name)
    m.var_info[vref.idx] = empty_info(dual_info)
    vref
end
function empty_info(::DualVariableInfo)
    return VariableInfo()
end

function get_constrain_ref(vref::BilevelVariableRef)
    model = vref.model
    ctr_ref = model.upper_var_to_lower_ctr_link[model.var_upper[vref.idx]]
    idx = -1
    for (ind, ref) in model.ctr_lower
        if ref == ctr_ref
            idx = ind
        end
    end
    @assert idx != -1
    return BilevelConstraintRef(model, idx)
end

function JuMP.dual(cref::BilevelConstraintRef)
    # Right now this code assumes there is no possibility for vectorized constraints
    if level(cref) == BilevelJuMP.LOWER_ONLY
        # Constraint index on the lower model
        con_lower_ref = cref.model.ctr_lower[cref.index]
        con_lower_idx = con_lower_ref.index
        # Dual variable associated with constraint index
        model_var_idxs = cref.model.lower_primal_dual_map.primal_con_dual_var[con_lower_idx]
        # Single bilevel model variable associated with the dual variable
        sblm_var_idxs = MOI.VariableIndex[]
        for vi in model_var_idxs
            push!(sblm_var_idxs, cref.model.lower_dual_to_sblm[vi])
        end
        # Solver variable associated withe the sblm model
        solver_var_idxs = MOI.VariableIndex[]
        for vi in sblm_var_idxs
            push!(solver_var_idxs, cref.model.sblm_to_solver[vi])
        end
        pre_duals = MOI.get(cref.model.solver, MOI.VariablePrimal(), solver_var_idxs)
        return JuMP.reshape_vector(
            pre_duals,
            JuMP.dual_shape(con_lower_ref.shape)
            )
    elseif level(cref) == BilevelJuMP.UPPER_ONLY
        m = cref.model
        con_upper_ref = cref.model.ctr_upper[cref.index]
        solver_ctr_idx = m.sblm_to_solver[m.upper_to_sblm[JuMP.index(con_upper_ref)]]
        pre_duals = MOI.get(cref.model.solver, MOI.ConstraintDual(), solver_ctr_idx)
        return JuMP.reshape_vector(
            pre_duals,
            JuMP.dual_shape(con_upper_ref.shape)
            )
    else
        error("Dual solutions of upper level constraints are not available. Either the solution method does nto porvide duals or or the solver failed to get one.")
    end
end

function JuMP.normalized_rhs(cref::BilevelConstraintRef)
    return JuMP.normalized_rhs(raw_ref(cref))
end


text:
unique_frame_name(frame::CartesianFrame3D) =
    Symbol("$(frame)_(#$(frame.id))")

function setgeometry!(vis::Visualizer, mechanism::Mechanism, source::AbstractGeometrySource=Skeleton(false, true))
    elements = visual_elements(mechanism, source)
    setgeometry!(vis, mechanism, elements)
end

function setgeometry!(vis::Visualizer, mechanism::Mechanism, elements::AbstractVector{<:VisualElement})
    batch(vis) do v
        delete!(v)
        addgeometry!(vis, mechanism, elements)
    end
    nothing
end

apply_scaling(geometry::HyperSphere, scale::Vec) = HyperEllipsoid(center(geometry), scale)
apply_scaling(geometry, scale) = error("Geometry $geometry has a non-unit scale: $scale, but I don't know how to transform it into a scaled geometry automatically.")

"""
The remote tree viewer interface has no notion of scale
(only a quaternion and translation), so we have to do some
monkey business to push the scaling into the geometry itself.

Currently, this only works for converting spheres into ellipsoids
and assumes the scaling is positive.
"""
function remove_scaling(geometry, transform)
    H = [transform_deriv(transform, Vec(0., 0, 0)) transform(Vec(0., 0, 0));
    Vec(0, 0, 0, 1)']
    scale = Vec(norm(H[1:3, 1]), norm(H[1:3, 2]), norm(H[1:3, 3]))
    if norm(scale - Vec(1, 1, 1)) > 1e-3
        H[:,1] ./= scale[1]
        H[:,2] ./= scale[2]
        H[:,3] ./= scale[3]
        tform = AffineMap(H[1:3, 1:3], transform(Vec(0., 0, 0)))
        return apply_scaling(geometry, scale), tform
    else
        return geometry, transform
    end
end

function addgeometry!(vis::Visualizer, mechanism::Mechanism, elements::AbstractVector{<:VisualElement})
    frame_to_elements = Dict{CartesianFrame3D, Vector{VisualElement}}()
    for element in elements
        push!(get!(() -> VisualElement[], frame_to_elements, element.frame), element)
    end
    body_names = Set{Symbol}()
    batch(vis) do v
        for body in bodies(mechanism)
            body_name = Symbol(body)
            if body_name in body_names
                error("Duplicate body name in mechanism")
            end
            push!(body_names, body_name)
            bodyvis = v[body_name]
            found_geometry = false
            for transform in rbd.frame_definitions(body)
                frame = transform.from
                framename = unique_frame_name(frame)
                if haskey(frame_to_elements, frame)
                    found_geometry = true
                    for element in frame_to_elements[frame]
                        geometry, transform = remove_scaling(element.geometry, element.transform)
                        geomdata = GeometryData(geometry, element.color, transform)
                        addgeometry!(bodyvis[framename], geomdata)
                    end
                    settransform!(bodyvis[framename], rbd.frame_definition(body, frame))
                end
            end
            if !found_geometry
                setgeometry!(bodyvis, [])
            end
        end
    end
    nothing
end

function addgeometry!(vis::Visualizer, mechanism::Mechanism, frame::CartesianFrame3D, geometry::Union{AbstractGeometry, AbstractMesh})
    addgeometry!(vis, mechanism, [VisualElement(frame, geometry, RGBA{Float32}(1, 0, 0, 0.5), IdentityTransformation())])
end

function addgeometry!(vis::Visualizer, mechanism::Mechanism, frame::CartesianFrame3D; scale=1.0)
    addgeometry!(vis, mechanism, frame, Triad(scale, false))
end

function addgeometry!(vis::Visualizer, mechanism::Mechanism, point::Point3D; radius=0.03)
    addgeometry!(vis, mechanism, point.frame, HyperSphere(Point{3, Float64}(point.v), radius))
end

function addgeometry!(vis::Visualizer, mechanism::Mechanism, frame::CartesianFrame3D, geomdata::GeometryData)
    addgeometry!(vis, mechanism, [VisualElement(frame, geomdata.geometry, convert(RGBA{Float32}, geomdata.color), geomdata.transform)])
end

settransform!(vis::Visualizer, tform::rbd.Transform3D) = settransform!(vis, convert(AffineMap, tform))

function settransform!(vis::Visualizer, state::MechanismState)
    batch(vis) do v
        for body in bodies(state.mechanism)
            body_name = Symbol(body)
            if body_name in keys(v.core.tree[v.path].children)
                settransform!(v[body_name], transform_to_root(state, body))
            end
        end
    end
end

"""
Construct a DrakeVisualizer.Visualizer for the given mechanism by constructing
simple geometries just from the structure of the kinematic tree. If
`show_intertias` is true, then also construct equivalent inertial ellipsoids
for every link.
"""
function Visualizer(mechanism::Mechanism, prefix=[:robot1];
                    show_inertias::Bool=false, randomize_colors::Bool=true)
    source = Skeleton(inertias=show_inertias, randomize_colors=randomize_colors)
    Visualizer(mechanism, source, prefix)
end

function Visualizer(mechanism::Mechanism, elements::AbstractVector{<:VisualElement},
                    prefix=[:robot1])
    vis = Visualizer()[prefix]
    setgeometry!(vis, mechanism, elements)
    vis
end

function Visualizer(mechanism::Mechanism, source::AbstractGeometrySource, prefix=[:robot1])
    elements = visual_elements(mechanism, source)
    vis = Visualizer()[prefix]
    setgeometry!(vis, mechanism, elements)
    vis
end

convert(::Type{AffineMap}, T::Transform3D) =
    AffineMap(rotation(T), translation(T))

inspect!(state::MechanismState, vis::Visualizer) = manipulate!(state) do state
    settransform!(vis, state)
end

function inspect!(state::MechanismState;
        show_inertias::Bool=false, randomize_colors::Bool=true)
    vis = Visualizer()[:robot1]
    elements = visual_elements(state.mechanism, Skeleton(show_inertias, randomize_colors))
    setgeometry!(vis, state.mechanism, elements)
    inspect!(state, vis)
end

inspect(mechanism::Mechanism, args...; kwargs...) =
    inspect!(MechanismState{Float64}(mechanism), args...; kwargs...)


text: include("../src/StructuralInheritance.jl")

max_stars_repo_name: WschW/StrucralInheritence

max_issues_repo_name: WschW/StrucralInheritence

max_forks_repo_name: WschW/StrucralInheritence


text:
function updateModel(vector)
    s = vector[1];
    i = vector[2];
    r = vector[3];

    updated_s = s - lambda*s*i*dt
    updated_i = i + lambda*s*i*dt - gam*i*dt
    updated_r = r + gam*i*dt

    return [updated_s updated_i updated_r]
end

# initial setup

lambda = 1/20000 # infection rate parameter
gam = 1/10       # recovery rate parameter
dt = 0.5         # length of time step in days (two days)
tfinal = 610;    # number of days
s0 = 10000.0     # initial S
i0 = 4.          # initial I
r0 = 0.          # initial R

nsteps = round(Int64, tfinal/dt)
resultvals = Array{Float64}(undef, nsteps+1, 3)
timevec = Array{Float64}(undef, nsteps+1)
resultvals[1,:] = [s0, i0, r0]
timevec[1] = 0.

# loop!

for step  = 1:nsteps
    resultvals[step+1, :] = updateModel(resultvals[step, :])
    timevec[step+1] = timevec[step] + dt
end

using Plots
gr()

plot(timevec, resultvals,  # we should of course at a minimum provide some labels
title  = "SIR Simulation",
xlabel = "Day",
ylabel = "Population",
label  = ["Susceptibles" "Infecteds" "Removeds"]
)

savefig("Assets/SIR.png")


text:
function build_so(i::String)
    # Remove any existing shared library
    run(`rm -f binary/lib$(i).so`)

    # Build a shared library
    run(`gcc -Wall -Werror -fpic -shared -o binary/lib$(i).so c/$(i).c`)

    return
end


text:
# conditional loop
i = 0
while i <= 5
    print(i)
    print(",")
    global i += 1
end
print(i)

max_stars_repo_name: elucian/julia

max_issues_repo_name: elucian/julia

max_forks_repo_name: elucian/julia


text:
using LASindex
using FileIO
using RegionTrees
using StaticArrays
using Test

workdir = dirname(@__FILE__)

# source: http://www.liblas.org/samples/ and lasindex
index = joinpath(workdir, "..", "sample", "Palm Beach Post Hurricane.lax")

# Header check
qt = load(index)

# Determine correct quadtree levels
@test LASindex.zlevels(0) == [0]
@test LASindex.zlevels(1) == [1]
@test LASindex.zlevels(4) == [4]
@test LASindex.zlevels(5) == [1,1]
@test LASindex.zlevels(20) == [4,4]
@test LASindex.zlevels(21) == [1,1,1]
@test LASindex.zlevels(22) == [2,1,1]
@test LASindex.zlevels(84) == [4,4,4]
@test LASindex.zlevels(85) == [1,1,1,1]

# Single integer indexing on RegionTrees cells
# based on Z order
root = Cell(SVector(0., 0), SVector(1., 1))
split!(root)
@test root[1] == root[1,1]
@test root[2] == root[2,1]
@test root[3] == root[1,2]
@test root[4] == root[2,2]
@test root[0] == root
@test_throws ErrorException root[5]

# Split cell based on given index level
# until cell with index is found, return it
root = Cell(SVector(0., 0), SVector(1., 1))
@test LASindex.quadtree!(root, 1).boundary.origin ≈ Cell(SVector(0., 0), SVector(0.5, 0.5)).boundary.origin
@test LASindex.quadtree!(root, 4).boundary.origin ≈ Cell(SVector(0.5, 0.5), SVector(0.5, 0.5)).boundary.origin
@test LASindex.quadtree!(root, 5).boundary.origin ≈ Cell(SVector(0., 0), SVector(0.25, 0.25)).boundary.origin
@test LASindex.quadtree!(root, 20).boundary.origin ≈ Cell(SVector(0.75, 0.75), SVector(0.25, 0.25)).boundary.origin
@test LASindex.quadtree!(root, 1).boundary.widths ≈ Cell(SVector(0., 0), SVector(0.5, 0.5)).boundary.widths
@test LASindex.quadtree!(root, 4).boundary.widths ≈ Cell(SVector(0.5, 0.5), SVector(0.5, 0.5)).boundary.widths
@test LASindex.quadtree!(root, 5).boundary.widths ≈ Cell(SVector(0., 0), SVector(0.25, 0.25)).boundary.widths
@test LASindex.quadtree!(root, 20).boundary.widths ≈ Cell(SVector(0.75, 0.75), SVector(0.25, 0.25)).boundary.widths

# Adding data to quadtree cell based on index
root = Cell(SVector(0., 0), SVector(1., 1), Vector{UnitRange{Int}}())
LASindex.quadtree!(root, 1, [2:5])
@test root[1,1].data == [2:5]

# Merge UnitRanges that do have overlap
r = Vector{UnitRange{Integer}}([1:5, 2:4, 12:13, 14:16])
@test LASindex.merge(r) == Vector{UnitRange{Integer}}([1:5, 12:16])
# and again to make sure it's not mutating
@test LASindex.merge(r) == Vector{UnitRange{Integer}}([1:5, 12:16])
r = Vector{UnitRange{Integer}}([1:5, 2:6, 12:13, 14:16, -1:2])
@test LASindex.merge(r) == Vector{UnitRange{Integer}}([-1:6, 12:16])
r = Vector{UnitRange{Integer}}()
@test LASindex.merge(r) == Vector{UnitRange{Integer}}()
r = Vector{UnitRange{Integer}}([1:5])
@test LASindex.merge(r) == Vector{UnitRange{Integer}}([1:5])

# Test spatial intersection
bbox = SVector(955000.0, 955002.0, 885000.0, 988800.0)  # small intersection
@test LASindex.intersect(qt, bbox) == Vector{UnitRange{Integer}}()
bbox = SVector(955000.0, 985002.0, 885000.0, 988800.0)  # fully overlaps
@test LASindex.intersect(qt, bbox) == Vector{UnitRange{Integer}}([1:1924631])


text:
# This file is mostly generated by `generate_builtins.jl` with some additional @static checks
# for recently added builtins.

function getargs(args, frame)
    nargs = length(args)-1  # skip f
    callargs = resize!(frame.framedata.callargs, nargs)
    for i = 1:nargs
        callargs[i] = @lookup(frame, args[i+1])
    end
    return callargs
end

const kwinvoke_name = isdefined(Core, Symbol("#kw##invoke")) ? Symbol("#kw##invoke") : Symbol("##invoke")
const kwinvoke_instance = getfield(Core, kwinvoke_name).instance

"""
    ret = maybe_evaluate_builtin(frame, call_expr, expand::Bool)

If `call_expr` is to a builtin function, evaluate it, returning the result inside a `Some` wrapper.
Otherwise, return `call_expr`.

If `expand` is true, `Core._apply` calls will be resolved as a call to the applied function.
"""
function maybe_evaluate_builtin(frame, call_expr, expand::Bool)
    # By having each call appearing statically in the "switch" block below,
    # each gets call-site optimized.
    args = call_expr.args
    nargs = length(args) - 1
    fex = args[1]
    if isa(fex, QuoteNode)
        f = fex.value
    else
        f = @lookup(frame, fex)
    end
    # Builtins and intrinsics have empty method tables. We can circumvent
    # a long "switch" check by looking for this.
    mt = typeof(f).name.mt
    if isa(mt, Core.MethodTable)
        isempty(mt) || return call_expr
    end
    # Builtins
    if f === <:
        if nargs == 2
            return Some{Any}(<:(@lookup(frame, args[2]), @lookup(frame, args[3])))
        else
            return Some{Any}(<:(getargs(args, frame)...))
        end
    elseif f === ===
        if nargs == 2
            return Some{Any}(===(@lookup(frame, args[2]), @lookup(frame, args[3])))
        else
            return Some{Any}(===(getargs(args, frame)...))
        end
    elseif f === Core._apply
        argswrapped = getargs(args, frame)
        if !expand
            return Some{Any}(Core._apply(argswrapped...))
        end
        new_expr = Expr(:call, argswrapped[1])
        popfirst!(argswrapped)
        argsflat = append_any(argswrapped...)
        for x in argsflat
            push!(new_expr.args, (isa(x, Symbol) || isa(x, Expr) || isa(x, QuoteNode)) ? QuoteNode(x) : x)
        end
        return new_expr
    elseif f === Core._apply_latest
        argswrapped = getargs(args, frame)
        if !expand
            return Some{Any}(Core._apply_latest(argswrapped...))
        end
        new_expr = Expr(:call, argswrapped[1])
        popfirst!(argswrapped)
        argsflat = append_any(argswrapped...)
        for x in argsflat
            push!(new_expr.args, (isa(x, Symbol) || isa(x, Expr) || isa(x, QuoteNode)) ? QuoteNode(x) : x)
        end
        return new_expr
    elseif @static isdefined(Core, :_apply_iterate) ? f === Core._apply_iterate : false
        argswrapped = getargs(args, frame)
        if !expand
            return Some{Any}(Core._apply_iterate(argswrapped...))
        end
        aw1 = argswrapped[1]::Function
        @assert aw1 === Core.iterate || aw1 === Core.Compiler.iterate || aw1 === Base.iterate "cannot handle `_apply_iterate` with non iterate as first argument, got $(aw1), $(typeof(aw1))"
        new_expr = Expr(:call, argswrapped[2])
        popfirst!(argswrapped) # pop the iterate
        popfirst!(argswrapped) # pop the function
        argsflat = append_any(argswrapped...)
        for x in argsflat
            push!(new_expr.args, (isa(x, Symbol) || isa(x, Expr) || isa(x, QuoteNode)) ? QuoteNode(x) : x)
        end
        return new_expr
    elseif f === Core._apply_pure
        return Some{Any}(Core._apply_pure(getargs(args, frame)...))
    elseif f === Core._expr
        return Some{Any}(Core._expr(getargs(args, frame)...))
    elseif @static isdefined(Core, :_typevar) ? f === Core._typevar : false
        if nargs == 3
            return Some{Any}(Core._typevar(@lookup(frame, args[2]), @lookup(frame, args[3]), @lookup(frame, args[4])))
        else
            return Some{Any}(Core._typevar(getargs(args, frame)...))
        end
    elseif f === Core.apply_type
        return Some{Any}(Core.apply_type(getargs(args, frame)...))
    elseif f === Core.arrayref
        if nargs == 0
            return Some{Any}(Core.arrayref())
        elseif nargs == 1
            return Some{Any}(Core.arrayref(@lookup(frame, args[2])))
        elseif nargs == 2
            return Some{Any}(Core.arrayref(@lookup(frame, args[2]), @lookup(frame, args[3])))
        elseif nargs == 3
            return Some{Any}(Core.arrayref(@lookup(frame, args[2]), @lookup(frame, args[3]), @lookup(frame, args[4])))
        elseif nargs == 4
            return Some{Any}(Core.arrayref(@lookup(frame, args[2]), @lookup(frame, args[3]), @lookup(frame, args[4]), @lookup(frame, args[5])))
        elseif nargs == 5
            return Some{Any}(Core.arrayref(@lookup(frame, args[2]), @lookup(frame, args[3]), @lookup(frame, args[4]), @lookup(frame, args[5]), @lookup(frame, args[6])))
        else
            return Some{Any}(Core.arrayref(getargs(args, frame)...))
        end
    elseif f === Core.arrayset
        if nargs == 0
            return Some{Any}(Core.arrayset())
        elseif nargs == 1
            return Some{Any}(Core.arrayset(@lookup(frame, args[2])))
        elseif nargs == 2
            return Some{Any}(Core.arrayset(@lookup(frame, args[2]), @lookup(frame, args[3])))
        elseif nargs == 3
            return Some{Any}(Core.arrayset(@lookup(frame, args[2]), @lookup(frame, args[3]), @lookup(frame, args[4])))
        elseif nargs == 4
            return Some{Any}(Core.arrayset(@lookup(frame, args[2]), @lookup(frame, args[3]), @lookup(frame, args[4]), @lookup(frame, args[5])))
        elseif nargs == 5
            return Some{Any}(Core.arrayset(@lookup(frame, args[2]), @lookup(frame, args[3]), @lookup(frame, args[4]), @lookup(frame, args[5]), @lookup(frame, args[6])))
        elseif nargs == 6
            return Some{Any}(Core.arrayset(@lookup(frame, args[2]), @lookup(frame, args[3]), @lookup(frame, args[4]), @lookup(frame, args[5]), @lookup(frame, args[6]), @lookup(frame, args[7])))
        else
            return Some{Any}(Core.arrayset(getargs(args, frame)...))
        end
    elseif f === Core.arraysize
        if nargs == 2
            return Some{Any}(Core.arraysize(@lookup(frame, args[2]), @lookup(frame, args[3])))
        else
            return Some{Any}(Core.arraysize(getargs(args, frame)...))
        end
    elseif @static isdefined(Core, :const_arrayref) ? f === Core.const_arrayref : false
        return Some{Any}(Core.const_arrayref(getargs(args, frame)...))
    elseif f === Core.sizeof
        if nargs == 1
            return Some{Any}(Core.sizeof(@lookup(frame, args[2])))
        else
            return Some{Any}(Core.sizeof(getargs(args, frame)...))
        end
    elseif f === Core.svec
        return Some{Any}(Core.svec(getargs(args, frame)...))
    elseif f === applicable
        return Some{Any}(applicable(getargs(args, frame)...))
    elseif f === fieldtype
        if nargs == 2
            return Some{Any}(fieldtype(@lookup(frame, args[2]), @lookup(frame, args[3])))
        elseif nargs == 3
            return Some{Any}(fieldtype(@lookup(frame, args[2]), @lookup(frame, args[3]), @lookup(frame, args[4])))
        else
            return Some{Any}(fieldtype(getargs(args, frame)...))
        end
    elseif f === getfield
        if nargs == 2
            return Some{Any}(getfield(@lookup(frame, args[2]), @lookup(frame, args[3])))
        elseif nargs == 3
            return Some{Any}(getfield(@lookup(frame, args[2]), @lookup(frame, args[3]), @lookup(frame, args[4])))
        else
            return Some{Any}(getfield(getargs(args, frame)...))
        end
    elseif f === ifelse
        if nargs == 3
            return Some{Any}(ifelse(@lookup(frame, args[2]), @lookup(frame, args[3]), @lookup(frame, args[4])))
        else
            return Some{Any}(ifelse(getargs(args, frame)...))
        end
    elseif f === invoke
            if !expand
                argswrapped = getargs(args, frame)
                return Some{Any}(invoke(argswrapped...))
            end
            return Expr(:call, invoke, args[2:end]...)
    elseif f === isa
        if nargs == 2
            return Some{Any}(isa(@lookup(frame, args[2]), @lookup(frame, args[3])))
        else
            return Some{Any}(isa(getargs(args, frame)...))
        end
    elseif f === isdefined
        if nargs == 1
            return Some{Any}(isdefined(@lookup(frame, args[2])))
        elseif nargs == 2
            return Some{Any}(isdefined(@lookup(frame, args[2]), @lookup(frame, args[3])))
        else
            return Some{Any}(isdefined(getargs(args, frame)...))
        end
    elseif f === nfields
        if nargs == 1
            return Some{Any}(nfields(@lookup(frame, args[2])))
        else
            return Some{Any}(nfields(getargs(args, frame)...))
        end
    elseif f === setfield!
        if nargs == 3
            return Some{Any}(setfield!(@lookup(frame, args[2]), @lookup(frame, args[3]), @lookup(frame, args[4])))
        else
            return Some{Any}(setfield!(getargs(args, frame)...))
        end
    elseif f === throw
        if nargs == 1
            return Some{Any}(throw(@lookup(frame, args[2])))
        else
            return Some{Any}(throw(getargs(args, frame)...))
        end
    elseif f === tuple
        return Some{Any}(ntupleany(i->@lookup(frame, args[i+1]), length(args)-1))
    elseif f === typeassert
        if nargs == 2
            return Some{Any}(typeassert(@lookup(frame, args[2]), @lookup(frame, args[3])))
        else
            return Some{Any}(typeassert(getargs(args, frame)...))
        end
    elseif f === typeof
        if nargs == 1
            return Some{Any}(typeof(@lookup(frame, args[2])))
        else
            return Some{Any}(typeof(getargs(args, frame)...))
        end
    # Intrinsics
    elseif f === Base.cglobal
        if nargs == 1
            return Some{Any}(Core.eval(moduleof(frame), call_expr))
        elseif nargs == 2
            call_expr = copy(call_expr)
            call_expr.args[3] = @lookup(frame, args[3])
            return Some{Any}(Core.eval(moduleof(frame), call_expr))
        end
    elseif f === Base.llvmcall
        return Some{Any}(Base.llvmcall(getargs(args, frame)...))
    end
    if isa(f, Core.IntrinsicFunction)
        cargs = getargs(args, frame)
        return Some{Any}(ccall(:jl_f_intrinsic_call, Any, (Any, Ptr{Any}, UInt32), f, cargs, length(cargs)))
    end
    if isa(f, getfield(Core, kwinvoke_name))
        return Some{Any}(kwinvoke_instance(getargs(args, frame)...))
    end
    return call_expr
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    VectorizeBridge{T,F,S,G} <: Bridges.Constraint.AbstractBridge

`VectorizeBridge` implements the following reformulations:

  * ``g(x) \\ge a`` into ``[g(x) - a] \\in \\mathbb{R}_+``
  * ``g(x) \\le a`` into ``[g(x) - a] \\in \\mathbb{R}_-``
  * ``g(x) == a`` into ``[g(x) - a] \\in \\{0\\}``

where `T` is the coefficient type of `g(x) - a`.

## Source node

`VectorizeBridge` supports:

  * `G` in [`MOI.GreaterThan{T}`](@ref)
  * `G` in [`MOI.LessThan{T}`](@ref)
  * `G` in [`MOI.EqualTo{T}`](@ref)

## Target nodes

`VectorizeBridge` creates:

  * `F` in `S`, where `S` is one of [`MOI.Nonnegatives`](@ref),
    [`MOI.Nonpositives`](@ref), [`MOI.Zeros`](@ref) depending on the type of the
    input set.
"""
mutable struct VectorizeBridge{T,F,S,G} <: AbstractBridge
    vector_constraint::MOI.ConstraintIndex{F,S}
    set_constant::T # constant in scalar set
end

const Vectorize{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{VectorizeBridge{T},OT}

function bridge_constraint(
    ::Type{VectorizeBridge{T,F,S,G}},
    model::MOI.ModelLike,
    scalar_f::G,
    set::MOI.Utilities.ScalarLinearSet{T},
) where {T,F,S,G}
    scalar_const = MOI.constant(scalar_f, T)
    if !iszero(scalar_const)
        throw(MOI.ScalarFunctionConstantNotZero{T,G,typeof(set)}(scalar_const))
    end
    vector_f = convert(F, scalar_f)
    set_const = MOI.constant(set)
    MOI.Utilities.operate_output_index!(-, T, 1, vector_f, set_const)
    vector_constraint = MOI.add_constraint(model, vector_f, S(1))
    return VectorizeBridge{T,F,S,G}(vector_constraint, set_const)
end

function MOI.supports_constraint(
    ::Type{VectorizeBridge{T}},
    ::Type{<:MOI.AbstractScalarFunction},
    ::Type{<:MOI.Utilities.ScalarLinearSet{T}},
) where {T}
    return true
end

function MOI.Bridges.added_constrained_variable_types(::Type{<:VectorizeBridge})
    return Tuple{Type}[]
end

function MOI.Bridges.added_constraint_types(
    ::Type{<:VectorizeBridge{T,F,S}},
) where {T,F,S}
    return Tuple{Type,Type}[(F, S)]
end

function concrete_bridge_type(
    ::Type{<:VectorizeBridge{T}},
    G::Type{<:MOI.AbstractScalarFunction},
    S::Type{<:MOI.Utilities.ScalarLinearSet{T}},
) where {T}
    H = MOI.Utilities.promote_operation(-, T, G, T)
    F = MOI.Utilities.promote_operation(vcat, T, H)
    return VectorizeBridge{T,F,MOI.Utilities.vector_set_type(S),G}
end

function MOI.get(
    ::VectorizeBridge{T,F,S},
    ::MOI.NumberOfConstraints{F,S},
)::Int64 where {T,F,S}
    return 1
end

function MOI.get(
    bridge::VectorizeBridge{T,F,S},
    ::MOI.ListOfConstraintIndices{F,S},
) where {T,F,S}
    return [bridge.vector_constraint]
end

function MOI.delete(model::MOI.ModelLike, bridge::VectorizeBridge)
    MOI.delete(model, bridge.vector_constraint)
    return
end

function MOI.supports(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimalStart,MOI.ConstraintDualStart},
    ::Type{VectorizeBridge{T,F,S,G}},
) where {T,F,S,G}
    return MOI.supports(model, attr, MOI.ConstraintIndex{F,S})
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimalStart,
    bridge::VectorizeBridge,
)
    x = MOI.get(model, attr, bridge.vector_constraint)
    @assert length(x) == 1
    return x[1] + bridge.set_constant
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimal,
    bridge::VectorizeBridge,
)
    x = MOI.get(model, attr, bridge.vector_constraint)
    @assert length(x) == 1
    if MOI.Utilities.is_ray(MOI.get(model, MOI.PrimalStatus(attr.result_index)))
        # If it is an infeasibility certificate, it is a ray and satisfies the
        # homogenized problem, see https://github.com/jump-dev/MathOptInterface.jl/issues/433
        return x[1]
    else
        # Otherwise, we need to add the set constant since the ConstraintPrimal
        # is defined as the value of the function and the set_constant was
        # removed from the original function
        return x[1] + bridge.set_constant
    end
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimalStart,
    bridge::VectorizeBridge,
    value,
)
    MOI.set(
        model,
        attr,
        bridge.vector_constraint,
        [value - bridge.set_constant],
    )
    return
end

function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintDual,MOI.ConstraintDualStart},
    bridge::VectorizeBridge,
)
    x = MOI.get(model, attr, bridge.vector_constraint)
    @assert length(x) == 1
    return x[1]
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintDualStart,
    bridge::VectorizeBridge,
    value,
)
    MOI.set(model, attr, bridge.vector_constraint, [value])
    return
end

function MOI.modify(
    model::MOI.ModelLike,
    bridge::VectorizeBridge,
    change::MOI.ScalarCoefficientChange,
)
    MOI.modify(
        model,
        bridge.vector_constraint,
        MOI.MultirowChange(change.variable, [(1, change.new_coefficient)]),
    )
    return
end

function MOI.set(
    model::MOI.ModelLike,
    ::MOI.ConstraintSet,
    bridge::VectorizeBridge,
    new_set::MOI.Utilities.ScalarLinearSet,
)
    bridge.set_constant = MOI.constant(new_set)
    MOI.modify(
        model,
        bridge.vector_constraint,
        MOI.VectorConstantChange([-bridge.set_constant]),
    )
    return
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintFunction,
    bridge::VectorizeBridge{T,F,S,G},
) where {T,F,S,G}
    f = MOI.Utilities.scalarize(
        MOI.get(model, attr, bridge.vector_constraint),
        true,
    )
    @assert length(f) == 1
    return convert(G, f[1])
end

function MOI.get(
    ::MOI.ModelLike,
    ::MOI.ConstraintSet,
    bridge::VectorizeBridge{T,F,S},
) where {T,F,S}
    return MOI.Utilities.scalar_set_type(S, T)(bridge.set_constant)
end


text:
@info "build started."
using Documenter
using PersistenceDiagrams

makedocs(;
    sitename="PersistenceDiagrams.jl",
    format=Documenter.HTML(
        # Use clean URLs, unless built as a "local" build;
        ;
        prettyurls=!("local" in ARGS),
        assets=["assets/favicon.ico"],
    ),
    pages=[
        "Home" => "index.md",
        "Basics" => "basics.md",
        "Distances and Matchings" => "distances.md",
        "Vectorization" => "vectorization.md",
        "MLJ Models" => "mlj.md",
    ],
    doctest=false, # Doctests are run as part of testing -- no need to run them twice.
)

deploydocs(; repo="github.com/mtsch/PersistenceDiagrams.jl.git")


text:
"""
    ReducedMatrix{I<:AbstractCell,C<:Chain}

Representation of the reduced part of the matrix. Is indexed by `I`, indexing returns views
into a `Chain` of type `C`.

!!! warning
    This type does no input validation. Doing things like adding columns that already exist
    is undefined.
"""
struct ReducedMatrix{I<:AbstractCell,F,S,C<:Chain{F,S}}
    column_index::Dict{I,Int}
    indices::Vector{Int}
    values::C

    function ReducedMatrix{I,F,S}() where {I,F,S}
        column_index = Dict{S,Int}()
        indices = Int[1]
        values = Chain{F,S}()
        return new{I,F,S,typeof(values)}(column_index, indices, values)
    end
end

function Base.sizehint!(matrix::ReducedMatrix, size)
    sizehint!(matrix.column_index, size)
    sizehint!(matrix.indices, size + 1)
    sizehint!(matrix.values, size)
    return matrix
end

function Base.setindex!(matrix::ReducedMatrix, column::Chain, simplex)
    @assert sign(simplex) == 1
    n = length(column)
    if n > 0
        append!(matrix.values, column)
        push!(matrix.indices, matrix.indices[end] + n)
        matrix.column_index[simplex] = length(matrix.indices) - 1
    end
    return column
end

Base.length(matrix::ReducedMatrix) = length(matrix.indices) - 1
function Base.haskey(matrix::ReducedMatrix{I}, simplex::I) where {I}
    return haskey(matrix.column_index, simplex)
end
function Base.haskey(matrix::ReducedMatrix{I}, element::AbstractChainElement{I}) where {I}
    return haskey(matrix.column_index, simplex(element))
end

function Base.getindex(matrix::ReducedMatrix{S}, element::AbstractChainElement{S}) where {S}
    return matrix[simplex(element)]
end
function Base.getindex(matrix::ReducedMatrix{S}, simplex::S) where {S}
    index = get(matrix.column_index, simplex, 0)
    if index == 0
        from = 0
        to = -1
    else
        from = matrix.indices[index]
        to = matrix.indices[index + 1] - 1
    end
    # Constructing directly is equivalent to Base.unsafe_view
    return SubArray(matrix.values, (from:to,))
end

"""
    bulk_insert!(matrix::ReducedMatrix, pairs)

Insert apparent pairs into reduced matrix. Columns are inserted so `matrix[τ] == [σ]`, where
`pairs` is a collection of tuples `(σ, τ)`.
"""
function bulk_insert!(matrix::ReducedMatrix{<:Any,F,S}, pairs) where {F,S}
    n = length(pairs)
    n_values = length(matrix.values)
    resize!(matrix.values, n_values + n)
    n_indices = length(matrix.indices)
    last_index = last(matrix.indices)
    resize!(matrix.indices, n_indices + n)
    @inbounds for (i, (σ, τ)) in enumerate(pairs)
        matrix.indices[n_indices + i] = last_index + i
        matrix.values[n_values + i] = ChainElement{S,F}(σ, sign(τ))
        matrix.column_index[abs(τ)] = n_indices + i - 1
    end
    return matrix
end
function bulk_insert!(matrix::ReducedMatrix, ::Tuple{})
    return matrix
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 0.8.6
#   kernelspec:
#     display_name: Julia 1.0.3
#     language: julia
#     name: julia-1.0
# ---

module ExactRecursiveExpressions
using Plots

#Get Final Size for SIR Model
@fastmath function getSIR_final_size(N,β,γ)
    final_size = zeros(N+1,1);
    final_size[2] = 1;
    for Z2 in 0:N
        @inbounds for Z1 in Z2+1:N-1
            p1 = 1 / ( 1 + γ/(β*(N-Z1)));
            final_size[Z1+2] = final_size[Z1+2] + final_size[Z1+1]*p1;
            final_size[Z1+1] = final_size[Z1+1]*(1-p1);
        end
    end
    return final_size;
end

#Sample Size
N = 20;
#Infected per-Capita Rate
β = 2/(N-1);
#Recover per-capita rate
γ = 1.0;

#Get Final Size for SIR Model
final_size = getSIR_final_size(N,β,γ);

#Get Time Final Size Function
@time final_size = getSIR_final_size(N,β,γ);

bar(0:N,final_size)

@fastmath function getSI4R_final_size(N::Int64,β::Float64,γ::Float64)
    Psi = (N+1)*(N+2)*(N+3)*(N+4)/24;
    p_vec = zeros(Float64,Int(Psi),1);
    p_vec[2]=1;
    @inbounds for Z5 in 0:N
        w::Int64 = Psi - (N-Z5+1)*(N-Z5+2)*(N-Z5+3)*(N-Z5+4)/24 + 1;
        @inbounds for Z4 in Z5:N
            a5 = 4*γ*(Z4-Z5);
            @inbounds for Z3 in Z4:N
                a4 = 4*γ*(Z3-Z4);
                @inbounds for Z2 in Z3:N;
                    a3 = 4*γ*(Z2-Z3);
                    @inbounds for Z1 in Z2:N
                    a1 = β*(N-Z1)*(Z1-Z5);
                    a2 = 4*γ*(Z1-Z2);
                    tot = a1+a2+a3+a4+a5;
                    if Z1-Z5 == 0
                        final_size[Z5+1] = p_vec[w];
                    end
                    if a1 > 0
                       p_vec[w+1] = p_vec[w+1]+ p_vec[w]*a1/tot;
                    end
                    if a2 > 0
                        p_vec[w+N-Z2] = p_vec[w+N-Z2]+ p_vec[w]*a2/tot;
                    end
                    if a3 > 0
                        place3::Int64 = (N-Z3)*(N-Z3+1)/2;
                        p_vec[w+place3] = p_vec[w+place3]+ p_vec[w]*a3/tot;
                    end
                    if a4 > 0
                        place4::Int64 = (N-Z4)*(N-Z4+1)*(N-Z4+2)/6;
                        p_vec[w+place4] = p_vec[w+place4] + p_vec[w]*a4/tot;
                    end
                    if a5 > 0
                        p_vec[w] = p_vec[w]*a5/tot;
                    end
                    w = w + 1;
                end
            end
        end
    end
end
    return final_size
end

#Get Final Size for SI(4)R model
final_size = getSI4R_final_size(N,β,γ);

#Get Run time of SI(4)R Final Size Model
@time final_size = getSI4R_final_size(N,β,γ);

bar(0:N,final_size)

end


text:
module ColumnwiseLUSolver

export ManyColumnLU, SingleColumnLU

using ..Mesh.Grids
using ..Mesh.Topologies
using ..DGmethods
using ..DGmethods: BalanceLaw, DGModel, num_state, num_diffusive
using ..LinearSolvers
const LS = LinearSolvers
using ..MPIStateArrays
using LinearAlgebra
using KernelAbstractions

include("ColumnwiseLUSolver_kernels.jl")

abstract type AbstractColumnLUSolver <: AbstractLinearSolver end

"""
    ManyColumnLU()

This solver is used for systems that are block diagonal where each block is
associated with a column of the mesh.  The systems are solved using a
non-pivoted LU factorization.
"""
struct ManyColumnLU <: AbstractColumnLUSolver end

"""
    SingleColumnLU()

This solver is used for systems that are block diagonal where each block is
associated with a column of the mesh.  Moreover, each block is assumed to be
the same.  The systems are solved using a non-pivoted LU factorization.
"""
struct SingleColumnLU <: AbstractColumnLUSolver end

struct ColumnwiseLU{F, AT}
    f::F
    A::AT
end

function LS.prefactorize(op, solver::AbstractColumnLUSolver, Q, args...)
    dg = op.f!

    # TODO: can we get away with just passing the grid?
    A = banded_matrix(
        op,
        dg,
        similar(Q),
        similar(Q),
        args...;
        single_column = typeof(solver) <: SingleColumnLU,
    )

    band_lu!(A, dg)

    ColumnwiseLU(dg, A)
end

function LS.linearsolve!(
    clu::ColumnwiseLU{F},
    ::AbstractColumnLUSolver,
    Q,
    Qrhs,
    args...,
) where {F <: DGModel}
    device = typeof(Q.data) <: Array ? CPU() : CUDA()

    dg = clu.f
    A = clu.A
    Q .= Qrhs

    band_forward!(Q, A, dg)
    band_back!(Q, A, dg)
end

"""
    band_lu!(A, dg::DGModel)

"""
function band_lu!(A, dg::DGModel)
    bl = dg.balancelaw
    grid = dg.grid
    topology = grid.topology
    @assert isstacked(topology)
    @assert typeof(dg.direction) <: VerticalDirection

    FT = eltype(A)
    device = typeof(A) <: Array ? CPU() : CUDA()

    nstate = num_state(bl, FT)
    N = polynomialorder(grid)
    Nq = N + 1
    Nqj = dimensionality(grid) == 2 ? 1 : Nq

    eband = num_diffusive(bl, FT) == 0 ? 1 : 2

    nrealelem = length(topology.elems)
    nvertelem = topology.stacksize
    nhorzelem = div(nrealelem, nvertelem)

    groupsize = (Nq, Nqj)
    ndrange = (nhorzelem * Nq, Nqj)

    if ndims(A) == 2
        # single column case
        #
        # TODO Would it be faster to copy the matrix to the host and factorize it
        # there?
        groupsize = (1, 1)
        ndrange = groupsize
        A = reshape(A, 1, 1, size(A)..., 1)
    end

    event = Event(device)
    event = band_lu_knl!(device, groupsize)(
        A,
        Val(Nq),
        Val(groupsize[1]),
        Val(groupsize[2]),
        Val(nstate),
        Val(nvertelem),
        Val(ndrange[end]),
        Val(eband);
        ndrange = ndrange,
        dependencies = (event,),
    )
    wait(device, event)
end

function band_forward!(Q, A, dg::DGModel)
    bl = dg.balancelaw
    grid = dg.grid
    topology = grid.topology
    @assert isstacked(topology)
    @assert typeof(dg.direction) <: VerticalDirection

    FT = eltype(A)
    device = typeof(Q.data) <: Array ? CPU() : CUDA()

    nstate = num_state(bl, FT)
    N = polynomialorder(grid)
    Nq = N + 1
    Nqj = dimensionality(grid) == 2 ? 1 : Nq

    eband = num_diffusive(bl, FT) == 0 ? 1 : 2

    nrealelem = length(topology.elems)
    nvertelem = topology.stacksize
    nhorzelem = div(nrealelem, nvertelem)

    event = Event(device)
    event = band_forward_knl!(device, (Nq, Nqj))(
        Q.data,
        A,
        Val(Nq),
        Val(Nqj),
        Val(nstate),
        Val(nvertelem),
        Val(nhorzelem),
        Val(eband);
        ndrange = (nhorzelem * Nq, Nqj),
        dependencies = (event,),
    )
    wait(device, event)
end

function band_back!(Q, A, dg::DGModel)
    bl = dg.balancelaw
    grid = dg.grid
    topology = grid.topology
    @assert isstacked(topology)
    @assert typeof(dg.direction) <: VerticalDirection

    FT = eltype(A)
    device = typeof(Q.data) <: Array ? CPU() : CUDA()

    nstate = num_state(bl, FT)
    N = polynomialorder(grid)
    Nq = N + 1
    Nqj = dimensionality(grid) == 2 ? 1 : Nq

    eband = num_diffusive(bl, FT) == 0 ? 1 : 2

    nrealelem = length(topology.elems)
    nvertelem = topology.stacksize
    nhorzelem = div(nrealelem, nvertelem)

    event = Event(device)
    event = band_back_knl!(device, (Nq, Nqj))(
        Q.data,
        A,
        Val(Nq),
        Val(Nqj),
        Val(nstate),
        Val(nvertelem),
        Val(nhorzelem),
        Val(eband);
        ndrange = (nhorzelem * Nq, Nqj),
        dependencies = (event,),
    )
    wait(device, event)
end

"""
    banded_matrix(dg::DGModel, [Q::MPIStateArray, dQ::MPIStateArray,
                  single_column=false])

Forms the banded matrices for each the column operator defined by the `DGModel`
dg.  If `single_column=false` then a banded matrix is stored for each column and
if `single_column=true` only the banded matrix associated with the first column
of the first element is stored. The bandwidth of the DG column banded matrix is
`p = q = (polynomialorder + 1) * nstate * nvertelem - 1` with `p` and `q` being
the upper and lower bandwidths.

The banded matrices are stored in the LAPACK band storage format
<https://www.netlib.org/lapack/lug/node124.html>.

The banded matrices are returned as an arrays where the array type matches that
of `Q`. If `single_column=false` then the returned array has 5 dimensions, which
are:
- first horizontal column index
- second horizontal column index
- band index (-q:p)
- vertical DOF index with state `s`, vertical DOF index `k`, and vertical
  element `ev` mapping to `s + nstate * (k - 1) + nstate * nvertelem * (ev - 1)`
- horizontal element index

If the `single_column=true` then the returned array has 2 dimensions which are
the band index and the vertical DOF index.
"""
function banded_matrix(
    dg::DGModel,
    Q::MPIStateArray = MPIStateArray(dg),
    dQ::MPIStateArray = MPIStateArray(dg);
    single_column = false,
)
    banded_matrix(
        (dQ, Q) -> dg(dQ, Q, nothing, 0; increment = false),
        dg,
        Q,
        dQ;
        single_column = single_column,
    )
end

"""
    banded_matrix(f!::Function, dg::DGModel,
                  Q::MPIStateArray = MPIStateArray(dg),
                  dQ::MPIStateArray = MPIStateArray(dg), args...;
                  single_column = false, args...)

Forms the banded matrices for each the column operator defined by the linear
operator `f!` which is assumed to have the same banded structure as the
`DGModel` dg.  If `single_column=false` then a banded matrix is stored for each
column and if `single_column=true` only the banded matrix associated with the
first column of the first element is stored. The bandwidth of the DG column
banded matrix is `p = q = (polynomialorder + 1) * nstate * nvertelem - 1` with
`p` and `q` being the upper and lower bandwidths.

The banded matrices are stored in the LAPACK band storage format
<https://www.netlib.org/lapack/lug/node124.html>.

The banded matrices are returned as an arrays where the array type matches that
of `Q`. If `single_column=false` then the returned array has 5 dimensions, which
are:
- first horizontal column index
- second horizontal column index
- band index (-q:p)
- vertical DOF index with state `s`, vertical DOF index `k`, and vertical
  element `ev` mapping to `s + nstate * (k - 1) + nstate * nvertelem * (ev - 1)`
- horizontal element index

If the `single_column=true` then the returned array has 2 dimensions which are
the band index and the vertical DOF index.

Here `args` are passed to `f!`.
"""
function banded_matrix(
    f!,
    dg::DGModel,
    Q::MPIStateArray = MPIStateArray(dg),
    dQ::MPIStateArray = MPIStateArray(dg),
    args...;
    single_column = false,
)
    bl = dg.balancelaw
    grid = dg.grid
    topology = grid.topology
    @assert isstacked(topology)
    @assert typeof(dg.direction) <: VerticalDirection

    FT = eltype(Q.data)
    device = typeof(Q.data) <: Array ? CPU() : CUDA()

    nstate = num_state(bl, FT)
    N = polynomialorder(grid)
    Nq = N + 1

    # p is lower bandwidth
    # q is upper bandwidth
    eband = num_diffusive(bl, FT) == 0 ? 1 : 2
    p = q = nstate * Nq * eband - 1

    nrealelem = length(topology.elems)
    nvertelem = topology.stacksize
    nhorzelem = div(nrealelem, nvertelem)

    dim = dimensionality(grid)

    Nqj = dim == 2 ? 1 : Nq

    # first horizontal DOF index
    # second horizontal DOF index
    # band index -q:p
    # vertical DOF index
    # horizontal element index
    A = if single_column
        similar(Q.data, p + q + 1, Nq * nstate * nvertelem)
    else
        similar(Q.data, Nq, Nqj, p + q + 1, Nq * nstate * nvertelem, nhorzelem)
    end
    fill!(A, zero(FT))

    # loop through all DOFs in a column and compute the matrix column
    for ev in 1:nvertelem
        for s in 1:nstate
            for k in 1:Nq
                # Set a single 1 per column and rest 0
                event = Event(device)
                event = knl_set_banded_data!(device, (Nq, Nqj, Nq))(
                    bl,
                    Val(dim),
                    Val(N),
                    Val(nvertelem),
                    Q.data,
                    k,
                    s,
                    ev,
                    1:nhorzelem,
                    1:nvertelem;
                    ndrange = (nvertelem * Nq, nhorzelem * Nqj, Nq),
                    dependencies = (event,),
                )
                wait(device, event)

                # Get the matrix column
                f!(dQ, Q, args...)

                # Store the banded matrix
                event = Event(device)
                event = knl_set_banded_matrix!(device, (Nq, Nqj, Nq))(
                    bl,
                    Val(dim),
                    Val(N),
                    Val(nvertelem),
                    Val(p),
                    Val(q),
                    Val(eband + 1),
                    A,
                    dQ.data,
                    k,
                    s,
                    ev,
                    1:nhorzelem,
                    (-eband):eband;
                    ndrange = ((2eband + 1) * Nq, nhorzelem * Nqj, Nq),
                    dependencies = (event,),
                )
                wait(device, event)
            end
        end
    end
    A
end

"""
    banded_matrix_vector_product!(dg::DGModel, A, dQ::MPIStateArray,
                                  Q::MPIStateArray)

Compute a matrix vector product `dQ = A * Q` where `A` is assumed to be a matrix
created using the `banded_matrix` function.

This function is primarily for testing purposes.
"""
function banded_matrix_vector_product!(
    dg::DGModel,
    A,
    dQ::MPIStateArray,
    Q::MPIStateArray,
)
    bl = dg.balancelaw
    grid = dg.grid
    topology = grid.topology
    @assert isstacked(topology)
    @assert typeof(dg.direction) <: VerticalDirection

    FT = eltype(Q.data)
    device = typeof(Q.data) <: Array ? CPU() : CUDA()

    eband = num_diffusive(bl, FT) == 0 ? 1 : 2
    nstate = num_state(bl, FT)
    N = polynomialorder(grid)
    Nq = N + 1
    p = q = nstate * Nq * eband - 1

    nrealelem = length(topology.elems)
    nvertelem = topology.stacksize
    nhorzelem = div(nrealelem, nvertelem)

    dim = dimensionality(grid)

    Nqj = dim == 2 ? 1 : Nq

    event = Event(device)
    event = knl_banded_matrix_vector_product!(device, (Nq, Nqj, Nq))(
        bl,
        Val(dim),
        Val(N),
        Val(nvertelem),
        Val(p),
        Val(q),
        dQ.data,
        A,
        Q.data,
        1:nhorzelem,
        1:nvertelem;
        ndrange = (nvertelem * Nq, nhorzelem * Nqj, Nq),
        dependencies = (event,),
    )
    wait(device, event)
end

end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
# Solve the n-queens problem to test optional constraints.

"Creates the constraints matrix for the n-queens problem."
function make_queens_matrix(n)
  slash_diagonal_for(row, column) = row + column - 1
  backslash_diagonal_for(row, column) = n + row - column

  constraints_for(row, column) = Set(Any[
    (:row, row),
    (:column, column),
    (:slash, slash_diagonal_for(row, column)),
    (:backslash, backslash_diagonal_for(row, column))
  ])

  CoverMatrix(
    Dict([(row, column) => constraints_for(row, column) for row in 1:n, column in 1:n]),
    optionals = union(
      Set(Any[(:slash, x) for x in 1:(2n-1)]),
      Set(Any[(:backslash, x) for x in 1:(2n-1)])
    )
  )
end

"See [OEIS sequence A000170](http://oeis.org/A000170)."
const oeis_A000170 = [1, 0, 0, 2, 10, 4, 40, 92, 352, 724]

let
  for (queen_count, position_count) in enumerate(oeis_A000170)
    @test exact_cover_count(make_queens_matrix(queen_count)) == position_count
  end
end


text:
## In-place testing
#======================================#

# function surface_velocity!(vel, r, horseshoes, Γs, U, Ω)
#     U_hat = -normalize(U)
#     for i in eachindex(horseshoes)
#         vel += @views induced_trailing_velocity(r, horseshoes[i], Γs[i], U_hat) - (U + Ω × r)
#     end

#     vel
# end

# function evaluate_aerodynamics!(R, horseshoes, Γs, U, Ω, ρ, speed)
#     forces = Vector{typeof(U)}(undef, length(horseshoes))

#     U_hat = U / speed
#     Ω_hat = Ω / speed
#     for i in eachindex(horseshoes)
#         # Residual
#         R[i] = @views velocity(horseshoes, Γs / speed, collocation_point(horseshoes[i]), horseshoe_normal(horseshoes[i]), U_hat, Ω_hat)

#         # Forces
#         V_trail   = @views surface_velocity(horseshoes[i], Γs, horseshoes, U, Ω)
#         forces[i] = @views kutta_joukowsky(ρ, Γs[i], V_trail, bound_leg_vector(horseshoes[i]))
#     end

#     forces
# end

# function new_surface_forces(Γs_comp, hs_comp, Γs, horseshoes, U, Ω, ρ, speed)
#     U_hat  = -U / speed
#     forces = similar(Γs_comp, typeof(U))
#     for i in eachindex(hs_comp)
#         r_i = @views bound_leg_center(hs_comp[i])
#         l_i = @views bound_leg_vector(hs_comp[i])
#         vel = zeros(eltype(r_i), 3)
#         vel_fs = -(U + Ω × r_i)
#         for j in eachindex(horseshoes)
#             vel += @views trailing_velocity(r_i, horseshoes[j], Γs[j], U_hat) + vel_fs
#         end
#         forces[i] = @views kutta_joukowsky(ρ, Γs_comp[i], vel, l_i)
#    end

#    forces
# end


text:
export Sampler, RandomSampler, GridSampler, ParemeterGenerator

abstract type Sampler end

abstract type RandomSampler <: Sampler end
abstract type GridSampler <: Sampler end

abstract type ParameterGenerator <: AbstractVector{Any} end


text:
@testset "Brownian Bridge" begin

using DiffEqNoiseProcess, DiffEqBase, Test, DataStructures, Random, DiffEqBase.EnsembleAnalysis

Random.seed!(100)
W = BrownianBridge(0.0,1.0,0.0,1.0,0.0,0.0)
prob = NoiseProblem(W,(0.0,1.0))
ensemble_prob = EnsembleProblem(prob)
@time sol = solve(ensemble_prob,dt=0.1,trajectories=100000)

# Spot check the mean and the variance
qs = 0:0.1:1
for i in 2:10
  q = qs[i]
  @test ≈(timestep_mean(sol,i),q,atol=1e-2)
  @test ≈(timestep_meanvar(sol,i)[2],(1-q)*q,atol=1e-2)
end
@test ≈(timestep_mean(sol,1)[1],0.0,atol=1e-16)
@test ≈(timestep_meanvar(sol,1)[2],0.0,atol=1e-16)
@test ≈(timestep_mean(sol,11)[1],1.0,atol=1e-16)
@test ≈(timestep_meanvar(sol,11)[2],0.0,atol=1e-16)

μ = 1.2
σ = 2.2
W = GeometricBrownianBridge(μ,σ,0.0,1.0,0.0,1.0,0.0,0.0)
prob = NoiseProblem(W,(0.0,1.0))
ensemble_prob = EnsembleProblem(prob)
@time sol = solve(ensemble_prob,dt=0.1,trajectories=100)

Random.seed!(100)
r = 100 # should be independent of the rate, so make it crazy
rate(u,p,t) = r
W = CompoundPoissonBridge(rate,0.0,1.0,0.0,1.0)
prob = NoiseProblem(W,(0.0,1.0))
ensemble_prob = EnsembleProblem(prob)
@time sol = solve(ensemble_prob,dt=0.1,trajectories=100000)

# Spot check the mean and the variance
qs = 0:0.1:1
for i in 2:10
  q = qs[i]
  # Mean and variance of binomial matches that of the Brownian bridge!
  @test ≈(timestep_mean(sol,i),q,atol=1e-2)
  @test ≈(timestep_meanvar(sol,i)[2],(1-q)*q,atol=1e-2)
end
@test ≈(timestep_mean(sol,1)[1],0.0,atol=1e-16)
@test ≈(timestep_meanvar(sol,1)[2],0.0,atol=1e-16)
@test ≈(timestep_mean(sol,11)[1],1.0,atol=1e-16)
@test ≈(timestep_meanvar(sol,11)[2],0.0,atol=1e-16)

end


text:
"""
    SpaceShuttle

This module implements a little bit about the space shuttle
"""
module SpaceShuttle

export launch, deploy_canadarm

"""
    launch()

Tell the _space shuttle_ to take off!

Does the following: \$\\sum_0^N x_i\$
"""
launch() = "3..2...1... launch!"

"""
    deploy_canadarm

This will deploy the Canada 🍁 arm
"""
deploy_canadarm() = "grabbing the satellite"

end # module


text:
@testset "fdiff" begin
    @testset "API" begin
        # fdiff! works the same as fdiff
        mat_in = rand(3, 3, 3)
        mat_out = similar(mat_in)
        fdiff!(mat_out, mat_in, dims = 2)
        @test mat_out == fdiff(mat_in, dims = 2)

        mat_in = rand(3, 3, 3)
        mat_out = similar(mat_in)
        fdiff!(mat_out, mat_in, dims = 3, rev=true)
        @test mat_out == fdiff(mat_in, dims = 3, rev=true)

        # default `dims` are only available when input is a vector
        v = rand(9)
        @test fdiff(v) == fdiff(v; dims=1)
        m = rand(3, 3)
        @test_throws UndefKeywordError fdiff(m)
    end

    @testset "NumericalTests" begin
        a = reshape(collect(1:9), 3, 3)
        b_fd_1 = [1 1 1; 1 1 1; -2 -2 -2]
        b_fd_2 = [3 3 -6; 3 3 -6; 3 3 -6]
        b_bd_1 = [-2 -2 -2; 1 1 1; 1 1 1]
        b_bd_2 = [-6 3 3; -6 3 3; -6 3 3]
        out = similar(a)

        @test fdiff(a, dims = 1) == b_fd_1
        @test fdiff(a, dims = 2) == b_fd_2
        @test fdiff(a, dims = 1, rev=true) == b_bd_1
        @test fdiff(a, dims = 2, rev=true) == b_bd_2
        fdiff!(out, a, dims = 1)
        @test out == b_fd_1
        fdiff!(out, a, dims = 2)
        @test out == b_fd_2
        fdiff!(out, a, dims = 1, rev=true)
        @test out == b_bd_1
        fdiff!(out, a, dims = 2, rev=true)
        @test out == b_bd_2

        @testset "Boundary" begin
            # These originally lives in Images and we use it to check the compatibility
            forwarddiffy(u::AbstractMatrix) = [u[2:end,:]; u[end:end,:]] - u
            forwarddiffx(u::AbstractMatrix) = [u[:,2:end] u[:,end:end]] - u
            backdiffy(u::AbstractMatrix) = u - [u[1:1,:]; u[1:end-1,:]]
            backdiffx(u::AbstractMatrix) = u - [u[:,1:1] u[:,1:end-1]]

            X = rand(3, 3)
            @test forwarddiffx(X) == fdiff(X, dims=2, boundary=:zero)
            @test forwarddiffy(X) == fdiff(X, dims=1, boundary=:zero)
            @test backdiffx(X) == fdiff(X, dims=2, rev=true, boundary=:zero)
            @test backdiffy(X) == fdiff(X, dims=1, rev=true, boundary=:zero)
        end

        if VERSION >= v"1.1"
            # dims keyword for diff only exists after Julia 1.1
            # check numerical results with base implementation
            drop_last_slice(X, dims) = collect(StackView(collect(eachslice(X, dims=dims))[1:end-1]..., dims=dims))
            function drop_last_slice(X::AbstractVector, dims)
                @assert dims==1
                X[1:end-1]
            end
            drop_first_slice(X, dims) = collect(StackView(collect(eachslice(X, dims=dims))[2:end]..., dims=dims))
            function drop_first_slice(X::AbstractVector, dims)
                @assert dims==1
                X[2:end]
            end

            for N in 1:3
                sz = ntuple(_->5, N)
                A = rand(sz...)
                A_out = similar(A)

                for dims = 1:N
                    out_base = diff(A; dims=dims)
                    out = fdiff(A; dims=dims)
                    @test out_base == drop_last_slice(out, dims)

                    out_base = .-reverse(diff(reverse(A; dims=dims); dims=dims); dims=dims)
                    out = fdiff(A; dims=dims, rev=true)
                    @test out_base == drop_first_slice(out, dims)
                end
            end
        end
    end

    @testset "OffsetArrays" begin
        A = OffsetArray(rand(3, 3), -1, -1)
        A_out = fdiff(A, dims=1)
        @test axes(A_out) == (0:2, 0:2)
        @test A_out.parent == fdiff(parent(A), dims=1)

        A = OffsetArray(rand(3, 3), -1, -1)
        A_out = fdiff(A, dims=1, rev=true)
        @test axes(A_out) == (0:2, 0:2)
        @test A_out.parent == fdiff(parent(A), dims=1, rev=true)
    end
end