text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Random

ioslength(io::IOBuffer) = (io.seekable ? io.size : bytesavailable(io))

bufcontents(io::Base.GenericIOBuffer) = unsafe_string(pointer(io.data), io.size)

@testset "Read/write empty IOBuffer" begin
    io = IOBuffer()
    @test eof(io)
    @test_throws EOFError read(io,UInt8)
    @test write(io,"abc") === 3
    @test isreadable(io)
    @test iswritable(io)
    @test isopen(io)
    @test ioslength(io) == 3
    @test position(io) == 3
    @test eof(io)
    seek(io, 0)
    @test read(io,UInt8) == convert(UInt8, 'a')
    a = Vector{UInt8}(undef, 2)
    @test read!(io, a) == a
    @test a == UInt8['b','c']
    @test bufcontents(io) == "abc"
    seek(io, 1)
    truncate(io, 2)
    @test position(io) == 1
    @test !eof(io)
    seekend(io)
    @test position(io) == 2
    truncate(io, 0)
    @test position(io) == 0
    truncate(io, 10)
    @test position(io) == 0
    @test all(io.data .== 0)
    @test write(io,Int16[1,2,3,4,5,6]) === 12
    seek(io, 2)
    truncate(io, 10)
    @test ioslength(io) == 10
    io.readable = false
    @test_throws ArgumentError read!(io,UInt8[0])
    truncate(io, 0)
    @test write(io,"boston\ncambridge\n") > 0
    @test String(take!(io)) == "boston\ncambridge\n"
    @test String(take!(io)) == ""
    @test write(io, Complex{Float64}(0)) === 16
    @test write(io, Rational{Int64}(1//2)) === 16
    close(io)
    @test_throws ArgumentError write(io,UInt8[0])
    @test_throws ArgumentError seek(io,0)
    @test eof(io)
end

@testset "Read/write readonly IOBuffer" begin
    io = IOBuffer("hamster\nguinea pig\nturtle")
    @test position(io) == 0
    @test readline(io) == "hamster"
    @test read(io, String) == "guinea pig\nturtle"
    @test_throws EOFError read(io,UInt8)
    seek(io,0)
    @test read(io,UInt8) == convert(UInt8, 'h')
    @test_throws ArgumentError truncate(io,0)
    @test_throws ArgumentError write(io,UInt8(0))
    @test_throws ArgumentError write(io,UInt8[0])
    @test String(take!(io)) == "hamster\nguinea pig\nturtle"
    @test String(take!(io)) == "hamster\nguinea pig\nturtle" #should be unchanged
    @test_throws ArgumentError Base.compact(io) # not writeable
    close(io)
end

@testset "PipeBuffer" begin
    io = PipeBuffer()
    @test_throws EOFError read(io,UInt8)
    @test write(io,"pancakes\nwaffles\nblueberries\n") > 0
    @test position(io) == 0
    @test readline(io) == "pancakes"
    Base.compact(io)
    @test readline(io) == "waffles"
    @test write(io,"whipped cream\n") > 0
    @test readline(io) == "blueberries"
    @test_throws ArgumentError seek(io,0)
    @test_throws ArgumentError truncate(io,0)
    @test readline(io) == "whipped cream"
    @test write(io,"pancakes\nwaffles\nblueberries\n") > 0
    @test readlines(io) == String["pancakes", "waffles", "blueberries"]
    write(io,"\n\r\n\n\r \n") > 0
    @test readlines(io, keep=true) == String["\n", "\r\n", "\n", "\r \n"]
    write(io,"\n\r\n\n\r \n") > 0
    @test readlines(io, keep=false) == String["", "", "", "\r "]
    @test write(io,"α\nβ\nγ\nδ") > 0
    @test readlines(io, keep=true) == String["α\n","β\n","γ\n","δ"]
    @test write(io,"α\nβ\nγ\nδ") > 0
    @test readlines(io, keep=false) == String["α", "β", "γ", "δ"]
    @test readlines(IOBuffer(""), keep=true) == []
    @test readlines(IOBuffer(""), keep=false) == []
    @test readlines(IOBuffer("first\nsecond"), keep=true) == String["first\n", "second"]
    @test readlines(IOBuffer("first\nsecond"), keep=false) == String["first", "second"]

    let fname = tempname()
        for dokeep in [true, false],
            endline in ["\n", "\r\n"],
            i in -5:5

            ref = ("1"^(2^17 - i)) * endline
            open(fname, "w") do io
                write(io, ref)
            end
            x = readlines(fname, keep = dokeep)
            if !dokeep
                ref = chomp(ref)
            end
            @test ref == x[1]
        end
        rm(fname)
    end

    Base.compact(io)
    @test position(io) == 0
    @test ioslength(io) == 0
    Base.ensureroom(io,50)
    @test position(io) == 0
    @test ioslength(io) == 0
    @test length(io.data) == 50
    Base.ensureroom(io,10)
    @test ioslength(io) == 0
    @test length(io.data) == 50
    io.maxsize = 75
    Base.ensureroom(io,100)
    @test ioslength(io) == 0
    @test length(io.data) == 75
    seekend(io)
    @test ioslength(io) == 0
    @test position(io) == 0
    write(io,zeros(UInt8,200))
    @test ioslength(io) == 75
    @test length(io.data) == 75
    write(io,1)
    @test ioslength(io) == 75
    @test length(io.data) == 75
    write(io,[1,2,3])
    @test ioslength(io) == 75
    @test length(io.data) == 75
    skip(io,1)
    @test write(io,UInt8(104)) === 1
    skip(io,3)
    @test write(io,b"apples") === 3
    skip(io,71)
    @test write(io,'y') === 1
    @test read(io, String) == "happy"
    @test eof(io)
    write(io,zeros(UInt8,73))
    write(io,'a')
    write(io,'b')
    write(io,'c')
    write(io,'d')
    write(io,'e')
    @test ioslength(io) == 75
    @test length(io.data) == 75
    @test position(io) == 0
    skip(io,72)
    @test String(take!(io)) == "\0ab"
    @test String(take!(io)) == ""

    # issues 4021
    print(io, true)
    close(io)
end

@testset "issue 5453" begin
    io = IOBuffer("abcdef")
    a = Vector{UInt8}(undef, 1024)
    @test_throws EOFError read!(io,a)
    @test eof(io)
end

@test isempty(readlines(IOBuffer(), keep=true))

@testset "issue #8193" begin
    io = IOBuffer("asdf")
    @test position(skip(io, -1)) == 0
    @test position(skip(io, 6)) == 4
    @test position(seek(io, -1)) == 0
    @test position(seek(io, 6)) == 4
end

@testset "issue #10658" begin
    io = IOBuffer("hello")
    @test position(skip(io, 4)) == 4
    @test position(skip(io, 10)) == 5
    @test position(skip(io, -2)) == 3
    @test position(skip(io, -3)) == 0
    @test position(skip(io, -3)) == 0
end

@testset "pr #11554" begin
    io  = IOBuffer(SubString("***αhelloworldω***", 4, 16))
    io2 = IOBuffer(Vector{UInt8}(b"goodnightmoon"), read=true, write=true)

    @test read(io, Char) == 'α'
    @test_throws ArgumentError write(io,"!")
    @test_throws ArgumentError write(io,'β')
    a = Vector{UInt8}(undef, 10)
    @test read!(io, a) === a
    @test String(a) == "helloworld"
    @test read(io, Char) == 'ω'
    @test_throws EOFError read(io,UInt8)
    skip(io, -3)
    @test read(io, String) == "dω"
    @test bufcontents(io) == "αhelloworldω"
    @test_throws ArgumentError write(io,"!")
    @test take!(io) == b"αhelloworldω"
    seek(io, 2)
    seekend(io2)
    write(io2, io)
    @test read(io, String) == ""
    @test read(io2, String) == ""
    @test String(take!(io)) == "αhelloworldω"
    seek(io2, 0)
    truncate(io2, io2.size - 2)
    @test read(io2, String) == "goodnightmoonhelloworld"
    seek(io2, 0)
    write(io2, io2)
    @test read(io2, String) == ""
    @test bufcontents(io2) == "goodnightmoonhelloworld"
end

# issue #11917
# (previous tests triggered this sometimes, but this should trigger nearly all the time)
let io = IOBuffer(maxsize=0)
   write(io, fill(0x01, 1048577))
end

@testset "BufferStream" begin
    bstream = Base.BufferStream()
    @test isopen(bstream)
    @test isreadable(bstream)
    @test iswritable(bstream)
    @test bytesavailable(bstream) == 0
    @test sprint(show, bstream) == "BufferStream() bytes waiting:$(bytesavailable(bstream.buffer)), isopen:true"
    a = rand(UInt8,10)
    write(bstream,a)
    @test !eof(bstream)
    @test flush(bstream) === nothing
    b = read(bstream,UInt8)
    @test a[1] == b
    b = read(bstream,UInt8)
    @test a[2] == b
    c = zeros(UInt8,8)
    @test bytesavailable(bstream) == 8
    @test !eof(bstream)
    read!(bstream,c)
    @test c == a[3:10]
    @test close(bstream) === nothing
    @test eof(bstream)
    @test bytesavailable(bstream) == 0
    flag = Ref{Bool}(false)
    event = Base.Event()
    bstream = Base.BufferStream()
    task = @async begin
        notify(event)
        read(bstream, 16)
        flag[] = true
    end
    wait(event)
    write(bstream, rand(UInt8, 16))
    wait(task)
    @test flag[] == true
end

@test flush(IOBuffer()) === nothing # should be a no-op

# pr #19461
let io = IOBuffer()
    @test Base.buffer_writes(io) === io
end

@testset "skipchars" begin
    io = IOBuffer("")
    @test eof(skipchars(isspace, io))

    io = IOBuffer("   ")
    @test eof(skipchars(isspace, io))

    io = IOBuffer("#    \n     ")
    @test eof(skipchars(isspace, io, linecomment='#'))

    io = IOBuffer("      text")
    skipchars(isspace, io)
    @test String(readavailable(io)) == "text"

    io = IOBuffer("   # comment \n    text")
    skipchars(isspace, io, linecomment='#')
    @test String(readavailable(io)) == "text"

    for char in ['@','߷','࿊','𐋺']
        io = IOBuffer("alphabeticalstuff$char")
        @test !eof(skipchars(isletter, io))
        @test read(io, Char) == char
    end
end

@testset "Test constructor with a generic type argument." begin
    io = IOBuffer(maxsize=Int16(10))
    @test io isa IOBuffer
    io = IOBuffer(maxsize=Int32(10))
    @test io isa IOBuffer
    io = IOBuffer(maxsize=Int64(10))
    @test io isa IOBuffer
end

@testset "# 25398 return value for write(::IO, ::IO)" begin
    ioa = IOBuffer()
    iob = IOBuffer("World")
    n = write(ioa, iob)
    @test String(take!(ioa)) == "World"
    @test n == 5
end

@testset "Base.compact" begin
    a = Base.GenericIOBuffer(UInt8[], true, true, false, true, typemax(Int))
    mark(a) # mark at position 0
    write(a, "Hello!")
    @test Base.compact(a) == nothing # because pointer > mark
    close(a)
    b = Base.GenericIOBuffer(UInt8[], true, true, false, true, typemax(Int))
    write(b, "Hello!")
    read(b)
    mark(b) # mark at position 6
    write(b, "Goodbye!") # now pointer is > mark but mark is > 0
    Base.compact(b)
    @test readline(b) == "Goodbye!"
    close(b)
end

@testset "peek(::GenericIOBuffer)" begin
    io = Base.GenericIOBuffer(UInt8[], true, true, false, true, typemax(Int))
    write(io, "こんにちは")
    @test peek(io) == 0xe3
    @test peek(io, Char) == 'こ'
    @test peek(io, Int32) == -476872221
    close(io)
end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
# ---
# title: 720. Longest Word in Dictionary
# id: problem720
# author: zhwang
# date: 2022-03-17
# difficulty: Easy
# categories: Hash Table, Trie
# link: <https://leetcode.com/problems/longest-word-in-dictionary/description/>
# hidden: true
# ---
#
# Given a list of strings `words` representing an English Dictionary, find the
# longest word in `words` that can be built one character at a time by other
# words in `words`. If there is more than one possible answer, return the
# longest word with the smallest lexicographical order.
#
# If there is no answer, return the empty string.
#
# **Example 1:**
#
#
#
#     Input:
#     words = ["w","wo","wor","worl", "world"]
#     Output: "world"
#     Explanation:
#     The word "world" can be built one character at a time by "w", "wo", "wor", and "worl".
#
#
# **Example 2:**
#
#
#
#     Input:
#     words = ["a", "banana", "app", "appl", "ap", "apply", "apple"]
#     Output: "apple"
#     Explanation:
#     Both "apply" and "apple" can be built from other words in the dictionary. However, "apple" is lexicographically smaller than "apply".
#
#
# **Note:**
#
# * All the strings in the input will only contain lowercase letters.
# * The length of `words` will be in the range `[1, 1000]`.
# * The length of `words[i]` will be in the range `[1, 30]`.
#
#
## @lc code=start
using LeetCode

## use Trie
longest_common_word(words::Vector{String}) = longest_common_word!(words)
function longest_common_word!(words::Vector{String})::String
    root, res = Dict{Char,Dict}(), ""
    for word in sort!(words; by=length)
        node, n = root, length(word)
        for (i, c) in enumerate(word)
            if haskey(node, c)
                node = node[c]
            elseif i == n
                node[c] = Dict{Char,Dict}() ## add new child
                res = word
            else
                break
            end
        end
    end
    return res
end

## @lc code=end

ext: jl


text:
include("experiment.jl")

returns, qtable = Experiment.run()
testhistory, qtable = Experiment.test(qtable)

open("states.csv", "w") do io
    for s in testhistory.states
        println(io, join(s, ","))
    end
end

open("actions.csv", "w") do io
    for a in testhistory.actions
        println(io, join(a, ","))
    end
end

open("rewards.csv", "w") do io
    for r in testhistory.rewards
        println(io, r)
    end
end

open("returns.csv", "w") do io
    for r in returns
        println(io, r)
    end
end

max_stars_repo_name: high-moctane/anylang-rl

max_issues_repo_name: high-moctane/anylang-rl

max_forks_repo_name: high-moctane/anylang-rl


text:
printstyled("Loading parameters ..."; color = :green);println(" ");
##################################################################

### Parameters
D1 = 3.2; #MAE of NiAl NiCp2
D2 = 4.2; #MAE of tip NiCp2
g1 = 1.75;
g2 = 1.65;
B = 0; # tesla T

### Constants
mu_B = 0.05788; # meV/T

println("Success loading parameters!");println(" ");
##################################################################


text:
module Lyceum

using Reexport
using Pkg

@reexport using LyceumBase
@reexport using LyceumMuJoCo
@reexport using LyceumAI
@reexport using MuJoCo
@reexport using Shapes

export LyceumBase, LyceumMuJoCo, LyceumAI, MuJoCo, Shapes

const LYCEUM_PACKAGES = [
    "LyceumBase",
    "LyceumMuJoCo",
    "LyceumAI",
    "MuJoCo",
    "Shapes",
]

# packages which are in the tomls but not aren't loaded (i.e. with using) by default
const UNUSED_LYCEUM_PACKAGES = [
    "LyceumMuJoCoViz" # so that `using Lyceum` works on headless machines
]

function pkgspecs(;ignore_versions::Bool = false, rev::Union{AbstractString, Nothing} = nothing)
    if !ignore_versions && rev !== nothing
        throw(ArgumentError("Cannot ignore versions and specify a revision at the same time"))
    end
    toml = Pkg.TOML.parsefile(joinpath(@__DIR__, "../Project.toml"))
    specs = Pkg.Types.PackageSpec[]
    for (name, uuid) in toml["deps"]
        if rev !== nothing
            spec = PackageSpec(name = name, uuid = uuid, rev = rev)
        elseif !ignore_versions && haskey(toml["compat"], name)
            v = toml["compat"][name]
            spec = PackageSpec(name = name, uuid = uuid, version = v)
        else
            spec = PackageSpec(name = name, uuid = uuid)
        end
        push!(specs, spec)
    end
    specs
end

function devall()
    lyceumpkgs = [LYCEUM_PACKAGES..., UNUSED_LYCEUM_PACKAGES...]
    specs = filter(s -> s.name in lyceumpkgs, pkgspecs(ignore_versions = true))
    Pkg.develop(specs)
end

function addall(;ignore_versions::Bool = false, rev::Union{AbstractString, Nothing} = nothing)
    lyceumpkgs = [LYCEUM_PACKAGES..., UNUSED_LYCEUM_PACKAGES...]
    specs = filter(s -> s.name in lyceumpkgs, pkgspecs(ignore_versions = ignore_versions, rev = rev))
    Pkg.add(specs)
end

end # module


text:
using JuMP
using SCS
using LinearAlgebra

###### Constants
length = 5

e = MathConstants.e
pi = MathConstants.pi

function get_A(n)
    A = zeros(Int32, n, n)
    println(size(A))
    for i in 1:1:size(A)[1] - 1
        A[i+1, i] = 1
    end
    return A
end

function get_B(n)
    B = zeros(Int32, n, 1)
    B[1, 1] = 1
    return B
end

function get_AB(n)
    A = get_A(n)
    B = get_B(n)
    return A, B
end

A, B = get_AB(length - 1)
show(stdout, MIME("text/plain"), A)
show(stdout, MIME("text/plain"), B)

###### Model
model = Model(SCS.Optimizer)
set_silent(model)

###### Variables
# @variable(model, P[1:length - 1, 1:length - 1], Symmetric)
@variable(model, X[1:length, 1:length], PSD)
@variable(model, r[1:length])
@variable(model, R[1:length])  # Fourier transoform of r

###### Constraints
for i in 1:1:size(R)[1]
    # @constraint(model, 1 <= r[i] <= 2)
    euler_fact = zeros(ComplexF64, 5, 1)
    for j in 1:1:size(euler_fact)[1]
        euler_fact[j] = e^(-im*2*pi/length*(j-1)*(i-1))
    end
    show(stdout, MIME("text/plain"), euler_fact)
    @constraint(model, R[i] == dot(r, euler_fact))
end

# # Up left
# @constraint(model, X[1:length-1, 1: length-1] .== P - transpose(A) * P * A)
# # Up right
# @constraint(model, X[1:length-1, length] .== transpose(r[1:length-1]) - transpose(A) * P * B)
# # Down left
# @constraint(model, X[length, 1:length-1] .== r[2:length] - transpose(B) * P * A)
# # Down right
# @constraint(model, X[length, length] == 2*r[1] - transpose(B) * P * B)

###### Objecitive Function
@objective(model, Min, sum(r) + tr(X))
JuMP.optimize!(model)

println("optimal value: ", objective_value(model))
# println("X: ")
show(stdout, MIME("text/plain"), JuMP.value.(X))
show(stdout, MIME("text/plain"), JuMP.value.(r))

# println(
#     "\nx: ", sqrt(value(X[1,1])),
#     "\ny: ", sqrt(value(X[2,2])),
#     "\nz: ", sqrt(value(X[3,3])))

max_stars_repo_name: troncosoae/opti-electric-project

max_issues_repo_name: troncosoae/opti-electric-project

max_forks_repo_name: troncosoae/opti-electric-project


text:
include("Testing.jl")
include("SPA.jl")

Wl = 2048
N = 10000

err = Testing.run_test(SPA.detect_pitch, Wl = Wl, N = N)
print(string(err*100, "% error"))

###########################################
using GLMakie, Makie.Colors
set_theme!(theme_black())
using WAV
using DSP
include("Utils.jl")
include("Testing.jl")
include("SPA.jl")
include("Segmenter.jl")

# Inputs and parameters
Wl = 2048

function get_song(file)
    y, Fs = wavread(string("Inputs/", file, ".wav"))
    low_note = Utils.getNote(1) #hz

    # Signal pre-processing
    y = [(y[i, 1] + y[i, 2]) / 2 for i = 1:length(y[:, 1])]
    function highpassfilter(signals, Fs, cutoff, order = 4)
        wdo = 2.0 * cutoff / Fs
        filth = digitalfilter(Highpass(wdo), Butterworth(order))
        filtfilt(filth, signals)
    end
    y = highpassfilter(y, Fs, low_note)
    return y, Fs
end

file = "Arpeggio"
y, Fs = get_song(file)
windows = Utils.getWindows(y, Wl)

#############
# ANALYSIS
#############
sound = []
block_length = Int(Wl * 0.125)
for i = 1:length(windows)
    window = Utils.hannWindow(windows[i])
    x_start = 1 + (i * Wl)
    x_end = (i * Wl) + Wl
    S = fft(window)[1:Int(Wl / 2)+1]
    push!(sound, Segmenter.Segment(:Pitch, x_start, x_end, SPA.detect_pitch(S, Fs)[1]))
end
sound = Segmenter.consolidate_segments(sound)

#############
# GRAPHING
#############
fig = Figure(resolution = (3840, 2160))
begin
    ax = Axis(fig[1, 1], xticks = MultiplesTicks(10, Wl, "f"))
    for segment in sound
        lines!(ax, Point2.(segment.x_start:segment.x_end, Utils.getNoteFromFreq.(float.(segment.value))), color = :green)
    end
end
display(fig)
# save(string("Outputs/", file, ".png"), fig)


text:
using PyPlot

function visPreprocessing(KS::KineticSolver, file::String)

    rho = readdlm(file, ',')
    #println("Final mass = $(sum(rho[:]))." )
    #nicewrite(string(KS.outputfolderprefix,"data/rhofinal.txt"),rho)

    if KS.ictype == "linesource"
        fig, ax = subplots(figsize = (10.5, 8), dpi = 100)
        minVal = 0.0
        maxVal = 0.5
        #	title("$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude), rot=$(KS.rotationmagnitude)", fontsize=20)
        pcolormesh(rho, cmap = "RdBu_r", vmin = minVal, vmax = maxVal)
        colorbar()
        PyPlot.savefig(string(KS.outputfolderprefix, "rho.png"))

        rhoana = readdlm("../vis/exactLineSource.txt", ',')
        #fig, ax = subplots(figsize=(10.5, 8), dpi=100)
        plot(rhoana[:, 1], rhoana[:, 2])
        rhong = rho[3:end-2, 3:end-2] # remove ghosts
        cuthori = rhong[Int64(KS.ny / 2), :]
        x = range(-1.5, stop = 1.5, length = KS.nx - 4)
        plot(x, cuthori)
        cutverti = rhong[:, Int64(KS.nx / 2)]
        plot(x, cutverti)
        cutdiag = [rhong[i, i] for i = 1:KS.nx-4]
        plot(sqrt(2) * x, cutdiag)
        cutdiag = [rhong[KS.ny-4+1-i, i] for i = 1:KS.nx-4]
        plot(sqrt(2) * x, cutdiag)
        title(
            "$(KS.Q.quadtype), nq=$(KS.nquadpoints), convolution=$(KS.convolutionmagnitude) ",
            fontsize = 30,
        )
        PyPlot.savefig(string(KS.outputfolderprefix, "rhocut.png"))

    elseif KS.ictype == "checkerboard"
        fig, ax = subplots(figsize = (10.5, 8), dpi = 100)
        minVal = -7.0
        maxVal = 0.0
        title(
            "$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude)",
            fontsize = 20,
        )
        pcolormesh(log10.(abs.(rho)), cmap = "RdBu_r", vmin = minVal, vmax = maxVal)
        colorbar()
        PyPlot.savefig(string(KS.outputfolderprefix, "rho_log10.png"))

    elseif KS.ictype == "beam"

        fig, ax = subplots(figsize = (10.5, 8), dpi = 100)
        minVal = -7.0
        maxVal = 0.0
        title(
            "$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude)",
            fontsize = 20,
        )
        pcolormesh(log10.(abs.(rho[3:end-2, 3:end-2])), cmap = "RdBu_r")#,vmin=minVal, vmax=maxVal)
        colorbar()
        PyPlot.savefig(string(KS.outputfolderprefix, "rho_log10.png"))

        fig, ax = subplots(figsize = (10.5, 8), dpi = 100)
        minVal = -7.0
        maxVal = 0.0
        title(
            "$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude)",
            fontsize = 20,
        )
        pcolormesh(((rho[3:end-2, 3:end-2])), cmap = "RdBu_r")#,vmin=minVal, vmax=maxVal)
        colorbar()
        PyPlot.savefig(string(KS.outputfolderprefix, "rho.png"))

    end
end

max_stars_repo_name: vavrines/CSN

max_issues_repo_name: vavrines/CSN

max_forks_repo_name: vavrines/CSN


text:
include("godunov.jl")
using .Godunov
using Plots
fig = Plots.plot()

"""
NN = 100 # no of cells
Nt = 64000000 # no of time steps
Np = 6400000 # how often to plot
Ns = 5000 # when to start plotting
"""
NN = 40 # no of cells
#Nt = 8000 # no of time steps
#Np = 1000 # how often to plot
#Ns = 1#5000 # when to start plotting
tmax = 1.0e3 # max time
dt = 1.0e2 # how often to plot
tp = 0.0e0#5000 # when to start plotting

r=Godunov.ring(NN)#,xi=0.4,xf=3.0)
Godunov.getAlphas!(r)
#for i in 1:10
for tt in 0.0:dt:tmax
    tx = r.t
    if tx ≥ tp
        Plots.plot!(r.X_[1:NN÷1],r.Σ_[1:NN÷1],label="t: $tx",legend=false,xlim=(0,3))#:bottomright)#,title="Σ")#;legend=false)
        display(fig)
        print("time: $tx\n")
        Plots.plot!(r.X_[1:NN÷1],r.Φ_[1:NN÷1],label=false,linestyle=:dash)#,label="t: $t",legend=:topleft)#,title="J")#;legend=false)
        display(fig)
    end
#    for j in 1:800
        Godunov.step!(r)#,dt)
#    end
    Godunov.tstep!(r,dt)
end
tf=r.t
print("Done at t= $tf\n")

max_stars_repo_name: diakosmos/razor

max_issues_repo_name: diakosmos/razor

max_forks_repo_name: diakosmos/razor


text:
using Test
using HMMBase
using Distributions
using Random

Random.seed!(2019)

hmms = [
    HMM([0.9 0.1; 0.1 0.9], [Normal(10, 1), Gamma(1, 1)]),
    HMM([0.9 0.1; 0.1 0.9], [Categorical([0.1, 0.2, 0.7]), Categorical([0.5, 0.5])]),
    HMM(
        [0.9 0.1; 0.1 0.9],
        [MvNormal([0.0, 0.0], [1.0, 1.0]), MvNormal([10.0, 10.0], [1.0, 1.0])],
    ),
]

@testset "Integration $(typeof(hmm))" for hmm in hmms, T in [0, 1, 1000]
    K = size(hmm, 1)

    # HMM API
    @test hmm !== copy(hmm)

    z, y = rand(hmm, T, seq = true)
    @test size(z, 1) == size(y, 1)
    @test size(y, 2) == size(hmm, 2)

    yp = rand(hmm, z)
    @test size(z, 1) == size(y, 1)
    @test size(y, 2) == size(hmm, 2)

    LL = loglikelihoods(hmm, y)
    @test size(LL) == (T, K)

    # Forward/Backward
    α1, logtot1 = forward(hmm, y)
    β1, logtot2 = backward(hmm, y)
    γ1 = posteriors(hmm, y)
    logtot3 = loglikelihood(hmm, y)

    @test size(α1) == size(β1) == size(γ1)
    @test logtot1 ≈ logtot2 ≈ logtot3

    # Viterbi
    zv1 = viterbi(hmm, y)
    @test size(zv1) == size(z)

    # MLE
    if T > 2
        hmm2, _ = fit_mle(hmm, y, maxiter = 1, display = :iter)
        @test size(hmm2) == size(hmm)
        @test typeof(hmm2) == typeof(hmm)

        hmm2, _ = fit_mle(hmm, y, init = nothing)
        @test size(hmm2) == size(hmm)
        @test typeof(hmm2) == typeof(hmm)

        hmm2, _ = fit_mle(hmm, y, init = :kmeans, robust = true)
        @test size(hmm2) == size(hmm)
        @test typeof(hmm2) == typeof(hmm)
    end
end


text:
j = 30
maxn = 70
ω = 1.0
ω₀ = 1.0
δ = 1.0

λ = 1.0
λᵪ= sqrt(ω * ω₀) / (1.0 + δ)

G = 2.0 * λ / (sqrt(2.0 * j) * ω)
logG = log(G)

function Index(n, m)
    return n * (2.0 * j + 1) + m + j + 1
end

function LogFactorial(n)
    if n <= 1
        return 0
    else
        return log(n) + LogFactorial(n - 1)
    end
end

function Minus(n)
    if n % 2 == 0
        return 1
    else
        return -1
    end
end

function fu(n1, n2)
    coef = 0.5 * (LogFactorial(n1) + LogFactorial(n2) - G*G)

    sum = 0
    for r = 0:min(n1, n2)
        sum += Minus(r) * exp(logG * (n1 + n2 - 2.0 * r) - LogFactorial(n1 - r) - LogFactorial(n2 - r) - LogFactorial(r))
    end

    return sum
end

function Fuf(np, mp, n, m)
    if m + 1 == mp
        return -0.5 * ω₀ * sqrt(j * (j + 1) - m * (m + 1)) * Minus(n) * fu(np + 1, n + 1)
    elseif m - 1 == mp
        return -0.5 * ω₀ * sqrt(j * (j + 1) - m * (m - 1)) * Minus(np) * fu(np + 1, n + 1)
    end

    return 0.0
end

max_stars_repo_name: PavelStransky/ClassicalDynamics

max_issues_repo_name: PavelStransky/ClassicalDynamics

max_forks_repo_name: PavelStransky/ClassicalDynamics


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

#     CurvedFace{RealT<:Real}
#
# Contains the data needed to represent a curved face with data points (x,y,z) as a Lagrange polynomial
# interpolant written in barycentric form at a given set of nodes.
struct CurvedFace{RealT<:Real}
  nodes               ::Vector{RealT}
  barycentric_weights ::Vector{RealT}
  coordinates         ::Array{RealT, 3} #[ndims, nnodes, nnodes]
end

# evalute the Gamma face interpolant at a particular point s = (s_1, s_2) and return the (x,y,z) coordinate
function evaluate_at(s, boundary_face::CurvedFace)

   @unpack nodes, barycentric_weights, coordinates = boundary_face

   x_coordinate_at_s_on_boundary_face = lagrange_interpolation_2d(s, nodes, view(coordinates, 1, :, :),
                                                                  barycentric_weights)
   y_coordinate_at_s_on_boundary_face = lagrange_interpolation_2d(s, nodes, view(coordinates, 2, :, :),
                                                                  barycentric_weights)
   z_coordinate_at_s_on_boundary_face = lagrange_interpolation_2d(s, nodes, view(coordinates, 3, :, :),
                                                                  barycentric_weights)

   return x_coordinate_at_s_on_boundary_face,
          y_coordinate_at_s_on_boundary_face,
          z_coordinate_at_s_on_boundary_face
end

# Calculate a 2D Lagrange interpolating polynomial in barycentric 2 form
# of a function f(x,y) at a given coordinate (x,y) for a given node distribution.
function lagrange_interpolation_2d(x, nodes, function_values, barycentric_weights)

  f_intermediate = zeros(eltype(function_values), length(nodes))
  for j in eachindex(nodes)
    f_intermediate[j] = lagrange_interpolation(x[2], nodes, view(function_values, j, :),
                                               barycentric_weights)
  end
  point_value = lagrange_interpolation(x[1], nodes, f_intermediate, barycentric_weights)

  return point_value
end

end # @muladd


text:
mutable struct Options
    width::Int
    height::Int
    negate::Bool
    simplify::Bool
end

width(opt::Options) = opt.width
height(opt::Options) = opt.height

function set_default_size(width, height)
    _options.width = width
    _options.height = height
end

function set_options(; kwargs...)
    for opt in fieldnames(Options)
        if opt in keys(kwargs)
            setfield!(_options, opt, kwargs[opt])
        end
    end
end

max_stars_repo_name: davidanthoff/ProfileVega

max_issues_repo_name: davidanthoff/ProfileVega

max_forks_repo_name: davidanthoff/ProfileVega


text:
#= -------------------------------------------------------------------------
# @file problema2_PID_transfer_functions.jl
#
# @date 12/09/17 12:54:35
# @author Martin Noblia
# @email mnoblia@disroot.org
#
# @brief
#
# @detail
#
#  Licence:
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
---------------------------------------------------------------------------=#
using ControlSystems


text:
using Test
using MARY_fRG.Basics

@testset "点相关的测试" begin
    pt1 = Point2D(1., 1.)
    println(pt1)
    pt2 = Point2D(2., 3.)
    pt3 = pt1 + pt2
    @test (pt3.y == 4) && (pt3.x == 3)
    pt4 = middle_point(pt1, pt3)
    @test (pt4.x == 2) && (pt4.y == 2.5)
    ang1 = absolute_angle(pt1)
    @test ang1 == pi / 4
    ang1 = absolute_angle(Point2D(-1., 1.))
    @test ang1 == 3pi / 4
    ang1 = absolute_angle(Point2D(1., -1.))
    @test ang1 == 7pi / 4
    pt5 = -pt1
    println(pt5)
    pt6 = pt5 - pt1
    println(pt6)
    pt7 = pt5 - 2pt1
    println(pt7)
end


text:
module WT
export  DiscreteWavelet,
        FilterWavelet,
        LSWavelet,
        OrthoFilter,
        GLS,
        wavelet

using ..Util
using Compat
import Compat: ASCIIString

# TYPE HIERARCHY

abstract DiscreteWavelet{T}
#abstract ContinuousWavelet{T}
# discrete transforms via filtering
abstract FilterWavelet{T} <: DiscreteWavelet{T}
# discrete transforms via lifting
abstract LSWavelet{T} <: DiscreteWavelet{T}
# all wavelet types
#typealias WaveletTransformType Union{DiscreteWavelet, ContinuousWavelet}

"""Get wavelet type name."""
function name(::DiscreteWavelet) end
"""Get wavelet filter length."""
function Base.length(::FilterWavelet) end

immutable FilterTransform end
immutable LiftingTransform end
"""Transform by filtering."""
const Filter = FilterTransform()
"""Transform by lifting."""
const Lifting = LiftingTransform()

# BOUNDARY TYPES

abstract WaveletBoundary
# periodic (default)
immutable PerBoundary <: WaveletBoundary end
# zero padding
#immutable ZPBoundary <: WaveletBoundary end
# constant padding
#immutable CPBoundary <: WaveletBoundary end
# and so on...

const Periodic = PerBoundary()
const DEFAULT_BOUNDARY = PerBoundary()

# WAVELET CLASSES

"""
The `WaveletClass` type has subtypes `OrthoWaveletClass`
and `BiOrthoWaveletClass`.

The `WT` module has for convenience constants defined named
as the class short name and optionally amended with a number
specifing the number of vanishing moments.

A class can also be explicitly constructed as e.g. `Daubechies{4}()`.

# Examples
`WT.db2`, `WT.haar`, `WT.cdf97`

**See also:** `WT.class`, `WT.name`, `WT.vanishingmoments`
"""
abstract WaveletClass
abstract OrthoWaveletClass <: WaveletClass
abstract BiOrthoWaveletClass <: WaveletClass

# Single classes
for (TYPE, NAMEBASE, MOMENTS) in (
        (:Haar, "haar", 1),
        (:Beylkin, "beyl", -1), # TODO moments
        (:Vaidyanathan, "vaid",-1), # TODO moments
        )
    @eval begin
        immutable $TYPE <: OrthoWaveletClass end
        class(::$TYPE) = $(string(TYPE))
        name(::$TYPE) = string($NAMEBASE)
        vanishingmoments(::$TYPE) = $MOMENTS
    end
    CONSTNAME = Symbol(NAMEBASE)
    @eval begin
        const $CONSTNAME = $TYPE()                  # type shortcut
    end
end

# Parameterized classes
for (TYPE, NAMEBASE, RANGE) in (
        (:Daubechies, "db", 1:10),
        (:Coiflet, "coif", 2:2:8),
        (:Symlet, "sym", 4:10),
        (:Battle, "batt", 2:2:6),
        )
    @eval begin
        immutable $TYPE{N} <: OrthoWaveletClass end
        class(::$TYPE) = $(string(TYPE))
        name{N}(::$TYPE{N}) = string($NAMEBASE,N)
        vanishingmoments{N}(::$TYPE{N}) = N
    end
    for NUM in RANGE
        CONSTNAME = Symbol(string(NAMEBASE, NUM))
        @eval begin
            const $CONSTNAME = $TYPE{$NUM}()        # type shortcut
        end
    end
end

# Parameterized BiOrtho classes
for (TYPE, NAMEBASE, RANGE1, RANGE2) in (
        (:CDF, "cdf", [9], [7]),
        )
    @eval begin
        immutable $TYPE{N1, N2} <: BiOrthoWaveletClass end
        class(::$TYPE) = $(string(TYPE))
        name{N1, N2}(::$TYPE{N1, N2}) = string($NAMEBASE,N1,"/",N2)
        vanishingmoments{N1, N2}(::$TYPE{N1, N2}) = (N1, N2)
    end
    for i in length(RANGE1)
        CONSTNAME = Symbol(string(NAMEBASE,RANGE1[i],RANGE2[i]))
        @eval begin
            const $CONSTNAME = $TYPE{$RANGE1[$i],$RANGE2[$i]}()        # type shortcut
        end
    end
end

# IMPLEMENTATIONS OF FilterWavelet

"""
Wavelet type for discrete orthogonal transforms by filtering.

**See also:** `GLS`, `wavelet`
"""
immutable OrthoFilter{T<:WaveletBoundary} <: FilterWavelet{T}
    qmf     ::Vector{Float64}        # quadrature mirror filter
    name    ::ASCIIString            # filter short name
end

@compat function (::Type{OrthoFilter}){WC<:OrthoWaveletClass, T<:WaveletBoundary}(w::WC, ::T=DEFAULT_BOUNDARY)
    name = WT.name(w)
    if WC <: Daubechies
        qmf = daubechies(vanishingmoments(w))
    else
        qmf = get(FILTERS, name, nothing)
        qmf == nothing && throw(ArgumentError("filter not found"))
    end
    # make sure it is normalized in l2-norm
    return OrthoFilter{T}(qmf./norm(qmf), name)
end

Base.length(f::OrthoFilter) = length(f.qmf)
qmf(f::OrthoFilter) = f.qmf
name(f::OrthoFilter) = f.name

"""Scale filter by scalar."""
function scale{T<:WaveletBoundary}(f::OrthoFilter{T}, a::Number)
    return OrthoFilter{T}(f.qmf.*a, f.name)
end

"""Quadrature mirror filter pair."""
function makeqmfpair(f::OrthoFilter, fw::Bool=true, T::Type=eltype(qmf(f)))
    scfilter, dcfilter = makereverseqmfpair(f, fw, T)
    return reverse(scfilter), reverse(dcfilter)
end

"""Reversed quadrature mirror filter pair."""
function makereverseqmfpair(f::OrthoFilter, fw::Bool=true, T::Type=eltype(qmf(f)))
    h = convert(Vector{T}, qmf(f))
    if fw
        scfilter = reverse(h)
        dcfilter = mirror(h)
    else
        scfilter = h
        dcfilter = reverse(mirror(h))
    end
    return scfilter, dcfilter
end

#immutable BiOrthoFilter{T<:WaveletBoundary} <: FilterWavelet{T}
#    qmf1::Vector{Float64}       # quadrature mirror filter 1
#    qmf2::Vector{Float64}       # quadrature mirror filter 2
#    name::String                # filter short name
#    BiOrthoFilter(qmf1, qmf2, name) = new(qmf1, qmf2, name)
#end

# IMPLEMENTATIONS OF LSWavelet

abstract StepType
immutable PredictStep <: StepType end
immutable UpdateStep <: StepType end
const Predict = PredictStep()
const Update = UpdateStep()

immutable LSStepParam{T<:Number}
    coef    ::Vector{T}        # lifting coefficients
    shift   ::Int              # + left shift, - right shift
end

immutable LSStep{T<:Number}
    param::LSStepParam{T}
    steptype::StepType
end

@compat function (::Type{LSStep}){T}(st::StepType, coef::Vector{T}, shift::Int)
    return LSStep{T}(LSStepParam{T}(coef, shift), st)
end

Base.length(s::LSStep) = length(s.param)
Base.length(s::LSStepParam) = length(s.coef)

"""
Wavelet type for discrete general (bi)orthogonal transforms
by using a lifting scheme.

**See also:** `OrthoFilter`, `wavelet`
"""
immutable GLS{T<:WaveletBoundary} <: LSWavelet{T}
    step    ::Vector{LSStep{Float64}}    # steps to be taken
    norm1   ::Float64           # normalization of scaling coefs.
    norm2   ::Float64           # normalization of detail coefs.
    name    ::ASCIIString       # name of scheme
end

@compat function (::Type{GLS}){WC<:WaveletClass, T<:WaveletBoundary}(w::WC, ::T=DEFAULT_BOUNDARY)
    name = WT.name(w)
    schemedef = get(SCHEMES, name, nothing)
    schemedef == nothing && throw(ArgumentError("scheme not found"))
    return GLS{T}(schemedef[1], schemedef[2], schemedef[3], name)
end

name(s::GLS) = s.name

# IMPLEMENTATIONS OF ContinuousWavelet

# ...

# TRANSFORM TYPE CONSTRUCTORS

"""
    wavelet(c[, t=WT.Filter][, boundary=WT.Periodic])

Construct wavelet type where `c` is a wavelet class,
`t` is the transformation type (`WT.Filter` or `WT.Lifting`),
and `boundary` is the type of boundary treatment.

# Examples
```julia
wavelet(WT.coif6)
wavelet(WT.db1, WT.Lifting)
```

**See also:** `WT.WaveletClass`
"""
function wavelet end

wavelet(c::WaveletClass, boundary::WaveletBoundary=DEFAULT_BOUNDARY) = wavelet(c, Filter, boundary)
wavelet(c::OrthoWaveletClass, t::FilterTransform, boundary::WaveletBoundary=DEFAULT_BOUNDARY) = OrthoFilter(c, boundary)
wavelet(c::WaveletClass, t::LiftingTransform, boundary::WaveletBoundary=DEFAULT_BOUNDARY) = GLS(c, boundary)

# ------------------------------------------------------------

# Compute filters from the Daubechies class
# N is the number of zeros at -1
function daubechies(N::Int)
    @assert N > 0
    # Create polynomial
    C = Array(Int, N)
    @inbounds for n = 0:N-1
        C[N-n] = binomial(N-1+n, n)
    end

    # Find roots in y domain (truncated binomial series; (1 - y)^{-N})
    Y = roots(C)

    # Find roots in z domain:
    # z + z^{-1} = 2 - 4*y
    # where y is a root from above
    Z = zeros(Complex128, 2*N-2)
    @inbounds for i = 1:N-1
        Yi = Y[i]
        d = 2*sqrt( Yi*Yi - Yi )
        y2 = 1 - 2*Yi
        Z[i] = y2 + d
        Z[i+N-1] = y2 -d
    end

    # Retain roots inside unit circle
    nr = 0  # count roots
    @inbounds for i = eachindex(Z)
        if abs(Z[i]) <= 1 + eps()
            nr += 1
        end
    end

    # Find coefficients of the polynomial
    # (1 + z)^N * \prod_i (z - z_i)
    R = Array(Complex128, N+nr)
    @inbounds for i = 1:N
        R[i] = -1
    end
    k = N
    @inbounds for i = eachindex(Z)
        if abs(Z[i]) <= 1 + eps()
            k += 1
            R[k] = Z[i]
        end
    end
    HH = vieta( R )

    # Normalize coefficients
    scale!(HH, 1/norm(HH))
    return real(HH)
end

# Compute roots of polynomial
# Input is a coefficient vector with highest powers first
function roots(C::AbstractVector)
    A = compan(C)
    return eigvals(A)
end

# Create companion matrix for a polynomial
# Input is a coefficient vector with highest powers first
function compan(C::AbstractVector)
    n = length(C)
    A = zeros(n-1, n-1)

    if n > 1
        @inbounds A[1,:] = -C[2:end] ./ C[1]
        @inbounds A[2:n:end] = 1
    end
    return A
end

# Vieta-like formula for computing polynomial coefficients from roots
# See
# http://www.mathworks.se/help/matlab/ref/poly.html
function vieta(R::AbstractVector)
    n = length( R )
    C = zeros(Complex128, n+1)
    C[1] = 1
    Ci::Complex128 = 0
    Cig::Complex128 = 0

    @inbounds for k = 1:n
        Ci = C[1]
        for i = 1:k
            Cig = C[i+1]
            C[i+1] = Cig - R[k] * Ci
            Ci = Cig
        end
    end
    return C
end

# scaling filters h (low pass)
# the number at end of a filter name is the
# number of vanishing moments of the mother wavelet function
# sources:
# http://statweb.stanford.edu/~wavelab/ (Orthogonal/MakeONFilter.m)
# http://www.mathworks.com/matlabcentral/fileexchange/5502-filter-coefficients-to-popular-wavelets
### https://github.com/nigma/pywt/blob/master/src/wavelets_coeffs.template.h
# name => qmf
const FILTERS = @compat Dict{ASCIIString,Vector{Float64}}(
# Haar filter
"haar" =>
[0.7071067811865475,0.7071067811865475]
,
# Daubechies filters, see daubechies()

# Coiflet filters
"coif2" =>
[-0.072732619513,0.337897662458,0.852572020212,0.384864846864,-0.072732619513,-0.015655728135]
,
"coif4" =>
[0.0163873364635998,-0.0414649367819558,-0.0673725547222826,0.3861100668229939,0.8127236354493977,0.4170051844236707,-0.0764885990786692,-0.0594344186467388,0.0236801719464464,0.0056114348194211,-0.0018232088707116,-0.0007205494453679]
,
"coif6" =>
[-0.0037935128644910,0.0077825964273254,0.0234526961418362,-0.0657719112818552,-0.0611233900026726,0.4051769024096150,0.7937772226256169,0.4284834763776168,-0.0717998216193117,-0.0823019271068856,0.0345550275730615,0.0158805448636158,-0.0090079761366615,-0.0025745176887502,0.0011175187708906,0.0004662169601129,-0.0000709833031381,-0.0000345997728362]
,
"coif8" =>
[0.0008923136685824,-0.0016294920126020,-0.0073461663276432,0.0160689439647787,0.0266823001560570,-0.0812666996808907,-0.0560773133167630,0.4153084070304910,0.7822389309206135,0.4343860564915321,-0.0666274742634348,-0.0962204420340021,0.0393344271233433,0.0250822618448678,-0.0152117315279485,-0.0056582866866115,0.0037514361572790,0.0012665619292991,-0.0005890207562444,-0.0002599745524878,0.0000623390344610,0.0000312298758654,-0.0000032596802369,-0.0000017849850031]
,
"coif10" =>
[-0.0002120808398259,0.0003585896879330,0.0021782363583355,-0.0041593587818186,-0.0101311175209033,0.0234081567882734,0.0281680289738655,-0.0919200105692549,-0.0520431631816557,0.4215662067346898,0.7742896037334738,0.4379916262173834,-0.0620359639693546,-0.1055742087143175,0.0412892087544753,0.0326835742705106,-0.0197617789446276,-0.0091642311634348,0.0067641854487565,0.0024333732129107,-0.0016628637021860,-0.0006381313431115,0.0003022595818445,0.0001405411497166,-0.0000413404322768,-0.0000213150268122,0.0000037346551755,0.0000020637618516,-0.0000001674428858,-0.0000000951765727]
,
# Symmlet filter
"sym4" =>
[0.0455703458960000,-0.0178247014420000,-0.1403176241790000,0.4212345342040000,1.1366582434079999,0.7037390686560000,-0.0419109651250000,-0.1071489014180000]
,
"sym5" =>
[0.0276321529580000,-0.0298424998690000,-0.2479513626130000,0.0234789231360000,0.8965816483800000,1.0230529668940000,0.2819906968540000,-0.0553441861170000,0.0417468644220000,0.0386547959550000]
,
"sym6" =>
[-0.0110318675090000,0.0024999220930000,0.0632505626600000,-0.0297837512990000,-0.1027249698620000,0.4779043713330000,1.1138927839260000,0.6944579729580000,-0.0683231215870000,-0.1668632154120000,0.0049366123720000,0.0217847003270000]
,
"sym7" =>
[0.0145213947620000,0.0056713426860000,-0.1524638718960000,-0.1980567068070000,0.4081839397250000,1.0857827098140000,0.7581626019640000,0.0246656594890000,-0.0700782912220000,0.0960147679360000,0.0431554525820000,-0.0178704316510000,-0.0014812259150000,0.0037926585340000]
,
"sym8" =>
[-0.0047834585120000,-0.0007666908960000,0.0448236230420000,0.0107586117510000,-0.2026486552860000,-0.0866536154060000,0.6807453471900000,1.0991066305370001,0.5153986703740000,-0.0734625087610000,-0.0384935212630000,0.0694904659110000,0.0053863887540000,-0.0211456865280000,-0.0004283943000000,0.0026727933930000]
,
"sym9" =>
[0.0019811937360000,0.0008765025390000,-0.0187693968360000,-0.0163033512260000,0.0427444336020000,0.0008251409290000,-0.0771721610970000,0.3376589236020000,1.0152597908320000,0.8730484073490000,0.0498828309590000,-0.2708937835030000,-0.0257864459300000,0.0877912515540000,0.0125288962420000,-0.0145155785530000,-0.0006691415090000,0.0015124873090000]
,
"sym10" =>
[-0.0006495898960000,0.0000806612040000,0.0064957283750000,-0.0011375353140000,-0.0287862319260000,0.0081528167990000,0.0707035675500000,-0.0452407722180000,-0.0502565400920000,0.5428130112130000,1.0882515305000000,0.6670713381540000,-0.1002402150310000,-0.2255589722340000,0.0164188694260000,0.0649509245790000,-0.0020723639230000,-0.0122206426300000,0.0001352450200000,0.0010891704470000]
,
# Battle-Lemarie filter
"batt2" =>
[-0.0000867523000000,-0.0001586010000000,0.0003617810000000,0.0006529220000000,-0.0015570100000000,-0.0027458800000000,0.0070644200000000,0.0120030000000000,-0.0367309000000000,-0.0488618000000000,0.2809310000000000,0.5781630000000000,0.2809310000000000,-0.0488618000000000,-0.0367309000000000,0.0120030000000000,0.0070644200000000,-0.0027458800000000,-0.0015570100000000,0.0006529220000000,0.0003617810000000,-0.0001586010000000,-0.0000867523000000]
,
"batt4" =>
[0.0001033070000000,-0.0001642640000000,-0.0002018180000000,0.0003267490000000,0.0003959460000000,-0.0006556200000000,-0.0007804680000000,0.0013308600000000,0.0015462400000000,-0.0027452900000000,-0.0030786300000000,0.0057993200000000,0.0061414300000000,-0.0127154000000000,-0.0121455000000000,0.0297468000000000,0.0226846000000000,-0.0778079000000000,-0.0354980000000000,0.3068300000000000,0.5417360000000000,0.3068300000000000,-0.0354980000000000,-0.0778079000000000,0.0226846000000000,0.0297468000000000,-0.0121455000000000,-0.0127154000000000,0.0061414300000000,0.0057993200000000,-0.0030786300000000,-0.0027452900000000,0.0015462400000000,0.0013308600000000,-0.0007804680000000,-0.0006556200000000,0.0003959460000000,0.0003267490000000,-0.0002018180000000,-0.0001642640000000,0.0001033070000000]
,
"batt6" =>
[0.0001011130000000,0.0001107090000000,-0.0001591680000000,-0.0001726850000000,0.0002514190000000,0.0002698420000000,-0.0003987590000000,-0.0004224850000000,0.0006355630000000,0.0006628360000000,-0.0010191200000000,-0.0010420700000000,0.0016465900000000,0.0016413200000000,-0.0026864600000000,-0.0025881600000000,0.0044400200000000,0.0040788200000000,-0.0074684800000000,-0.0063988600000000,0.0128754000000000,0.0099063500000000,-0.0229951000000000,-0.0148537000000000,0.0433544000000000,0.0208414000000000,-0.0914068000000000,-0.0261771000000000,0.3128690000000000,0.5283740000000000,0.3128690000000000,-0.0261771000000000,-0.0914068000000000,0.0208414000000000,0.0433544000000000,-0.0148537000000000,-0.0229951000000000,0.0099063500000000,0.0128754000000000,-0.0063988600000000,-0.0074684800000000,0.0040788200000000,0.0044400200000000,-0.0025881600000000,-0.0026864600000000,0.0016413200000000,0.0016465900000000,-0.0010420700000000,-0.0010191200000000,0.0006628360000000,0.0006355630000000,-0.0004224850000000,-0.0003987590000000,0.0002698420000000,0.0002514190000000,-0.0001726850000000,-0.0001591680000000,0.0001107090000000,0.0001011130000000]
,
# Beylkin filter
"beyl" =>
[0.0993057653740000,0.4242153608130000,0.6998252140570000,0.4497182511490000,-0.1109275983480000,-0.2644972314460000,0.0269003088040000,0.1555387318770000,-0.0175207462670000,-0.0885436306230000,0.0196798660440000,0.0429163872740000,-0.0174604086960000,-0.0143658079690000,0.0100404118450000,0.0014842347820000,-0.0027360316260000,0.0006404853290000]
,
# Vaidyanathan filter
"vaid" =>
[-0.0000629061180000,0.0003436319050000,-0.0004539566200000,-0.0009448971360000,0.0028438345470000,0.0007081375040000,-0.0088391034090000,0.0031538470560000,0.0196872150100000,-0.0148534480050000,-0.0354703986070000,0.0387426192930000,0.0558925236910000,-0.0777097509020000,-0.0839288843660000,0.1319716614170000,0.1350842271290000,-0.1944504717660000,-0.2634948024880000,0.2016121617750000,0.6356010598720000,0.5727977932110000,0.2501841295050000,0.0457993341110000]

)

# biortho filters
# name => (qmf1,qmf2)
const BIFILTERS = @compat Dict{ASCIIString,NTuple{2,Vector{Float64}}}(
# test
"test" =>
([0.7071067811865475,0.7071067811865475],
 [0.7071067811865475,0.7071067811865475])

)

# in matlab (d,p) -> (predict, update)
const SCHEMES = @compat Dict{ASCIIString,NTuple{3}}(
# cdf 5/3 -> bior 2.2, cdf 9/7 -> bior 4.4
# Cohen-Daubechies-Feauveau [Do Quan & Yo-Sung Ho. Optimized median lifting scheme for lossy image compression.]
"cdf9/7" => ([  LSStep(Update,  [1.0,1.0]*1.5861343420604, 0),
                LSStep(Predict, [1.0,1.0]*0.05298011857291494, 1),
                LSStep(Update,  [1.0,1.0]*-0.882911075531393, 0),
                LSStep(Predict, [1.0,1.0]*-0.44350685204384654, 1)],
                1.1496043988603355,
                0.8698644516247099)
,

# Haar, same as db1
"haar" => ([LSStep(Predict, [-1.0], 0),
            LSStep(Update,  [0.5], 0)],
            0.7071067811865475,
            1.4142135623730951)
,
# Daubechies
"db1" => ([ LSStep(Predict, [-1.0], 0),
            LSStep(Update,  [0.5], 0)],
            0.7071067811865475,
            1.4142135623730951)
,
"db2" => ([ LSStep(Predict, [-1.7320508075688772], 0),
            LSStep(Update,  [-0.0669872981077807,0.4330127018922193], 1),
            LSStep(Predict, [1.0], -1)],
            0.5176380902050414,
            1.9318516525781364)

)

end


text:
#=
Directly bridging SV in GT -> LT triggers
=#
MOI.supports_constraint(::Optimizer{T},
                        ::Type{<:Union{SV,SAF{T},SQF{T}}},
                        ::Type{<:Union{LT{T},GT{T},ET{T},IT{T}}},
                        ) where T <: Real = true

# No automatic bridge from ZO to Integers exist. Closest is semi-integer to
# ZO + 2x SAF.
MOI.supports_constraint(::Optimizer,
                        ::Type{SV},
                        ::Type{<:Union{ZO, MOI.Integer}}) = true

MOI.supports_constraint(::Optimizer,
                        ::Type{<:Union{VECVAR}},
                        ::Type{<:Union{MOI.SOS1, MOI.SOS2, SOC_CONE, EXP_CONE, POW_CONE{T}, PSD_CONE}},
                        ) where T <: Real = true

# Using SingleVariable or ScalarAffineFunction as the objective greatly
# simplifies objective bound cuts and access to objective values. Using
# SAF prevents the introduce of a new auxilliary variable if the user
# defines a SAF objective type.
MOI.supports(::Optimizer{T},
             ::Union{MOI.ObjectiveSense,
                     MOI.ObjectiveFunction{SAF{T}}}
                     ) where T <: Real = true

function MOI.copy_to(model::Optimizer, src::MOI.ModelLike; copy_names = false)
    return MOIU.default_copy_to(model, src, copy_names)
end

MOI.get(m::Optimizer, ::MOI.ObjectiveFunctionType) = MOI.get(m._model._input_model, MOI.ObjectiveFunctionType())
MOI.get(m::Optimizer, ::MOI.SolverName) = "EAGO: Easy Advanced Global Optimization"
MOI.get(m::Optimizer, ::MOI.TerminationStatus) = m._termination_status_code
function MOI.get(m::Optimizer, v::MOI.PrimalStatus)
    if v.N > 1
        MOI.check_result_index_bounds(m, v)
    end
    m._primal_status_code
end
MOI.get(m::Optimizer{T}, ::MOI.SolveTime) where T <: Real = m._run_time
MOI.get(m::Optimizer, ::MOI.NodeCount) = m._node_count
MOI.get(m::Optimizer, ::MOI.ResultCount) = (m._primal_status_code === MOI.FEASIBLE_POINT) ? 1 : 0

MOI.supports(::Optimizer, ::MOI.Silent) = true
function MOI.set(m::Optimizer, ::MOI.Silent, b::Bool)
    m._is_silent = b
    return
end
MOI.get(m::Optimizer, ::MOI.Silent) = m._is_silent

MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
MOI.get(m::Optimizer{T}, ::MOI.TimeLimitSec)  where T <: Real = isfinite(m._time_limit) ? m._time_limit : nothing
MOI.set(m::Optimizer{T}, ::MOI.TimeLimitSec, ::Nothing) where T <: Real  = (m._time_limit = typemax(T);   return)
MOI.set(m::Optimizer{T}, ::MOI.TimeLimitSec, v::T)      where T <: Real = (m._time_limit = v; return)

function MOI.set(m::Optimizer, p::MOI.RawParameter, value)
    if !(p.name isa String) && !(p.name isa Symbol)
        error("EAGO only supports raw parameters with Symbol or String names.")
    end
    setfield!(m, Symbol(p.name), value)

    return nothing
end

for attr in (MOI.ListOfVariableIndices, MOI.ListOfConstraints)
    @eval function MOI.get(m::Optimizer, ::$attr)
        MOI.get(m._model._input_model, $attr())
    end
end

for attr in (MOI.ListOfConstraintAttributesSet, MOI.ListOfConstraintIndices)
    @eval function MOI.get(d::Optimizer, ::$attr{F,S}) where {F,S}
        return MOI.get(d._model._input_model, $attr{F,S}())
    end
end

for attr in (MOI.ConstraintFunction, MOI.ConstraintSet)
    for S in (MOI.LessThan,MOI.GreaterThan,MOI.EqualTo,MOI.Interval)
        @eval function MOI.set(d::Optimizer, ::$attr, ci::CI{SV,$S{T}}, v) where T<:Real
            return MOI.set(d._model._input_model, $attr(), ci, v)
        end
        @eval function MOI.get(d::Optimizer, ::$attr, ci::CI{SV,$S{T}}) where T<:Real
            return MOI.get(d._model._input_model, $attr(), ci)
        end
        for F in (MOI.ScalarAffineFunction,MOI.ScalarQuadraticFunction)
            @eval function MOI.set(d::Optimizer, ::$attr, ci::CI{$F{T},$S{T}}, v) where T<:Real
                return MOI.set(d._model._input_model, $attr(), ci, v)
            end
            @eval function MOI.get(d::Optimizer, ::$attr, ci::CI{$F{T},$S{T}}) where T<:Real
                return MOI.get(d._model._input_model, $attr(), ci)
            end
        end
    end
    @eval function MOI.get(d::Optimizer, ::$attr, ci::CI{SV,S}) where  S <: Union{ZO, MOI.Integer}
        return MOI.get(d._model._input_model, $attr(), ci)
    end
    @eval function MOI.get(d::Optimizer, ::$attr, ci::CI{F,S}) where {F <: Union{VECVAR},
                                                                      S <: Union{MOI.SOS1, MOI.SOS2, SOC_CONE, EXP_CONE, POW_CONE, PSD_CONE}}
        return MOI.get(d._mode._input_modell, $attr(), ci)
    end
    @eval function MOI.set(d::Optimizer, ::$attr, ci::CI{SV,S}, v) where  S <: Union{ZO, MOI.Integer}
        return MOI.set(d._model._input_model, $attr(), ci, v)
    end
    @eval function MOI.set(d::Optimizer, ::$attr, ci::CI{F,S}) where {F <: Union{VECVAR},
                                                                      S <: Union{MOI.SOS1, MOI.SOS2, SOC_CONE, EXP_CONE, POW_CONE, PSD_CONE}}
        return MOI.set(d._model._input_model, $attr(), ci, v)
    end
end

function MOI.get(m::Optimizer{T}, v::MOI.VariablePrimal, vi::MOI.VariableIndex) where T <: Real
    MOI.check_result_index_bounds(m, v)
    m._solution[vi.value]
end
MOI.get(m::Optimizer{T}, p::MOI.VariablePrimal, vi::Vector{MOI.VariableIndex}) where T <: Real = MOI.get.(m, p, vi)
function MOI.get(m::Optimizer{T}, v::MOI.ConstraintPrimal, ci::MOI.ConstraintIndex{SV, S}) where {S <: Union{LT,GT,ET,IT,ZO,MOI.Integer}, T}
    MOI.check_result_index_bounds(m, v)
    return m._solution[ci.value]
end

function MOI.get(m::Optimizer{T}, v::MOI.ConstraintPrimal, ci::MOI.ConstraintIndex{F, S}) where {F <: Union{SAF, SQF}, S <: Union{LT,GT,ET,IT}, T}
    MOI.check_result_index_bounds(m, v)
    return m._constraint_primal[ci]
end
function MOI.get(m::Optimizer{T}, v::MOI.ConstraintPrimal, ci::MOI.ConstraintIndex{F, S}) where {F <: Union{VECVAR}, S <: Union{MOI.SOS1, MOI.SOS2, SOC_CONE, EXP_CONE, POW_CONE, PSD_CONE}, T}
    MOI.check_result_index_bounds(m, v)
    return m._constraint_primal[ci]
end
MOI.get(opt::Optimizer, a::MOI.ConstraintPrimal, ci::Vector{MOI.ConstraintIndex}) = MOI.get.(opt, a, ci)

const EAGO_OPTIMIZER_ATTRIBUTES = Symbol[:relaxed_optimizer, :lp_optimizer, :mip_optimizer,
                                         :socp_optimizer, :semidefinite_optimizer,
                                         :nlp_optimizer, :minlp_optimizer,
                                         :relaxed_optimizer_kwargs, :upper_optimizer,
                                         :enable_optimize_hook, :ext, :ext_type, :_parameters]
const EAGO_MODEL_STRUCT_ATTRIBUTES = Symbol[:_stack, :_log, :_current_node, :_working_problem, :_model]
const EAGO_MODEL_NOT_STRUCT_ATTRIBUTES = setdiff(fieldnames(Optimizer), union(EAGO_OPTIMIZER_ATTRIBUTES,
                                                                              EAGO_MODEL_STRUCT_ATTRIBUTES))
const EAGO_MODEL_EITHER_ATTRIBUTE = union(EAGO_MODEL_STRUCT_ATTRIBUTES, EAGO_MODEL_NOT_STRUCT_ATTRIBUTES)

function MOI.empty!(m::Optimizer{T}) where {N,T<:Real,S}
    m._model = InputModel{T}()
    m._solver = GlobalOptimizer{}(T)
    m._problem_type = UNCLASSIFIED

    m._termination_status_code = MOI.OPTIMIZE_NOT_CALLED
    m._primal_status_code      = MOI.OTHER_RESULT_STATUS

    m._node_count      = 0
    m._objective_value = typemax(T)
    m._objective_bound = typemin(T)
    m._relative_gap = typemax(T)

    m._start_time   = zero(T)
    m._time_limit   = typemax(T)
    m._time_left    = 1000*one(T)
    m._run_time     = zero(T)
    m._parse_time   = zero(T)
    return
end

function MOI.is_empty(m::Optimizer{T}) where T<:Real

    is_empty_flag = MOI.is_empty(m._model)
    is_empty_flag &= MOI.is_empty(m._solver)
    is_empty_flag &= (m._problem_type == UNCLASSIFIED)

    is_empty_flag &= (m._termination_status_code == MOI.OPTIMIZE_NOT_CALLED)
    is_empty_flag &= (m._primal_status_code      == MOI.OTHER_RESULT_STATUS)

    is_empty_flag &= (m._node_count      == 0)
    is_empty_flag &= (m._objective_value == typemax(T))
    is_empty_flag &= (m._objective_bound == typemin(T))
    is_empty_flag &= (m._relative_gap    == typemax(T))

    is_empty_flag &= (m._start_time == zero(T))
    is_empty_flag &= (m._time_limit == typemax(T))
    is_empty_flag &= (m._time_left  == 1000*one(T))
    is_empty_flag &= (m._run_time   == zero(T))
    is_empty_flag &= (m._parse_time == zero(T))

    return is_empty_flag
end

function MOI.get(m::Optimizer{T}, v::MOI.ObjectiveValue) where T<:Real
    MOI.check_result_index_bounds(m, v)
    m._objective_value
end
MOI.get(m::Optimizer{T}, ::MOI.ObjectiveBound) where T<:Real = m._objective_bound
MOI.get(m::Optimizer{T}, ::MOI.NumberOfVariables) where T<:Real = MOI.get(m._model._input_model, MOI.NumberOfVariables())

function MOI.get(m::Optimizer{T}, ::MOI.RelativeGap) where T<:Real
    LBD = m._objective_value
    UBD = m._objective_bound
    if m._model._optimization_sense === MOI.MAX_SENSE
        gap = abs(UBD - LBD)/abs(LBD)
    else
        gap = abs(UBD - LBD)/abs(UBD)
    end
    return gap
end

function MOI.get(m::Optimizer{T}, p::MOI.RawParameter) where T<:Real
    if p.name isa String
        psym = Symbol(p.name)
    elseif p.name isa Symbol
        psym = p.name
    else
        error("EAGO only supports raw parameters with Symbol or String names.")
    end

    if psym in fieldnames(EAGOParameters)
        return getfield(m._parameters, psym)
    end
    return getfield(m, psym)
end

MOI.add_variable(d::Optimizer) = MOI.add_variable(d._model._input_model)

for S in (LT,GT,ET,IT)
    for F in (SAF,SQF)
        @eval function MOI.add_constraint(d::Optimizer{T}, f::$F{T}, s::$S{T}) where T<:Real
            MOI.add_constraint(d._model._input_model, f, s)
        end
    end
end

for S in (LT,GT,ET,IT)
    @eval function MOI.add_constraint(d::Optimizer{T}, f::SV, s::$S{T}) where T<:Real
        MOI.add_constraint(d._model._input_model, f, s)
    end
end
function MOI.add_constraint(d::Optimizer{T}, f::SV, s::S) where {S<:Union{ZO,MOI.Integer},T<:Real}
    MOI.add_constraint(d._model._input_model, f, s)
end

function MOI.add_constraint(d::Optimizer, f::F, s::S) where {F<:Union{VECVAR},
                                                             S<:Union{MOI.SOS1, MOI.SOS2, SOC_CONE, EXP_CONE, POW_CONE, PSD_CONE}}
    MOI.add_constraint(d._model._input_model, f, s)
end

function MOI.set(d::Optimizer, ::MOI.ObjectiveFunction{SAF{T}}, func::SAF{T}) where T<:Real
    MOI.set(d._model._input_model, MOI.ObjectiveFunction{SAF{T}}(), func)
end

MOI.supports(d::Optimizer, ::MOI.NLPBlock) = true
function MOI.set(d::Optimizer, ::MOI.NLPBlock, nlp_data)
    d._model._nlp_data = nlp_data
    return
end
function MOI.set(d::Optimizer, ::MOI.NLPBlock, ::Nothing)
    d._model._nlp_data = nothing
    return
end

function MOI.set(d::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    MOI.set(d._model._input_model, MOI.ObjectiveSense(), sense)
end


text:
function reset_seed(seed)
    Random.seed!(seed)
end

function basetype_string(set::LazySet)
    return string(basetype(set))
end

function flatten_dyn(dyn)
    return dyn.A(), dyn.b()
end

function flatten_dyn(dyn::AffDyn)
    return dyn.A, dyn.b
end

function flatten_Interval(I)
    I = convert(Interval, I)
    return min(I), max(I)
end

function flatten_BallInf(B)
    B = convert(BallInf, B)
    return B.center, B.radius
end

function flatten_Hyperrectangle(H)
    H = convert(Hyperrectangle, H)
    return H.center, H.radius
end

function flatten_HalfSpace(hs)
    hs = convert(HalfSpace, hs)
    return hs.a, hs.b
end

function flatten_HPolyhedron(P)
    P = convert(HPolyhedron, P)
    return P.constraints
end

function flatten_VPolygon(P)
    P = convert(VPolygon, P)
    return P.vertices
end

function successor(b::AbstractVector, t, x0)
    result = x0 + b * t
    return result
end

function successor(A::AbstractMatrix, t, x0)
    result = exp(A * t) * x0
    return result
end

function successor(dyn::AffDyn, t, x0)
    if iszero(dyn.A)
        return successor(dyn.b, t, x0)
    elseif iszero(dyn.b)
        return successor(dyn.A, t, x0)
    end
    # affine solution that only works for invertible matrix
#    try
#        inv_A = inv(dyn.A)
#        id = Diagonal(ones(n))
#        result = exp(dyn.A * t) * x0 + (exp(dyn.A * t) - id) * inv_A * dyn.b
#        return result
#    catch e
#        # A not invertible
#        # return nothing
#    end
    # make affine system linear by extending the dimension
    A_high = [dyn.A dyn.b; zeros(length(dyn.b) + 1)']
    x1_high = successor(A_high, t, vcat(x0, 1.0))
    x1 = x1_high[1:end-1]
    return x1
end

function successor(b::AbstractVector, time, P0::LazySet)
    offset = b * time
    P1 = translate(P0, offset)
    return P1
end

function successor(A::AbstractMatrix, time, P0::LazySet)
    B = exp(A * time)
    P1 = linear_map(B, P0)
    return P1
end

function successor(dyn::AffDyn, time, P0::LazySet)
    # make affine system linear by extending the dimension
    n = length(dyn.b)
    A_high = [dyn.A dyn.b; zeros(n + 1)']
    if n == 1
        # 1D case can be solved more efficiently and with numerical stability
        P0_itv = convert(Interval, P0)
        l = low(P0_itv)[1]
        u = high(P0_itv)[1]
        P0_high = Hyperrectangle(low=[l, 1.0], high=[u, 1.0])
        P1_high = successor(A_high, time, P0_high)
        P1 = LazySets.project(P1_high, 1:n)
        P1 = convert(Interval, P1)
    else
        one = Interval(1.0 - 1e-5, 1.0 + 1e-5)
        # one = Singleton([1.0])
        P0_high = CartesianProduct(P0, one)
        P1_high = successor(A_high, time, P0_high)
        P1_high = LazySets.normalize(P1_high)
        P1 = LazySets.project(P1_high, 1:n)
        P1 = LazySets.normalize(P1)
    end
    return P1
end

function is_similar_flow(P0, dyn, P1, time; return_set::Bool=false)
    Q1 = successor(dyn, time, P0)
    if return_set
        return !isdisjoint(P1, Q1), Q1
    else
        return !isdisjoint(P1, Q1)
    end
end

function contains(X::LazySet, point::AbstractVector)
    return point ∈ X
end

function polytope_hull(X, Y)
    return convex_hull(X, Y)
end

function get_ivp(b::AbstractVector, X0)
    problem = @ivp(x' = b, x(0) ∈ X0)
    return problem
end

function get_ivp(A::AbstractMatrix, X0)
    problem = @ivp(x' = Ax, x(0) ∈ X0)
    return problem
end

function get_ivp(dyn::AffDyn, X0)
    A = dyn.A
    b = dyn.b
    problem = @ivp(x' = Ax + b, x(0) ∈ X0)
    return problem
end

function compute_reach_tube_hull(dyn, X0::LazySet, time_horizon,
                                 time_step::Number=1e-3, oa=Hyperrectangle)
    problem = get_ivp(dyn, X0)
    if USE_REACHABILITY
        sets = compute_reach_tube_hull_Reachability(problem, time_horizon,
            time_step)
    else
        n = dim(X0)
        sets = compute_reach_tube_hull_ReachabilityAnalysis(problem,
            time_horizon, time_step, n)
    end
    hull = overapproximate(ConvexHullArray(sets), oa)
    return hull
end

function compute_reach_tube_hull_Reachability(problem, time_horizon, time_step)
    options = Options(:T => time_horizon)
    algorithm = BFFPSV18(:δ => time_step)
    solution = Reachability.solve(problem, options; op=algorithm)
    flowpipe = solution.flowpipes[1]
    sets = [Reachability.set(rs) for rs in flowpipe.reachsets]
    return sets
end

function compute_reach_tube_hull_ReachabilityAnalysis(problem, time_horizon,
                                                      time_step, n)
    algorithm = BOX(δ=time_step)
    solution = ReachabilityAnalysis.solve(problem, T=time_horizon, algorithm)
    flowpipe = solution.F
    sets = [ReachabilityAnalysis.set(rs) for rs in ReachabilityAnalysis.array(flowpipe)]
    return sets
end

function getInterval(lo, hi)
    return Interval(lo, hi)
end

function getBallInf(center, radius)
    return BallInf(center, radius)
end

function getHyperrectangle(center, radius)
    return Hyperrectangle(center, radius)
end

function getHalfSpace(a, b)
    return HalfSpace(a, b)
end

function getHPolyhedron(constraints)
    return HPolyhedron(constraints)
end

function getUniverse(dimension)
    return Universe(dimension)
end

max_stars_repo_name: HySynth/HySynth

max_issues_repo_name: HySynth/HySynth

max_forks_repo_name: HySynth/HySynth


text:
using Ferrite: Cell

"""
    abstract type StiffnessTopOptProblem{dim, T} <: AbstractTopOptProblem end

An abstract stiffness topology optimization problem. All subtypes must have the following fields:
- `ch`: a `Ferrite.ConstraintHandler` struct
- `metadata`: Metadata having various cell-node-dof relationships
- `black`: a `BitVector` of length equal to the number of elements where `black[e]` is 1 iff the `e`^th element must be part of the final design
- `white`:  a `BitVector` of length equal to the number of elements where `white[e]` is 1 iff the `e`^th element must not be part of the final design
- `varind`: an `AbstractVector{Int}` of length equal to the number of elements where `varind[e]` gives the index of the decision variable corresponding to element `e`. Because some elements can be fixed to be black or white, not every element has a decision variable associated.
"""
abstract type StiffnessTopOptProblem{dim,T} <: AbstractTopOptProblem end

# Fallbacks
getdim(::StiffnessTopOptProblem{dim,T}) where {dim,T} = dim
floattype(::StiffnessTopOptProblem{dim,T}) where {dim,T} = T
getE(p::StiffnessTopOptProblem) = p.E
getν(p::StiffnessTopOptProblem) = p.ν
getgeomorder(p::StiffnessTopOptProblem) = nnodespercell(p) == 9 ? 2 : 1
getdensity(::StiffnessTopOptProblem{dim,T}) where {dim,T} = T(0)
getmetadata(p::StiffnessTopOptProblem) = p.metadata
getdh(p::StiffnessTopOptProblem) = p.ch.dh
getcloaddict(p::StiffnessTopOptProblem{dim,T}) where {dim,T} = Dict{String,Vector{T}}()
getpressuredict(p::StiffnessTopOptProblem{dim,T}) where {dim,T} = Dict{String,T}()
getfacesets(p::StiffnessTopOptProblem{dim,T}) where {dim,T} = Dict{String,Tuple{Int,T}}()
Ferrite.getncells(problem::StiffnessTopOptProblem) = Ferrite.getncells(getdh(problem).grid)

"""
```
///**********************************
///*                                *
///*                                * |
///*                                * |
///********************************** v

@params struct PointLoadCantilever{dim, T, N, M} <: StiffnessTopOptProblem{dim, T}
    rect_grid::RectilinearGrid{dim, T, N, M}
    E::T
    ν::T
    ch::ConstraintHandler{<:DofHandler{dim, <:Cell{dim,N,M}, T}, T}
    force::T
    force_dof::Integer
    black::AbstractVector
    white::AbstractVector
    varind::AbstractVector{Int}
    metadata::Metadata
end
```

- `dim`: dimension of the problem
- `T`: number type for computations and coordinates
- `N`: number of nodes in a cell of the grid
- `M`: number of faces in a cell of the grid
- `rect_grid`: a RectilinearGrid struct
- `E`: Young's modulus
- `ν`: Poisson's ration
- `force`: force at the center right of the cantilever beam (positive is downward)
- `force_dof`: dof number at which the force is applied
- `ch`: a `Ferrite.ConstraintHandler` struct
- `metadata`: Metadata having various cell-node-dof relationships
- `black`: a `BitVector` of length equal to the number of elements where `black[e]` is 1 iff the `e`^th element must be part of the final design
- `white`:  a `BitVector` of length equal to the number of elements where `white[e]` is 1 iff the `e`^th element must not be part of the final design
- `varind`: an `AbstractVector{Int}` of length equal to the number of elements where `varind[e]` gives the index of the decision variable corresponding to element `e`. Because some elements can be fixed to be black or white, not every element has a decision variable associated.
"""
@params struct PointLoadCantilever{dim,T,N,M} <: StiffnessTopOptProblem{dim,T}
    rect_grid::RectilinearGrid{dim,T,N,M}
    E::T
    ν::T
    ch::ConstraintHandler{<:DofHandler{dim,<:Cell{dim,N,M},T},T}
    force::T
    force_dof::Integer
    black::AbstractVector
    white::AbstractVector
    varind::AbstractVector{Int}
    metadata::Metadata
end
function Base.show(::IO, ::MIME{Symbol("text/plain")}, ::PointLoadCantilever)
    return println("TopOpt point load cantilever beam problem")
end

"""
    PointLoadCantilever(::Type{Val{CellType}}, nels::NTuple{dim,Int}, sizes::NTuple{dim}, E, ν, force) where {dim, CellType}

- `dim`: dimension of the problem
- `E`: Young's modulus
- `ν`: Poisson's ration
- `force`: force at the center right of the cantilever beam (positive is downward)
- `nels`: number of elements in each direction, a 2-tuple for 2D problems and a 3-tuple for 3D problems
- `sizes`: the size of each element in each direction, a 2-tuple for 2D problems and a 3-tuple for 3D problems
- `CellType`: can be either `:Linear` or `:Quadratic` to determine the order of the geometric and field basis functions and element type. Only isoparametric elements are supported for now.

Example:
```
nels = (60,20);
sizes = (1.0,1.0);
E = 1.0;
ν = 0.3;
force = 1.0;

# Linear elements and linear basis functions
celltype = :Linear

# Quadratic elements and quadratic basis functions
#celltype = :Quadratic

problem = PointLoadCantilever(Val{celltype}, nels, sizes, E, ν, force)
```
"""
function PointLoadCantilever(
    ::Type{Val{CellType}},
    nels::NTuple{dim,Int},
    sizes::NTuple{dim},
    E=1.0,
    ν=0.3,
    force=1.0,
) where {dim,CellType}
    iseven(nels[2]) && (length(nels) < 3 || iseven(nels[3])) ||
        throw("Grid does not have an even number of elements along the y and/or z axes.")

    _T = promote_type(eltype(sizes), typeof(E), typeof(ν), typeof(force))
    if _T <: Integer
        T = Float64
    else
        T = _T
    end
    if CellType === :Linear || dim === 3
        rect_grid = RectilinearGrid(Val{:Linear}, nels, T.(sizes))
    else
        rect_grid = RectilinearGrid(Val{:Quadratic}, nels, T.(sizes))
    end

    if haskey(rect_grid.grid.facesets, "fixed_all")
        pop!(rect_grid.grid.facesets, "fixed_all")
    end
    #addfaceset!(rect_grid.grid, "fixed_all", x -> left(rect_grid, x));
    addnodeset!(rect_grid.grid, "fixed_all", x -> left(rect_grid, x))

    if haskey(rect_grid.grid.nodesets, "down_force")
        pop!(rect_grid.grid.nodesets, "down_force")
    end
    addnodeset!(
        rect_grid.grid, "down_force", x -> right(rect_grid, x) && middley(rect_grid, x)
    )

    # Create displacement field u
    dh = DofHandler(rect_grid.grid)
    if CellType === :Linear || dim === 3
        push!(dh, :u, dim) # Add a displacement field
    else
        ip = Lagrange{2,RefCube,2}()
        push!(dh, :u, dim, ip) # Add a displacement field
    end
    close!(dh)

    ch = ConstraintHandler(dh)

    #dbc = Dirichlet(:u, getfaceset(rect_grid.grid, "fixed_all"), (x,t) -> zeros(T, dim), collect(1:dim))
    dbc = Dirichlet(
        :u, getnodeset(rect_grid.grid, "fixed_all"), (x, t) -> zeros(T, dim), collect(1:dim)
    )
    add!(ch, dbc)
    close!(ch)
    t = T(0)
    update!(ch, t)

    metadata = Metadata(dh)

    fnode = Tuple(getnodeset(rect_grid.grid, "down_force"))[1]
    node_dofs = metadata.node_dofs
    force_dof = node_dofs[2, fnode]

    N = nnodespercell(rect_grid)
    M = nfacespercell(rect_grid)

    black, white = find_black_and_white(dh)
    varind = find_varind(black, white)

    return PointLoadCantilever(
        rect_grid, E, ν, ch, force, force_dof, black, white, varind, metadata
    )
end

"""
```
 |
 |
 v
O*********************************
O*                               *
O*                               *
O*                               *
O*********************************
                                 O

struct HalfMBB{dim, T, N, M} <: StiffnessTopOptProblem{dim, T}
    rect_grid::RectilinearGrid{dim, T, N, M}
    E::T
    ν::T
    ch::ConstraintHandler{<:DofHandler{dim, <:Cell{dim,N,M}, T}, T}
    force::T
    force_dof::Integer
    black::AbstractVector
    white::AbstractVector
    varind::AbstractVector{Int}
    metadata::Metadata
end
```

- `dim`: dimension of the problem
- `T`: number type for computations and coordinates
- `N`: number of nodes in a cell of the grid
- `M`: number of faces in a cell of the grid
- `rect_grid`: a RectilinearGrid struct
- `E`: Young's modulus
- `ν`: Poisson's ration
- `force`: force at the top left of half the MBB (positive is downward)
- `force_dof`: dof number at which the force is applied
- `ch`: a `Ferrite.ConstraintHandler` struct
- `metadata`: Metadata having various cell-node-dof relationships
- `black`: a `BitVector` of length equal to the number of elements where `black[e]` is 1 iff the `e`^th element must be part of the final design
- `white`:  a `BitVector` of length equal to the number of elements where `white[e]` is 1 iff the `e`^th element must not be part of the final design
- `varind`: an `AbstractVector{Int}` of length equal to the number of elements where `varind[e]` gives the index of the decision variable corresponding to element `e`. Because some elements can be fixed to be black or white, not every element has a decision variable associated.
"""
@params struct HalfMBB{dim,T,N,M} <: StiffnessTopOptProblem{dim,T}
    rect_grid::RectilinearGrid{dim,T,N,M}
    E::T
    ν::T
    ch::ConstraintHandler{<:DofHandler{dim,<:Cell{dim,N,M},T},T}
    force::T
    force_dof::Integer
    black::AbstractVector
    white::AbstractVector
    varind::AbstractVector{Int}
    metadata::Metadata
end
function Base.show(::IO, ::MIME{Symbol("text/plain")}, ::HalfMBB)
    return println("TopOpt half MBB problem")
end

"""
    HalfMBB(::Type{Val{CellType}}, nels::NTuple{dim,Int}, sizes::NTuple{dim}, E, ν, force) where {dim, CellType}

- `dim`: dimension of the problem
- `E`: Young's modulus
- `ν`: Poisson's ration
- `force`: force at the top left of half the MBB (positive is downward)
- `nels`: number of elements in each direction, a 2-tuple for 2D problems and a 3-tuple for 3D problems
- `sizes`: the size of each element in each direction, a 2-tuple for 2D problems and a 3-tuple for 3D problems
- `CellType`: can be either `:Linear` or `:Quadratic` to determine the order of the geometric and field basis functions and element type. Only isoparametric elements are supported for now.

Example:
```
nels = (60,20);
sizes = (1.0,1.0);
E = 1.0;
ν = 0.3;
force = -1.0;

# Linear elements and linear basis functions
celltype = :Linear

# Quadratic elements and quadratic basis functions
#celltype = :Quadratic

problem = HalfMBB(Val{celltype}, nels, sizes, E, ν, force)
```
"""
function HalfMBB(
    ::Type{Val{CellType}},
    nels::NTuple{dim,Int},
    sizes::NTuple{dim},
    E=1.0,
    ν=0.3,
    force=1.0,
) where {dim,CellType}
    _T = promote_type(eltype(sizes), typeof(E), typeof(ν), typeof(force))
    if _T <: Integer
        T = Float64
    else
        T = _T
    end
    if CellType === :Linear || dim === 3
        rect_grid = RectilinearGrid(Val{:Linear}, nels, T.(sizes))
    else
        rect_grid = RectilinearGrid(Val{:Quadratic}, nels, T.(sizes))
    end

    if haskey(rect_grid.grid.facesets, "fixed_u1")
        pop!(rect_grid.grid.facesets, "fixed_u1")
    end
    #addfaceset!(rect_grid.grid, "fixed_u1", x -> left(rect_grid, x));
    addnodeset!(rect_grid.grid, "fixed_u1", x -> left(rect_grid, x))

    if haskey(rect_grid.grid.nodesets, "fixed_u2")
        pop!(rect_grid.grid.nodesets, "fixed_u2")
    end
    addnodeset!(
        rect_grid.grid, "fixed_u2", x -> bottom(rect_grid, x) && right(rect_grid, x)
    )

    if haskey(rect_grid.grid.nodesets, "down_force")
        pop!(rect_grid.grid.nodesets, "down_force")
    end
    addnodeset!(rect_grid.grid, "down_force", x -> top(rect_grid, x) && left(rect_grid, x))

    # Create displacement field u
    dh = DofHandler(rect_grid.grid)
    if CellType === :Linear || dim === 3
        push!(dh, :u, dim)
    else
        ip = Lagrange{2,RefCube,2}()
        push!(dh, :u, dim, ip)
    end
    close!(dh)

    ch = ConstraintHandler(dh)
    #dbc1 = Dirichlet(:u, getfaceset(rect_grid.grid, "fixed_u1"), (x,t)->T[0], [1])
    dbc1 = Dirichlet(:u, getnodeset(rect_grid.grid, "fixed_u1"), (x, t) -> T[0], [1])
    add!(ch, dbc1)
    dbc2 = Dirichlet(:u, getnodeset(rect_grid.grid, "fixed_u2"), (x, t) -> T[0], [2])
    add!(ch, dbc2)
    close!(ch)

    t = T(0)
    update!(ch, t)

    metadata = Metadata(dh)

    fnode = Tuple(getnodeset(rect_grid.grid, "down_force"))[1]
    node_dofs = metadata.node_dofs
    force_dof = node_dofs[2, fnode]

    N = nnodespercell(rect_grid)
    M = nfacespercell(rect_grid)

    black, white = find_black_and_white(dh)
    varind = find_varind(black, white)

    return HalfMBB(rect_grid, E, ν, ch, force, force_dof, black, white, varind, metadata)
end

nnodespercell(p::Union{PointLoadCantilever,HalfMBB}) = nnodespercell(p.rect_grid)
function getcloaddict(p::Union{PointLoadCantilever{dim,T},HalfMBB{dim,T}}) where {dim,T}
    f = T[0, -p.force, 0]
    fnode = Tuple(getnodeset(p.rect_grid.grid, "down_force"))[1]
    return Dict{Int,Vector{T}}(fnode => f)
end

"""
```
////////////
............
.          .
.          .
.          .
.          .
.          ......................
.                               .
.                               .
.                               . |
................................. v
                                force

struct LBeam{T, N, M} <: StiffnessTopOptProblem{2, T}
    E::T
    ν::T
    ch::ConstraintHandler{<:DofHandler{2, <:Cell{2,N,M}, T}, T}
    force::T
    force_dof::Integer
    black::AbstractVector
    white::AbstractVector
    varind::AbstractVector{Int}
    metadata::Metadata
end
```

- `T`: number type for computations and coordinates
- `N`: number of nodes in a cell of the grid
- `M`: number of faces in a cell of the grid
- `E`: Young's modulus
- `ν`: Poisson's ration
- `force`: force at the center right of the cantilever beam (positive is downward)
- `force_dof`: dof number at which the force is applied
- `ch`: a `Ferrite.ConstraintHandler` struct
- `metadata`: Metadata having various cell-node-dof relationships
- `black`: a `BitVector` of length equal to the number of elements where `black[e]` is 1 iff the `e`^th element must be part of the final design
- `white`:  a `BitVector` of length equal to the number of elements where `white[e]` is 1 iff the `e`^th element must not be part of the final design
- `varind`: an `AbstractVector{Int}` of length equal to the number of elements where `varind[e]` gives the index of the decision variable corresponding to element `e`. Because some elements can be fixed to be black or white, not every element has a decision variable associated.
"""
@params struct LBeam{T,N,M} <: StiffnessTopOptProblem{2,T}
    E::T
    ν::T
    ch::ConstraintHandler{<:DofHandler{2,<:Cell{2,N,M},T},T}
    force::T
    force_dof::Integer
    black::AbstractVector
    white::AbstractVector
    varind::AbstractVector{Int}
    metadata::Metadata
end
Base.show(::IO, ::MIME{Symbol("text/plain")}, ::LBeam) = println("TopOpt L-beam problem")

"""
    LBeam(::Type{Val{CellType}}, ::Type{T}=Float64; length = 100, height = 100, upperslab = 50, lowerslab = 50, E = 1.0, ν = 0.3, force = 1.0) where {T, CellType}

- `T`: number type for computations and coordinates
- `E`: Young's modulus
- `ν`: Poisson's ration
- `force`: force at the center right of the cantilever beam (positive is downward)
- `length`, `height`, `upperslab` and `lowerslab` are explained in [`LGrid`](@ref).
- `CellType`: can be either `:Linear` or `:Quadratic` to determine the order of the geometric and field basis functions and element type. Only isoparametric elements are supported for now.

Example:
```
E = 1.0;
ν = 0.3;
force = 1.0;

# Linear elements and linear basis functions
celltype = :Linear

# Quadratic elements and quadratic basis functions
#celltype = :Quadratic

problem = LBeam(Val{celltype}, E = E, ν = ν, force = force)
```
"""
function LBeam(
    ::Type{Val{CellType}},
    ::Type{T}=Float64;
    length=100,
    height=100,
    upperslab=50,
    lowerslab=50,
    E=1.0,
    ν=0.3,
    force=1.0,
) where {T,CellType}
    # Create displacement field u
    grid = LGrid(
        Val{CellType},
        T;
        length=length,
        height=height,
        upperslab=upperslab,
        lowerslab=lowerslab,
    )

    dh = DofHandler(grid)
    if CellType === :Linear
        push!(dh, :u, 2)
    else
        ip = Lagrange{2,RefCube,2}()
        push!(dh, :u, 2, ip)
    end
    close!(dh)

    ch = ConstraintHandler(dh)
    dbc = Dirichlet(:u, getfaceset(grid, "top"), (x, t) -> T[0, 0], [1, 2])
    add!(ch, dbc)
    close!(ch)

    t = T(0)
    update!(ch, t)

    metadata = Metadata(dh)

    fnode = Tuple(getnodeset(grid, "load"))[1]
    node_dofs = metadata.node_dofs
    force_dof = node_dofs[2, fnode]

    black, white = find_black_and_white(dh)
    varind = find_varind(black, white)

    TInds = typeof(varind)
    TMeta = typeof(metadata)
    return LBeam(E, ν, ch, force, force_dof, black, white, varind, metadata)
end

function boundingbox(grid::Ferrite.Grid{dim}) where {dim}
    xmin1 = minimum(n -> n.x[1], grid.nodes)
    xmax1 = maximum(n -> n.x[1], grid.nodes)
    xmin2 = minimum(n -> n.x[2], grid.nodes)
    xmax2 = maximum(n -> n.x[2], grid.nodes)
    if dim === 2
        return ((xmin1, xmin2), (xmax1, xmax2))
    else
        xmin3 = minimum(n -> n.x[3], grid.nodes)
        xmax3 = maximum(n -> n.x[3], grid.nodes)
        return ((xmin1, xmin2, xmin3), (xmax1, xmax2, xmax3))
    end
end

function RectilinearTopology(b, topology=ones(getncells(getdh(b).grid)))
    bb = boundingbox(getdh(b).grid)
    go = getgeomorder(b)
    nels = Int.(round.(bb[2] .- bb[1]))
    dim = length(nels)
    if go === 1
        rectgrid = generate_grid(Quadrilateral, nels, Vec{dim}(bb[1]), Vec{dim}(bb[2]))
    elseif go === 2
        rectgrid = generate_grid(
            QuadraticQuadrilateral, nels, Vec{dim}(bb[1]), Vec{dim}(bb[2])
        )
    else
        throw("Unsupported geometry.")
    end
    new_topology = zeros(prod(nels))
    for (i, cell) in enumerate(CellIterator(getdh(b)))
        sub = Int.(round.((cell.coords[1]...,))) .+ (1, 1)
        ind = LinearIndices(nels)[sub...]
        new_topology[ind] = topology[i]
    end
    return copy(reshape(new_topology, nels)')
end

nnodespercell(p::LBeam{T,N}) where {T,N} = N
getdim(::LBeam) = 2
function getcloaddict(p::LBeam{T}) where {T}
    f = T[0, -p.force]
    fnode = Tuple(getnodeset(getdh(p).grid, "load"))[1]
    return Dict{Int,Vector{T}}(fnode => f)
end

"""
```
                                                               1

                                                              OOO
                                                              ...
                                                              . .
                                                           4  . .
                                30                            . .
/ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . <-
/ .                                                                 . <- 2 f
/ .    3                                                            . <-
/ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . <-
                                                              ^^^
                                                              |||
                                                              1 f

struct TieBeam{T, N, M} <: StiffnessTopOptProblem{2, T}
    E::T
    ν::T
    force::T
    ch::ConstraintHandler{<:DofHandler{2, N, T, M}, T}
    black::AbstractVector
    white::AbstractVector
    varind::AbstractVector{Int}
    metadata::Metadata
end
```

- `T`: number type for computations and coordinates
- `N`: number of nodes in a cell of the grid
- `M`: number of faces in a cell of the grid
- `E`: Young's modulus
- `ν`: Poisson's ration
- `force`: force at the center right of the cantilever beam (positive is downward)
- `ch`: a `Ferrite.ConstraintHandler` struct
- `metadata`: Metadata having various cell-node-dof relationships
- `black`: a `BitVector` of length equal to the number of elements where `black[e]` is 1 iff the `e`^th element must be part of the final design
- `white`:  a `BitVector` of length equal to the number of elements where `white[e]` is 1 iff the `e`^th element must not be part of the final design
- `varind`: an `AbstractVector{Int}` of length equal to the number of elements where `varind[e]` gives the index of the decision variable corresponding to element `e`. Because some elements can be fixed to be black or white, not every element has a decision variable associated.
"""
@params struct TieBeam{T,N,M} <: StiffnessTopOptProblem{2,T}
    E::T
    ν::T
    force::T
    ch::ConstraintHandler{<:DofHandler{2,<:Cell{2,N,M},T},T}
    black::AbstractVector
    white::AbstractVector
    varind::AbstractVector{Int}
    metadata::Metadata
end
function Base.show(::IO, ::MIME{Symbol("text/plain")}, ::TieBeam)
    return println("TopOpt tie-beam problem")
end

"""
    TieBeam(::Type{Val{CellType}}, ::Type{T} = Float64, refine = 1, force = T(1); E = T(1), ν = T(0.3)) where {T, CellType}

- `T`: number type for computations and coordinates
- `E`: Young's modulus
- `ν`: Poisson's ration
- `force`: force at the center right of the cantilever beam (positive is downward)
- `refine`: an integer value of 1 or greater that specifies the mesh refinement extent. A value of 1 gives the standard tie-beam problem in literature.
- `CellType`: can be either `:Linear` or `:Quadratic` to determine the order of the geometric and field basis functions and element type. Only isoparametric elements are supported for now.
"""
function TieBeam(
    ::Type{Val{CellType}}, ::Type{T}=Float64, refine=1, force=T(1); E=T(1), ν=T(0.3)
) where {T,CellType}
    grid = TieBeamGrid(Val{CellType}, T, refine)
    dh = DofHandler(grid)
    if CellType === :Linear
        push!(dh, :u, 2)
    else
        ip = Lagrange{2,RefCube,2}()
        push!(dh, :u, 2, ip)
    end
    close!(dh)

    ch = ConstraintHandler(dh)
    dbc1 = Dirichlet(:u, getfaceset(grid, "leftfixed"), (x, t) -> T[0, 0], [1, 2])
    add!(ch, dbc1)
    dbc2 = Dirichlet(:u, getfaceset(grid, "toproller"), (x, t) -> T[0], [2])
    add!(ch, dbc2)
    close!(ch)

    t = T(0)
    update!(ch, t)

    metadata = Metadata(dh)

    black, white = find_black_and_white(dh)
    varind = find_varind(black, white)

    return TieBeam(E, ν, force, ch, black, white, varind, metadata)
end

getdim(::TieBeam) = 2
nnodespercell(::TieBeam{T,N}) where {T,N} = N
function getpressuredict(p::TieBeam{T}) where {T}
    return Dict{String,T}("rightload" => 2 * p.force, "bottomload" => -p.force)
end
getfacesets(p::TieBeam) = getdh(p).grid.facesets

"""
```
    ******************************
    * Pin1         F1 ―――>       *
    *  o            |            *
    *               v            *
    *  Pin2                      *
    *    o     F2 ―――>           *
    *           |                *
    *           v                *
    ******************************

    RayProblem(nels, pins, loads)

Constructs an instance of the type `RayProblem` that is a 2D beam with:
 - Number of elements `nels`, e.g. `(60, 20)` where each element is a 1 x 1 square,
 - Pinned locations `pins` where each pinned location is a `Vector` of length 2, e.g. `[[1, 18], [2, 8]]` indicating the locations of the pins, and
 - Loads specified in `loads` where `loads` is a dictionary mapping the location of each load to its vector value, e.g. `Dict([10, 18] => [1.0, -1.0], [5, 5] => [1.0, -1.0])` which defines a load of `[1.0, -1.0]` at the point located at `[10, 18]` and a similar load at the point located at `[5, 5]`.
```
"""
@params struct RayProblem{T,N,M} <: StiffnessTopOptProblem{2,T}
    rect_grid::RectilinearGrid{2,T,N,M}
    E::T
    ν::T
    ch::ConstraintHandler{<:DofHandler{2,<:Cell{2,N,M},T},T}
    loads::Dict
    black::AbstractVector
    white::AbstractVector
    varind::AbstractVector{Int}
    metadata::Metadata
end
function RayProblem(
    nels::NTuple{2,Int}, pins::Vector{<:Vector}, loads::Dict{<:Vector,<:Vector}
)
    T = Float64
    rect_grid = RectilinearGrid(Val{:Linear}, nels, (1.0, 1.0))
    dim = length(nels)

    for (i, pin) in enumerate(pins)
        if haskey(rect_grid.grid.nodesets, "fixed$i")
            pop!(rect_grid.grid.nodesets, "fixed$i")
        end
        addnodeset!(rect_grid.grid, "fixed$i", x -> x ≈ pin)
    end
    for (i, k) in enumerate(keys(loads))
        if haskey(rect_grid.grid.nodesets, "force$i")
            pop!(rect_grid.grid.nodesets, "force$i")
        end
        addnodeset!(rect_grid.grid, "force$i", x -> x ≈ k)
    end

    # Create displacement field u
    dh = DofHandler(rect_grid.grid)
    ip = Lagrange{2,RefCube,1}()
    push!(dh, :u, dim, ip) # Add a displacement field
    close!(dh)

    ch = ConstraintHandler(dh)

    for i in 1:length(pins)
        dbc = Dirichlet(
            :u,
            getnodeset(rect_grid.grid, "fixed$i"),
            (x, t) -> zeros(T, dim),
            collect(1:dim),
        )
        add!(ch, dbc)
    end
    close!(ch)
    t = T(0)
    update!(ch, t)

    metadata = Metadata(dh)
    black, white = find_black_and_white(dh)
    varind = find_varind(black, white)

    loadsdict = Dict{Int,Vector{Float64}}(
        map(enumerate(keys(loads))) do (i, k)
            fnode = Tuple(getnodeset(rect_grid.grid, "force$i"))[1]
            (fnode => loads[k])
        end,
    )

    return RayProblem(rect_grid, 1.0, 0.3, ch, loadsdict, black, white, varind, metadata)
end
nnodespercell(p::RayProblem) = nnodespercell(p.rect_grid)
getcloaddict(p::RayProblem) = p.loads