text:
using HttpServer
using WebSockets
using JSON

include("depparse.jl")

type Token
    i::Int
    j::Int
    label::String
end

function Token(i::Int, j::Int)
    #h = i + j
    #id = Int(rem(h,length(posconf))) + 1
    #label = posconf[id][1]
    Token(i, j, "label")
end

function readconf(path::String)
    lines = open(readlines, path)
    map(lines) do line
        items = split(chomp(line), '\t')
        (items[1], items[2])
    end
end

const filepath = dirname(@__FILE__)
const clients = Dict()
#const posconf = readconf("$(filepath)/pos.conf")
#const entityconf = readconf("$(filepath)/entity.conf")

function create_error(msg::String)
    dict = Dict(
        "error": msg,
    )
    JSON.json(dict)
end

wsh = WebSocketHandler() do req, client
    println("Client: $(client.id) is connected.")
    while true
        #println("Request from $(client.id) recieved.")
        msg = String(read(client))
        if length(msg) > 2000
            write(client, create_error("Error: the mssage length exceeds 2000."))
            continue
        end

        json = try
            JSON.parse(msg)
        catch
            write(client, create_error("Error: invalid json format: $(msg)"))
            continue
        end
        if !haskey(json, "text")
            write(client, create_error("Error: key: text does not exist. Check the JSON format."))
            continue
        end
        text = json["text"]

        # line
        lines = Vector{Int}[]
        start = 0
        for i = 1:length(text)
            char = text[i]
            if char == '\n'
                start > 0 && push!(lines,[start-1,i-2])
                start = 0
            else
                start == 0 && (start = i)
            end
        end
        start > 0 && push!(lines,[start-1,length(text)-1])

        # token
        spans = Dict()
        id = 1
        start = 0
        for i = 1:length(text)
            char = text[i]
            if char == ' ' || char == '\n'
                selecter = rand()
                label = selecter<0.5 ? "token" : "superLongToken"
                color = selecter<0.5 ? "#84b62b" : "#FBFF5E"
                start > 0 && (spans["span-$id"] = [start-1,i-2,label,color])
                id += 1
                start = 0
            else
                start == 0 && (start = i)
            end
        end
        start > 0 && (spans["span-$id"] = [start-1,length(text)-1,"token","#84b62b"])

        # relation
        spankeys = collect(keys(spans))
        rels = Dict()
        id = 1
        for i = 1:length(spans)/3
            ids = rand(1:length(spans), 2)
            rels["rel-$id"] = ["one-way", spankeys[ids[1]], spankeys[ids[2]], "label-$id", "#84b62b"]
            id += 1
        end

        res = Dict("line"=>lines, "span"=>spans, "relation" => rels)
        println(JSON.json(res))
        write(client, JSON.json(res))
    end
end

onepage = readstring(joinpath(dirname(@__FILE__), "index.html"))
httph = HttpHandler() do req::Request, res::Response
    Response(onepage)
end
server = Server(httph, wsh)
run(server, 3000)

max_stars_repo_name: hshindo/react-nlp

max_issues_repo_name: hshindo/react-nlp

max_forks_repo_name: hshindo/react-nlp


text:
# This file is part of Kpax3. License is MIT.
RecipesBase.@userplot Kpax3PlotTrace

RecipesBase.@recipe function plottrace(g::Kpax3PlotTrace;
                                       maxlag=200,
                                       M=20000,
                                       main="")
  if (length(g.args) != 1) ||
     (typeof(g.args[1]) != Vector{Float64})
      throw(KInputError("One argument required. Got: $(typeof(g.args))"))
  end

  entropy = g.args[1]

  nsim = length(entropy)

  # plot the last M points
  if M < nsim
    xtr = max(1, nsim - M):nsim
    ytr = entropy[xtr]
  else
    xtr = 1:nsim
    ytr = copy(entropy)
  end

  plot_title = if main != ""
    ac = StatsBase.autocov(entropy, 0:maxlag)

    variid = ac[1] / nsim
    mcvar = imsevar(ac, nsim)
    mcse = sqrt(mcvar)
    eff = min(nsim, ess(variid, mcvar, nsim))

    effstr = if eff <= 100000
      Printf.@sprintf("%d", eff)
    else
      Printf.@sprintf("%.3e", eff)
    end

    sestr = if 0.001 <= mcse <= 1000
      Printf.@sprintf("%.3f", mcse)
    else
      Printf.@sprintf("%.3e", mcse)
    end

    string(main, " (mcse = ", sestr, ", ESS = ", effstr, ")")
  else
    ""
  end

  RecipesBase.@series begin
    seriestype := :path
    linecolor := :black
    background_color := :white
    label := ""
    legend := :none
    title := plot_title
    window_title := ""
    xlims := (minimum(xtr), maximum(xtr))
    ylims := (max(0.0, minimum(ytr) - 0.05), maximum(ytr) + 0.05)

    html_output_format --> :svg
    size --> (800, 600)
    xlabel --> "Iteration"
    ylabel --> "Entropy"

    xtr, ytr
  end
end

RecipesBase.@userplot Kpax3PlotDensity

RecipesBase.@recipe function plotdensity(g::Kpax3PlotDensity;
                                         maxlag=200,
                                         main="")
  if (length(g.args) != 1) ||
     (typeof(g.args[1]) != Vector{Float64})
      throw(KInputError("One argument required. Got: $(typeof(g.args))"))
  end

  entropy = g.args[1]

  nsim = length(entropy)

  plot_title = if main != ""
    ac = StatsBase.autocov(entropy, 0:maxlag)

    variid = ac[1] / nsim
    mcvar = imsevar(ac, nsim)
    mcse = sqrt(mcvar)
    eff = min(nsim, ess(variid, mcvar, nsim))

    effstr = if eff <= 100000
      Printf.@sprintf("%d", eff)
    else
      Printf.@sprintf("%.3e", eff)
    end

    sestr = if 0.001 <= mcse <= 1000
      Printf.@sprintf("%.3f", mcse)
    else
      Printf.@sprintf("%.3e", mcse)
    end

    string(main, " (mcse = ", sestr, ", ESS = ", effstr, ")")
  else
    ""
  end

  kd = KernelDensity.kde(entropy, npoints=256)

  RecipesBase.@series begin
    seriestype := :path
    linecolor := :black
    background_color := :white
    label := ""
    legend := :none
    title := plot_title
    window_title := ""
    xlims := (max(0.0, minimum(entropy) - 0.05), maximum(entropy) + 0.05)

    html_output_format --> :svg
    size --> (800, 600)
    xlabel --> "Entropy"
    ylabel --> "Density"

    kd.x, kd.density
  end
end

RecipesBase.@userplot Kpax3PlotJump

RecipesBase.@recipe function plotjump(g::Kpax3PlotJump;
                                      main="")
  if (length(g.args) != 1) ||
     (typeof(g.args[1]) != Vector{Float64})
      throw(KInputError("One argument required. Got: $(typeof(g.args))"))
  end

  avgd = g.args[1]

  maxlag = length(avgd)

  # Exponential smoothing for estimating the long term mean
  asy = avgd[1]
  for t in 2:maxlag
    asy = 0.6 * asy + 0.4 * avgd[t]
  end

  RecipesBase.@series begin
    seriestype := :scatter
    markercolor := :black
    label := ""
    legend := :none

    1:maxlag, avgd
  end

  RecipesBase.@series begin
    seriestype := :hline
    linestyle := :dash
    linecolor := :black
    background_color := :white
    label := ""
    legend := :none
    title := main
    window_title := ""
    xlims := (0, maxlag + 1)

    html_output_format --> :svg
    size --> (800, 600)
    xlabel --> "Lag"
    ylabel --> "Average distance"

    [asy]
  end
end

ext: jl


text:
using ConstrainedDynamics

# Parameters
ex = [1.;0.;0.]
ey = [0.;1.;0.]
ez = [0.;0.;1.]

l1 = 1.0
x, z = .1, .1
b1 = Box(x, l1, l1, l1 * 1, color = RGBA(1., 1., 0.))
b2 = Box(x, l1, z, l1 * 1, color = RGBA(1., 0., 0.))

vert11 = [0;l1 / 2;l1 / 2]
vert12 = -vert11

vert21 = [0;l1 / 2;0.]
vert22 = -vert21

# Links
origin = Origin{Float64}()
link1 = Body(b1)
link2 = Body(b2)

# Constraints
# joint0to1 = EqualityConstraint(OriginConnection(origin,link1))
# joint0to1 = EqualityConstraint(Fixed(origin,link1,zeros(3),zeros(3)))
joint0to1 = EqualityConstraint(Revolute(origin, link1, zeros(3), zeros(3), ex))
joint1to2 = EqualityConstraint(Cylindrical(link1, link2, vert12, vert21, -ey))

links = [link1; link2]
constraints = [joint0to1;joint1to2]
shapes = [b1;b2]

function control!(mechanism, t)
    setForce!(mechanism, joint1to2, [[.1];nothing])
    return
end

mech = Mechanism(origin, links, constraints, shapes = shapes, g = 0., tend = 10.)
setPosition!(mech,link1,q = Quaternion(RotX(pi / 2)))
setPosition!(mech,link1,link2,p1 = vert12,p2 = vert21)
# setVelocity!(mech,link1,ω=[2.;0;0])
# setVelocity!(mech,link1,link2,p1=vert12,p2=vert21)

# setForce!([[0.1];nothing],joint1to2,mech)

simulate!(mech,control!,save = true)
visualize!(mech)


text:
"""
A set of tools to deal with GPS track (find intersection with line, interpolate time of intersection...).
"""
module TrackTools

using Dates

const gamma = 1.0
const G = [gamma; gamma]

"""
    TrackPoint(time::T, x::Float64, y::Float64) where {T <: Dates.AbstractDateTime}

A point on a track.
"""
struct TrackPoint{T<:Dates.AbstractDateTime}
    time::T
    x::Float64
    y::Float64
end

"""
    matrix_line(x1, y1, x2, y2)

Compute A matrix parameters (ie parameters α and β) of a line between 2 points (x1, y1) and (x2, y2).

Equation of line (with `𝛾 = 1`)

``
α \\cdot x + β \\cdot y = 𝛾
``

with 2 points

``
\\begin{cases}
α \\cdot x_1 + β \\cdot y_1 = 𝛾 \\
α \\cdot x_2 + β \\cdot y_1 = 𝛾
\\end{cases}
``

with matrix notation

``
P \\cdot A = G

A = P^-1 \\cdot G
``

where

``
G =
\\begin{bmatrix}
𝛾 \\
𝛾
\\end{bmatrix}
``

and

``
P =
\\begin{bmatrix}
x_1 & y_1 \\
x_2 & y_2
\\end{bmatrix}
``

and

``
A =
\\begin{bmatrix}
α \\
β
\\end{bmatrix}
``
"""
function matrix_line(x1, y1, x2, y2)
    P = [x1 y1; x2 y2]
    return inv(P) * G
end

"""
    hasintersect(x1, y1, x2, y2, xI, yI) -> Bool

Return `true` if (`xI`, `yI`) is inside a rectangle defined by (`x1`, `y1`) and (`x2`, `y2`).
"""
function hasintersect(x1, y1, x2, y2, xI, yI)
    isinxrange(x1, x2, xI) = min(x1, x2) <= xI <= max(x1, x2)
    isinyrange(y1, y2, yI) = min(y1, y2) <= yI <= max(y1, y2)
    return isinxrange(x1, x2, xI) && isinyrange(y1, y2, yI)
end

"""
    interpolate_time(x1, t1::T, x2, t2::T, xI) -> T where {T<:Dates.AbstractDateTime}

Interpolate time of intersection (with a millisecond resolution) given time and position on an x-axis (may also be y)

x(t) = (x2 - x1) / (t2 - t1) * (t - t1) + x1

Δt = t2 - t1

t - t1 = (x - x1) / (x2 - x1) * Δt

t = (x - x1) / (x2 - x1) * Δt + t1
"""
function interpolate_time(x1, t1::T, x2, t2::T, xI) where {T<:Dates.AbstractDateTime}
    Δt = Dates.value(t2 - t1)
    return Millisecond(round(Int, Δt * (xI - x1) / (x2 - x1))) + t1
end

"""
    interpolate_time(x1, y1, t1, x2, y2, t2, xI, yI) -> Dates.AbstractDateTime

Interpolate time of intersection (with a millisecond resolution) given time and axis position where difference was the most significative.
"""
function interpolate_time(x1, y1, t1::T, x2, y2, t2::T, xI, yI) where {T<:Dates.AbstractDateTime}
    if abs(x2 - x1) >= abs(y2 - y1)
        return interpolate_time(x1, t1, x2, t2, xI)
    else
        return interpolate_time(y1, t1, y2, t2, yI)
    end
end

"""
    interpolate_position(x1, y1, t1, x2, y2, t2, t) -> (x, y)

Interpolate position (x, y) at time t given 2 positions (x1, y1) and (x2, y2) at time t1 and t2.

x(t) = (x2 - x1) / (t2 - t1) * (t - t1) + x1

"""
function interpolate_position(x1, y1, t1::T, x2, y2, t2::T, t::T) where {T<:Dates.AbstractDateTime}
    x = (x2 - x1) * (t - t1) / (t2 - t1) + x1
    y = (y2 - y1) * (t - t1) / (t2 - t1) + y1
    return x, y
end

"""
    calc_intersect_coordinates(x1, y1, x2, y2, AL) -> (xI, yI)

Calculate intersect coordinates given 2 positions coordinates (`x1`, `y1`) and (`x2`, `y2`)
and matrix line parameters `AL`.

``
\begin{cases}
α \\cdot x + β \\cdot y = 𝛾
α_L \\cdot x + β_L \\cdot y = 𝛾
\end{cases}
``

or with matrices

``
T \\cdot X = G
X = T^-1 \\cdot G
``

with

``
T =
\begin{bmatrix}
α & β \\
α_L & β_L
\end{bmatrix}
``

and

``
X =
\begin{bmatrix}
xI \\
yI
\end{bmatrix}
``

"""
function calc_intersect_coordinates(x1, y1, x2, y2, AL)
    A = matrix_line(x1, y1, x2, y2)
    T = [transpose(A); transpose(AL)]
    X = inv(T) * G
    xI, yI = X
    return xI, yI
end

"""
    find_intersections(x1L, y1L, x2L, y2L, points) -> TrackPoint[]

Find all intersections between track and a line defined by 2 points (`x1L`, `y1L`) and (`x2L`, `y2L`).
Tracks points are defined using a generator of tuples `(t, x, y)`` named `points`
where t is time (`DateTime` or `ZonedDateTime`) and position (`x`, `y`)
"""
function find_intersections(x1L, y1L, x2L, y2L, points)
    AL = matrix_line(x1L, y1L, x2L, y2L)
    x1, y1, x2, y2 = zero([x1L, y1L, x2L, y2L])
    ((t0, _, _), _) = iterate(points)
    t1 = t2 = t0
    I = TrackPoint(t1, x1, y1)
    lst_I = TrackPoint[]
    for (i, (t2, x2, y2)) in enumerate(points)
        if i > 1
            if (x1 != x2) || (y1 != y2)
                (xI, yI) = calc_intersect_coordinates(x1, y1, x2, y2, AL)
                if hasintersect(x1, y1, x2, y2, xI, yI)
                    tI = interpolate_time(x1, y1, t1, x2, y2, t2, xI, yI)
                    I = TrackPoint(tI, x1, y1)
                    push!(lst_I, I)
                end
            end
        end
        (x1, y1, t1) = (x2, y2, t2)
    end
    return lst_I
end

end # module


text:
using GalacticCMAEvolutionStrategy, GalacticOptim
using Test

@testset "GalacticCMAEvolutionStrategy.jl" begin
    rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
    x0 = zeros(2)
    _p = [1.0, 100.0]
    l1 = rosenbrock(x0, _p)
    f = OptimizationFunction(rosenbrock)
    prob = OptimizationProblem(f, x0, _p, lb=[-1.0, -1.0], ub=[0.8, 0.8])
    sol = solve(prob, CMAEvolutionStrategyOpt())
    @test 10 * sol.minimum < l1
end


text:
using WorkflowTest
using Test

@testset "WorkflowTest.jl" begin
    # Write your tests here.
end


text:
@testset "898.bitwise-ors-of-subarrays.jl" begin
    @test subarray_bitwise_ors([1,1,2]) == 3
    @test subarray_bitwise_ors([1,2,4]) == 6
    @test subarray_bitwise_ors([0]) == 1
end


text:
pretty_print(x::Signal) = string(value(x))
pretty_print(x::Vector) = string("[", first(x), "...", last(x), "]")
pretty_print(x::Matrix{T}) where {T} = string("Matrix{", T, "}")
pretty_print(x::DataType) = string(x.name.name)
pretty_print(x) = string(x)

print_type(x::Function) = string(x.env.name)
print_type(x::Symbol) = string(":", x)
function print_type(x::Type{T}) where T<:Tuple
    "($(join(map(print_type, x.parameters), ", ")))"
end
function print_type(x::TypeVar)
    string(x.name)
end
print_type(x::TypeName) = print_type(x.name)
print_type(x::UnionAll) = sprint(io-> show(io, x))
print_type(x) = print_type(x.name)
function print_type(x::Union)
    "Union{$(join(map(print_type, Base.uniontypes(x)), " "))}"
end

s_length(x::Symbol) = length(string(x))

const default_kwargs_docs = Dict(
    :position => "must be `Vector{Point}` or `Grid` or `nothing`",
    :position_x => "must be `Vector{Float}` or `Float` or `nothing`",
    :position_y => "must be `Vector{Float}` or `Float` or `nothing`",
    :position_z => "must be `Vector{Float}` or `Float` or `nothing`",

    :scale => "must be `Vector{Vec}` or `Vec`",
    :scale_x => "must be `Vector{Float}` or `Float` or `nothing`",
    :scale_y => "must be `Vector{Float}` or `Float` or `nothing`",
    :scale_z => "must be `Vector{Float}` or `Float` or `nothing`",
    :color => "Must be a Colorant",

    :model => "Must be`Mat4f0`. Scales, rotates and positions the object",
    :preferred_camera => "Preferred camera",
    :light => "Light to illuminate the model",
    :shading => "Accepts true or false. When true, uses blinn phong shading, otherwise flat"
)
function get_docs(arg1, style=:default)
    s = Style{style}()
    data = _default(arg1, s, Dict{Symbol, Any}(
        :shader => nothing # don't let _default create a shader, a window is necessary for that
    ))
    doc = get(data, :doc_string, Dict())
    targets = get(data, :gl_convert_targets, Dict())
    # we don't want to have these in the docs, so remove them
    delete!(data, :doc_string)
    delete!(data, :gl_convert_targets)

    maxlen = maximum(map(s_length, keys(data)))
    io = IOBuffer()
    for (k, value) in data
        spaces = maxlen - s_length(k)
        docstring = get(doc, k, get(default_kwargs_docs, k, ""))
        if docstring != ""
            docstring = docstring
        end
        t = string(get(targets, k, ""))
        if t != ""
            t = "can be gpu type: "*pretty_print(t)
        end
        println(io, " ", k, " "^spaces, " = ", pretty_print(value))
        println(io, docstring)
        println(io, t)
        println(io)
        println(io, "---")
        println(io)
    end
    seekstart(io)
    Markdown.parse(io)
end

function get_doc(dict, sig)
    for (k, v) in dict
        k == sig && return v.text[1]
    end
    ""
end
_tuple(x) = (x,)
_tuple(x::Core.SimpleVector) = tuple(x...)
_tuple(x::Tuple) = x
_tuple(x::Type{T}) where {T<:Tuple} = tuple(x.parameters...)

function all_docs(io = STDOUT)
    method_table = methods(GLVisualize._default)
    metadict = Docs.meta(GLVisualize)
    default_docs = metadict[Base.Docs.Binding(GLVisualize, :_default)].docs

    for method in method_table
        println(io, method)
        println(io)
        # sig_types = method.sig
        # if isa(sig_types, UnionAll)
        #     sig_types = Base.unwrap_unionall(sig_types)
        # end
        # sig = (sig_types.parameters[2:end]...,)
        # world = typemax(UInt)
        # x = first(Base._methods(GLVisualize._default, Tuple{sig...}, -1, world))
        # tvars = _tuple(x)[2:end]
        # sparam = if !isempty(tvars)
        #     string("{", join(map(x-> string(x), tvars), ", "), "}")
        # else
        #     ""
        # end
        # docstr = get_doc(default_docs, Tuple{sig...})
        # arg1 = print_type(sig[1])
        # arg2 = if isa(sig[2],  UnionAll) #
        #     ""
        # else
        #     string(", :", sig[2].parameters[1])
        # end
        # print(io, docstr)
        # println(io, "visualize", sparam, "(", arg1, arg2, ")")
        # println(io)
    end

end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

@testset "graph.disjointunion" begin

    G = pathgraph(5)
    H = pathgraph(5)
    I = cyclegraph(5)
    disjointunion!(G, H, I)
    @test nodecount(G) == 15
    @test edgecount(G) == 13
    (g, h, i) = connectedcomponents(G)
    @test length(g) == 5
    @test length(h) == 5
    @test length(i) == 5

    gsub = nodesubgraph(G, Set(1:5))
    U = disjointunion(gsub, H, I)
    @test nodecount(U) == 15
    @test edgecount(U) == 13
    @test length(connectedcomponents(U)) == 3
    @test getunionnode(U, 3, 1) == 11
    @test getunionedge(U, 3, 1) == 9
    @test getsourcenode(U, 15).source == 3
    @test getsourcenode(U, 15).sourcekey == 5
    @test getsourceedge(U, 13).source == 3
    @test getsourceedge(U, 13).sourcekey == 5

    k5 = completegraph(5)
    k5sub = nodesubgraph(k5, Set(1:4))
    U = disjointunion(k5sub, k5sub)
    @test nodecount(U) == 8
    @test edgecount(U) == 12
    @test length(connectedcomponents(U)) == 2

end # graph.disjointunion


text:
# Autogenerated wrapper script for Netgen_jll for armv6l-linux-musleabihf-cxx11
export inetcomp, netcomp, netgen, ntk2adl, ntk2xnf

JLLWrappers.@generate_wrapper_header("Netgen")
JLLWrappers.@declare_executable_product(inetcomp)
JLLWrappers.@declare_executable_product(netcomp)
JLLWrappers.@declare_executable_product(netgen)
JLLWrappers.@declare_executable_product(ntk2adl)
JLLWrappers.@declare_executable_product(ntk2xnf)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        inetcomp,
        "bin/inetcomp",
    )

    JLLWrappers.@init_executable_product(
        netcomp,
        "bin/netcomp",
    )

    JLLWrappers.@init_executable_product(
        netgen,
        "bin/netgen",
    )

    JLLWrappers.@init_executable_product(
        ntk2adl,
        "bin/ntk2adl",
    )

    JLLWrappers.@init_executable_product(
        ntk2xnf,
        "bin/ntk2xnf",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module RigidBodySimTest

using Compat
using Compat.Test
using Compat.Random
using Compat.LinearAlgebra

using RigidBodySim

using RigidBodyDynamics

using DiffEqCallbacks: DiscreteCallback
using DiffEqBase: add_tstop!
using OrdinaryDiffEq: Rodas4P

using MechanismGeometries
using MeshCatMechanisms
using InteractBase: observe
using Blink: Window

if VERSION < v"0.7-"
    const seed! = srand
else
    import Random: seed!
end

function dynamics_allocations(dynamics::Dynamics, state::MechanismState) # introduce function barrier
    x = Vector(state)
    ẋ = similar(x)
    p = nothing
    t = 3.
    dynamics(ẋ, x, p, t)
    allocs = @allocated dynamics(ẋ, x, p, t)
end

@testset "Dynamics" begin
    seed!(134)
    mechanism = rand_tree_mechanism(Float64, [Revolute{Float64} for i = 1 : 30]...)
    dynamics = Dynamics(mechanism)
    state = MechanismState(mechanism)
    rand!(state)
    @test dynamics_allocations(dynamics, state) <= 80
end

@testset "compare to simulate" begin
    seed!(1)
    urdf = joinpath(@__DIR__, "urdf", "Acrobot.urdf")
    mechanism = parse_urdf(Float64, urdf)
    state = MechanismState(mechanism)
    rand!(state)
    x0 = Vector(state)
    final_time = 5.
    problem = ODEProblem(Dynamics(mechanism), state, (0., final_time))
    sol = solve(problem, Vern7(), abs_tol = 1e-10, dt = 0.05)
    Compat.copyto!(state, x0)
    ts, qs, vs = RigidBodyDynamics.simulate(state, final_time)
    @test [qs[end]; vs[end]] ≈ sol[end] atol = 1e-2
end

@testset "renormalization callback" begin
    mechanism = rand_tree_mechanism(Float64, QuaternionFloating{Float64})
    floatingjoint = first(joints(mechanism))
    state = MechanismState(mechanism)

    rand!(configuration(state))
    @test !RigidBodyDynamics.is_configuration_normalized(floatingjoint, configuration(state, floatingjoint))

    problem = ODEProblem(Dynamics(mechanism), state, (0., 1e-3))
    sol = solve(problem, Vern7(), dt = 1e-4, callback = configuration_renormalizer(state))

    Compat.copyto!(state, sol[end])
    @test RigidBodyDynamics.is_configuration_normalized(floatingjoint, configuration(state, floatingjoint))
end

@testset "RealtimeRateLimiter" begin
    du = [0; 0]
    u0 = [0.; 0.]
    dynamics = (u, p, t) -> eltype(u).(du)
    tmin = 10.1
    tmax = 12.2
    for max_rate in [2.0, 0.5]
        prob = ODEProblem(dynamics, u0, (tmin, tmax))
        rate_limiter = RealtimeRateLimiter(max_rate = max_rate)
        sol = solve(prob, Tsit5(); callback = rate_limiter)
        soltime = @elapsed solve(prob, Tsit5(); callback = rate_limiter)
        expected = (tmax - tmin) / max_rate
        @show soltime
        @show expected
        @test soltime ≈ expected atol = 0.5
    end
end

@testset "PeriodicController" begin
    urdf = joinpath(@__DIR__, "urdf", "Acrobot.urdf")
    mechanism = parse_urdf(Float64, urdf)
    state = MechanismState(mechanism)
    controltimes = Float64[]
    initialize = (c, u, t, integrator) -> empty!(controltimes)
    τ = similar(velocity(state))
    Δt = 0.25

    make_controller = function ()
        PeriodicController(τ, Δt, function (τ, t, state)
            push!(controltimes, t)
            τ[1] = sin(t)
            τ[2] = cos(t)
        end; initialize = initialize)
    end

    controller = make_controller()
    final_time = 25.3
    problem = ODEProblem(Dynamics(mechanism, controller), state, (0., final_time))

    # ensure that controller gets called at appropriate times:
    sol = solve(problem, Tsit5(), abs_tol = 1e-10, dt = 0.05)
    @test controltimes == collect(0. : Δt : final_time - rem(final_time, Δt))

    # ensure that we can solve the same problem again without errors and with a different integrator
    empty!(controltimes)
    sol = solve(problem, Vern7(), abs_tol = 1e-10, dt = 0.05)
    @test controltimes == collect(0. : Δt : final_time - rem(final_time, Δt))

    # issue #60
    empty!(controltimes)
    problem60 = ODEProblem(Dynamics(mechanism, (τ, t, state) -> controller(τ, t, state)), state, (0., final_time))
    @test_throws RigidBodySim.Control.PeriodicControlFailure solve(problem60, Vern7(), abs_tol = 1e-10, dt = 0.05)
    controller = controller = make_controller()
    problem60 = ODEProblem(Dynamics(mechanism, (τ, t, state) -> controller(τ, t, state)), state, (0., final_time))
    @test_throws RigidBodySim.Control.PeriodicControlFailure solve(problem60, Vern7(), abs_tol = 1e-10, dt = 0.05)
    problem60_fixed = ODEProblem(Dynamics(mechanism, (τ, t, state) -> controller(τ, t, state)), state, (0., final_time),
        callback = PeriodicCallback(controller))
    sol60 = solve(problem60_fixed, Vern7(), abs_tol = 1e-10, dt = 0.05)
    @test controltimes == collect(0. : Δt : final_time - rem(final_time, Δt))
    @test sol60.t == sol.t
    @test sol60.u == sol.u
end

@testset "visualizer callbacks" begin
    mechanism = rand_tree_mechanism(Float64, [Revolute{Float64} for i = 1 : 30]...)
    state = MechanismState(mechanism)
    gui = GUI(mechanism, MeshCatMechanisms.Skeleton(randomize_colors = true, inertias = true))
    vis = gui.visualizer
    controls = gui.controls
    open(vis, Window())
    open(controls, Window())
    open(gui)
    dt = 1e-4

    set_configuration!(vis, configuration(state))
    @test Vector(vis.state) == Vector(state)

    tfinal = 0.5
    problem = ODEProblem(Dynamics(mechanism), state, (0., tfinal), callback=CallbackSet(gui))

    # Simulate without interaction
    sol = solve(problem, RK4(), adaptive = false, dt = dt)
    @test sol.t[end] == tfinal
    @test Vector(vis.state) == sol.u[end]

    # Test simulation controls
    tfinal = 100.0
    problem = ODEProblem(Dynamics(mechanism), state, (0., tfinal), callback=CallbackSet(gui))

    # Simulate for 3 seconds (wall time) and then send a termination command
    @async (sleep(3.); observe(controls.terminate)[] += 1)
    sol = solve(problem, RK4(), adaptive = false, dt = dt)
    @test sol.t[end] > 2 * dt
    @test sol.t[end] < tfinal
    println("last(sol.t) after early termination 1: $(last(sol.t))")

    # Rinse and repeat with the same ODEProblem (make sure that we don't terminate straight away)
    observe(controls.terminate)[] += 1
    sleep(0.1)
    @async (sleep(3.); observe(controls.terminate)[] += 1)
    sol = solve(problem, RK4(), adaptive = false, dt = dt)
    @test sol.t[end] > 2 * dt
    @test sol.t[end] < tfinal
    println("last(sol.t) after early termination 2: $(last(sol.t))")

    # Pause and unpause a short simulation, make sure that the simulation takes longer than without pausing
    problem = ODEProblem(Dynamics(mechanism), state, (0., 1.), callback=CallbackSet(gui))
    tpause = 0.5
    pausecondition = Condition()
    havepaused = Ref(false)
    condition = (u, t, integrator) -> !havepaused[] && t >= tpause
    action = function (integrator)
        observe(controls.pause)[] += 1
        notify(pausecondition)
        havepaused[] = true
    end
    initialize = (c, t, u, integrator) -> (havepaused[] = false; add_tstop!(integrator, tpause))
    pauser = DiscreteCallback(condition, action, save_positions=(false, false), initialize = initialize)
    integrator = init(problem, RK4(), adaptive = false, dt = dt, callback=pauser)
    @async begin
        wait(pausecondition)
        yield() # wait for message to reach command handler callback
        integrator_time = integrator.t
        @test integrator_time < tpause + 0.1 # it takes a while for the pause message to reach command handler callback
        sleep(2.)
        @test integrator.t == integrator_time # make sure simulation remains paused
        havepaused[] = true
        observe(controls.pause)[] += 1 # unpause
    end
    solve!(integrator)
    @test havepaused[]

    # Simulate for 3 seconds wall time, then pause, and then terminate a second later to make sure terminating works while paused
    tfinal = 100.0
    problem = ODEProblem(Dynamics(mechanism), state, (0., tfinal), callback=CallbackSet(gui))
    @async (sleep(3.); observe(controls.pause)[] += 1)
    @async (sleep(4.); observe(controls.terminate)[] += 1)
    sol = solve(problem, RK4(), adaptive = false, dt = dt)
    @test sol.t[end] > 2 * dt
    @test sol.t[end] < tfinal
    println("last(sol.t) after early termination 3: $(last(sol.t))")
end

@testset "ODESolution animation" begin
    mechanism = rand_tree_mechanism(Float64, [Revolute{Float64} for i = 1 : 30]...)
    state = MechanismState(mechanism)
    vis = MechanismVisualizer(mechanism, MeshCatMechanisms.Skeleton(randomize_colors = true, inertias = true))
    open(vis, Window())
    wait(vis)

    final_time = 5.
    problem = ODEProblem(Dynamics(mechanism), state, (0., final_time))
    sol = solve(problem, Vern7(), abs_tol = 1e-10, dt = 0.05)

    # regular playback
    realtime_rate = 2.
    setanimation!(vis, sol, realtime_rate = realtime_rate, max_fps = 60.)
    sleep(final_time / realtime_rate)

    # broken:
    #=
    setanimation!(vis, sol, realtime_rate = 1000.)
    realtime_rate = 2.
    elapsed = @elapsed setanimation!(vis, sol, realtime_rate = realtime_rate, max_fps = 60.)
    @test elapsed ≈ final_time / realtime_rate atol = 0.1

    # premature termination
    termination_time = 1.5
    @async (sleep(termination_time); send_terminate_message())
    elapsed = @elapsed animate(vis, state, sol, realtime_rate = realtime_rate)
    @test elapsed ≈ termination_time atol = 0.1

    # pause and unpause
    pause_time = 1.0
    unpause_time = 3.5
    @async (sleep(pause_time); send_pause_message())
    @async (sleep(unpause_time); send_pause_message())
    elapsed = @elapsed animate(vis, state, sol, realtime_rate = realtime_rate)
    @show elapsed
    @test elapsed ≈ final_time / realtime_rate + (unpause_time - pause_time) atol = 0.2 # higher atol because of pause poll int

    # pause and terminate
    @async (sleep(pause_time); send_pause_message())
    @async (sleep(termination_time); send_terminate_message())
    elapsed = @elapsed animate(vis, state, sol, realtime_rate = realtime_rate)
    @show elapsed
    @test elapsed ≈ termination_time atol = 0.2 # higher atol because of pause poll int
    =#
end

@testset "Stiff integrator" begin
    urdf = joinpath(@__DIR__, "urdf", "Acrobot.urdf")
    mechanism = parse_urdf(Float64, urdf)
    state = MechanismState(mechanism)
    seed!(1)
    rand!(state)
    x0 = Vector(state)
    final_time = 1.
    problem = ODEProblem(Dynamics(mechanism), state, (0., final_time))

    sol_nonstiff = solve(problem, Vern7(), abs_tol = 1e-10, dt = 0.05)
    sol_stiff = solve(problem, Rodas4P(), abs_tol = 1e-10, dt = 0.05)
    @test last(sol_nonstiff) ≈ last(sol_stiff) atol = 1e-2
end

# notebooks
include("test_notebooks.jl")

end


text:
using ITensors
using Test

using ChainRulesCore: rrule_via_ad

include("utils/chainrulestestutils.jl")

using Zygote: ZygoteRuleConfig, gradient

@testset "ChainRules rrules: Ops" begin
  s = siteinds("S=1/2", 4)

  x = 2.4
  V = randomITensor(s[1], s[2])

  f = function (x)
    y = ITensor(Op("Ry", 1; θ=x), s)
    return y[1, 1]
  end
  args = (x,)
  test_rrule(
    ZygoteRuleConfig(),
    f,
    args...;
    rrule_f=rrule_via_ad,
    check_inferred=false,
    rtol=1.0e-7,
    atol=1.0e-7,
  )

  f = function (x)
    y = exp(ITensor(Op("Ry", 1; θ=x), s))
    return y[1, 1]
  end
  args = (x,)
  test_rrule(
    ZygoteRuleConfig(),
    f,
    args...;
    rrule_f=rrule_via_ad,
    check_inferred=false,
    rtol=1.0e-7,
    atol=1.0e-7,
  )

  f = function (x)
    y = Op("Ry", 1; θ=x) + Op("Ry", 1; θ=x)
    return y[1].params.θ
  end
  args = (x,)
  test_rrule(
    ZygoteRuleConfig(),
    f,
    args...;
    rrule_f=rrule_via_ad,
    check_inferred=false,
    rtol=1.0e-7,
    atol=1.0e-7,
  )

  f = function (x)
    y = ITensor(Op("Ry", 1; θ=x) + Op("Ry", 1; θ=x), s)
    return y[1, 1]
  end
  args = (x,)
  test_rrule(
    ZygoteRuleConfig(),
    f,
    args...;
    rrule_f=rrule_via_ad,
    check_inferred=false,
    rtol=1.0e-7,
    atol=1.0e-7,
  )

  f = function (x)
    y = exp(ITensor(Op("Ry", 1; θ=x) + Op("Ry", 1; θ=x), s))
    return y[1, 1]
  end
  args = (x,)
  test_rrule(
    ZygoteRuleConfig(),
    f,
    args...;
    rrule_f=rrule_via_ad,
    check_inferred=false,
    rtol=1.0e-7,
    atol=1.0e-7,
  )

  f = function (x)
    y = ITensor(exp(Op("Ry", 1; θ=x) + Op("Ry", 1; θ=x)), s)
    return y[1, 1]
  end
  args = (x,)
  test_rrule(
    ZygoteRuleConfig(),
    f,
    args...;
    rrule_f=rrule_via_ad,
    check_inferred=false,
    rtol=1.0e-7,
    atol=1.0e-7,
  )

  f = function (x)
    y = ITensor(2 * Op("Ry", 1; θ=x), s)
    return y[1, 1]
  end
  args = (x,)
  test_rrule(
    ZygoteRuleConfig(),
    f,
    args...;
    rrule_f=rrule_via_ad,
    check_inferred=false,
    rtol=1.0e-7,
    atol=1.0e-7,
  )

  f = function (x)
    y = ITensor(2 * (Op("Ry", 1; θ=x) + Op("Ry", 1; θ=x)), s)
    return y[1, 1]
  end
  args = (x,)
  test_rrule(
    ZygoteRuleConfig(),
    f,
    args...;
    rrule_f=rrule_via_ad,
    check_inferred=false,
    rtol=1.0e-7,
    atol=1.0e-7,
  )

  f = function (x)
    y = ITensor(Op("Ry", 1; θ=x) * Op("Ry", 2; θ=x), s)
    return y[1, 1]
  end
  args = (x,)
  test_rrule(
    ZygoteRuleConfig(),
    f,
    args...;
    rrule_f=rrule_via_ad,
    check_inferred=false,
    rtol=1.0e-7,
    atol=1.0e-7,
  )

  f = function (x)
    y = ITensor(exp(-x * Op("X", 1) * Op("X", 2)), s)
    return norm(y)
  end
  args = (x,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)

  f = function (x)
    y = exp(-x * Op("X", 1) * Op("X", 2))
    y *= exp(-x * Op("X", 1) * Op("X", 2))
    U = Prod{ITensor}(y, s)
    return norm(U(V))
  end
  args = (x,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)

  f = function (x)
    y = exp(-x * (Op("X", 1) + Op("Z", 1) + Op("Z", 1)); alg=Trotter{1}(1))
    U = Prod{ITensor}(y, s)
    return norm(U(V))
  end
  args = (x,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)

  ## ## XXX: Error in vcat!
  ## f = function (x)
  ##   y = -x * (Op("X", 1) * Op("X", 2) + Op("Z", 1) * Op("Z", 2))
  ##   U = ITensor(y, s)
  ##   return norm(U * V)
  ## end
  ##
  ## ## XXX: Error in vcat!
  ## f = function (x)
  ##   y = exp(-x * (Op("X", 1) * Op("X", 2) + Op("Z", 1) * Op("Z", 2)); alg=Trotter{1}(1))
  ##   U = ITensor(y, s)
  ##   return norm(U * V)
  ## end
end


text:
function loopLayout!(loop, nTurns, n, prevX, nxtX, prevY, side; yOffset = 0.1, expand=0)
	yOffset=0.1
    rem = (n - nTurns*2) % nTurns
	turnSize = div((n-rem-nTurns*2), nTurns)
    lftleg = 1:div(nTurns*2,2)
    rgtleg = (n-div(nTurns*2,2)+1):n
    turns = []
    position = lftleg[end] + 1

    for i in 1:nTurns
        if rem != 0
            push!(turns, position:position+turnSize)
            position += turnSize+1
            rem -= 1
        else
            push!(turns, position:position+turnSize-1)
            position += turnSize
        end
    end

    if nTurns < 5
        dX = (nxtX-prevX)/(nTurns-1)
        if side == 1
            expand = 1
        end
    else
        dX = (nxtX-prevX)/5
    end
    dY = 0.1

    ## position legs
	for i in lftleg
        loop[i, :x] = prevX
        loop[i, :y] = prevY + side*(yOffset + i * dY)
	end

    loop[rgtleg, :x] .= nxtX

    ## check for loop expand
    if expand > 0
        dX2 = 0.2
        for i in 0:expand
            if expand-i == expand
                loop[collect(lftleg)[end-i], :y] = loop[collect(lftleg)[end]-(expand+1), :y] + side*dY
                loop[collect(lftleg)[end-i], :x] -= (expand-i) * dX2
            else
                loop[collect(lftleg)[end-i], :y] = loop[collect(lftleg)[end]-(expand+1), :y]
                loop[collect(lftleg)[end-i], :x] -= (expand-(i-1)) * dX2
            end

            if i == 1
                loop[collect(rgtleg)[i], :x] = nxtX + expand*dX2
            end
            loop[collect(rgtleg)[i+2], :x] = nxtX + (expand-i)*dX2
        end

        prevX -= expand * dX2
        nxtX += expand * dX2
        dX = (nxtX-prevX)/(nTurns-1)
    end

    loop[rgtleg, :y] = reverse(loop[lftleg, :y])

    ## position first turn
    for i in turns[1]
        loop[i, :x] = prevX
        loop[i, :y] = prevY + side*(yOffset + (i-expand)*dY)
    end

    for i in 2:length(turns)
        loop[turns[i], :x] .= prevX + dX*(i-1)
        if i%2 == 0
            loop[turns[i], :y] = reverse(loop[turns[1],:y])[end-length(turns[i])+1:end]
        else
            loop[turns[i], :y] = loop[turns[1],:y][1:length(turns[i])]
        end
    end

    xOffset = (nxtX - prevX)/((nTurns-1)*2)
    for i in 2:length(turns)
        if loop[turns[i][1], :y] != loop[turns[i-1][end], :y]
            loop[turns[i-1][end], :x] += xOffset
            loop[turns[i-1][end], :y] += side*0.5*dY*(-1^(i+1))
        else
            loop[turns[i-1][end], :x] += xOffset/2
            loop[turns[i][1], :x] -= xOffset/2
        end
    end
end


text:
import DataFrames

"""
"""
function transform_columns! end

function transform_columns!(
        df::DataFrames.AbstractDataFrame,
        f::Function,
        column_name::Symbol,
        )::Nothing
    old_column = df[column_name]
    DataFrames.deletecols!(df, column_name,)
    new_column = fix_type(f(old_column))
    df[column_name] = new_column
    return nothing
end

function transform_columns!(
        df::DataFrames.AbstractDataFrame,
        f::Function,
        column_names::AbstractArray{Symbol},
        )::Nothing
    for column_name in column_names
        transform_columns!(df, f, column_name)
    end
    return nothing
end


text:
"""
    hamming(a::DNASequence, b::DNASequence)

Compute number of mismatches between two DNA sequences of equal length
"""
function Distances.hamming(a::DNASequence, b::DNASequence)
    n = length(a)
    (n != length(b)) && throw(DimensionMismatch("Lengths of sequences must be equal"))
    mismatches = 0
    for i in 1:n
        if a[i] != b[i]
            mismatches += 1
        end
    end
    mismatches
end

"""
    gen_barcode_mapping(barcodes::DataFrame, len::Int)

Generate a dictionary mapping of sequence of type Bio.Seq.DNASequence to their
string id from `barcodes`. Truncate all barcodes to length `len`. The final
result is not guaranteed to be unique due to the truncation step.
"""
function gen_barcode_mapping(barcodes::DataFrame, len::Int)
    mapping = Dict{DNASequence, String}()
    for row in eachrow(barcodes)
        mapping[DNASequence(row[:sequence][1:len])] = row[:id]
    end
    mapping
end

_try_get(x) = x
_try_get(x::Nullable) = get(x)


text:
################################################################################
#
#  Types and memory management for arf_interval
#
################################################################################

struct arf_interval_struct
    a::arf_t
    b::arf_t

    arf_interval_struct() = new(arf(), arf())
    arf_interval_struct(a::arf_t, b::arf_t) = new(a, b)
end

const arf_interval_t = arf_interval_struct
const arf_interval_ptr = Vector{arf_interval_t}
const arf_interval_srcptr = arf_interval_ptr

arf_interval() = arf_interval_struct()
arf_interval(a::arf_t, b::arf_t) = arf_interval_struct(a, b)

################################################################################
#
#  Conversions
#
################################################################################

convert(::Type{arf_interval_struct}, x) = arf_interval(convert(arf_struct, x[1]), convert(arf_struct, x[2]))
convert(::Type{arf_interval_struct}, x::arf_interval_struct) = x

convert(::Type{Tuple{T,S}}, x::arf_interval_struct) where {T,S}   = (convert(T, x.a), convert(S, x.b))
convert(::Type{Vector{T}}, x::arf_interval_struct) where {T}      = [convert(T, x.a), convert(T, x.b)]

convert(::Type{Tuple{Fslong,Fslong}}, x::arf_interval_struct)     = (convert(Fslong, x.a, ARF_RND_DOWN), convert(Fslong, x.b, ARF_RND_UP))
convert(::Type{Vector{Fslong}}, x::arf_interval_struct)           = [convert(Fslong, x.a, ARF_RND_DOWN), convert(Fslong, x.b, ARF_RND_UP)]

convert(::Type{Tuple{Cdouble,Cdouble}}, x::arf_interval_struct)   = (convert(Cdouble, x.a, ARF_RND_DOWN), convert(Cdouble, x.b, ARF_RND_UP))
convert(::Type{Vector{Cdouble}}, x::arf_interval_struct)          = [convert(Cdouble, x.a, ARF_RND_DOWN), convert(Cdouble, x.b, ARF_RND_UP)]

convert(::Type{Tuple{BigFloat,BigFloat}}, x::arf_interval_struct) = (convert(BigFloat, x.a, MPFR_RNDD), convert(BigFloat, x.b, MPFR_RNDU))
convert(::Type{Vector{BigFloat}}, x::arf_interval_struct)         = [convert(BigFloat, x.a, MPFR_RNDD), convert(BigFloat, x.b, MPFR_RNDU)]

convert(::Type{Tuple{BigInt,BigInt}}, x::arf_interval_struct)     = (floor(BigInt, convert(BigFloat, x.a, MPFR_RNDD)), ceil(BigInt, convert(BigFloat, x.b, MPFR_RNDU)))
convert(::Type{Vector{BigInt}}, x::arf_interval_struct)           = [floor(BigInt, convert(BigFloat, x.a, MPFR_RNDD)), ceil(BigInt, convert(BigFloat, x.b, MPFR_RNDU))]

################################################################################
#
#  String I/O
#
################################################################################

arf_interval_printd(v::arf_interval_t, n) = begin
    print("[")
    arf_printd(v.a, n)
    print(", ")
    arf_printd(v.b, n)
    print("]")
end

show(io::IO, x::arf_interval_t) = begin
    prec = 256
    d = ceil(Int, prec * 0.30102999566398119521)
    arf_interval_printd(x, d)
end

################################################################################
#
#  Unsafe setting
#
################################################################################

for typeofv in (arf_interval_t, Ptr{arf_interval_t})
    @eval begin
        arf_interval_set(v::($typeofv), u::arf_interval_t)  = (arf_set(v.a, u.a);  arf_set(v.b, u.b))
        arf_interval_swap(v::($typeofv), u::arf_interval_t) = (arf_swap(v.a, u.a); arf_swap(v.b, u.b))
    end
end

arf_interval_get_arb(x::arb_t, v::arf_interval_t, prec) = arb_set_interval_arf(x, v.a, v.b, prec)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "MariaDB_Connector_C"
version = v"3.1.6"

# Collection of sources required to build MariaDB_Connector_C
sources = [
    "https://github.com/mariadb-corporation/mariadb-connector-c.git" =>
    "63df45ce3df3fbc04d8fab9bceb77f9d1cccd4aa",
    "./bundled"
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/mariadb-*/

if [[ "${target}" == *-mingw* ]]; then
    for p in ../patches/*.patch; do
        atomic_patch -p1 "${p}"
    done
    export CFLAGS="-std=c99"
    # Minimum version of Windows supported by MariaDB is 7,
    # see tables in https://docs.microsoft.com/en-us/windows/win32/winprog/using-the-windows-headers
    if [[ "${nbits}" == 32 ]]; then
        export CFLAGS="${CFLAGS} -D_WIN32_WINNT=0x0601"
    elif [[ "${nbits}" == 64 ]]; then
        export CFLAGS="${CFLAGS} -DNTDDI_VERSION=0x06010000"
    fi
elif [[ "${target}" == *-apple-* ]]; then
    # Nuke header files in the sysroot
    # to avoid they'll be picked up by CMake
    rm -rf /opt/${target}/${target}/sys-root/usr/include/openssl
    rm -rf /opt/${target}/${target}/sys-root/usr/include/{iconv.h,libcharset.h,localcharset.h}

    SYMBS_DEFS=()
    for sym in iconv iconv_close iconv_open; do
        SYMBS_DEFS+=(-D${sym}=lib${sym})
    done
    export CFLAGS="${SYMBS_DEFS[@]}"
fi

mkdir build && cd build
cmake .. -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DWITH_MYSQLCOMPAT=OFF \
    -DWITH_EXTERNAL_ZLIB=ON \
    -DZLIB_FOUND=ON \
    -DZLIB_INCLUDE_DIR=${prefix}/include \
    -DZLIB_LIBRARY=${libdir}/libz.${dlext} \
    -DOPENSSL_FOUND=ON \
    -DOPENSSL_CRYPTO_LIBRARY=${libdir}/libcrypto.${dlext} \
    -DOPENSSL_SSL_LIBRARY=${libdir}/libssl.${dlext} \
    -DICONV_LIBRARIES=${libdir}/libiconv.${dlext} \
    -DICONV_INCLUDE_DIR=${prefix}/include
make -j${nproc}
make install
install_license ../COPYING.LIB
"""

# MariaDB doesn't support non *86 platforms with Musl, see
# https://gitlab.alpinelinux.org/alpine/aports/commit/ad7c54bd9b6d8e60e88a313a92c83d083f196db8
platforms = supported_platforms(; exclude = [Platform("aarch64", "linux"; libc="musl"),
                                             Platform("armv7l", "linux"; libc="musl")])

# The products that we will ensure are always built
products = [
    LibraryProduct("libmariadb", :libmariadb, ["lib/mariadb"]),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "LibCURL_jll",
    "Libiconv_jll",
    "OpenSSL_jll",
    "Zlib_jll",
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using TermoQuimica
using Test

@testset "TermoQuimica.jl" begin
    # Write your tests here.
end

max_stars_repo_name: EmilioAlvizo/TermoQuimica

max_issues_repo_name: EmilioAlvizo/TermoQuimica

max_forks_repo_name: EmilioAlvizo/TermoQuimica


text:
using RayTracing
using Gridap
using Test

jsonfile = joinpath(@__DIR__,"../demo/pincell.json")
model = DiscreteModelFromFile(jsonfile)

@testset "Principal tests" begin

    tg = TrackGenerator(model, 8, 0.02)
    trace!(tg)
    segmentize!(tg)

    @testset "Tracing" begin
        @test tg.n_total_tracks == 420
        @test tg.n_tracks_x == [31, 74, 74, 31]
        @test tg.n_tracks_y == [74, 31, 31, 74]
        @test tg.n_tracks == [105, 105, 105, 105]
    end

    @testset "Azimuthal quadradure" begin
        @test isequal(RayTracing.nazim(tg.azimuthal_quadrature), 8)
        @test isequal(RayTracing.nazim2(tg.azimuthal_quadrature), 4)
        @test isequal(RayTracing.nazim4(tg.azimuthal_quadrature), 2)
        @test isapprox(tg.azimuthal_quadrature.δ, 0.02)
        @test all(isapprox.(tg.azimuthal_quadrature.δs, 0.01994243696980254))
        @test tg.azimuthal_quadrature.ϕs ≈ [0.39670866289121387, 1.1740876639036828, 1.9675049896861103, 2.7448839906985794]
    end

    @testset "Entry and exit points" begin
        for track in tg.tracks_by_uid
            @test isapprox(track.p, track.segments[begin].p)
            @test isapprox(track.q, track.segments[end].q)
        end
    end

    @testset "Track length" begin
        for track in tg.tracks_by_uid
            l1 = track.ℓ
            l2 = sum(RayTracing.ℓ.(track.segments))
            @test isapprox(l1, l2)
        end
    end
end

@testset "Reflection tests" begin

    bcs = BoundaryConditions(top=Vaccum, bottom=Reflective, left=Reflective, right=Vaccum)
    tg = TrackGenerator(model, 4, 0.8; bcs=bcs)
    trace!(tg)

    @testset "Track 1" begin
        track = tg.tracks_by_uid[1]
        next_track_fwd = track.next_track_fwd
        next_track_bwd = track.next_track_bwd
        DirNextTrackFwd = RayTracing.dir_next_track_fwd(track)
        DirNextTrackBwd = RayTracing.dir_next_track_bwd(track)
        @test RayTracing.bc_fwd(track) == Vaccum
        @test RayTracing.bc_bwd(track) == Reflective
        @test next_track_fwd.uid == 7
        @test next_track_bwd.uid == 6
        @test DirNextTrackFwd == RayTracing.Forward
        @test DirNextTrackBwd == RayTracing.Forward
    end

    @testset "Track 2" begin
        track = tg.tracks_by_uid[2]
        next_track_fwd = track.next_track_fwd
        next_track_bwd = track.next_track_bwd
        DirNextTrackFwd = RayTracing.dir_next_track_fwd(track)
        DirNextTrackBwd = RayTracing.dir_next_track_bwd(track)
        @test RayTracing.bc_fwd(track) == Vaccum
        @test RayTracing.bc_bwd(track) == Reflective
        @test next_track_fwd.uid == 8
        @test next_track_bwd.uid == 5
        @test DirNextTrackFwd == RayTracing.Forward
        @test DirNextTrackBwd == RayTracing.Forward
    end

    @testset "Track 3" begin
        track = tg.tracks_by_uid[3]
        next_track_fwd = track.next_track_fwd
        next_track_bwd = track.next_track_bwd
        DirNextTrackFwd = RayTracing.dir_next_track_fwd(track)
        DirNextTrackBwd = RayTracing.dir_next_track_bwd(track)
        @test RayTracing.bc_fwd(track) == Vaccum
        @test RayTracing.bc_bwd(track) == Reflective
        @test next_track_fwd.uid == 8
        @test next_track_bwd.uid == 5
        @test DirNextTrackFwd == RayTracing.Backward
        @test DirNextTrackBwd == RayTracing.Backward
    end

    @testset "Track 4" begin
        track = tg.tracks_by_uid[4]
        next_track_fwd = track.next_track_fwd
        next_track_bwd = track.next_track_bwd
        DirNextTrackFwd = RayTracing.dir_next_track_fwd(track)
        DirNextTrackBwd = RayTracing.dir_next_track_bwd(track)
        @test RayTracing.bc_fwd(track) == Vaccum
        @test RayTracing.bc_bwd(track) == Reflective
        @test next_track_fwd.uid == 7
        @test next_track_bwd.uid == 6
        @test DirNextTrackFwd == RayTracing.Backward
        @test DirNextTrackBwd == RayTracing.Backward
    end

    @testset "Track 5" begin
        track = tg.tracks_by_uid[5]
        next_track_fwd = track.next_track_fwd
        next_track_bwd = track.next_track_bwd
        DirNextTrackFwd = RayTracing.dir_next_track_fwd(track)
        DirNextTrackBwd = RayTracing.dir_next_track_bwd(track)
        @test RayTracing.bc_fwd(track) == Reflective
        @test RayTracing.bc_bwd(track) == Reflective
        @test next_track_fwd.uid == 3
        @test next_track_bwd.uid == 2
        @test DirNextTrackFwd == RayTracing.Forward
        @test DirNextTrackBwd == RayTracing.Forward
    end

    @testset "Track 6" begin
        track = tg.tracks_by_uid[6]
        next_track_fwd = track.next_track_fwd
        next_track_bwd = track.next_track_bwd
        DirNextTrackFwd = RayTracing.dir_next_track_fwd(track)
        DirNextTrackBwd = RayTracing.dir_next_track_bwd(track)
        @test RayTracing.bc_fwd(track) == Reflective
        @test RayTracing.bc_bwd(track) == Reflective
        @test next_track_fwd.uid == 4
        @test next_track_bwd.uid == 1
        @test DirNextTrackFwd == RayTracing.Forward
        @test DirNextTrackBwd == RayTracing.Forward
    end

    @testset "Track 7" begin
        track = tg.tracks_by_uid[7]
        next_track_fwd = track.next_track_fwd
        next_track_bwd = track.next_track_bwd
        DirNextTrackFwd = RayTracing.dir_next_track_fwd(track)
        DirNextTrackBwd = RayTracing.dir_next_track_bwd(track)
        @test RayTracing.bc_fwd(track) == Vaccum
        @test RayTracing.bc_bwd(track) == Vaccum
        @test next_track_fwd.uid == 4
        @test next_track_bwd.uid == 1
        @test DirNextTrackFwd == RayTracing.Backward
        @test DirNextTrackBwd == RayTracing.Backward
    end

    @testset "Track 8" begin
        track = tg.tracks_by_uid[8]
        next_track_fwd = track.next_track_fwd
        next_track_bwd = track.next_track_bwd
        DirNextTrackFwd = RayTracing.dir_next_track_fwd(track)
        DirNextTrackBwd = RayTracing.dir_next_track_bwd(track)
        @test RayTracing.bc_fwd(track) == Vaccum
        @test RayTracing.bc_bwd(track) == Vaccum
        @test next_track_fwd.uid == 3
        @test next_track_bwd.uid == 2
        @test DirNextTrackFwd == RayTracing.Backward
        @test DirNextTrackBwd == RayTracing.Backward
    end

    # a more complicated case
    tg = TrackGenerator(model, 8, 0.8; bcs=bcs)
    trace!(tg)

    @testset "Track 1" begin
        track = tg.tracks_by_uid[1]
        next_track_fwd = track.next_track_fwd
        next_track_bwd = track.next_track_bwd
        DirNextTrackFwd = RayTracing.dir_next_track_fwd(track)
        DirNextTrackBwd = RayTracing.dir_next_track_bwd(track)
        @test RayTracing.bc_fwd(track) == Vaccum
        @test RayTracing.bc_bwd(track) == Reflective
        @test next_track_fwd.uid == 11
        @test next_track_bwd.uid == 10
        @test DirNextTrackFwd == RayTracing.Forward
        @test DirNextTrackBwd == RayTracing.Forward
    end

    @testset "Track 2" begin
        track = tg.tracks_by_uid[2]
        next_track_fwd = track.next_track_fwd
        next_track_bwd = track.next_track_bwd
        DirNextTrackFwd = RayTracing.dir_next_track_fwd(track)
        DirNextTrackBwd = RayTracing.dir_next_track_bwd(track)
        @test RayTracing.bc_fwd(track) == Vaccum
        @test RayTracing.bc_bwd(track) == Reflective
        @test next_track_fwd.uid == 12
        @test next_track_bwd.uid == 10
        @test DirNextTrackFwd == RayTracing.Forward
        @test DirNextTrackBwd == RayTracing.Backward
    end

    @testset "Track 3" begin
        track = tg.tracks_by_uid[3]
        next_track_fwd = track.next_track_fwd
        next_track_bwd = track.next_track_bwd
        DirNextTrackFwd = RayTracing.dir_next_track_fwd(track)
        DirNextTrackBwd = RayTracing.dir_next_track_bwd(track)
        @test RayTracing.bc_fwd(track) == Vaccum
        @test RayTracing.bc_bwd(track) == Reflective
        @test next_track_fwd.uid == 12
        @test next_track_bwd.uid == 11
        @test DirNextTrackFwd == RayTracing.Backward
        @test DirNextTrackBwd == RayTracing.Backward
    end

    @testset "Track 4" begin
        track = tg.tracks_by_uid[4]
        next_track_fwd = track.next_track_fwd
        next_track_bwd = track.next_track_bwd
        DirNextTrackFwd = RayTracing.dir_next_track_fwd(track)
        DirNextTrackBwd = RayTracing.dir_next_track_bwd(track)
        @test RayTracing.bc_fwd(track) == Vaccum
        @test RayTracing.bc_bwd(track) == Reflective
        @test next_track_fwd.uid == 9
        @test next_track_bwd.uid == 8
        @test DirNextTrackFwd == RayTracing.Forward
        @test DirNextTrackBwd == RayTracing.Forward
    end

    @testset "Track 5" begin
        track = tg.tracks_by_uid[5]
        next_track_fwd = track.next_track_fwd
        next_track_bwd = track.next_track_bwd
        DirNextTrackFwd = RayTracing.dir_next_track_fwd(track)
        DirNextTrackBwd = RayTracing.dir_next_track_bwd(track)
        @test RayTracing.bc_fwd(track) == Vaccum
        @test RayTracing.bc_bwd(track) == Reflective
        @test next_track_fwd.uid == 9
        @test next_track_bwd.uid == 7
        @test DirNextTrackFwd == RayTracing.Backward
        @test DirNextTrackBwd == RayTracing.Forward
    end

    @testset "Track 6" begin
        track = tg.tracks_by_uid[6]
        next_track_fwd = track.next_track_fwd
        next_track_bwd = track.next_track_bwd
        DirNextTrackFwd = RayTracing.dir_next_track_fwd(track)
        DirNextTrackBwd = RayTracing.dir_next_track_bwd(track)
        @test RayTracing.bc_fwd(track) == Vaccum
        @test RayTracing.bc_bwd(track) == Reflective
        @test next_track_fwd.uid == 8
        @test next_track_bwd.uid == 7
        @test DirNextTrackFwd == RayTracing.Backward
        @test DirNextTrackBwd == RayTracing.Backward
    end

    @testset "Track 7" begin
        track = tg.tracks_by_uid[7]
        next_track_fwd = track.next_track_fwd
        next_track_bwd = track.next_track_bwd
        DirNextTrackFwd = RayTracing.dir_next_track_fwd(track)
        DirNextTrackBwd = RayTracing.dir_next_track_bwd(track)
        @test RayTracing.bc_fwd(track) == Reflective
        @test RayTracing.bc_bwd(track) == Reflective
        @test next_track_fwd.uid == 6
        @test next_track_bwd.uid == 5
        @test DirNextTrackFwd == RayTracing.Forward
        @test DirNextTrackBwd == RayTracing.Forward
    end

    @testset "Track 8" begin
        track = tg.tracks_by_uid[8]
        next_track_fwd = track.next_track_fwd
        next_track_bwd = track.next_track_bwd
        DirNextTrackFwd = RayTracing.dir_next_track_fwd(track)
        DirNextTrackBwd = RayTracing.dir_next_track_bwd(track)
        @test RayTracing.bc_fwd(track) == Vaccum
        @test RayTracing.bc_bwd(track) == Reflective
        @test next_track_fwd.uid == 6
        @test next_track_bwd.uid == 4
        @test DirNextTrackFwd == RayTracing.Backward
        @test DirNextTrackBwd == RayTracing.Forward
    end

    @testset "Track 9" begin
        track = tg.tracks_by_uid[9]
        next_track_fwd = track.next_track_fwd
        next_track_bwd = track.next_track_bwd
        DirNextTrackFwd = RayTracing.dir_next_track_fwd(track)
        DirNextTrackBwd = RayTracing.dir_next_track_bwd(track)
        @test RayTracing.bc_fwd(track) == Vaccum
        @test RayTracing.bc_bwd(track) == Vaccum
        @test next_track_fwd.uid == 5
        @test next_track_bwd.uid == 4
        @test DirNextTrackFwd == RayTracing.Backward
        @test DirNextTrackBwd == RayTracing.Backward
    end

    @testset "Track 10" begin
        track = tg.tracks_by_uid[10]
        next_track_fwd = track.next_track_fwd
        next_track_bwd = track.next_track_bwd
        DirNextTrackFwd = RayTracing.dir_next_track_fwd(track)
        DirNextTrackBwd = RayTracing.dir_next_track_bwd(track)
        @test RayTracing.bc_fwd(track) == Reflective
        @test RayTracing.bc_bwd(track) == Reflective
        @test next_track_fwd.uid == 2
        @test next_track_bwd.uid == 1
        @test DirNextTrackFwd == RayTracing.Forward
        @test DirNextTrackBwd == RayTracing.Forward
    end

    @testset "Track 11" begin
        track = tg.tracks_by_uid[11]
        next_track_fwd = track.next_track_fwd
        next_track_bwd = track.next_track_bwd
        DirNextTrackFwd = RayTracing.dir_next_track_fwd(track)
        DirNextTrackBwd = RayTracing.dir_next_track_bwd(track)
        @test RayTracing.bc_fwd(track) == Reflective
        @test RayTracing.bc_bwd(track) == Vaccum
        @test next_track_fwd.uid == 3
        @test next_track_bwd.uid == 1
        @test DirNextTrackFwd == RayTracing.Forward
        @test DirNextTrackBwd == RayTracing.Backward
    end

    @testset "Track 12" begin
        track = tg.tracks_by_uid[12]
        next_track_fwd = track.next_track_fwd
        next_track_bwd = track.next_track_bwd
        DirNextTrackFwd = RayTracing.dir_next_track_fwd(track)
        DirNextTrackBwd = RayTracing.dir_next_track_bwd(track)
        @test RayTracing.bc_fwd(track) == Vaccum
        @test RayTracing.bc_bwd(track) == Vaccum
        @test next_track_fwd.uid == 3
        @test next_track_bwd.uid == 2
        @test DirNextTrackFwd == RayTracing.Backward
        @test DirNextTrackBwd == RayTracing.Backward
    end
end


text:
#!/usr/bin/env julia --startup-file=no
# Script for converting BATSRUS unstructured binary outputs to VTK formats.
# Example usage:
# julia -p 8 convert2vtk.jl
#
# Currently the relative location of the output files to be processed is "IO2".
#
# Hongyang Zhou, hyzhou@umich.edu

using Distributed
using Pkg
if "VisAna" ∈ keys(Pkg.installed())
   @everywhere using Batsrus
else
   @warn "Batsrus.jl not installed. Install the package by `Pkg.add("Batsrus")`"
end
@everywhere using Glob

dir = "IO2"

# Concatenate files
# Before deleting the .T
# you need to make sure all the .tec files have been concatenated!
headers = glob("*.T", dir)
@sync @distributed for header in headers
   name = basename(header)[1:end-2]
   println("header=$(header)")
   nameIn = glob(name*"*.tec", dir)
   nameOut= joinpath(dir, name*".dat")
   run(pipeline(`cat $(nameIn)`, nameOut))
   rm(header)
   rm.(nameIn)
end

filenamesIn = "*.dat"
filenames = Vector{String}(undef,0)
filesfound = glob(filenamesIn, dir)
filenames = vcat(filenames, filesfound)
# Do not work on files that have already been converted
# This won't work if new files are generated with exactly the same name!
filenames = [fname for fname in filenames if ~isfile(fname[1:end-3]*"vtu")]

@sync @distributed for outname in filenames
   println("filename=$(outname)")
   head, data, connectivity = readtecdata(outname)
   convertVTK(head, data, connectivity, outname[1:end-4])
end

# Choose whether or not to delete the original *.dat files.
#datFiles = glob("*.dat", dir)
# Parallel:
#@sync @distributed for file in datFiles
#   rm(file)
#end
# Serial (might be faster if I/O bound):
#rm.(datfiles)


text:
using JuMP
using GLPK

G = [
0 100 30  0  0;
0   0 20  0  0;
0   0  0 10 60;
0  15  0  0 50;
0   0  0  0  0
]

n = size(G)[1]

shortest_path = Model(with_optimizer(GLPK.Optimizer))

@variable(shortest_path, x[1:n,1:n], Bin)
@constraint(shortest_path, [i = 1:n, j = 1:n; G[i,j] == 0], x[i,j] == 0) # Arcs with zero cost are not a part of the path as they do no exist
@constraint(shortest_path, [i = 1:n; i != 1 && i != 2], sum(x[i,:]) == sum(x[:,i])) # Flow conservation constraint
@constraint(shortest_path, sum(x[1,:]) - sum(x[:,1]) == 1) # Flow coming out of source = 1
@constraint(shortest_path, sum(x[2,:]) - sum(x[:,2]) == -1) # Flowing coming out of destination = -1 i.e. Flow entering destination = 1
@objective(shortest_path, Min, sum(G .* x))

optimize!(shortest_path)
@show objective_value(shortest_path);
@show value.(x);

G = [
6 4 5 0;
0 3 6 0;
5 0 4 3;
7 5 5 5;
]

n = size(G)[1]

assignment = Model(with_optimizer(GLPK.Optimizer))
@variable(assignment, y[1:n,1:n], Bin)
@constraint(assignment, [i = 1:n], sum(y[:,i]) == 1) # One person can only be assigned to one object
@constraint(assignment, [j = 1:n], sum(y[j,:]) == 1) # One object can only be assigned to one person
@objective(assignment, Max, sum(G .* y))

optimize!(assignment)
@show objective_value(assignment);
@show value.(y);

G = [
0 3 2 2 0 0 0 0
0 0 0 0 5 1 0 0
0 0 0 0 1 3 1 0
0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 4
0 0 0 0 0 0 0 2
0 0 0 0 0 0 0 4
0 0 0 0 0 0 0 0
]

n = size(G)[1]

max_flow = Model(with_optimizer(GLPK.Optimizer))

@variable(max_flow, f[1:n,1:n] >= 0)
@constraint(max_flow, [i = 1:n, j = 1:n], f[i,j] <= G[i,j]) # Capacity constraints
@constraint(max_flow, [i = 1:n; i != 1 && i != 8], sum(f[i,:]) == sum(f[:,i])) # Flow conservation contraints
@objective(max_flow, Max, sum(f[1, :]))

optimize!(max_flow)
@show objective_value(max_flow);
@show value.(f);


text:
function labels!(ec::EChart, series::Int; show::Bool = true, position::String = "insideTop")

	ec.series[series].label = Label(show = show, position = position)

    return ec

end

function labels!(ec::EChart; show::Bool = true, position::String = "insideTop")

	for n in 1:length(ec.series)
		labels!(ec, n, show = show, position = position)
	end

	return ec

end


text:
using Graphs
using SimpleGraphs

function swap_edge!(g::GenericGraph,e_o::Edge,newSrc::Int,newDest::Int)
    # swap edge
    eind = edge_index(e_o)
    fwd_edge = Edge(eind,newSrc,newDest)
    rev_edge = Edge(eind,newDest,newSrc)
    g.edges[eind] = fwd_edge
    e_o_prime = Edge(eind,target(e_o),source(e_o))

    # we are moving the entire edge delete the old one
    inc_idx_f = findfirst(g.finclist[source(e_o)].==e_o,true)
    if inc_idx_f == 0
        # delete old edge at target
        inc_idx_f = findfirst(g.finclist[source(e_o)].==e_o_prime,true)
        deleteat!(g.finclist[source(e_o)],inc_idx_f)

        # delete old edge at source
        inc_idx_r = findfirst(g.finclist[target(e_o)].==e_o,true)
        deleteat!(g.finclist[source(e_o)],inc_idx_r)
    else
        # delete old edge at target
        deleteat!(g.finclist[source(e_o)],inc_idx_f)
        # delete old edge at target
        inc_idx_r = findfirst(g.finclist[target(e_o)].==e_o_prime,true)
        deleteat!(g.finclist[target(e_o)],inc_idx_r)
    end
    # add forward edge
    push!(g.finclist[newSrc],fwd_edge)

    # add reverse edge
    push!(g.finclist[newDest], rev_edge)
end
function rewire(g::GenericGraph,p::Float64;self_loops=false,parallel_edges=false,keep_vertex_source=false)
    # Inputs:
    #   g: graph of type GenericGraph
    #   p: probability of rewiring an edge
    #   self_loops: bool indicating if self loops are allowed
    #   parallel_edges: bool indicating whether multiple edges can be between two vertices
    #   keep_vertex_source: whether an edge is constrained to keep its original source
    num_e = num_edges(g)
    num_v = num_vertices(g)
    edg_list = collect(edges(g))
    max_its = 3000
    for i in 1:num_e
        r = rand()
        if r < p
            e = edg_list[i]
            if keep_vertex_source
                # keep the same source to preserve vertex edges already connected to that vertex
                source_v = source(e)
            else
                # choose random source
                source_v = rand(collect(1:num_v))
            end

            if !parallel_edges
                neighbors = collect(Graphs.out_neighbors(source_v,g))
            end
            # rewire the edge
            ok = false
            ok_sl = false
            ok_pl = false
            cand = -1
            its = 0
            while !ok && its < max_its
                its += 1
                if its > max_its
                    println("exceeded its")
                end
                cand = rand(collect(1:num_v))
                if !self_loops
                    if cand != source_v
                        ok_sl = true
                    else
                        ok_sl = false
                    end
                else
                    ok_sl = true
                end
                if !parallel_edges
                    if cand in neighbors
                        ok_pl = false
                    else
                        ok_pl = true
                    end
                else
                    ok_pl = true
                end
                ok = ok_sl && ok_pl
            end
            # println(collect(out_neighbors(source_v,g)))
            swap_edge!(g,e,source_v,cand)
            # println(collect(out_neighbors(source_v,g)))
        else
            # don't rewire this edge, move on
            continue
        end
    end
    return g
end

function rewire_old(g::SimpleGraphs.SimpleGraph,p::Float64)
    max_its = 10
    for e in elist(g)
        s = e[1]
        t = e[2]
        if rand() < p
            vl = vlist(g)
            source_v = s
            dest_v = t
            its = 0
            while (has(g,source_v,dest_v) || source_v == dest_v) && its < max_its
                its += 1
                source_v = rand(vl)
                dest_v = rand(vl)
            end
            if its == max_its
                println("too many tries")
            end
        else
            continue
        end
        delete!(g,s,t)
        add!(g,source_v,dest_v)
    end
    return g
end

function rewire(g::SimpleGraphs.SimpleGraph,p::Float64)
    e_lst = elist(g)
    # choose which will be rewired
    num_e = NE(g)
    num_swaps = round(Int,p*num_e)
    if isodd(num_swaps)
        # make it even
        num_swaps -= 1
    end
    rewire_list = randperm(num_e)[1:num_swaps]
    # take edges two at a time and rewire sources to source and targets to targets
    for i = 1:2:num_swaps
        s1 = e_lst[rewire_list[i]][1]
        t1 = e_lst[rewire_list[i]][2]
        s2 = e_lst[rewire_list[i+1]][1]
        t2 = e_lst[rewire_list[i+1]][2]
        delete!(g,s1,t1)
        delete!(g,s2,t2)
        add!(g,s1,s2)
        add!(g,t1,t2)
    end
    if NE(g) != num_e
        e_diff = num_e - NE(g)
        for i = 1:e_diff
            ns = 1
            nt = 1
            while ns == nt || has(g,ns,nt)
                ns = rand(1:NV(g))
                nt = rand(1:NV(g))
            end
            add!(g,ns,nt)
        end
    end
    return g
end


text:
using KernelAbstractions: @kernel, @index, Event, MultiEvent
using OffsetArrays: OffsetArray

import Oceananigans.BoundaryConditions:
    fill_halo_regions!,
    fill_west_halo!, fill_east_halo!, fill_south_halo!,
    fill_north_halo!, fill_bottom_halo!, fill_top_halo!

#####
##### MPI tags for halo communication BCs
#####

sides  = (:west, :east, :south, :north, :top, :bottom)
side_id = Dict(side => n for (n, side) in enumerate(sides))

opposite_side = Dict(
    :west => :east, :east => :west,
    :south => :north, :north => :south,
    :bottom => :top, :top => :bottom
)

# Define functions that return unique send and recv MPI tags for each side.
# It's an integer where
#   digit 1: the side
#   digits 2-4: the "from" rank
#   digits 5-7: the "to" rank

RANK_DIGITS = 3

for side in sides
    side_str = string(side)
    send_tag_fn_name = Symbol("$(side)_send_tag")
    recv_tag_fn_name = Symbol("$(side)_recv_tag")
    @eval begin
        function $send_tag_fn_name(local_rank, rank_to_send_to)
            from_digits = string(local_rank, pad=RANK_DIGITS)
            to_digits = string(rank_to_send_to, pad=RANK_DIGITS)
            side_digit = string(side_id[Symbol($side_str)])
            return parse(Int, from_digits * to_digits * side_digit)
        end

        function $recv_tag_fn_name(local_rank, rank_to_recv_from)
            from_digits = string(rank_to_recv_from, pad=RANK_DIGITS)
            to_digits = string(local_rank, pad=RANK_DIGITS)
            side_digit = string(side_id[opposite_side[Symbol($side_str)]])
            return parse(Int, from_digits * to_digits * side_digit)
        end
    end
end

#####
##### Filling halos for halo communication boundary conditions
#####

function fill_halo_regions!(c::OffsetArray, bcs, arch::MultiArch, grid, c_location, args...; kwargs...)

    barrier = Event(device(child_architecture(arch)))

    x_events_requests = fill_west_and_east_halos!(c, bcs.west, bcs.east, arch, barrier, grid, c_location, args...; kwargs...)
    y_events_requests = fill_south_and_north_halos!(c, bcs.south, bcs.north, arch, barrier, grid, c_location, args...; kwargs...)
    z_events_requests = fill_bottom_and_top_halos!(c, bcs.bottom, bcs.top, arch, barrier, grid, c_location, args...; kwargs...)

    events_and_requests = [x_events_requests..., y_events_requests..., z_events_requests...]

    mpi_requests = filter(e -> e isa MPI.Request, events_and_requests) |> Array{MPI.Request}
    # Length check needed until this PR is merged: https://github.com/JuliaParallel/MPI.jl/pull/458
    length(mpi_requests) > 0 && MPI.Waitall!(mpi_requests)

    events = filter(e -> e isa Event, events_and_requests)
    wait(device(child_architecture(arch)), MultiEvent(Tuple(events)))

    return nothing
end

#####
##### fill_west_and_east_halos!   }
##### fill_south_and_north_halos! } for non-communicating boundary conditions (fallback)
##### fill_bottom_and_top_halos!  }
#####

for (side, opposite_side) in zip([:west, :south, :bottom], [:east, :north, :top])
    fill_both_halos! = Symbol("fill_$(side)_and_$(opposite_side)_halos!")
    fill_side_halo! = Symbol("fill_$(side)_halo!")
    fill_opposite_side_halo! = Symbol("fill_$(opposite_side)_halo!")

    @eval begin
        function $fill_both_halos!(c, bc_side, bc_opposite_side, arch, barrier, grid, args...; kwargs...)
            event_side = $fill_side_halo!(c, bc_side, child_architecture(arch), barrier, grid, args...; kwargs...)
            event_opposite_side = $fill_opposite_side_halo!(c, bc_opposite_side, child_architecture(arch), barrier, grid, args...)
            return event_side, event_opposite_side
        end
    end
end

#####
##### fill_west_and_east_halos!   }
##### fill_south_and_north_halos! } for when both halos are communicative
##### fill_bottom_and_top_halos!  }
#####

for (side, opposite_side) in zip([:west, :south, :bottom], [:east, :north, :top])
    fill_both_halos! = Symbol("fill_$(side)_and_$(opposite_side)_halos!")
    send_side_halo = Symbol("send_$(side)_halo")
    send_opposite_side_halo = Symbol("send_$(opposite_side)_halo")
    recv_and_fill_side_halo! = Symbol("recv_and_fill_$(side)_halo!")
    recv_and_fill_opposite_side_halo! = Symbol("recv_and_fill_$(opposite_side)_halo!")

    @eval begin
        function $fill_both_halos!(c, bc_side::HaloCommunicationBC, bc_opposite_side::HaloCommunicationBC,
                                   arch, barrier, grid, c_location, args...; kwargs...)

            @assert bc_side.condition.from == bc_opposite_side.condition.from  # Extra protection in case of bugs
            local_rank = bc_side.condition.from

            send_req1 = $send_side_halo(c, grid, c_location, local_rank, bc_side.condition.to)
            send_req2 = $send_opposite_side_halo(c, grid, c_location, local_rank, bc_opposite_side.condition.to)

            recv_req1 = $recv_and_fill_side_halo!(c, grid, c_location, local_rank, bc_side.condition.to)
            recv_req2 = $recv_and_fill_opposite_side_halo!(c, grid, c_location, local_rank, bc_opposite_side.condition.to)

            return send_req1, send_req2, recv_req1, recv_req2
        end
    end
end

#####
##### Sending halos
#####

for side in sides
    side_str = string(side)
    send_side_halo = Symbol("send_$(side)_halo")
    underlying_side_boundary = Symbol("underlying_$(side)_boundary")
    side_send_tag = Symbol("$(side)_send_tag")

    @eval begin
        function $send_side_halo(c, grid, c_location, local_rank, rank_to_send_to)
            send_buffer = $underlying_side_boundary(c, grid, c_location)
            send_tag = $side_send_tag(local_rank, rank_to_send_to)

            @debug "Sending " * $side_str * " halo: local_rank=$local_rank, rank_to_send_to=$rank_to_send_to, send_tag=$send_tag"
            send_req = MPI.Isend(send_buffer, rank_to_send_to, send_tag, MPI.COMM_WORLD)

            return send_req
        end
    end
end

#####
##### Receiving and filling halos (buffer is a view so it gets filled upon receive)
#####

for side in sides
    side_str = string(side)
    recv_and_fill_side_halo! = Symbol("recv_and_fill_$(side)_halo!")
    underlying_side_halo = Symbol("underlying_$(side)_halo")
    side_recv_tag = Symbol("$(side)_recv_tag")

    @eval begin
        function $recv_and_fill_side_halo!(c, grid, c_location, local_rank, rank_to_recv_from)
            recv_buffer = $underlying_side_halo(c, grid, c_location)
            recv_tag = $side_recv_tag(local_rank, rank_to_recv_from)

            @debug "Receiving " * $side_str * " halo: local_rank=$local_rank, rank_to_recv_from=$rank_to_recv_from, recv_tag=$recv_tag"
            recv_req = MPI.Irecv!(recv_buffer, rank_to_recv_from, recv_tag, MPI.COMM_WORLD)

            return recv_req
        end
    end
end


text:
# 1D linear convection problem

# Packages and inclusions
using ModelingToolkit, MethodOfLines, DiffEqBase, LinearAlgebra, Test, DomainSets
#using Plots
# Tests

@testset "Test 00: Dt(u(t,x)) ~ -Dx(u(t,x))" begin
    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)

    # 1D PDE and boundary conditions
    eq = Dt(u(t, x)) ~ -Dx(u(t, x))

    asf(x) = (0.5 / (0.2 * sqrt(2.0 * 3.1415))) * exp(-(x - 1.0)^2 / (2.0 * 0.2^2))
    bcs = [u(0, x) ~ asf(x),
        u(t, 0) ~ u(t, 2)]

    # Space and time domains
    domains = [t ∈ Interval(0.0, 2.0),
        x ∈ Interval(0.0, 2.0)]

    # PDE system
    @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)])

    # Method of lines discretization
    dx = 2 / 80
    order = 1
    discretization = MOLFiniteDifference([x => dx], t)
    # explicitly specify upwind order
    discretization_upwind = MOLFiniteDifference([x => dx], t; upwind_order=order)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys, discretization)
    prob_upwind = discretize(pdesys, discretization_upwind)

    # Solve ODE problem
    using OrdinaryDiffEq
    sol = solve(prob, Euler(), dt=0.025, saveat=0.1)
    sol_upwind = solve(prob_upwind, Euler(), dt=0.025, saveat=0.1)

    x_interval = infimum(domains[2].domain)+dx:dx:supremum(domains[2].domain)
    u = asf.(x_interval)
    # Plot and save results

    # Test
    t_f = size(sol, 3)

    x_sol = x_interval
    t_f = size(sol, 3)
    exact = u
    # plot(x_sol, u, seriestype = :scatter,label="Analytic solution")
    # plot!(x_sol, sol.u[end], label="Numeric solution")
    # plot!(x_sol, u.-sol.u[end], label="Differential Error")

    # savefig("plots/MOL_Linear_Convection_Test00.png")

    @test sol.u[end] ≈ u atol = 0.1
    @test sol_upwind.u[end] ≈ u atol = 0.1
end

@testset "Test 00a: Dt(u(t,x)) ~ Dx(u(t,x))" begin
    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)

    # 1D PDE and boundary conditions
    eq = Dt(u(t, x)) ~ Dx(u(t, x))

    asf(x) = (0.5 / (0.2 * sqrt(2.0 * 3.1415))) * exp(-(x - 1.0)^2 / (2.0 * 0.2^2))
    bcs = [u(0, x) ~ asf(x),
        u(t, 0) ~ u(t, 2)]

    # Space and time domains
    domains = [t ∈ Interval(0.0, 2.0),
        x ∈ Interval(0.0, 2.0)]

    # PDE system
    @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)])

    # Method of lines discretization
    dx = 2 / 80
    order = 1
    discretization = MOLFiniteDifference([x => dx], t)
    # explicitly specify upwind order
    discretization_upwind = MOLFiniteDifference([x => dx], t; upwind_order=order)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys, discretization)
    prob_upwind = discretize(pdesys, discretization_upwind)

    # Solve ODE problem
    using OrdinaryDiffEq
    sol = solve(prob, Euler(), dt=0.025, saveat=0.1)
    sol_upwind = solve(prob_upwind, Euler(), dt=0.025, saveat=0.1)

    x_interval = infimum(domains[2].domain)+dx:dx:supremum(domains[2].domain)
    u = asf.(x_interval)
    # Plot and save results

    # Test
    t_f = size(sol, 3)

    x_sol = x_interval
    t_f = size(sol, 3)
    exact = u
    # plot(x_sol, u, seriestype = :scatter,label="Analytic solution")
    # plot!(x_sol, sol.u[end], label="Numeric solution")
    # plot!(x_sol, u.-sol.u[end], label="Differential Error")

    # savefig("plots/MOL_Linear_Convection_Test00.png")

    @test sol.u[end] ≈ u atol = 0.1
    @test sol_upwind.u[end] ≈ u atol = 0.1
end

@testset "Test 00b: Dt(u(t,x)) - Dx(u(t,x)) ~ 0" begin
    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)

    # 1D PDE and boundary conditions
    eq = Dt(u(t, x)) - Dx(u(t, x)) ~ 0

    asf(x) = (0.5 / (0.2 * sqrt(2.0 * 3.1415))) * exp(-(x - 1.0)^2 / (2.0 * 0.2^2))
    bcs = [u(0, x) ~ asf(x),
        u(t, 0) ~ u(t, 2)]

    # Space and time domains
    domains = [t ∈ Interval(0.0, 2.0),
        x ∈ Interval(0.0, 2.0)]

    # PDE system
    @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)])

    # Method of lines discretization
    dx = 2 / 80
    order = 1
    discretization = MOLFiniteDifference([x => dx], t)
    # explicitly specify upwind order
    discretization_upwind = MOLFiniteDifference([x => dx], t; upwind_order=order)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys, discretization)
    prob_upwind = discretize(pdesys, discretization_upwind)

    # Solve ODE problem
    using OrdinaryDiffEq
    sol = solve(prob, Euler(), dt=0.025, saveat=0.1)
    sol_upwind = solve(prob_upwind, Euler(), dt=0.025, saveat=0.1)

    x_interval = infimum(domains[2].domain)+dx:dx:supremum(domains[2].domain)
    u = asf.(x_interval)
    # Plot and save results

    # Test
    t_f = size(sol, 3)

    x_sol = x_interval
    t_f = size(sol, 3)
    exact = u
    # plot(x_sol, u, seriestype = :scatter,label="Analytic solution")
    # plot!(x_sol, sol.u[end], label="Numeric solution")
    # plot!(x_sol, u.-sol.u[end], label="Differential Error")

    # savefig("plots/MOL_Linear_Convection_Test00.png")

    @test sol.u[end] ≈ u atol = 0.1
    @test sol_upwind.u[end] ≈ u atol = 0.1
end

@testset "Test 00c: Dt(u(t,x)) + Dx(u(t,x)) ~ 0" begin
    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)

    # 1D PDE and boundary conditions
    eq = Dt(u(t, x)) + Dx(u(t, x)) ~ 0

    asf(x) = (0.5 / (0.2 * sqrt(2.0 * 3.1415))) * exp(-(x - 1.0)^2 / (2.0 * 0.2^2))
    bcs = [u(0, x) ~ asf(x),
        u(t, 0) ~ u(t, 2)]

    # Space and time domains
    domains = [t ∈ Interval(0.0, 2.0),
        x ∈ Interval(0.0, 2.0)]

    # PDE system
    @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)])

    # Method of lines discretization
    dx = 2 / 80
    order = 1
    discretization = MOLFiniteDifference([x => dx], t)
    # explicitly specify upwind order
    discretization_upwind = MOLFiniteDifference([x => dx], t; upwind_order=order)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys, discretization)
    prob_upwind = discretize(pdesys, discretization_upwind)

    # Solve ODE problem
    using OrdinaryDiffEq
    sol = solve(prob, Euler(), dt=0.025, saveat=0.1)
    sol_upwind = solve(prob_upwind, Euler(), dt=0.025, saveat=0.1)

    x_interval = infimum(domains[2].domain)+dx:dx:supremum(domains[2].domain)
    u = asf.(x_interval)
    # Plot and save results

    # Test
    t_f = size(sol, 3)

    x_sol = x_interval
    t_f = size(sol, 3)
    exact = u
    # plot(x_sol, u, seriestype = :scatter,label="Analytic solution")
    # plot!(x_sol, sol.u[end], label="Numeric solution")
    # plot!(x_sol, u.-sol.u[end], label="Differential Error")

    # savefig("plots/MOL_Linear_Convection_Test00.png")

    @test sol.u[end] ≈ u atol = 0.1
    @test sol_upwind.u[end] ≈ u atol = 0.1
end

@testset "Test 01: Dt(u(t,x)) ~ -Dx(u(t,x)) + 0.001" begin

    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)

    # 1D PDE and boundary conditions
    eq = Dt(u(t, x)) ~ -Dx(u(t, x)) + 0.001
    asf(x) = (0.5 / (0.2 * sqrt(2.0 * 3.1415))) * exp(-(x - 1.0)^2 / (2.0 * 0.2^2))
    bcs = [u(0, x) ~ asf(x),
        u(t, 0) ~ u(t, 2)]

    # Space and time domains
    domains = [t ∈ Interval(0.0, 2.0),
        x ∈ Interval(0.0, 2.0)]

    # PDE system
    @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)])

    # Method of lines discretization
    dx = 2 / 80
    order = 1
    discretization = MOLFiniteDifference([x => dx], t)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys, discretization)

    # Solve ODE problem
    using OrdinaryDiffEq
    sol = solve(prob, Euler(), dt=0.025, saveat=0.1)

    # Test
    x_interval = infimum(domains[2].domain)+dx:dx:supremum(domains[2].domain)
    u = asf.(x_interval)
    # Plot and save results

    # Test
    t_f = size(sol, 3)

    x_sol = x_interval
    t_f = size(sol, 3)
    exact = u
    # plot(x_sol, u, seriestype = :scatter,label="Analytic solution")
    # plot!(x_sol, sol.u[end], label="Numeric solution")
    # plot!(x_sol, u.-sol.u[end], label="Differential Error")

    # savefig("plots/MOL_Linear_Convection_Test01.png")

    @test sol.u[end] ≈ u atol = 0.1

end

@testset "Test 02: Dt(u(t,x)) ~ -v*Dx(u(t,x))" begin
    # Parameters, variables, and derivatives
    @parameters t x v
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)

    v = 1.0

    # 1D PDE and boundary conditions
    eq = Dt(u(t, x)) ~ -v * Dx(u(t, x))
    asf(x) = (0.5 / (0.2 * sqrt(2.0 * 3.1415))) * exp(-(x - 1.0)^2 / (2.0 * 0.2^2))
    bcs = [u(0, x) ~ asf(x),
        u(t, 0) ~ u(t, 2)]

    # Space and time domains
    domains = [t ∈ Interval(0.0, 2.0),
        x ∈ Interval(0.0, 2.0)]

    # PDE system
    @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)])

    # Method of lines discretization
    dx = 2 / 80
    order = 1
    discretization = MOLFiniteDifference([x => dx], t, upwind_order=1)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys, discretization)

    # Solve ODE problemdoes this mean that - doesn't seem to affect things
    using OrdinaryDiffEq
    sol = solve(prob, Euler(), dt=0.025, saveat=0.1)

    # Plot and save results
    # using Plots
    # plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,5]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,6]))
    # savefig("MOL_1D_Linear_Convection_Test02.png")

    # Test
    x_interval = infimum(domains[2].domain)+dx:dx:supremum(domains[2].domain)
    u = asf.(x_interval)
    # Plot and save results

    # Test
    t_f = size(sol, 3)

    x_sol = x_interval
    t_f = size(sol, 3)
    exact = u
    # plot(x_sol, u, seriestype = :scatter,label="Analytic solution")
    # plot!(x_sol, sol.u[end], label="Numeric solution")
    # plot!(x_sol, u.-sol.u[end], label="Differential Error")

    # savefig("plots/MOL_1D_Linear_Convection_Test02.png")

    # savefig("plots/MOL_Linear_Convection_Test02.png")
    # plot()
    # anim = @animate for i in eachindex(sol.t)
    #     plot!(x_sol, sol.u[i], label="Numeric solution at $(sol.t[i])")
    # end
    # gif(anim, "plots/MOL_Linear_Convection_Test02.gif", fps = 5)

    @test sol.u[end] ≈ u atol = 0.1
end
@test_broken begin#@testset "Test 03: Dt(u(t,x)) ~ -Dx(v(t,x))*u(t,x)-v(t,x)*Dx(u(t,x)) with v(t,x)=1" begin
    # Parameters, variables, and derivatives
    @parameters t x
    @variables v(..) u(..)
    Dt = Differential(t)
    Dx = Differential(x)

    # 1D PDE and boundary conditions
    eq = [Dt(u(t, x)) ~ -(Dx(v(t, x)) * u(t, x) + v(t, x) * Dx(u(t, x))),
        v(t, x) ~ 1.0]
    asf(x) = (0.5 / (0.2 * sqrt(2.0 * 3.1415))) * exp(-(x - 1.0)^2 / (2.0 * 0.2^2))
    bcs = [u(0, x) ~ asf(x),
        u(t, 0) ~ u(t, 2),
        v(0, x) ~ 1.0,
        v(t, 0) ~ 1.0,
        v(t, 2) ~ 1.0
    ]

    # Space and time domains

    # Space and time domains
    domains = [t ∈ Interval(0.0, 2.0),
        x ∈ Interval(0.0, 2.0)]

    # PDE system
    @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x), v(t, x)])

    # Method of lines discretization
    dx = 2 / 80
    order = 1
    discretization = MOLFiniteDifference([x => dx], t)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys, discretization)

    # Solve ODE problem
    using OrdinaryDiffEq
    sol = solve(prob, Euler(), dt=0.025, saveat=0.1)

    #Plot and save results
    # using Plots
    # plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,5]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,6]))
    # savefig("MOL_1D_Linear_Convection_Test03.png")

    # Test
    x_interval = infimum(domains[2].domain)+dx:dx:supremum(domains[2].domain)-dx
    u = asf.(x_interval)
    # Plot and save results

    # Test
    t_f = size(sol, 3)

    x_sol = x_interval
    t_f = size(sol, 3)
    exact = u
    # plot(x_sol, u, seriestype = :scatter,label="Analytic solution")
    # plot!(x_sol, sol.u[end], label="Numeric solution")
    # plot!(x_sol, u.-sol.u[end], label="Differential Error")

    # savefig("plots/MOL_Linear_Convection_Test03.png")

    @test sol.u[end] ≈ u atol = 0.1
end

@testset "Test 04: Dt(u(t,x)) ~ -Dx(v(t,x))*u(t,x)-v(t,x)*Dx(u(t,x)) with v(t,x)=0.999 + 0.001 * t * x " begin
    # Parameters, variables, and derivatives
    @parameters t x
    @variables v(..) u(..)
    Dt = Differential(t)
    Dx = Differential(x)

    # 1D PDE and boundary conditions
    eq = [Dt(u(t, x)) ~ -Dx(v(t, x)) * u(t, x) - v(t, x) * Dx(u(t, x)),
        v(t, x) ~ 0.999 + 0.001 * t * x]
    asf(x) = (0.5 / (0.2 * sqrt(2.0 * 3.1415))) * exp(-(x - 1.0)^2 / (2.0 * 0.2^2))

    bcs = [u(0, x) ~ asf(x),
        u(t, 0) ~ u(t, 2),
        v(0, x) ~ 0.999,
        v(t, 0) ~ 0.999,
        v(t, 2) ~ 0.9999 + 0.0001 * t * 2.0]

    # Space and time domains
    domains = [t ∈ Interval(0.0, 2.0),
        x ∈ Interval(0.0, 2.0)]

    # PDE system
    @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x), v(t, x)])

    # Method of lines discretization
    dx = 2 / 80
    order = 1
    discretization = MOLFiniteDifference([x => dx], t)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys, discretization)

    # Solve ODE problem
    using OrdinaryDiffEq
    sol = solve(prob, Euler(), dt=0.025, saveat=0.1)

    # Plot and save results
    # using Plots
    # plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,5]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,6]))
    # savefig("MOL_1D_Linear_Convection_Test04.png")

    # Test
    x_interval = infimum(domains[2].domain)+dx:dx:supremum(domains[2].domain)
    u = asf.(x_interval)
    # Plot and save results

    # Test
    t_f = size(sol, 3)

    x_sol = x_interval
    t_f = size(sol, 3)
    exact = u
    # plot(x_sol, u, seriestype = :scatter,label="Analytic solution")
    # plot!(x_sol, sol.u[end], label="Numeric solution")
    # plot!(x_sol, u.-sol.u[end], label="Differential Error")

    # savefig("plots/MOL_Linear_Convection_Test04.png")

    # plot()
    # anim = @animate for i in eachindex(sol.t)
    #     plot!(x_sol, sol.u[i], label="Numeric solution at $(sol.t[i])")
    # end
    # gif(anim, "plots/MOL_Linear_Convection_Test04.gif", fps = 5)

    @test sol.u[end] ≈ u atol = 0.1
end

@test_broken begin #@testset "Test 05: Dt(u(t,x)) ~ -Dx(v(t,x)*u(t,x)) with v(t,x)=0.999 + 0.001 * t * x " begin
    # Parameters, variables, and derivatives
    @parameters t x
    @variables v(..) u(..)
    Dt = Differential(t)
    Dx = Differential(x)

    # 1D PDE and boundary conditions
    eq = [Dt(u(t, x)) ~ -Dx(v(t, x) * u(t, x)),
        v(t, x) ~ 0.999 + 0.001 * t * x]
    asf(x) = (0.5 / (0.2 * sqrt(2.0 * 3.1415))) * exp(-(x - 1.0)^2 / (2.0 * 0.2^2))

    bcs = [u(0, x) ~ asf(x),
        u(t, 0) ~ u(t, 2),
        v(0, x) ~ 0.999,
        v(t, 0) ~ 0.999,
        v(t, 2) ~ 0.9999 + 0.0001 * t * 2.0]

    # Space and time domains
    domains = [t ∈ Interval(0.0, 2.0),
        x ∈ Interval(0.0, 2.0)]

    # PDE system
    @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x), v(t, x)])

    # Method of lines discretization
    dx = 2 / 80
    order = 1
    discretization = MOLFiniteDifference([x => dx], t)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys, discretization)

    # Solve ODE problem
    using OrdinaryDiffEq
    sol = solve(prob, Euler(), dt=0.025, saveat=0.1)

    # Plot and save results
    # using Plots
    # plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,5]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,6]))
    # savefig("MOL_1D_Linear_Convection_Test04.png")

    # Test
    x_interval = infimum(domains[2].domain)+dx:dx:supremum(domains[2].domain)
    u = asf.(x_interval)
    # Plot and save results

    # Test
    t_f = size(sol, 3)

    x_sol = x_interval
    t_f = size(sol, 3)
    exact = u
    # plot(x_sol, u, seriestype = :scatter,label="Analytic solution")
    # plot!(x_sol, sol.u[end], label="Numeric solution")
    # plot!(x_sol, u.-sol.u[end], label="Differential Error")

    # savefig("plots/MOL_Linear_Convection_Test04.png")

    # plot()
    # anim = @animate for i in eachindex(sol.t)
    #     plot!(x_sol, sol.u[i], label="Numeric solution at $(sol.t[i])")
    # end
    # gif(anim, "plots/MOL_Linear_Convection_Test04.gif", fps = 5)

    @test sol.u[end] ≈ u atol = 0.1
end


text:
@doc doc"Abstract Boolean Types: {{true},{false}.{true,false}}" ->
immutable AbstractBool <: Domain{Bool}
  v::UInt8
  AbstractBool(v::UInt8) = (@assert v == 0x1 || v == 0x2 || v== 0x3; new(v))
end

const t = AbstractBool(0x1)
const f = AbstractBool(0x2)
const tf = AbstractBool(0x3)

promote_rule(::Type{Bool}, ::Type{AbstractBool}) = AbstractBool
convert(::Type{AbstractBool}, b::Bool) = if b t else f end

## AbstractBool Set Operations
## ===========================
subsumes(x::AbstractBool, y::AbstractBool) = x === tf || x === y
subsumes(x::AbstractBool, y::Bool) = subsumes(x,convert(AbstractBool, y))

isintersect(x::AbstractBool, y::AbstractBool) = !((x === t && y === f) || (x === f && y === t))
isintersect(x::AbstractBool, y::Bool) = isintersect(x,convert(AbstractBool, y))
isintersect(x::Bool, y::AbstractBool) = isintersect(convert(AbstractBool, x),y)

isrelational(::Type{AbstractBool}) = false
isabstract(c::Type{Bool}, a::Type{AbstractBool}) = true

isequal(x::AbstractBool, y::AbstractBool) = x.v == y.v
domaineq(x::AbstractBool, y::AbstractBool) = x.v == y.v

⊔(a::AbstractBool) = a
⊔(a::AbstractBool, b::AbstractBool) = a === b ? a : tf
⊔(a::Bool, b::AbstractBool) = ⊔(convert(AbstractBool,a),b)
⊔(a::AbstractBool, b::Bool) = ⊔(a,convert(AbstractBool,b))
⊔(a::Bool, b::Bool) = a === b ? convert(AbstractBool,a) : tf

## =========================
## Lifted Boolean Arithmetic

function !(b::AbstractBool)
  if b === t
    f
  elseif b === f
    t
  elseif b === tf
    tf
  end
end

(==)(x::AbstractBool, y::AbstractBool) =
  x === tf || y === tf ? tf : x === t && y === t || x === f && y === f

function (==)(x::AbstractBool, y::AbstractBool)
  if x === tf || y === tf tf
  elseif x === t && y === t t
  elseif x === f && y === f t
  else f
  end
end
(==)(x::AbstractBool, y::Bool) = (==)(promote(x,y)...)
(==)(y::Bool,x::AbstractBool) = (==)(promote(y,x)...)

function (|)(x::AbstractBool, y::AbstractBool)
  if x === t || y === t t
  elseif x === tf || y === tf tf
  else f
  end
end
|(x::AbstractBool, y::Bool) = |(x,convert(AbstractBool,y))
|(y::Bool, x::AbstractBool) = |(convert(AbstractBool,y), x)

function (&)(x::AbstractBool, y::AbstractBool)
  if x === f || y === f f
  elseif x === tf || y === tf tf
  else t
  end
end

(&)(x::AbstractBool, y::Bool) = x & convert(AbstractBool, y)
(&)(y::Bool, x::AbstractBool) = convert(AbstractBool, y) & x

# When condition is TF we need to evaluate both branches
# and merge with ⊔
function ifelse(c::AbstractBool, x, y)
  if c === t
    x
  elseif c === f
    y
  elseif c === tf
    ⊔(x,y)
  end
end

## Printing
## ========
string(x::AbstractBool) = ["{true}","{false}","{true,false}"][x.v]
print(io::IO, x::AbstractBool) = print(io, string(x))
show(io::IO, x::AbstractBool) = print(io, string(x))
showcompact(io::IO, x::AbstractBool) = print(io, string(x))