text:
# Rational Quadratic ARD Covariance Function

"""
    RQArd <: StationaryARD{WeightedSqEuclidean}

ARD Rational Quadratic kernel (covariance)
```math
k(x,x') = σ²(1 + (x - x')ᵀL⁻²(x - x')/(2α))^{-α}
```
with length scale ``ℓ = (ℓ₁, ℓ₂, …)``, signal standard deviation ``σ``, and shape parameter
``α`` where ``L = diag(ℓ₁, ℓ₂, …)``.
"""
mutable struct RQArd <: StationaryARD{WeightedSqEuclidean}
    "Inverse squared length scale"
    iℓ2::Vector{Float64}
    "Signal variance"
    σ2::Float64
    "Shape parameter"
    α::Float64
    "Priors for kernel parameters"
    priors::Array

    """
        RQArd(ll::Vector{Float64}, lσ::Float64, lα::Float64)

    Create `RQArd` with length scale `exp.(ll)`, signal standard deviation `exp(lσ)`, and
    shape parameter `exp(lα)`.
    """
    RQArd(ll::Vector{Float64}, lσ::Float64, lα::Float64) =
        new(exp.(-2 .* ll), exp(2 * lσ), exp(lα), [])
end

function set_params!(rq::RQArd, hyp::VecF64)
    length(hyp) == num_params(rq) || throw(ArgumentError("RQArd kernel has $(num_params(rq_ard)) parameters"))
    @views @. rq.iℓ2 = exp(-2 * hyp[1:(end-2)])
    rq.σ2 = exp(2 * hyp[end-1])
    rq.α = exp(hyp[end])
end

get_params(rq::RQArd) = [-log.(rq.iℓ2) / 2; log(rq.σ2) / 2; log(rq.α)]
get_param_names(rq::RQArd) = [get_param_names(rq.iℓ2, :ll); :lσ; :lα]
num_params(rq::RQArd) = length(rq.iℓ2) + 2

Statistics.cov(rq::RQArd,r::Float64) = rq.σ2*(1+0.5*r/rq.α)^(-rq.α)

@inline dk_dll(rq::RQArd, r::Float64, wdiffp::Float64) =
    rq.σ2 * wdiffp * (1 + r / (2 * rq.α))^(-rq.α - 1)
@inline function dk_dlα(rq::RQArd, r::Float64)
    part = (1 + r / (2 * rq.α))
    return rq.σ2 * part^(-rq.α) * (r / (2 * part) - rq.α * log(part))
end
@inline function dKij_dθp(rq::RQArd, X::MatF64, i::Int, j::Int, p::Int, dim::Int)
    if p <= dim
        return dk_dll(rq, distij(metric(rq),X,i,j,dim), distijk(metric(rq),X,i,j,p))
    elseif p==dim+1
        return dk_dlσ(rq, distij(metric(rq),X,i,j,dim))
    else
        return dk_dlα(rq, distij(metric(rq),X,i,j,dim))
    end
end
@inline function dKij_dθp(rq::RQArd, X::MatF64, data::StationaryARDData, i::Int, j::Int, p::Int, dim::Int)
    return dKij_dθp(rq,X,i,j,p,dim)
end


text:
@testset "normal wrench" begin
    for k_quad_rule = (1, 2)
        p_pos = (0.1, 0.2)

        ### Box properties
        box_rad = 0.05
        i_prop_rigid     = InertiaProperties(400.0, d=0.09)
        c_prop_compliant = ContactProperties(Ē=1.0e9)  # , d=box_rad)

        ### Create mechanism and temporary structure
        mech_scen = MechanismScenario() # n_quad_rule=k_quad_rule)

        name_plane = "plane"
        eM_plane = eMesh_half_plane()
        add_contact!(mech_scen, name_plane, as_tet_eMesh(eM_plane), c_prop=c_prop_compliant)

        name_box = "box"
        eM_box = as_tri_eMesh(eMesh_box(box_rad))
        transform!(eM_box, SVector{3,Float64}(0.0, 0.0, box_rad))
        body_box, joint_box = add_body_contact!(mech_scen, name_box, eM_box, i_prop=i_prop_rigid)

        m_id_plane = find_mesh_id(mech_scen, name_plane)
        m_id_box = find_mesh_id(mech_scen, name_box)
        add_friction_bristle!(mech_scen, m_id_box, m_id_plane, μd=0.3, χ=0.6, k̄=1.0e6, τ=0.03, n_quad_rule=k_quad_rule)
        finalize!(mech_scen)

        pene = 0.1 * box_rad
        set_state_spq!(mech_scen, joint_box, trans=SVector{3,Float64}(p_pos..., -pene))

        x = get_state(mech_scen)
        calcXd(x, mech_scen)

        b = mech_scen.float.bodyBodyCache
        wrench = normal_wrench(b)
        wrench = as_static_vector(wrench)
        check = b.Ē * pene / 1.0 * box_rad^2*4
        f3 = SVector(0.0, 0.0, check)
        a3 = cross(SVector(p_pos..., 0.0), f3)
        check = -vcat(a3, f3)
        @test wrench ≈ check

        # wrench, p_center = normal_wrench_patch_center(b)
        # wrench = as_static_vector(wrench)
        # check = b.Ē * pene / 1.0 * box_rad^2*4
        # @test check ≈ wrench[6]
        # @test all((p_pos .* 0.999) .< p_center.v[1:2] .< (p_pos .* 1.001))
    end
end

ext: jl


text:
module CURAND

using ..CuArrays
using ..CuArrays: libcurand, @retry_reclaim

using CUDAapi

using CUDAdrv
using CUDAdrv: CUstream

using CUDAnative

using CEnum

# core library
include("libcurand_common.jl")
include("error.jl")
include("libcurand.jl")

# low-level wrappers
include("wrappers.jl")

# high-level integrations
include("random.jl")

# thread cache for task-local library handles
const thread_generators = Vector{Union{Nothing,RNG}}()

function generator()
    tid = Threads.threadid()
    if @inbounds thread_generators[tid] === nothing
        ctx = context()
        thread_generators[tid] = get!(task_local_storage(), (:CURAND, ctx)) do
            rng = RNG()
            Random.seed!(rng)
            rng
        end
    end
    @inbounds thread_generators[tid]
end

function __init__()
    resize!(thread_generators, Threads.nthreads())
    fill!(thread_generators, nothing)

    CUDAnative.atcontextswitch() do tid, ctx
        thread_generators[tid] = nothing
    end

    CUDAnative.attaskswitch() do tid, task
        thread_generators[tid] = nothing
    end
end

end

const seed! = CURAND.seed!
const rand = CURAND.rand
const randn = CURAND.randn
const rand_logn = CURAND.rand_logn
const rand_poisson = CURAND.rand_poisson

@deprecate curand CuArrays.rand
@deprecate curandn CuArrays.randn


text:
using Markdown
using Test

using Hyperscript

using HTM

const r = Hyperscript.render

@testset "HTM.jl" begin
    # Warning: some test cases may not represent supported usage.
    @testset "Features" begin
        @testset "Spread attributes" for attrs in ([:class => "fruit"], Dict(:class => "fruit"),
                                                   ["class" => "fruit"], Dict("class" => "fruit"))
            @test htm"

" |> r == "

"
        end

        @testset "Styles" begin
            @testset "Interpolations" for style in ([:background => "orange"], Dict(:background => "orange"),
                                                    ["background" => "orange"], Dict("background" => "orange"))
                @test htm"pineapple" |> r == "pineapple"
            end
            @test htm"pineapple" |> r == "pineapple"
            @test_broken htm"pineapple" |> r == "pineapple"
        end

        @testset "Classes" begin
            @testset "Interpolations" for class in ("fruit sour ", ["fruit", "sour", "sour"], Set(["fruit", "sour", "sour"]))
                # TODO: remove space at the end
                @test htm"

" |> r == "

"
            end
            @test htm"

" |> r == "

"
        end

        @testset "Self-closing tags" begin
            @test htm"
" |> r == "

"
            @test htm"<circle />" |> r == "<circle />"
        end

        @testset "Multiple root elements (fragments)" begin
            @test htm"

" == [htm"
", htm"
"]
        end

        @testset "Short circuit rendering" begin
            @test htm"
$(true || \"🍍\")
" |> r == "

"
            @test htm"
$(false || \"🍍\")
" |> r == "
🍍
"
            @test htm"$(true || \"🍍\")" === true
            @test htm"$(false || \"🍍\")" == "🍍"
        end

        @testset "Boolean attributes" begin
            @test htm"
" |> r == "

"
            @test htm"
" |> r == "

"
            @test htm"
" |> r == "

"
        end

        @testset "HTML optional quotes" begin
            @test htm"

" |> r == "

"

            @testset "URLs" begin
                @test htm""
                @test htm"" |> r == ""
            end
        end

        @testset "Components" begin
            struct Fruit
                name::String
                emoji::Char
            end
            @test htm"$(Fruit(\"pineapple\", '🍍'))" == Fruit("pineapple", '🍍')
            Base.show(io::IO, m::MIME"text/html", 🍎::Fruit) = show(io, m, htm"$(🍎.name):
$(🍎.emoji)
")
            Base.show(io::IO, ::MIME"text/plain", 🍎::Fruit) = show(io, MIME("text/html"), 🍎)
            @test_broken htm"
$(Fruit(\"pineapple\", '🍍'))
" |> r == "
pineapple:
🍍

"
        end

        @testset "Generic end-tags" begin
            @test htm"
🍍<//>" |> r == "
🍍
"
        end

        @testset "HTML-style comments" begin
            @test htm"" === nothing
            @test htm"

" |> r == "

"
            @test htm"

" |> r == "

"
            @test htm"

" |> r == "

"
        end

        @testset "Callbacks" begin
            @test_broken htm"<button onclick=$(() -> pineapples += 1)>🍍</button>" |> r == "<button onclick=\"pineapples += 1\">🍍</button>"
        end

        # TODO: can tag be empty? See <https://pt-br.reactjs.org/docs/fragments.html#short-syntax> for a usage.
    end

    @testset "Common HTML elements" begin
        # The Pareto slice (>80%) of <https://www.advancedwebranking.com/html/>
        @test htm"" |> r == ""
        @test htm"" |> r == ""
        @test htm"" |> r == ""
        @test htm"<title/>" |> r == "<title></title>"
        @test htm"<meta/>" |> r == "<meta />"
        @test htm"
" |> r == "

"
        @test htm"" |> r == ""
        @test htm"<link/>" |> r == "<link />"
        @test htm"" |> r == ""
        @test htm"
" |> r == "

"
        @test htm"
" |> r == "

"
        @test htm"
" |> r == "

"
        @test htm""

        # Surprisingly common specifics
        @test_broken htm"""""" |> r == """"""
        @test htm"" |> r == """"""
        @test htm"" |> r == """"""
        @test htm"
My favorite color is blue red!
" |> r == "
My favorite color is bluered!
"

        @test htm"<source src=horse.ogg type='audio/ogg'/>" |> r == """<source src="horse.ogg" type="audio/ogg" />"""
        @test htm"<template>

Flower

        </template>" |> r == """<template>
Flower
</template>"""
    end

    @testset "Common SVG elements" begin
        # TODO: make this more complete
        @test htm"""<svg width=100 height=100>
            <circle cx=50 cy=50 r=40 stroke=green stroke-width=4 fill=yellow/>
        </svg>""" |> r == """<svg height="100" width="100"><circle cy="50" stroke-width="4" stroke="green" r="40" fill="yellow" cx="50" /></svg>"""
    end

    @testset "Return types" begin
        @test htm"
🍍
" isa Hyperscript.Node
        @test htm"🍍" isa String
        @test htm"" === nothing
        @test htm"

" isa Vector{Hyperscript.Node{Hyperscript.HTMLSVG}}
        @test htm"🍍
" isa Vector{Any}
    end

    @testset "Whitespace" begin
        @test htm"""

                🍍

        """ |> r == "
\n        🍍\n
"
    end

    @testset "Edge cases" begin
        @test htm"" |> r == ""
    end

    @testset "Interpolations" begin
        @testset "Variables" begin
            @testset "As children" begin
                child = "🍍"
                @test htm"
$(child)
" |> r == "
🍍
"
            end

            @testset "As attributes" begin
                @testset "As keys" for key in ("class", :class, nothing, true, missing, 1, 1.0, [1, 2, 3], (1, 2, 3), "class" => "fruit")
                    @test_throws MethodError htm"

"
                end

                @testset "As values" begin
                    value = "fruit"
                    @test htm"

" |> r == "

"
                end
            end

            @testset "As tags" begin
                tag = "div"
                @testset "Matching end-tag" begin
                    @testset "Complete tags" begin
                        @test htm"<$(tag)></$(tag)>" |> r == "

"
                    end

                    @testset "Partial tags" begin
                        @test htm"<h$(tag)></h$(tag)>" |> r == "<hdiv></hdiv>"
                    end
                end

                @testset "Universal end-tag" begin
                    @testset "Complete tags" begin
                        @test htm"<$(tag)><//>" |> r == "

"
                    end

                    @testset "Partial tags" begin
                        @test htm"<h$(tag)><//>" |> r == "<hdiv></hdiv>"
                    end
                end
            end
        end

        @testset "Literals" begin
            @testset "As children" begin
                @test htm"
$(nothing)
" |> r == "

"
                @test htm"
$(missing)
" |> r == "
missing
"
                @test htm"
$(1)
" |> r == "
1
"
                @test htm"
$(1.0)
" |> r == "
1.0
"
                @test htm"
$(true)
" |> r == "

"
                @test htm"
$(:symbol)
" |> r == "
symbol
"
                @test htm"
$(\"string\")
" |> r == "
string
"
                @test htm"
$([1, 2, 3])
" |> r == "
123
"
                @test htm"
$((1, 2, 3))
" |> r == "
123
"
                @test htm"
$(\"class\" => \"fruit\")
" |> r == "
"class" => "fruit"
"
                @test htm"
$(Dict(\"class\" => \"fruit\"))
" |> r == "
Dict("class" => "fruit")
"

                @testset "Exotic objects" begin
                    @test htm"
$(md\"# 🍍\")
" |> r == "

🍍
\n

"
                    @test htm"
$(html\"

\")
" |> r == "

"
                    @test htm"
$(HTML(\"

\"))
" |> r == "

"
                end
            end

            @testset "As attributes" begin
                @testset "As keys" begin
                    @test_throws MethodError htm"
 \"fruit\"))=pineapple>
" |> r == "

"
                end

                @testset "As values" begin
                    @test htm"

" |> r == "

"
                    @test_throws TypeError htm"

" |> r == "

"
                    @test htm"

" |> r == "

"
                    @test htm"

" |> r == "

"
                    @test htm"

" |> r == "

"
                    @test htm"

" |> r == "

"
                    @test htm"

" |> r == "

"
                    @test htm"
 \"fruit\")>
" |> r == "

"
                    @test htm"
 \"fruit\"))>
" |> r == "

"
                end
            end

            @testset "As tags" begin
                @testset "Matching end-tag" begin
                    @testset "Complete tags" begin
                        @test_throws MethodError htm"<$(nothing)></$(nothing)>" |> r == "<nothing></nothing>"
                        @test_throws MethodError htm"<$(missing)></$(missing)>" |> r == "<missing></missing>"
                        @test_throws MethodError htm"<$(1)></$(1)>" |> r == "<1></1>"
                        @test_throws MethodError htm"<$(1.0)></$(1.0)>" |> r == "<1.0></1.0>"
                        @test_throws MethodError htm"<$(true)></$(true)>" |> r == "<true></true>"
                        @test_throws MethodError htm"<$(:symbol)></$(:symbol)>" |> r == "<symbol></symbol>"
                        @test htm"<$(\"string\")></$(\"string\")>" |> r == "<string></string>"
                        @test htm"<$([1, 2, 3])></$([1, 2, 3])>" |> r == "<123></123>"
                        @test_throws MethodError htm"<$(\"class\" => \"fruit\")></$(\"class\" => \"fruit\")>" |> r == "<"class" => "fruit"></"class" => "fruit">"
                        @test_throws MethodError htm"<$(Dict(\"class\" => \"fruit\"))></$(Dict(\"class\" => \"fruit\"))>" |> r == "<Dict("class" => "fruit")></Dict("class" => "fruit")>"
                    end

                    @testset "Partial tags" begin
                        @test htm"<h$(nothing)></h$(nothing)>" |> r == "<hnothing></hnothing>"
                        @test htm"<h$(missing)></h$(missing)>" |> r == "<hmissing></hmissing>"
                        @test htm"<h$(1)></h$(1)>" |> r == "

"
                        @test htm"<h$(1.0)></h$(1.0)>" |> r == "

"
                        @test htm"<h$(true)></h$(true)>" |> r == "<htrue></htrue>"
                        @test htm"<h$(:symbol)></h$(:symbol)>" |> r == "<hsymbol></hsymbol>"
                        @test htm"<h$(\"string\")></h$(\"string\")>" |> r == "<hstring></hstring>"
                        @test htm"<h$([1, 2, 3])></h$([1, 2, 3])>" |> r == "<h123></h123>"
                        @test htm"<h$(\"class\" => \"fruit\")></h$(\"class\" => \"fruit\")>" |> r == "<h"class" => "fruit"></h"class" => "fruit">"
                        @test htm"<h$(Dict(\"class\" => \"fruit\"))></h$(Dict(\"class\" => \"fruit\"))>" |> r == "<hDict("class" => "fruit")></hDict("class" => "fruit")>"
                    end
                end

                @testset "Universal end-tag" begin
                    @testset "Complete tags" begin
                        @test_throws MethodError htm"<$(nothing)><//>" |> r == "<nothing></nothing>"
                        @test_throws MethodError htm"<$(missing)><//>" |> r == "<missing></missing>"
                        @test_throws MethodError htm"<$(1)><//>" |> r == "<1></1>"
                        @test_throws MethodError htm"<$(1.0)><//>" |> r == "<1.0></1.0>"
                        @test_throws MethodError htm"<$(true)><//>" |> r == "<true></true>"
                        @test_throws MethodError htm"<$(:symbol)><//>" |> r == "<symbol></symbol>"
                        @test htm"<$(\"string\")><//>" |> r == "<string></string>"
                        @test htm"<$([1, 2, 3])><//>" |> r == "<123></123>"
                        @test_throws MethodError htm"<$(\"class\" => \"fruit\")><//>" |> r == "<"class" => "fruit"></"class" => "fruit">"
                        @test_throws MethodError htm"<$(Dict(\"class\" => \"fruit\"))><//>" |> r == "<Dict("class" => "fruit")></Dict("class" => "fruit")>"
                    end

                    @testset "Partial tags" begin
                        @test htm"<h$(nothing)><//>" |> r == "<hnothing></hnothing>"
                        @test htm"<h$(missing)><//>" |> r == "<hmissing></hmissing>"
                        @test htm"<h$(1)><//>" |> r == "

"
                        @test htm"<h$(1.0)><//>" |> r == "

"
                        @test htm"<h$(true)><//>" |> r == "<htrue></htrue>"
                        @test htm"<h$(:symbol)><//>" |> r == "<hsymbol></hsymbol>"
                        @test htm"<h$(\"string\")><//>" |> r == "<hstring></hstring>"
                        @test htm"<h$([1, 2, 3])><//>" |> r == "<h123></h123>"
                        @test htm"<h$(\"class\" => \"fruit\")><//>" |> r == "<h"class" => "fruit"></h"class" => "fruit">"
                        @test htm"<h$(Dict(\"class\" => \"fruit\"))><//>" |> r == "<hDict("class" => "fruit")></hDict("class" => "fruit")>"
                    end
                end
            end
        end
    end

    @testset "Escaped characters`" begin
        @test htm"
Look, Ma, $(\"automatic escaping\")!
" |> r == "
Look, Ma, <em>automatic escaping</em>!
"

        @testset "Escaped dollar signs`" begin
            @test htm"
\$notvar
" |> r == raw"
$notvar
"
            @test htm"
\$(notvar)
" |> r == raw"
$(notvar)
"

            @test htm"<\$notvar />" |> r == raw"<$notvar></$notvar>"
            @test_broken htm"
" |> r == raw"""

"""

            @test htm"<\$(notvar) />" |> r == raw"<$(notvar)></$(notvar)>"
            @test_broken htm"
" |> r == raw"""

"""

            @test_broken htm"
" |> r == raw"""

"""
            @test_broken htm"
" |> r == raw"""

"""

            @test_broken htm"
" |> r == raw"""

"""
            @test_broken htm"
" |> r == raw"""

"""

            @test_broken htm"
" |> r == raw"""

"""
            @test_broken htm"
" |> r == raw"""

"""

            @test_broken htm"
" |> r == raw"""

"""
            @test_broken htm"
" |> r == raw"""

"""
        end

        @testset "HTML entity escape trick" begin
            @test htm"
🍍 🍌
" |> r == "
🍍 🍌
"
            @test htm"
🍍 🍌
" |> r == "
🍍&nbsp;🍌
"
            @test htm"
🍍$(html\" \")🍌
" |> r == "
🍍 🍌
"
            @test htm"
🍍$(HTML(\" \"))🍌
" |> r == "
🍍 🍌
"
        end
    end

    @testset "Create elements" begin
        @test HTM.create_element("div", []) |> r == "

"
        @test HTM.create_element("div", [], []) |> r == "

"
        @test HTM.create_element("div", [], "Hi!") |> r == "
Hi!
"
        @test HTM.create_element("div", [], ["Hi!"]) |> r == "
Hi!
"
        @test HTM.create_element("div", ["class" => "fruit"]) |> r == "

"
        @test HTM.create_element("div", ["class" => "fruit"], "Hi!") |> r == "
Hi!
"
        @test HTM.create_element("div", ["class" => "fruit"], "Hi ", "there!") |> r == "
Hi there!
"

        @test HTM.create_element("button", ["class" => "fruit", "disabled" => nothing]) |> r == "<button class=\"fruit\" disabled></button>"
        @test HTM.create_element("button", ["class" => "fruit", "disabled" => nothing], "Click me") |> r == "<button class=\"fruit\" disabled>Click me</button>"

        @test HTM.create_element("circle", ["fill" => "orange"]) |> r == "<circle fill=\"orange\" />"
    end

    @testset "Internal representation" begin
        @test HTM.parsenode(IOBuffer("
")) == HTM.Node(["div"])
        @test HTM.parsenode(IOBuffer("
Hi!
")) == HTM.Node(["div"], [], ["Hi!"])
        @test HTM.parsenode(IOBuffer("
")) == HTM.Node(["div"], [:class => ["fruit"]])
        @test HTM.parsenode(IOBuffer("
Hi!
")) == HTM.Node(["div"], [:class => ["fruit"]], ["Hi!"])
        @test HTM.parsenode(IOBuffer("
Hi there!
")) == HTM.Node(["div"], [:class => ["fruit"]], ["Hi there!"])

        @test HTM.parsenode(IOBuffer("<button class=fruit disabled />")) == HTM.Node(["button"], [:class => ["fruit"], :disabled => [raw"$(true)"]])
        @test HTM.parsenode(IOBuffer("<button class=fruit disabled>Click me</button>")) == HTM.Node(["button"], [:class => ["fruit"], :disabled => [raw"$(true)"]], ["Click me"])

        @test HTM.parsenode(IOBuffer("<circle fill=orange />")) == HTM.Node(["circle"], [:fill => ["orange"]])

        @testset "Stress tests" begin
            quotes = ("", '"', '\'')
            separators = ("", ' ', '\n', '\t', "  ")

            @testset "Basic syntax" begin
                voiddiv = HTM.parsenode(IOBuffer("
"))
                litclass = HTM.parsenode(IOBuffer("
"))
                offlitclass = HTM.parsenode(IOBuffer("
"))

                let c = "Hi there!"
                    nonvoiddiv = HTM.parsenode(IOBuffer("
$(c)
"))
                    fruitnonvoiddiv = HTM.parsenode(IOBuffer("
$(c)
"))
                    draggablefruitnonvoiddiv = HTM.parsenode(IOBuffer("
$(c)
"))

                    @testset "Separator `$(s)`" for s in separators
                        @test HTM.parsenode(IOBuffer("
")) == voiddiv

                        @test HTM.parsenode(IOBuffer("

")) == voiddiv

                        @test HTM.parsenode(IOBuffer("
$(c)
")) == nonvoiddiv

                        @testset "Quote `$(q)`" for q in quotes
                            @test HTM.parsenode(IOBuffer("
")) == litclass
                            @test HTM.parsenode(IOBuffer("
")) == offlitclass

                            @test HTM.parsenode(IOBuffer("

")) == litclass
                            @test HTM.parsenode(IOBuffer("

")) == offlitclass

                            @test HTM.parsenode(IOBuffer("
$(c)
")) == fruitnonvoiddiv
                            @test HTM.parsenode(IOBuffer("
$(c)
")) == draggablefruitnonvoiddiv
                        end
                    end
                end
            end

            @testset "Dollar signs" begin
                @test HTM.parsenode(IOBuffer("
\$notvar
")) == HTM.parsenode(IOBuffer(raw"
$notvar
"))
                @test HTM.parsenode(IOBuffer("
\$(notvar)
")) == HTM.parsenode(IOBuffer(raw"
$(notvar)
"))

                @testset "Separator `$(s)`" for s in separators
                    @test HTM.parsenode(IOBuffer("<\$notvar$(s)/>")) == HTM.parsenode(IOBuffer(raw"<$notvar></$notvar>"))
                    @test HTM.parsenode(IOBuffer("
")) == HTM.parsenode(IOBuffer(raw"""

"""))

                    @test HTM.parsenode(IOBuffer("<\$(notvar)$(s)/>")) == HTM.parsenode(IOBuffer(raw"<$(notvar)></$(notvar)>"))
                    @test HTM.parsenode(IOBuffer("
")) == HTM.parsenode(IOBuffer(raw"""

"""))

                    @testset "Quote `$(q)`" for q in quotes
                        @test HTM.parsenode(IOBuffer("
")) == HTM.parsenode(IOBuffer(raw"""

"""))
                        @test HTM.parsenode(IOBuffer("
")) == HTM.parsenode(IOBuffer(raw"""

"""))

                        @test HTM.parsenode(IOBuffer("
")) == HTM.parsenode(IOBuffer(raw"""

"""))
                        @test HTM.parsenode(IOBuffer("
")) == HTM.parsenode(IOBuffer(raw"""

"""))
                    end
                end
            end
        end
    end
end


text:
using RequiredKeywordArguments
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@rka f(;y) = y
@test f(y=2) == 2
try
    f()
catch e
    @test isa(e, RequiredKeywordArgumentError)
end

@rka g(x; y) = x*y
@test g(2,y=3) == 6
try
    g(2)
catch e
    @test isa(e, RequiredKeywordArgumentError)
end

@rka h(;x=1, y) = x*y
@test h(y=3) == 3
@test h(y=3, x=2) == 6
try
    h(x=2)
catch e
    @test isa(e, RequiredKeywordArgumentError)
    Base.showerror(Base.STDOUT, e)
end

k(;y=RequiredKeywordArgument("y", "I can't compute `k` without `y`.")) = 3*y
try
    k()
catch e
    @test isa(e, RequiredKeywordArgumentError)
    @test e.reason == "I can't compute `k` without `y`."
    Base.showerror(Base.STDOUT, e)
end
l(;y=RequiredKeywordArgument("y", reason="I can't compute `l` without `y`.")) = 4*y
try
    l()
catch e
    @test isa(e, RequiredKeywordArgumentError)
    @test e.reason == "I can't compute `l` without `y`."
end
m(;y=RKA("y", reason="I can't compute `m` without `y`.")) = 5*y
try
    m()
catch e
    @test isa(e, RequiredKeywordArgumentError)
    @test e.reason == "I can't compute `m` without `y`."
end


text:
using FastRounding

using Test

const SysFloat = Union{Float64, Float32}

a = 4/3
b = 1/7
c = -b

function testrounding(op::Function, a::T, b::T, mode::RoundingMode) where T<: SysFloat
    S = typeof(op(a,b))

    c = setrounding(BigFloat, mode) do
            S(op(big(a), big(b)))
    end
    return c
end

@testset "Test Rounding" begin
    @testset "a + b" begin
        @test testrounding(+, a, b, RoundUp) == add_round(a, b, RoundUp)
        @test testrounding(+, a, b, RoundDown) == add_round(a, b, RoundDown)
        @test testrounding(+, a, b, RoundToZero) == add_round(a, b, RoundToZero)
        @test testrounding(+, a, b, RoundNearest) == add_round(a, b, RoundNearest)
    end
    @testset "a - b" begin
        @test testrounding(-, a, b, RoundUp) == sub_round(a, b, RoundUp)
        @test testrounding(-, a, b, RoundDown) == sub_round(a, b, RoundDown)
        @test testrounding(-, a, b, RoundToZero) == sub_round(a, b, RoundToZero)
        @test testrounding(-, a, b, RoundNearest) == sub_round(a, b, RoundNearest)
    end
    @testset "a * b" begin
        @test testrounding(*, a, b, RoundUp) == mul_round(a, b, RoundUp)
        @test testrounding(*, a, b, RoundDown) == mul_round(a, b, RoundDown)
        @test testrounding(*, a, b, RoundToZero) == mul_round(a, b, RoundToZero)
        @test testrounding(*, a, b, RoundNearest) == mul_round(a, b, RoundNearest)
    end
    @testset "a / b" begin
        @test testrounding(/, a, b, RoundUp) == div_round(a, b, RoundUp)
        @test testrounding(/, a, b, RoundDown) == div_round(a, b, RoundDown)
        @test testrounding(/, a, b, RoundToZero) == div_round(a, b, RoundToZero)
        @test testrounding(/, a, b, RoundNearest) == div_round(a, b, RoundNearest)
    end
    @testset "a + c" begin
        @test testrounding(+, a, c, RoundUp) == add_round(a, c, RoundUp)
        @test testrounding(+, a, c, RoundDown) == add_round(a, c, RoundDown)
        @test testrounding(+, a, c, RoundToZero) == add_round(a, c, RoundToZero)
        @test testrounding(+, a, c, RoundNearest) == add_round(a, c, RoundNearest)
    end
    @testset "a - c" begin
        @test testrounding(-, a, c, RoundUp) == sub_round(a, c, RoundUp)
        @test testrounding(-, a, c, RoundDown) == sub_round(a, c, RoundDown)
        @test testrounding(-, a, c, RoundToZero) == sub_round(a, c, RoundToZero)
        @test testrounding(-, a, c, RoundNearest) == sub_round(a, c, RoundNearest)
    end
    @testset "a * c" begin
        @test testrounding(*, a, c, RoundUp) == mul_round(a, c, RoundUp)
        @test testrounding(*, a, c, RoundDown) == mul_round(a, c, RoundDown)
        @test testrounding(*, a, c, RoundToZero) == mul_round(a, c, RoundToZero)
        @test testrounding(*, a, c, RoundNearest) == mul_round(a, c, RoundNearest)
    end
    @testset "a / c" begin
        @test testrounding(/, a, c, RoundUp) == div_round(a, c, RoundUp)
        @test testrounding(/, a, c, RoundDown) == div_round(a, c, RoundDown)
        @test testrounding(/, a, c, RoundToZero) == div_round(a, c, RoundToZero)
        @test testrounding(/, a, c, RoundNearest) == div_round(a, c, RoundNearest)
    end
    @testset "1.0 + 0.0" begin
        @test testrounding(+, 1.0, 0.0, RoundUp) == 1.0
        @test testrounding(+, 1.0, 0.0, RoundDown) == 1.0
        @test testrounding(+, 1.0, 0.0, RoundToZero) == 1.0
        @test testrounding(+, 1.0, 0.0, RoundNearest) == 1.0
    end
    @testset "1.0 - 0.0" begin
        @test testrounding(-, 1.0, 0.0, RoundUp) == 1.0
        @test testrounding(-, 1.0, 0.0, RoundDown) == 1.0
        @test testrounding(-, 1.0, 0.0, RoundToZero) == 1.0
        @test testrounding(-, 1.0, 0.0, RoundNearest) == 1.0
    end
    @testset "1.0 * 1.0" begin
        @test testrounding(*, 1.0, 1.0, RoundUp) == 1.0
        @test testrounding(*, 1.0, 1.0, RoundDown) == 1.0
        @test testrounding(*, 1.0, 1.0, RoundToZero) == 1.0
        @test testrounding(*, 1.0, 1.0, RoundNearest) == 1.0
    end
    @testset "1.0 / 1.0" begin
        @test testrounding(/, 1.0, 1.0, RoundUp) == 1.0
        @test testrounding(/, 1.0, 1.0, RoundDown) == 1.0
        @test testrounding(/, 1.0, 1.0, RoundToZero) == 1.0
        @test testrounding(/, 1.0, 1.0, RoundNearest) == 1.0
    end
end


text:
"""
  Compute the nth Pernicious number
"""

using BenchmarkTools
using Printf
using Primes

# Get the value of of an integer from the command line.
n, = size(ARGS)
if n < 1
    println("Usage: test_pernicious_numbers.jl max_num")
    println("       --->      Please specify a positive integer.")
    exit()
end

max_num = parse(Int64, ARGS[1])

ispernicious(n::Int64) = isprime(count_ones(n))

nextpernicious(n::Int64) = begin n += 1; while !ispernicious(n) n += 1 end; return n end

function perniciouses(n::Int64)
    rst = Vector{Int64}(undef, n)
    old_p = 3
    for i in 2:n
        new_p = nextpernicious(old_p)
        old_p = new_p
    end
    return old_p
end

perniciouses(a::Integer, b::Integer) = filter(ispernicious, a:b)

println("Pernicious number of order: ", max_num)

@btime perniciouses(max_num)

#println("Pernicious number of order ", max_num, ": ", perniciouses(max_num))


text:
using MPAlign
using Test
using Test
using LightGraphs
using LinearAlgebra

@testset "MPAlign" begin
    n = 20
    λ = 3
    s = 0.4

    data_check = [[2.554459, -1.346828, 1.013629, -0.615002, -0.975230, 1.207213, -0.051287, -2.375779, -0.975230, -0.200246, 1.101704, -0.472741, -1.331169, -0.848312, -0.051287, -2.375779, -0.167027, 1.528282, -2.375779, -2.375779],
    [-0.184265, 0.524227, 0.093951, 0.178841, -0.287554, 0.047793, -0.366399, 0.178349, -0.287554, -0.352418, 0.222273, 0.623926, 0.130005, -0.239924, -0.366399, 0.178349, 0.387049, 0.790944, 0.178349, 0.178349],
    [-1.494844, -0.004521, -0.126428, 0.594263, 0.984577, 0.414440, 0.989069, 0.178349, 0.984577, 0.781970, -0.220506, 0.003317, 0.733688, 0.935675, 0.989069, 0.178349, 0.558987, -0.650504, 0.178349, 0.178349],
    [-2.375779, 0.689174, -0.843302, 0.178349, 0.178349, -0.843302, -0.332477, 1.200000, 0.178349, -0.332477, -0.843302, 0.178349, 0.689174, 0.178349, -0.332477, 1.200000, 0.178349, -0.843302, 1.200000, 1.200000],
    [-2.375779, 0.689174, -0.843302, 0.178349, 0.178349, -0.843302, -0.332477, 1.200000, 0.178349, -0.332477, -0.843302, 0.178349, 0.689174, 0.178349, -0.332477, 1.200000, 0.178349, -0.843302, 1.200000, 1.200000],
    [2.227597, -0.957353, 1.440599, -0.678833, -1.093163, 1.075827, -0.318501, -1.864954, -1.093163, -0.418694, 1.500382, -0.265845, -1.100465, -0.924522, -0.318501, -1.864954, 0.002392, 1.700918, -1.864954, -1.864954],
    [1.457630, -0.601845, 1.013347, 0.487322, -0.073449, 0.755498, 0.264249, -1.354128, -0.073449, 0.369071, 1.168702, -0.053885, -0.724662, -0.073463, 0.264249, -1.354128, 0.316029, 1.425065, -1.354128, -1.354128],
    [-0.885897, 0.247145, -0.001501, 0.322286, 0.486080, 0.467487, 0.354896, 0.178349, 0.486080, 0.211582, -0.004608, 0.285872, 0.593118, 0.339161, 0.354896, 0.178349, 0.628307, 0.091947, 0.178349, 0.178349],
    [-0.954578, -0.284592, 0.163629, 0.827259, 1.029978, 0.494080, 1.198952, -0.332477, 1.029978, 1.065663, 0.074967, -0.122873, 0.296726, 0.991634, 1.198952, -0.332477, 0.506253, -0.449332, -0.332477, -0.332477],
    [-1.793645, -0.087575, -0.514993, 0.815151, 1.238781, 0.044013, 1.139419, 0.178349, 1.238781, 1.010177, -0.601054, -0.085596, 0.616506, 1.167244, 1.139419, 0.178349, 0.279952, -1.109414, 0.178349, 0.178349],
    [-1.494844, -0.004521, -0.126428, 0.594263, 0.984577, 0.414440, 0.989069, 0.178349, 0.984577, 0.781970, -0.220506, 0.003317, 0.733688, 0.935675, 0.989069, 0.178349, 0.558987, -0.650504, 0.178349, 0.178349],
    [-1.323959, 0.646442, -0.438503, 0.204155, 0.155202, -0.247447, -0.168178, 0.689174, 0.155202, -0.210661, -0.396930, 0.398289, 0.568804, 0.099775, -0.168178, 0.689174, 0.356128, -0.110561, 0.689174, 0.689174],
    [-1.278386, 0.673417, -0.346934, 0.221034, 0.034597, -0.313825, -0.211020, 0.689174, 0.034597, -0.224443, -0.303950, 0.422651, 0.489912, 0.073833, -0.211020, 0.689174, 0.318855, -0.095885, 0.689174, 0.689174],
    [-0.516603, -0.199114, 0.753929, 0.719172, 0.606255, 0.502954, 0.765130, -0.332477, 0.606255, 0.799976, 0.656983, 0.020333, 0.140557, 0.637789, 0.765130, -0.332477, 0.587032, 0.095781, -0.332477, -0.332477],
    [-0.524234, 0.398365, 0.117482, 0.174511, -0.047988, 0.265530, 0.009865, 0.178349, -0.047988, -0.084587, 0.187006, 0.613759, 0.326588, 0.051653, 0.009865, 0.178349, 0.455814, 0.471236, 0.178349, 0.178349],
    [-2.126634, 0.287073, -0.633868, 0.787181, 0.894609, -0.360548, 0.570338, 0.689174, 0.894609, 0.521494, -0.687224, 0.088912, 0.726697, 0.898115, 0.570338, 0.689174, 0.279251, -1.008489, 0.689174, 0.689174],
    [-1.091078, -0.374393, -0.032651, 0.750047, 1.150923, 0.380182, 1.243659, -0.332477, 1.150923, 1.194621, -0.150130, -0.165855, 0.365714, 1.006261, 1.243659, -0.332477, 0.334640, -0.776113, -0.332477, -0.332477],
    [0.458638, -0.379748, 0.894856, 0.402349, 0.220159, 1.223579, 0.596386, -0.843302, 0.220159, 0.488143, 0.903540, 0.001057, -0.109334, 0.247049, 0.596386, -0.843302, 0.646117, 0.882290, -0.843302, -0.843302],
    [-1.614477, -0.040615, -0.129549, 0.705434, 1.083739, 0.099925, 0.960101, 0.178349, 1.083739, 0.845376, -0.242586, -0.045018, 0.603338, 0.974986, 0.960101, 0.178349, 0.424062, -0.764559, 0.178349, 0.178349],
    [-2.375779, 0.689174, -0.843302, 0.178349, 0.178349, -0.843302, -0.332477, 1.200000, 0.178349, -0.332477, -0.843302, 0.178349, 0.689174, 0.178349, -0.332477, 1.200000, 0.178349, -0.843302, 1.200000, 1.200000]]

    data_check2 = [[2.341680, -1.365904, 0.499407, -1.071828, -1.350649, 0.502508, -0.749434, -2.375779, -1.350649, -0.749434, 0.619942, -0.578132, -1.634333, -1.350649, -0.749434, -2.375779, -0.704786, 1.158504, -2.375779, -2.375779],
    [-0.261587, 0.519571, -0.014854, 0.123443, -0.244960, 0.022747, -0.380084, 0.178349, -0.244960, -0.380084, 0.128128, 0.586818, 0.138806, -0.244960, -0.380084, 0.178349, 0.367401, 0.741968, 0.178349, 0.178349],
    [-1.501218, 0.004633, -0.325959, 0.191606, 0.378259, -0.323899, 0.236592, 0.178349, 0.378259, 0.236592, -0.403680, -0.092076, 0.237215, 0.378259, 0.236592, 0.178349, -0.001149, -0.822246, 0.178349, 0.178349],
    [-2.375779, 0.689174, -0.843302, 0.178349, 0.178349, -0.843302, -0.332477, 1.200000, 0.178349, -0.332477, -0.843302, 0.178349, 0.689174, 0.178349, -0.332477, 1.200000, 0.178349, -0.843302, 1.200000, 1.200000],
    [-2.375779, 0.689174, -0.843302, 0.178349, 0.178349, -0.843302, -0.332477, 1.200000, 0.178349, -0.332477, -0.843302, 0.178349, 0.689174, 0.178349, -0.332477, 1.200000, 0.178349, -0.843302, 1.200000, 1.200000],
    [2.018585, -0.969421, 0.598961, -0.822562, -1.070465, 0.581264, -0.564170, -1.864954, -1.070465, -0.564170, 0.697033, -0.342649, -1.215544, -1.070465, -0.564170, -1.864954, -0.408870, 1.138457, -1.864954, -1.864954],
    [1.250451, -0.601541, 0.369529, -0.311584, -0.620295, 0.394170, -0.301323, -1.354128, -0.620295, -0.301323, 0.490080, -0.104822, -0.883136, -0.620295, -0.301323, -1.354128, -0.225144, 0.998711, -1.354128, -1.354128],
    [-0.831924, 0.279415, -0.096843, 0.174277, 0.090550, -0.089423, -0.030350, 0.178349, 0.090550, -0.030350, -0.079310, 0.249045, 0.219000, 0.090550, -0.030350, 0.178349, 0.222042, -0.055388, 0.178349, 0.178349],
    [-1.052070, -0.291206, -0.198370, 0.151353, 0.282696, -0.184892, 0.304574, -0.332477, 0.282696, 0.304574, -0.268838, -0.225338, -0.111559, 0.282696, 0.304574, -0.332477, -0.149511, -0.676936, -0.332477, -0.332477],
    [-1.847771, -0.089896, -0.523658, 0.277113, 0.494442, -0.501378, 0.311430, 0.178349, 0.494442, 0.311430, -0.611657, -0.203051, 0.193373, 0.494442, 0.311430, 0.178349, -0.120488, -1.101600, 0.178349, 0.178349],
    [-1.501218, 0.004633, -0.325959, 0.191606, 0.378259, -0.323899, 0.236592, 0.178349, 0.378259, 0.236592, -0.403680, -0.092076, 0.237215, 0.378259, 0.236592, 0.178349, -0.001149, -0.822246, 0.178349, 0.178349],
    [-1.315154, 0.652493, -0.352207, 0.208426, 0.037262, -0.345487, -0.264559, 0.689174, 0.037262, -0.264559, -0.300586, 0.397119, 0.468357, 0.037262, -0.264559, 0.689174, 0.331796, -0.085008, 0.689174, 0.689174],
    [-1.315154, 0.652493, -0.352207, 0.208426, 0.037262, -0.345487, -0.264559, 0.689174, 0.037262, -0.264559, -0.300586, 0.397119, 0.468357, 0.037262, -0.264559, 0.689174, 0.331796, -0.085008, 0.689174, 0.689174],
    [-0.706605, -0.196349, -0.037332, 0.104717, 0.173243, -0.040016, 0.220251, -0.332477, 0.173243, 0.220251, -0.085279, -0.118141, -0.092777, 0.173243, 0.220251, -0.332477, -0.047669, -0.362737, -0.332477, -0.332477],
    [-0.564017, 0.397119, -0.066722, 0.154819, -0.071520, -0.032644, -0.202631, 0.178349, -0.071520, -0.202631, 0.018228, 0.475945, 0.182361, -0.071520, -0.202631, 0.178349, 0.300224, 0.336306, 0.178349, 0.178349],
    [-2.175616, 0.283536, -0.703951, 0.359891, 0.461001, -0.683502, 0.086742, 0.689174, 0.461001, 0.086742, -0.750165, -0.009726, 0.458666, 0.461001, 0.086742, 0.689174, 0.028317, -1.024733, 0.689174, 0.689174],
    [-1.219802, -0.366465, -0.246839, 0.114520, 0.333389, -0.239310, 0.359074, -0.332477, 0.333389, 0.359074, -0.346562, -0.318640, -0.101530, 0.333389, 0.359074, -0.332477, -0.206998, -0.882055, -0.332477, -0.332477],
    [0.423066, -0.362742, 0.308715, -0.153433, -0.258866, 0.306032, -0.035438, -0.843302, -0.258866, -0.035438, 0.333311, -0.054647, -0.446516, -0.258866, -0.035438, -0.843302, -0.063714, 0.393139, -0.843302, -0.843302],
    [-1.667888, -0.033709, -0.427301, 0.251748, 0.440922, -0.415132, 0.276257, 0.178349, 0.440922, 0.276257, -0.506857, -0.143111, 0.211547, 0.440922, 0.276257, 0.178349, -0.058401, -0.943435, 0.178349, 0.178349],
    [-2.375779, 0.689174, -0.843302, 0.178349, 0.178349, -0.843302, -0.332477, 1.200000, 0.178349, -0.332477, -0.843302, 0.178349, 0.689174, 0.178349, -0.332477, 1.200000, 0.178349, -0.843302, 1.200000, 1.200000]
    ]

    G1, G2 = read_graphs(n,"graphs_20.txt")
    PG = Pair_IO(G1,G2,n,λ,s)
    M_data = create_matrix_check(data_check2)
    @testset "message passing" begin
    @test check_data(PG,M_data,1)
    M_data = create_matrix_check(data_check)
    @test check_data(PG,M_data,4)
    end
    @testset "overlap function" begin
    ov1, _, v1, _ = eval_M(PG,Matrix(I,n,n))
    @test ov1 == n
    @test v1 == n
    end
    @testset "matched edges" begin
    m1, m2 = match_edges(PG,1:n,1:n)
    @test m1 == m2
    @test m2 == ne(intersect(PG.N1.G,PG.N2.G))
    end
    n = 100
    λ = 3
    s = 1#0.55
    PG = Pair_ER(n,λ,s)
    m1, m2 = match_edges(PG,PG.perm_G1G2,PG.perm_G2G1)
    @test m1 == m2
    @test m1 == ne(PG.Ginter)
    @testset "parameter computation" begin
    delta = 0
    for t in 0:0.1:1
        a,b = MPAlign.compute_ab(n,λ,t)
        delta += (a+(1-a)*b^2 - λ*s/n)^2+((1-a)*b*(1-b) - λ*(1-s)/n)^2
    end
    @test delta < 0.01
    end
end


text:
import DSP: resample
import SeisIO: resample, resample!
export bandpass, bandpass!, bandstop, bandstop!, lowpass, lowpass!
export highpass, highpass!, taper, taper!, envelope
export resample, resample!

"""
   bandpass!(A,freqmin,freqmax,fs,corners=4,zerophase=true)

Butterworth-Bandpass Filter.

Filter data `A` from `freqmin` to `freqmax` using `corners` corners.

# Arguments
- `A::AbstractArray`: Data to filter
- `freqmin::Real`: Pass band low corner frequency.
- `freqmax::Real`: Pass band high corner frequency.
- `fs::Real`: Sampling rate in Hz.
- `fs::Int`: Filter corners / order.
- `zerophase::Bool`: If True, apply filter once forwards and once backwards.
This results in twice the filter order but zero phase shift in
the resulting filtered trace.
"""
function bandpass!(A::AbstractArray{<:AbstractFloat},
                   freqmin::Real, freqmax::Real, fs::Real;
                   corners::Int=4, zerophase::Bool=true)
   T = eltype(A)
   fe = T(0.5 * fs)
   low = T(freqmin / fe)
   high = T(freqmax / fe)

   # warn if above Nyquist frequency
   if high - oneunit(high) > -1e-6
       @warn "Selected high corner frequency ($freqmax) of bandpass is at or
       above Nyquist ($fe). Applying a high-pass instead."
       highpass!(A,freqmin,fs,corners=corners,zerophase=zerophase)
       return nothing
   end

   # throw error if low above Nyquist frequency
   if low >= 1
       throw(ArgumentError("Selected low corner frequency is above Nyquist."))
   end

   # create filter
   responsetype = Bandpass(T(freqmin), T(freqmax); fs=fs)
   designmethod = Butterworth(T,corners)

   # use gpu-specific kernel if on the GPU
   if isa(A,AbstractGPUArray)
       gpufilter!(A,responsetype,designmethod)
       return nothing
   end

   # filter if on the CPU
   if zerophase
       A[:,:] .= filtfilt(digitalfilter(responsetype, designmethod), @view(A[:,:]))
   else
       A[:,:] .= filt(digitalfilter(responsetype, designmethod), @view(A[:,:]))
   end

   return nothing
end
bandpass(A::AbstractArray{<:AbstractFloat},freqmin::Real,
         freqmax::Real, fs::Real; corners::Int=4,zerophase::Bool=true) =
         (U = deepcopy(A);bandpass!(U,freqmin,freqmax, fs, corners=corners,
          zerophase=zerophase);return U)
bandpass!(R::RawData,freqmin::Real,freqmax::Real;
          corners::Int=4,zerophase::Bool=true) = (bandpass!(R.x,freqmin,freqmax,
          R.fs,corners=corners,zerophase=zerophase);setfield!(R,:freqmin,Float64(freqmin));
          setfield!(R,:freqmax,Float64(min(freqmax,R.fs/2)));return nothing)
bandpass(R::RawData,freqmin::Real,freqmax::Real;
        corners::Int=4,zerophase::Bool=true) = (U = deepcopy(R);bandpass!(U.x,
        freqmin,freqmax,U.fs,corners=corners,zerophase=zerophase);
        setfield!(U,:freqmin,Float64(freqmin));
        setfield!(U,:freqmax,Float64(min(freqmax,U.fs/2)));return U)
bandpass!(C::CorrData,freqmin::Real,freqmax::Real;
        corners::Int=4,zerophase::Bool=true) = (bandpass!(C.corr,freqmin,freqmax,
        C.fs,corners=corners,zerophase=zerophase);setfield!(C,:freqmin,Float64(freqmin));
        setfield!(C,:freqmax,Float64(min(freqmax,C.fs/2)));return nothing)
bandpass(C::CorrData,freqmin::Real,freqmax::Real;
        corners::Int=4,zerophase::Bool=true) = (U = deepcopy(C);bandpass!(U.corr,
        freqmin,freqmax,U.fs,corners=corners,zerophase=zerophase);
        setfield!(U,:freqmin,Float64(freqmin));
        setfield!(U,:freqmax,Float64(min(freqmax,U.fs/2)));return U)
bandpass!(C::SeisChannel, freqmin::Real, freqmax::Real;
                   corners::Int=4, zerophase::Bool=true) = filtfilt!(C,
                   fl=Float64(freqmin),fh=Float64(freqmax),np=corners,rt="Bandpass")
bandpass(C::SeisChannel,freqmin::Real, freqmax::Real; corners::Int=4,
         zerophase::Bool=true) = (U = deepcopy(C);bandpass!(U,freqmin,freqmax,
                                  corners=corners,zerophase=zerophase);return U)
bandpass!(S::SeisData, freqmin::Real, freqmax::Real;
                   corners::Int=4, zerophase::Bool=true) = filtfilt!(S,
                   fl=Float64(freqmin),fh=Float64(freqmax),np=corners,rt="Bandpass")
bandpass(S::SeisData, freqmin::Real, freqmax::Real; corners::Int=4,
         zerophase::Bool=true) = filtfilt(S,fl=Float64(freqmin),fh=Float64(freqmax),np=corners,rt="Bandpass")
bandpass!(N::NodalData,freqmin::Real,freqmax::Real;
         corners::Int=4,zerophase::Bool=true) = (bandpass!(N.data,freqmin,freqmax,
         N.fs[1],corners=corners,zerophase=zerophase);return nothing)
bandpass(N::NodalData,freqmin::Real,freqmax::Real;
       corners::Int=4,zerophase::Bool=true) = (U = deepcopy(N);bandpass!(N.data,
       freqmin,freqmax,N.fs[1],corners=corners,zerophase=zerophase);
       return U)

"""
  bandstop!(A,freqmin,freqmax,fs,corners=4,zerophase=true)

Butterworth-Bandstop Filter.

Filter data `A` removing data between frequencies `freqmin` to `freqmax` using
`corners` corners.

# Arguments
- `A::AbstractArray`: Data to filter
- `freqmin::Real`: Stop band low corner frequency.
- `freqmax::Real`: Stop band high corner frequency.
- `fs::Real`: Sampling rate in Hz.
- `fs::Int`: Filter corners / order.
- `zerophase::Bool`: If True, apply filter once forwards and once backwards.
This results in twice the filter order but zero phase shift in
the resulting filtered trace.
"""
function bandstop!(A::AbstractArray{<:AbstractFloat},
                 freqmin::Real,freqmax::Real,fs::Real;
                 corners::Int=4, zerophase::Bool=true)
    T = eltype(A)
    fe = T(0.5 * fs)
    low = T(freqmin / fe)
    high = T(freqmax / fe)

    # throw error if high above Nyquist frequency
    if high >= 1
        throw(ArgumentError("Selected high corner frequency $freqmax Hz is above Nyquist $(fs/2) Hz."))
    end

    # throw error if low above Nyquist frequency
    if low >= 1
        throw(ArgumentError("Selected low corner frequency $freqmin Hz is above Nyquist $(fs/2) Hz."))
    end

    # create filter
    responsetype = Bandstop(T(freqmin), T(freqmax); fs=fs)
    designmethod = Butterworth(T,corners)

    # use gpu-specific kernel if on the GPU
    if isa(A,AbstractGPUArray)
        gpufilter!(A,responsetype,designmethod)
        return nothing
    end

    # filter if on the CPU
    if zerophase
        A[:,:] .= filtfilt(digitalfilter(responsetype, designmethod), @view(A[:,:]))
    else
        A[:,:] .= filt(digitalfilter(responsetype, designmethod), @view(A[:,:]))
    end

    return nothing
end
bandstop(A::AbstractArray{<:AbstractFloat},freqmin::Real,
      freqmax::Real, fs::Real; corners::Int=4,zerophase::Bool=true) =
      (U = deepcopy(A);bandstop!(U,freqmin,freqmax,fs,corners=corners,
      zerophase=zerophase);return U)
bandstop!(R::RawData,freqmin::Real,freqmax::Real;
 corners::Int=4,zerophase::Bool=true) = bandstop!(R.x,freqmin,freqmax,
 R.fs,corners=corners,zerophase=zerophase)
bandstop(R::RawData,freqmin::Real,freqmax::Real;
corners::Int=4,zerophase::Bool=true) = (U = deepcopy(R);bandstop!(U.x,
freqmin,freqmax,U.fs,corners=corners,zerophase=zerophase);return U)
bandstop!(C::CorrData,freqmin::Real,freqmax::Real;
 corners::Int=4,zerophase::Bool=true) = bandstop!(C.corr,freqmin,freqmax,
 C.fs,corners=corners,zerophase=zerophase)
bandstop(C::CorrData,freqmin::Real,freqmax::Real;
corners::Int=4,zerophase::Bool=true) = (U = deepcopy(C);bandstop!(U.corr,
freqmin,freqmax,U.fs,corners=corners,zerophase=zerophase);return U)
bandstop!(C::SeisChannel, freqmin::Real, freqmax::Real;
                   corners::Int=4, zerophase::Bool=true) = filtfilt!(C,
                   fl=Float64(freqmin),fh=Float64(freqmax),np=corners,rt="Bandstop")
bandstop(C::SeisChannel,freqmin::Real, freqmax::Real; corners::Int=4,
         zerophase::Bool=true) = filtfilt(C,
         fl=Float64(freqmin),fh=Float64(freqmax),np=corners,rt="Bandstop")
bandstop!(S::SeisData, freqmin::Real, freqmax::Real;
                   corners::Int=4, zerophase::Bool=true) = filtfilt!(S,
                   fl=Float64(freqmin),fh=Float64(freqmax),np=corners,rt="Bandstop")
bandstop(S::SeisData,freqmin::Real, freqmax::Real; corners::Int=4,
         zerophase::Bool=true) = filtfilt(S,
         fl=Float64(freqmin),fh=Float64(freqmax),np=corners,rt="Bandstop")
bandstop!(N::NodalData,freqmin::Real,freqmax::Real;
    corners::Int=4,zerophase::Bool=true) = bandstop!(N.data,freqmin,freqmax,
    N.fs[1],corners=corners,zerophase=zerophase)
bandstop(N::NodalData,freqmin::Real,freqmax::Real;
    corners::Int=4,zerophase::Bool=true) = (U = deepcopy(N);bandstop!(N.data,
    freqmin,freqmax,U.fs[1],corners=corners,zerophase=zerophase);return U)

"""
lowpass(A,freq,fs,corners=4,zerophase=true)

Butterworth-Lowpass Filter.

Filter data `A` over certain frequency `freq` using `corners` corners.

# Arguments
- `A::AbstractArray`: Data to filter
- `freq::Real`: Filter corner frequency.
- `fs::Real`: Sampling rate in Hz.
- `fs::Int`: Filter corners / order.
- `zerophase::Bool`: If True, apply filter once forwards and once backwards.
This results in twice the filter order but zero phase shift in
the resulting filtered trace.
"""
function lowpass!(A::AbstractArray{<:AbstractFloat},freq::Real,
                fs::Real; corners::Int=4, zerophase::Bool=true)
    fe = 0.5 * fs
    f = freq / fe
    T = eltype(A)

    # warn if above Nyquist frequency
    if f >= 1
        @warn """Selected corner frequency ($freq) is
        above Nyquist ($fe). No filter applied."""
        return nothing
    end

    # create filter
    responsetype = Lowpass(T(freq); fs=fs)
    designmethod = Butterworth(T,corners)

    # use gpu-specific kernel if on the GPU
    if isa(A,AbstractGPUArray)
        gpufilter!(A,responsetype,designmethod)
        return nothing
    end

    # filter if on the CPU
    if zerophase
        A[:,:] .= filtfilt(digitalfilter(responsetype, designmethod), @view(A[:,:]))
    else
        A[:,:] .= filt(digitalfilter(responsetype, designmethod), @view(A[:,:]))
    end
    return nothing
end
lowpass(A::AbstractArray{<:AbstractFloat},freq::Real, fs::Real;
      corners::Int=4,zerophase::Bool=true) = (U = deepcopy(A);
      lowpass!(U,freq,fs,corners=corners,zerophase=zerophase);return U)
lowpass!(R::RawData,freq::Real; corners::Int=4,
       zerophase::Bool=true) = (lowpass!(R.x,freq,R.fs,corners=corners,
       zerophase=zerophase);setfield!(R,:freqmax,Float64(min(freq,R.fs/2)));
       return nothing)
lowpass(R::RawData,freq::Real; corners::Int=4,
       zerophase::Bool=true) = (U = deepcopy(R);lowpass!(U.x,freq,U.fs,
       corners=corners,zerophase=zerophase);
       setfield!(U,:freqmax,Float64(min(freq,U.fs/2)));return U)
lowpass!(C::CorrData,freq::Real; corners::Int=4,
      zerophase::Bool=true) = (lowpass!(C.corr,freq,C.fs,corners=corners,
      zerophase=zerophase);setfield!(C,:freqmax,Float64(min(freq,C.fs/2)));
      return nothing)
lowpass(C::CorrData,freq::Real; corners::Int=4,
      zerophase::Bool=true) = (U = deepcopy(C);lowpass!(U.corr,freq,U.fs,
      corners=corners,zerophase=zerophase);
      setfield!(U,:freqmax,Float64(min(freq,U.fs/2)));return U)
lowpass!(C::SeisChannel, freq::Real;corners::Int=4, zerophase::Bool=true) =
    filtfilt!(C,fh=Float64(freq),np=corners,rt="Lowpass")
lowpass(C::SeisChannel,freq::Real; corners::Int=4,zerophase::Bool=true) =
    filtfilt(C,fh=Float64(freq),np=corners,rt="Lowpass")
lowpass!(S::SeisData, freq::Real;corners::Int=4, zerophase::Bool=true) =
    filtfilt!(S,fh=Float64(freq),np=corners,rt="Lowpass")
lowpass(S::SeisData,freq::Real; corners::Int=4,zerophase::Bool=true) =
    filtfilt(S,fh=Float64(freq),np=corners,rt="Lowpass")
lowpass!(N::NodalData,freq::Real; corners::Int=4,
    zerophase::Bool=true) = (lowpass!(N.data,freq,N.fs[1],corners=corners,
    zerophase=zerophase);return nothing)
lowpass(N::NodalData,freq::Real; corners::Int=4,
    zerophase::Bool=true) = (U = deepcopy(N);lowpass!(N.data,freq,U.fs[1],
    corners=corners,zerophase=zerophase);return U)

"""
   highpass(A,freq,fs,corners=4,zerophase=true)

Butterworth-Highpass Filter.

Filter data `A` removing data below certain frequency `freq` using `corners` corners.

# Arguments
- `A::AbstractArray`: Data to filter
- `freq::Real`: Filter corner frequency.
- `fs::Real`: Sampling rate in Hz.
- `fs::Int`: Filter corners / order.
- `zerophase::Bool`: If True, apply filter once forwards and once backwards.
This results in twice the filter order but zero phase shift in
the resulting filtered trace.
"""
function highpass!(A::AbstractArray{<:AbstractFloat},freq::Real,
                  fs::Real; corners::Int=4, zerophase::Bool=true)
    fe = 0.5 * fs
    f = freq / fe
    T = eltype(A)

    # warn if above Nyquist frequency
    if f >= 1
        throw(ArgumentError("Selected low corner frequency $freq Hz is above Nyquist $(fs/2) Hz."))
    end

    # create filter
    responsetype = Highpass(T(freq); fs=fs)
    designmethod = Butterworth(T,corners)

    # use gpu-specific kernel if on the GPU
    if isa(A,AbstractGPUArray)
        gpufilter!(A,responsetype,designmethod)
        return nothing
    end

    # filter if on the CPU
    if zerophase
        A[:,:] .= filtfilt(digitalfilter(responsetype, designmethod), @view(A[:,:]))
    else
        A[:,:] .= filt(digitalfilter(responsetype, designmethod), @view(A[:,:]))
    end
    return nothing
end
highpass(A::AbstractArray{<:AbstractFloat},freq::Real,fs::Real;
        corners::Int=4,zerophase::Bool=true) = (U = deepcopy(A);
        highpass!(U,freq,fs,corners=corners,zerophase=zerophase);return U)
highpass!(R::RawData,freq::Real; corners::Int=4,
        zerophase::Bool=true) = (highpass!(R.x,freq,R.fs,corners=corners,
        zerophase=zerophase);setfield!(R,:freqmin,Float64(freq));return nothing)
highpass(R::RawData,freq::Real; corners::Int=4,
        zerophase::Bool=true) = (U = deepcopy(R);highpass!(U.x,freq,U.fs,
        corners=corners,zerophase=zerophase);setfield!(U,:freqmin,Float64(freq));
        return U)
highpass!(C::CorrData,freq::Real; corners::Int=4,
        zerophase::Bool=true) = (highpass!(C.corr,freq,C.fs,corners=corners,
        zerophase=zerophase);setfield!(C,:freqmin,Float64(freq));return nothing)
highpass(C::CorrData,freq::Real; corners::Int=4,
        zerophase::Bool=true) = (U = deepcopy(C);highpass!(U.corr,freq,U.fs,
        corners=corners,zerophase=zerophase);setfield!(U,:freqmin,Float64(freq));
        return U)
highpass!(C::SeisChannel, freq::Real;corners::Int=4, zerophase::Bool=true) =
    filtfilt!(C,fl=Float64(freq),np=corners,rt="Highpass")
highpass(C::SeisChannel,freq::Real; corners::Int=4,zerophase::Bool=true) =
    filtfilt(C,fl=Float64(freq),np=corners,rt="Highpass")
highpass!(S::SeisData, freq::Real;corners::Int=4, zerophase::Bool=true) =
    filtfilt!(S,fl=Float64(freq),np=corners,rt="Highpass")
highpass(S::SeisData,freq::Real; corners::Int=4,zerophase::Bool=true) =
    filtfilt(S,fl=Float64(freq),np=corners,rt="Highpass")
highpass!(N::NodalData,freq::Real; corners::Int=4,
    zerophase::Bool=true) = (highpass!(N.data,freq,N.fs[1],corners=corners,
    zerophase=zerophase);return nothing)
highpass(N::NodalData,freq::Real; corners::Int=4,
    zerophase::Bool=true) = (U = deepcopy(N);highpass!(N.data,freq,U.fs[1],
    corners=corners,zerophase=zerophase);return U)

"""
    envelope(A)

Envelope of a function.

Computes the upper and lower envelopes of the given function.

# Arguments
- `A::AbstractArray`: Data to make envelope of.
"""
function envelope(A::AbstractArray)
    Amean = mean(A,dims=1)
    Acentered = A .- Amean
    env = abs.(hilbert(Acentered))
    upper = env .+ Amean
    lower = -env .+ Amean
    return upper, lower
end

"""
    gpufilter!(A,responsetype,designmethod)

Apply filter to array `A` on the GPU.

# Arguments
- `A::AbstractGPUArray`: Array on the GPU to filter
- `responsetype::FilterType`: DSP.jl filter representation
- `designmethod::ZeroPoleGain`: Filter representation in terms of zeros `z`, poles `p`, and
gain `k`.
"""
function gpufilter!(A::AbstractGPUArray,responsetype::FilterType,designmethod::ZeroPoleGain)
    T = Float64
    N = size(A,1)
    tf = convert(PolynomialRatio, digitalfilter(responsetype, designmethod))

    # get filter coefficients in Float64 to avoid divide by zero
    b = coefb(tf)
    a = coefa(tf)
    newb = CUDA.zeros(T,N)
    newa = CUDA.zeros(T,N)
    copyto!(newb,b)
    copyto!(newa,a)
    bafft = rfft(newb) ./ rfft(newa)

    # convert back to Float32 for speed
    bafft = map(x->convert(ComplexF32,x), bafft)

    # apply filter
    A .= irfft(rfft(A,1) .* bafft, N, 1)
    reverse!(A,dims=1)
    A .= irfft(rfft(A,1) .* bafft, N, 1)
    reverse!(A,dims=1)
    return nothing
end

"""
    resample!(R, fs)

Resample data in `RawData` R to new sampling rate `fs`.

# Arguments
- `R::RawData`: RawData to resample.
- `fs::Real`: New sampling frequency.
"""
function resample!(R::RawData,fs::Real)
    @assert fs > 0 "New fs must be greater than 0!"
    if R.fs == fs
        return nothing
    end
    T = eltype(R.x)
    rate = T(fs / R.fs)
    R.x = resample_kernel(R.x, rate)
    R.fs = Float64(fs)
    if (R.freqmax <  fs / 2) && (rate < 1.0)
        R.freqmax = Float64(fs / 2)
    end
    return nothing
end
resample(R::RawData,fs::Real) = (U = deepcopy(R); resample!(U,fs); return U)

"""
    resample!(C, fs)

Resample data in CorrData` C to new sampling rate `fs`.

# Arguments
- `C::CorrData`: CorrData to resample.
- `fs::Real`: New sampling frequency.
"""
function resample!(C::CorrData,fs::Real)
    @assert fs > 0 "New fs must be greater than 0!"
    if C.fs == fs
        return nothing
    end
    T = eltype(C.corr)
    rate = T(fs / C.fs)
    C.corr = resample_kernel(C.corr, rate)
    C.fs = Float64(fs)
    if (C.freqmax <  fs / 2) && (rate < 1.0)
        C.freqmax = Float64(fs / 2)
    end
    return nothing
end
resample(C::CorrData,fs::Real) = (U = deepcopy(C); resample!(U,fs); return U)

function resample_kernel(x::AbstractMatrix,rate::Real)
    T = eltype(x)
    h = resample_filter(rate)
    self = FIRFilter(h, rate)

    # Get delay, in # of samples at the output rate, caused by filtering processes
    τ = timedelay(self)

    # Use setphase! to
    #   a) adjust the input samples to skip over before producing and output (integer part of τ)
    #   b) set the ϕ index of the PFB (fractional part of τ)
    setphase!(self, τ)

    # calculate number of zeros
    Nrows, Ncols = size(x)
    outLen = ceil(Int, Nrows*rate)
    reqInlen = inputlength(self, outLen)
    reqZerosLen = reqInlen - Nrows
    xPadded     = [x; zeros(T, reqZerosLen, Ncols)]

    # fill empty matrix
    out = zeros(T,outLen,Ncols)
    xout = filt(deepcopy(self), xPadded[:,1])
    xlength = length(xout)
    out[1:xlength,1] .= xout
    @inbounds for ii = 2:Ncols
        out[1:xlength,ii] .= filt(deepcopy(self), xPadded[:,ii])
    end
    return out
end


text:
using HashFunctions, Test

# from http://www.isthe.com/chongo/src/fnv/test_fnv.c
@testset "FNV1a 32" begin
    hf = FNV1A(32)
    @test calc(hf, "")        == 0x811c9dc5
    @test calc(hf, "a")       == 0xe40c292c
    @test calc(hf, "b")       == 0xe70c2de5
    @test calc(hf, "c")       == 0xe60c2c52
    @test calc(hf, "d")       == 0xe10c2473
    @test calc(hf, "e")       == 0xe00c22e0
    @test calc(hf, "f")       == 0xe30c2799
    @test calc(hf, "fo")      == 0x6222e842
    @test calc(hf, "foo")     == 0xa9f37ed7
    @test calc(hf, "foob")    == 0x3f5076ef
    @test calc(hf, "fooba")   == 0x39aaa18a
    @test calc(hf, "foobar")  == 0xbf9cf968
end

# from http://www.isthe.com/chongo/src/fnv/test_fnv.c
@testset "FNV1 32" begin
    hf = FNV1(32)
    @test calc(hf, "")        == 0x811c9dc5
    @test calc(hf, "a")       == 0x050c5d7e
    @test calc(hf, "b")       == 0x050c5d7d
    @test calc(hf, "c")       == 0x050c5d7c
    @test calc(hf, "d")       == 0x050c5d7b
    @test calc(hf, "e")       == 0x050c5d7a
    @test calc(hf, "f")       == 0x050c5d79
    @test calc(hf, "fo")      == 0x6b772514
    @test calc(hf, "foo")     == 0x408f5e13
    @test calc(hf, "foob")    == 0xb4b1178b
    @test calc(hf, "fooba")   == 0xfdc80fb0
    @test calc(hf, "foobar")  == 0x31f0b262
end

@testset "calcwithoffset equals calc on full string" begin
    for HT in [FNV1, FNV1A]
        hf = HT(32)
        @test calc(hf, "ba", 1, 1) == calc(hf, "a")
        @test calc(hf, "ab", 0, 1) == calc(hf, "a")
        @test calc(hf, "foobar", 0, 1) == calc(hf, "f")
        @test calc(hf, "foobar", 1, 3) == calc(hf, "oob")
        @test calc(hf, "foobar", 2, 2) == calc(hf, "ob")
        @test calc(hf, "foobar", 3, 2) == calc(hf, "ba")
        @test calc(hf, "foobar", 4, 2) == calc(hf, "ar")
        @test_throws AssertionError calc(hf, "foobar", 4, 3)
    end
end


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

function GitIndex(repo::GitRepo)
    idx_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_repository_index, :libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}), idx_ptr_ptr, repo.ptr)
    return GitIndex(repo, idx_ptr_ptr[])
end

function read!(idx::GitIndex, force::Bool = false)
    @check ccall((:git_index_read, :libgit2), Cint, (Ptr{Void}, Cint), idx.ptr, Cint(force))
    return idx
end

function write!(idx::GitIndex)
    @check ccall((:git_index_write, :libgit2), Cint, (Ptr{Void},), idx.ptr)
    return idx
end

function write_tree!(idx::GitIndex)
    oid_ptr = Ref(GitHash())
    @check ccall((:git_index_write_tree, :libgit2), Cint,
                 (Ptr{GitHash}, Ptr{Void}), oid_ptr, idx.ptr)
    return oid_ptr[]
end

function repository(idx::GitIndex)
    if isnull(idx.nrepo)
        throw(GitError(Error.Index, Error.ENOTFOUND, "Index does not have an owning repository."))
    else
        return Base.get(idx.nrepo)
    end
end

"""
    LibGit2.read_tree!(idx::GitIndex, tree::GitTree)
    LibGit2.read_tree!(idx::GitIndex, treehash::AbstractGitHash)

Read the tree `tree` (or the tree pointed to by `treehash` in the repository owned by
`idx`) into the index `idx`. The current index contents will be replaced.
"""
function read_tree!(idx::GitIndex, tree::GitTree)
    @check ccall((:git_index_read_tree, :libgit2), Cint,
                 (Ptr{Void}, Ptr{Void}), idx.ptr, tree.ptr)
end
read_tree!(idx::GitIndex, hash::AbstractGitHash) =
    read_tree!(idx, GitTree(repository(idx), hash))

function add!{T<:AbstractString}(idx::GitIndex, files::T...;
             flags::Cuint = Consts.INDEX_ADD_DEFAULT)
    @check ccall((:git_index_add_all, :libgit2), Cint,
                 (Ptr{Void}, Ptr{StrArrayStruct}, Cuint, Ptr{Void}, Ptr{Void}),
                 idx.ptr, collect(files), flags, C_NULL, C_NULL)
end

function update!{T<:AbstractString}(idx::GitIndex, files::T...)
    @check ccall((:git_index_update_all, :libgit2), Cint,
                 (Ptr{Void}, Ptr{StrArrayStruct}, Ptr{Void}, Ptr{Void}),
                 idx.ptr, collect(files), C_NULL, C_NULL)
end

function remove!{T<:AbstractString}(idx::GitIndex, files::T...)
    @check ccall((:git_index_remove_all, :libgit2), Cint,
                 (Ptr{Void}, Ptr{StrArrayStruct}, Ptr{Void}, Ptr{Void}),
                 idx.ptr, collect(files), C_NULL, C_NULL)
end

function add!{T<:AbstractString}(repo::GitRepo, files::T...;
             flags::Cuint = Consts.INDEX_ADD_DEFAULT)
    with(GitIndex, repo) do idx
        add!(idx, files..., flags = flags)
        write!(idx)
    end
    return
end

function update!{T<:AbstractString}(repo::GitRepo, files::T...)
    with(GitIndex, repo) do idx
        update!(idx, files...)
        write!(idx)
    end
    return
end

function remove!{T<:AbstractString}(repo::GitRepo, files::T...)
    with(GitIndex, repo) do idx
        remove!(idx, files...)
        write!(idx)
    end
    return
end

function read!(repo::GitRepo, force::Bool = false)
    with(GitIndex, repo) do idx
        read!(idx, force)
    end
    return
end

function Base.count(idx::GitIndex)
    return ccall((:git_index_entrycount, :libgit2), Csize_t, (Ptr{Void},), idx.ptr)
end

function Base.getindex(idx::GitIndex, i::Integer)
    ie_ptr = ccall((:git_index_get_byindex, :libgit2),
                   Ptr{IndexEntry},
                   (Ptr{Void}, Csize_t), idx.ptr, i-1)
    ie_ptr == C_NULL && return nothing
    return unsafe_load(ie_ptr)
end

function Base.find(path::String, idx::GitIndex)
    pos_ref = Ref{Csize_t}(0)
    ret = ccall((:git_index_find, :libgit2), Cint,
                  (Ref{Csize_t}, Ptr{Void}, Cstring), pos_ref, idx.ptr, path)
    ret == Cint(Error.ENOTFOUND) && return Nullable{Csize_t}()
    return Nullable(pos_ref[]+1)
end

stage(ie::IndexEntry) = ccall((:git_index_entry_stage, :libgit2), Cint, (Ptr{IndexEntry},), Ref(ie))

function Base.show(io::IO, idx::GitIndex)
    println(io, "GitIndex:\nRepository: ", repository(idx), "\nNumber of elements: ", count(idx))
end


text:
function Pump(; name, P = 1.0)
    @named oneport = OnePort()
    @unpack p = oneport
    ps = @parameters P=P
    eqs = [
        P ~ p
    ]
    extend(ODESystem(eqs, t, [], ps; name), oneport)
end


text:
mutable struct ImplicitRHS_Scalar{F,uType,tType,P} <: Function
  f::F
  tmp::uType
  a::tType
  t::tType
  dt::tType
  p::P
end

function (p::ImplicitRHS_Scalar)(resid,u)
  resid[1] = first(u) - p.tmp - p.a*first(p.f(first(u),p.p,p.t+p.dt))
end

mutable struct ImplicitRHS{F,uType,tType,DiffCacheType,P} <: Function
  f::F
  tmp::uType
  a::tType
  t::tType
  dt::tType
  dual_cache::DiffCacheType
  p::P
end

function (p::ImplicitRHS)(resid,u)
  _du1 = get_du(p.dual_cache, eltype(u))
  du1 = reinterpret(eltype(u),_du1)
  p.f(du1,u,p.p,p.t+p.dt)
  @.. resid = u - p.tmp - p.a*du1
end

function initialize!(integrator,
                     cache::Union{GenericImplicitEulerConstantCache,GenericTrapezoidConstantCache})
  cache.uhold[1] = integrator.uprev
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
  integrator.destats.nf += 1

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::GenericImplicitEulerConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack uhold,rhs,nl_rhs = cache
  alg = typeof(integrator.alg) <: CompositeAlgorithm ? integrator.alg.algs[integrator.cache.current] : integrator.alg
  rhs.tmp = uprev

  if integrator.success_iter > 0 && !integrator.reeval_fsal && alg.extrapolant == :interpolant
    uhold[1] = current_extrapolant(t+dt,integrator)
  elseif alg.extrapolant == :linear
    uhold[1] = uprev + dt*integrator.fsalfirst
  else # :constant
    uhold[1] = uprev
  end

  rhs.t = t
  rhs.dt = dt
  rhs.a = dt
  nlres = alg.nlsolve(nl_rhs,uhold)
  uhold[1] = nlres[1]
  integrator.fsallast = f(uhold[1],p,t+dt)
  integrator.destats.nf += 1
  u = uhold[1]

  if integrator.opts.adaptive && integrator.success_iter > 0
    # local truncation error (LTE) bound by dt^2/2*max|y''(t)|
    # use 2nd divided differences (DD) a la SPICE and Shampine

    # TODO: check numerical stability
    uprev2 = integrator.uprev2
    tprev = integrator.tprev

    dt1 = dt*(t+dt-tprev)
    dt2 = (t-tprev)*(t+dt-tprev)
    c = 7/12 # default correction factor in SPICE (LTE overestimated by DD)
    r = c*dt^2 # by mean value theorem 2nd DD equals y''(s)/2 for some s

    tmp = r*abs((u - uprev)/dt1 - (uprev - uprev2)/dt2)
    atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  else
    integrator.EEst = 1
  end

  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator,
                     cache::Union{GenericImplicitEulerCache,GenericTrapezoidCache})
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t)
  integrator.destats.nf += 1

  integrator.kshortsize = 2
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::GenericImplicitEulerCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack dual_cache,k,nl_rhs,rhs,tmp,atmp = cache
  alg = typeof(integrator.alg) <: CompositeAlgorithm ? integrator.alg.algs[integrator.cache.current] : integrator.alg
  copyto!(tmp,uprev)

  if integrator.success_iter > 0 && !integrator.reeval_fsal && alg.extrapolant == :interpolant
    current_extrapolant!(u,t+dt,integrator)
  elseif alg.extrapolant == :linear
    @.. u = uprev + dt*integrator.fsalfirst
  else
    copyto!(u,uprev)
  end

  rhs.t = t
  rhs.dt = dt
  rhs.a = dt
  nlres = alg.nlsolve(nl_rhs,u)
  copyto!(u,nlres)

  if integrator.opts.adaptive && integrator.success_iter > 0
    # local truncation error (LTE) bound by dt^2/2*max|y''(t)|
    # use 2nd divided differences (DD) a la SPICE and Shampine

    # TODO: check numerical stability
    uprev2 = integrator.uprev2
    tprev = integrator.tprev

    dt1 = dt*(t+dt-tprev)
    dt2 = (t-tprev)*(t+dt-tprev)
    c = 7/12 # default correction factor in SPICE (LTE overestimated by DD)
    r = c*dt^2 # by mean value theorem 2nd DD equals y''(s)/2 for some s

    @.. tmp = r*abs((u - uprev)/dt1 - (uprev - uprev2)/dt2)
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  else
    integrator.EEst = 1
  end

  f(k, u, p, t+dt)
  integrator.destats.nf += 1
end

function initialize!(integrator, cache::GenericTrapezoidConstantCache)
  cache.uhold[1] = integrator.uprev
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
  integrator.destats.nf += 1
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::GenericTrapezoidConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack uhold,rhs,nl_rhs = cache
  alg = typeof(integrator.alg) <: CompositeAlgorithm ? integrator.alg.algs[integrator.cache.current] : integrator.alg
  rhs.tmp = first(uprev) + (dt/2)*first(integrator.fsalfirst)

  if integrator.success_iter > 0 && !integrator.reeval_fsal && alg.extrapolant == :interpolant
    uhold[1] = current_extrapolant(t+dt,integrator)
  elseif alg.extrapolant == :linear
    uhold[1] = uprev + dt*integrator.fsalfirst
  else # :constant
    uhold[1] = uprev
  end

  rhs.t = t
  rhs.dt = dt
  rhs.a = dt/2
  nlres = alg.nlsolve(nl_rhs,uhold)
  uhold[1] = nlres[1]
  integrator.fsallast = f(uhold[1],p,t+dt)
  integrator.destats.nf += 1
  u = uhold[1]

  if integrator.opts.adaptive
    if integrator.iter > 2
      # local truncation error (LTE) bound by dt^3/12*max|y'''(t)|
      # use 3rd divided differences (DD) a la SPICE and Shampine

      # TODO: check numerical stability
      uprev2 = integrator.uprev2
      tprev = integrator.tprev
      uprev3 = cache.uprev3
      tprev2 = cache.tprev2

      dt1 = dt*(t+dt-tprev)
      dt2 = (t-tprev)*(t+dt-tprev)
      dt3 = (t-tprev)*(t-tprev2)
      dt4 = (tprev-tprev2)*(t-tprev2)
      dt5 = t+dt-tprev2
      c = 7/12 # default correction factor in SPICE (LTE overestimated by DD)
      r = c*dt^3/2 # by mean value theorem 3rd DD equals y'''(s)/6 for some s

      # tmp = r*abs(((u - uprev)/dt1 - (uprev - uprev2)/dt2) - ((uprev - uprev2)/dt3 - (uprev2 - uprev3)/dt4)/dt5)
      DD31 = (u - uprev)/dt1 - (uprev - uprev2)/dt2
      DD30 = (uprev - uprev2)/dt3 - (uprev2 - uprev3)/dt4
      tmp = r*abs((DD31 - DD30)/dt5)
      atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
      if integrator.EEst <= 1
        cache.uprev3 = uprev2
        cache.tprev2 = tprev
      end
    elseif integrator.success_iter > 0
      integrator.EEst = 1
      cache.uprev3 = integrator.uprev2
      cache.tprev2 = integrator.tprev
    else
      integrator.EEst = 1
    end
  end

  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator, cache::GenericTrapezoidCache)
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t)
  integrator.destats.nf += 1
  integrator.kshortsize = 2
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::GenericTrapezoidCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack dual_cache,k,rhs,nl_rhs,tmp,atmp = cache
  alg = typeof(integrator.alg) <: CompositeAlgorithm ? integrator.alg.algs[integrator.cache.current] : integrator.alg
  tmp .= uprev .+ (dt/2).*integrator.fsalfirst

  if integrator.success_iter > 0 && !integrator.reeval_fsal && alg.extrapolant == :interpolant
    current_extrapolant!(u,t+dt,integrator)
  elseif alg.extrapolant == :linear
    @.. u = uprev + dt*integrator.fsalfirst
  else
    copyto!(u,uprev)
  end

  # copyto!(rhs.fsalfirst,fsalfirst) Implicitly done by pointers: fsalfirst === fsalfirst == rhs.fsalfirst
  rhs.t = t
  rhs.dt = dt
  rhs.a = dt/2
  nlres = alg.nlsolve(nl_rhs,u)
  copyto!(u,nlres)

  if integrator.opts.adaptive
    if integrator.iter > 2
      # local truncation error (LTE) bound by dt^3/12*max|y'''(t)|
      # use 3rd divided differences (DD) a la SPICE and Shampine

      # TODO: check numerical stability
      uprev2 = integrator.uprev2
      tprev = integrator.tprev
      uprev3 = cache.uprev3
      tprev2 = cache.tprev2

      dt1 = dt*(t+dt-tprev)
      dt2 = (t-tprev)*(t+dt-tprev)
      dt3 = (t-tprev)*(t-tprev2)
      dt4 = (tprev-tprev2)*(t-tprev2)
      dt5 = t+dt-tprev2
      c = 7/12 # default correction factor in SPICE (LTE overestimated by DD)
      r = c*dt^3/2 # by mean value theorem 3rd DD equals y'''(s)/6 for some s

      # @.. tmp = r*abs(((u - uprev)/dt1 - (uprev - uprev2)/dt2) - ((uprev - uprev2)/dt3 - (uprev2 - uprev3)/dt4)/dt5)
      @inbounds for i in eachindex(u)
        DD31 = (u[i] - uprev[i])/dt1 - (uprev[i] - uprev2[i])/dt2
        DD30 = (uprev[i] - uprev2[i])/dt3 - (uprev2[i] - uprev3[i])/dt4
        tmp[i] = r*abs((DD31 - DD30)/dt5)
      end
      calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
      if integrator.EEst <= 1
        copyto!(cache.uprev3,uprev2)
        cache.tprev2 = tprev
      end
    elseif integrator.success_iter > 0
      integrator.EEst = 1
      copyto!(cache.uprev3,integrator.uprev2)
      cache.tprev2 = integrator.tprev
    else
      integrator.EEst = 1
    end
  end

  f(k, u, p, t+dt)
  integrator.destats.nf += 1
end


text:
"""
Altro.jl
"""

module Altro

using BenchmarkTools
using Crayons
using DocStringExtensions
using ForwardDiff
using Interpolations
using LinearAlgebra
using Logging
using Parameters
using RecipesBase
using SparseArrays
using SolverLogging
using StaticArrays
using Statistics
using UnsafeArrays

include("utils.jl")
include("model.jl")
include("integration.jl")
include("costfunctions.jl")
include("objective.jl")
include("constraints.jl")
include("constraint_list.jl")
include("convals.jl")
include("problem.jl")
include("solvers.jl")
include("solver_opts.jl")
include("ilqr/ilqr.jl")
include("ilqr/backwardpass.jl")
include("ilqr/rollout.jl")
include("ilqr/ilqr_methods.jl")
include("al/al.jl")
include("al/al_methods.jl")
include("al/al_objective.jl")
include("pn/pn.jl")
include("pn/pn_methods.jl")
include("altro/altro.jl")

export
    # solver
    ALTROSolver,
    SolverStats, SolverOptions,
    solve!, benchmark_solve!, states, controls,
    # constraints
    ConstraintList, add_constraint!, ConstraintParams, STATE, CONTROL, EQUALITY, INEQUALITY,
    GoalConstraint, BoundConstraint, NormConstraint,
    # problem
    Problem,
    # model
    AbstractModel,
    Explicit, RK4, RK3, RK2, Euler,
    dynamics, discrete_dynamics, discrete_dynamics!, discrete_jacobian!,
    # cost functions
    Objective, LQRObjective
end # module


text:
# This file was generated, do not modify it. # hide
LDA = @load LDA pkg=MultivariateStats
using Distances

classif = machine(LDA(dist=CosineDist()), X3, y)
fit!(classif, rows=train)
ŷ = predict_mode(classif, rows=test)

accuracy(ŷ, y[test]) |> r3

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/MortarContact2D.jl/blob/master/LICENSE

# Construct discrete coupling operator between two elements and eliminate
# boundary conditions from global matrix assembly. Thesis, p. 73.

using FEMBase, MortarContact2D, Test

X = Dict(
    1 => [1.0, 2.0],
    2 => [3.0, 2.0],
    3 => [2.0, 2.0],
    4 => [0.0, 2.0],
    5 => [3.0, 4.0],
    6 => [1.0, 4.0],
    7 => [0.0, 0.0],
    8 => [2.0, 0.0])

slave = Element(Seg2, [1, 2])
master = Element(Seg2, [3, 4])
elements = [slave, master]
update!(elements, "geometry", X)

problem = Problem(Mortar2D, "test interface", 2, "displacement")
add_slave_elements!(problem, [slave])
add_master_elements!(problem, [master])
assemble!(problem, 0.0)

K_SS = K_MM = [
 180.0   81.0  -126.0   27.0
  81.0  180.0   -27.0   36.0
-126.0  -27.0   180.0  -81.0
  27.0   36.0   -81.0  180.0]
M_SS = M_MM = [
 40.0   0.0  20.0   0.0
  0.0  40.0   0.0  20.0
 20.0   0.0  40.0   0.0
  0.0  20.0   0.0  40.0]
f_S = [-27.0, -81.0, 81.0, -135.0]
f_M = [ 0.0, 0.0, 0.0, 0.0]
K = [K_SS zeros(4,4); zeros(4,4) K_MM]
M = [M_SS zeros(4,4); zeros(4,4) M_MM]
f = [f_S; f_M]

eliminate_boundary_conditions!(problem, K, M, f)
K_expected = [
  441.0  -81.0  -279.0   81.0
  -81.0  684.0    81.0  -36.0
 -279.0   81.0   333.0  -81.0
   81.0  -36.0   -81.0  252.0]
f_expected = [108.0, -243.0, -54.0, 27.0]

@test isapprox(K[5:8,5:8], K_expected)
@test isapprox(f[5:8], f_expected)


text:
using ApproxFunOrthogonalPolynomials, SpecialFunctions, LinearAlgebra, Test
import DomainSets
import DomainSets: UnionDomain
import ApproxFunBase: space, SpaceOperator,
                    testspace, testbandedoperator, testraggedbelowoperator, testcalculus, testtransforms,
                    testfunctional

@testset "Spaces" begin
    @testset "evaluate" begin
        @test evaluate([1],Legendre(),0)==evaluate([1.0],Legendre(),0)
    end
    @testset "ChebyshevDirichlet" begin
        testtransforms(ChebyshevDirichlet{1,1}())

        @test Fun(exp,ChebyshevDirichlet{1,1})(.1) ≈ exp(.1)
        @test Fun(Fun(exp,ChebyshevDirichlet{1,1}),Ultraspherical(1))(.1) ≈ exp(.1)

        @test domain(ChebyshevDirichlet{1,1}()) === ChebyshevInterval()
        @test Fun(Fun(exp,ChebyshevDirichlet{1,1}(Segment(-1,1))),Ultraspherical(1))(.1) ≈ exp(.1)

        d = ChebyshevInterval()
        sp = ChebyshevDirichlet{1,1}(d)
        B=Dirichlet(sp)
        D=Derivative(d)
        L=D^2+I
        @test B[1:2,1:4] ≈ [1 -1 0 0; 1 1 0 0]

        @test csc(2)sin(1 - 0.1)  ≈ ([Dirichlet(d);L]\[[1.,0.],0.])(0.1)
        @test csc(2)sin(1 - 0.1)  ≈ ([B;L]\[[1.,0.],0.])(0.1)

        @test norm(([B;L]\[[1.,0],0])-([Dirichlet(d);L]\[[1.,0],0])) <10eps()
    end

    @testset "PiecewiseSpace" begin
        x = Fun(identity,UnionDomain(-1..0, 0..1))
        sp=space(x)
        testtransforms(sp;minpoints=2)

        D = Derivative(sp)
        testbandedoperator(D)

        B = [Dirichlet(sp); continuity(sp,0:1)]
        u = [B; D^2] \ [[1,0],zeros(2),0];
        @time u2 = [Dirichlet();Derivative(Chebyshev())^2] \ [[1.,0],0]
        @test u(0.) ≈ u2(0.)

        x = Fun(identity, UnionDomain(-10..0, 0..1, 1..15))
        sp=space(x)
        D=Derivative(sp)
        B=Dirichlet(sp)

        @time  u=[B;
            continuity(sp,0:1);
            D^2-x]\[[airyai(-10.),0],zeros(4),0];

        @test u(0.) ≈ airyai(0.)

        s=Fun(sin,-2..2)|>abs
        c=Fun(cos,-2..2)|>abs
        sc=Fun(x -> abs(sin(x))+abs(cos(x)), ContinuousSpace(PiecewiseSegment([-2,(-π/2),0,π/2,2])))
        @test norm(sc-(c+s))<100eps()
    end

    @testset "max/min creates breakpoints" begin
        x = Fun()
        g = 4*(x-0.2)
        f = max(-1,g)
        f2 = min(f,1)

        @test norm(max(x,x)-x)<100eps()
        @test norm(min(x,x)-x)<100eps()

        f3=Fun(x->x<-0.05 ? -1.0 : (x<0.45 ? 4*(x-.2) : 1), (-1..1) \ [-0.05,0.45])
        @test norm(f2.(range(-1,stop=1,length=10))-f3.(range(-1,stop=1,length=10))) < 4eps()
    end

    @testset "Complex piecewise" begin
        x = Fun(identity, Segment(im,0) ∪ Segment(0,1))
        @test x(0.5) ≈ 0.5
        @test x(0.5im) ≈ 0.5im

        @test Fun(Fun(1.0),space(x))(0.5) == 1.0
        @test Fun(Fun(1.0),space(x))(0.5im) == 1.0

        @test (x+1)(0.5) ≈ 1.5
        @test (x-1)(0.5) ≈ -0.5
        @test (1-x)(0.5) ≈ 0.5
    end

    @testset "Equivalent spaces" begin
        @test norm(Fun(cos,Chebyshev)-Fun(cos,Jacobi(-0.5,-0.5)))<100eps()
    end

    @testset "ContinuousSpace" begin
        d = PiecewiseSegment(1.,2.,3.,4.)
        S = ContinuousSpace(d)
        testtransforms(S; minpoints=3, invertibletransform=false)

        D = Derivative(S)
        testbandedoperator(D)

        A = [ldirichlet(S); D-I]
        testraggedbelowoperator(A)
        @time u = [ldirichlet(S); D-I] \ [exp(1.); 0]

        @test u(1.1) ≈ exp(1.1)
        @test u(3.4) ≈ exp(3.4)
        @test last(u) ≈ exp(4)

        d=PiecewiseSegment(0,1.,1+im,im,0.)
        s=Space(d)

        # narrow down bug
        @test s isa ContinuousSpace
        @test ApproxFunBase.rangetype(s) == Float64
        cs=ApproxFunBase.canonicalspace(s)
        @test ApproxFunBase.rangetype(cs) == Float64

        @test conversion_type(s,cs) == s

        @test Fun(exp,d)(.1) ≈ exp(.1)
    end

    @testset "Jacobi recurrence bug" begin
        S=Jacobi(0.5,-0.5)
        f=Fun(exp,S)
        @test f(0.1) ≈ exp(0.1)
    end

    @testset "Jacobi evaluation bug" begin
        S=Jacobi(0.5,-0.5)
        f=Fun(exp,S)
        B = ldirichlet(S)
        testfunctional(B)
        @test B*f ≈ exp(-1)
    end

    @testset "Hermite" begin
        f=Fun(x->x+x^2,Hermite())
        @test f(1.) ≈ 2.
        @test values(f) ≈ points(f)+points(f).^2

        D = Derivative(Hermite())
        @time testbandedoperator(D)

        g = D*f
        @test g(1.) ≈ 3.
    end

    @testset "Arc exp" begin
        z=Fun(identity,Arc(0.,.1,0.,π/2))
        @test exp(z)(0.1exp(0.2im)) ≈ exp(0.1exp(0.2im))
    end

    @testset "Line" begin
        f=Fun(z->2exp(z^2),Line(0.,π/2))
        @test f(1.1im) ≈ 2exp(-1.1^2)
    end

    @testset "Exp Curve" begin
        a = 1+10*im; b = 2-6*im
        d = IntervalCurve(Fun(x->1+a*x+x^2+b*x^3))

        x=Fun(d)

        @test exp(x)(1+a*0.1+0.1^2+b*0.1^3) ≈ exp(1+a*0.1+0.1^2+b*0.1^3)
    end

    @testset "ChebyshevDirichlet multiplication" begin
        S=ChebyshevDirichlet()
        x=Fun()
        @test norm((ApproxFunOrthogonalPolynomials.Recurrence(S)*Fun(exp,S)-Fun(x->x*exp(x),S)).coefficients) < 100eps()
        @test norm((x*Fun(exp,S)-Fun(x->x*exp(x),S)).coefficients) < 100eps()

        @test Integral(S,1)*Fun(S,[1.,2.,3.]) ≈ integrate(Fun(Fun(S,[1.,2.,3.]),Chebyshev()))
    end

    @testset "QuotientSpace" begin
        S = Chebyshev()
        for B in (Dirichlet(S),
                    Neumann(S),
                    vcat(bvp(S,4)...),
                    [Dirichlet(S); Neumann(S); Evaluation(S, -1, 3); Evaluation(S, 1, 3)])
            QS = QuotientSpace(B)
            A = Conversion(QS, S)

            @test norm(B[1:size(B,1),1:10]*A[1:10,1:10-size(B,1)]) < 10000eps()
        end
    end

    @testset "Fast Ultraspherical*Chebyshev" begin
        f = Fun(exp, Ultraspherical(1))
        g = Fun(exp)

        @test ApproxFunBase.hasfasttransform(f)
        @test ApproxFunBase.pointscompatible(f,g)
        @test ApproxFunBase.hasfasttransformtimes(f,g)
        @test ApproxFunBase.transformtimes(f,g)(0.1) ≈ exp(0.2)
    end

    @testset "PointSpace" begin
        δ = Fun(ApproxFunBase.PointSpace([2.0]),[1.0])
        f = Fun(x->cos(50x)) + δ
        g = Fun(x->cos(50x),Ultraspherical(1)) + δ
        @test domain(f) == domain(g)

        @test (f*g)(0.1) ≈ cos(50*0.1)^2
        @test (f*g)(2.0) ≈ 1

        @test exp(f)(0.1) ≈ exp(cos(50*0.1))
        @test exp(f)(2.0) ≈ exp(1)

        @test iszero(sign(-Fun(zero)))
        @test iszero(sign(Fun(zero)))
        @test iszero(abs(-Fun(zero)))
        @test iszero(abs(Fun(zero)))
        @test angle(-Fun(zero)) ≈ Fun(π)
        @test iszero(angle(Fun(zero)))

        @test sign(f)(0.1) ≈ sign(cos(50*0.1))
        @test sign(f)(2.0) ≈ 1
        @test abs(f)(0.1) ≈ abs(cos(50*0.1))
        @test abs(f)(2.0) ≈ 1
        @test angle(f)(0.1) ≈ angle(cos(50*0.1))
        @test angle(f)(2.0) ≈ 0
    end

    @testset "Jump Locations" begin
        x = Fun(UnionDomain(0..1, 2..3))
        @test length(jumplocations(sign(x))) == 0

        x = Fun(Chebyshev(-1..1))
        @test length(jumplocations(x)) == 0
        @test all(jumplocations(sign(x) + sign(x+0.2)) .≈ [-0.2, 0])
    end

    @testset "SumSpace Conversion" begin
        H = ApproxFunBase.HeavisideSpace([-1.0,0.0,1.0])
        C = ContinuousSpace(ApproxFunBase.PiecewiseSegment([-1.0,0,1]))
        S = H + C
        P = Ultraspherical(1,Segment(-1.0,0.0)) ∪ Ultraspherical(1,Segment(0.0,1.0))
        f = Fun(S, randn(100))
        @test f(0.1) ≈ Fun(f, P)(0.1)

        @test Conversion(S,P)[1,1] == 1.0
        @test Conversion(S,P)[1,2] == 0.5
        @test rangespace(Conversion(S,P)) == Ultraspherical(1,Segment(-1.0,0.0)) ∪ Ultraspherical(1,Segment(0.0,1.0))
        @test domainspace(Conversion(S,P)) ==
        ApproxFunBase.SumSpace(ApproxFunBase.HeavisideSpace([-1.0,0.0,1.0]),
            ContinuousSpace(ApproxFunBase.PiecewiseSegment([-1.0,0,1])))

        D = ApproxFunBase.DiracSpace([-1.0,0.0,1.0])
        S2 = ApproxFunBase.SumSpace(D , P)
        f = Fun(S, randn(100))
        @test (Conversion(S,S2) * f)(0.1) ≈ f(0.1)
    end

    @testset "real line integral" begin
        f = x -> x^2/(x^4+1)

        z₁,z₂,z₃,z₄ = exp(im*π/4), exp(3im*π/4), exp(5im*π/4), exp(7im*π/4)

        res₁ = z₁^2 / ((z₁ - z₂)*(z₁ - z₃)*(z₁ - z₄) )
        res₂ = z₂^2 / ((z₂ - z₁)*(z₂ - z₃)*(z₂ - z₄) )

        @test_skip sum(Fun(f, Line())) ≈ 2π*im*(res₁ + res₂)
    end

    @testset "Piecewise Legendre" begin
        f = Fun(x -> exp(-40(x-0.1)^2), Legendre(0..1), 1000)
        @test f(0.2) ≈ exp(-40(0.1)^2)

        f = Fun(x -> exp(-40(x-0.1)^2), Legendre(0..1))
        @test f(0.2) ≈ exp(-40(0.1)^2)

        sp = Legendre(Segment(0 , -1)) ⊕ Legendre(0 .. 1)
        @time f = Fun(x->sign(x)*exp(-40(x-0.1)^2), sp)
        @test f(0.2) ≈ exp(-40(0.1)^2)
        @test f(-0.2) ≈ -exp(-40(0.3)^2)
    end

    @testset "Ambiguous evaluate" begin
        f = Fun(1.0, Chebyshev(NaN..NaN))
        @test f(0.1) ≈ 1.0
    end

    @testset "x + abs(x) (#637)" begin
        x = Fun()
        g = abs(x) + x
        @test ncoefficients(g) == 3
        @test g(0.1) ≈ 0.2
        @test g(-0.2) ≈ 0.0
    end

    @testset "remove point" begin
        f = Fun(identity, (-1..1) \ 0)
        @test f(0.1) ≈ 0.1
        @test f(-0.1) ≈ -0.1
    end
end


text:
module TestExamples2dP4est

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "p4est_2d_dgsem")

@testset "P4estMesh" begin
  @trixi_testset "elixir_advection_basic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"),
      l2   = [9.14468177884088e-6],
      linf = [6.437440532947036e-5])
  end

  @trixi_testset "elixir_advection_nonconforming.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonconforming.jl"),
      l2   = [2.58174252995893e-5],
      linf = [2.65115939055204e-4])
  end

  @trixi_testset "elixir_advection_nonconforming_unstructured.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonconforming_unstructured.jl"),
      l2   = [3.038384623519386e-3],
      linf = [6.324792487776842e-2])
  end

  @trixi_testset "elixir_advection_amr_solution_independent.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_solution_independent.jl"),
      l2   = [7.945073295691892e-5],
      linf = [0.0007454287896710293])
  end

  @trixi_testset "elixir_advection_amr_unstructured_flag.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_flag.jl"),
      l2   = [1.29561551e-03],
      linf = [3.92061383e-02])
  end

  @trixi_testset "elixir_advection_restart.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
      l2   = [6.398955192910044e-6],
      linf = [3.474337336717426e-5])
  end

  @trixi_testset "elixir_euler_source_terms_nonperiodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic.jl"),
      l2   = [2.3653424742684444e-6, 2.1388875095440695e-6, 2.1388875095548492e-6, 6.010896863397195e-6],
      linf = [1.4080465931654018e-5, 1.7579850587257084e-5, 1.7579850592586155e-5, 5.956893531156027e-5])
  end

  @trixi_testset "elixir_euler_source_terms_nonperiodic_unstructured.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic_unstructured.jl"),
      l2   = [0.005689496253354319, 0.004522481295923261, 0.004522481295922983, 0.009971628336802528],
      linf = [0.05125433503504517, 0.05343803272241532, 0.053438032722404216, 0.09032097668196482])
  end

  @trixi_testset "elixir_euler_free_stream.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"),
      l2   = [2.063350241405049e-15, 1.8571016296925367e-14, 3.1769447886391905e-14, 1.4104095258528071e-14],
      linf = [1.9539925233402755e-14, 2e-12, 4.8e-12, 4e-12],
      atol = 2.0e-12, # required to make CI tests pass on macOS
    )
  end

  @trixi_testset "elixir_eulergravity_convergence.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_convergence.jl"),
      l2   = [0.0002456507718735638, 0.0003374236900407321, 0.00033742369004074064, 0.0007198398465056179],
      linf = [0.0013653498212873316, 0.0018996803803537077, 0.0018996803803517093, 0.004384818684607161],
      tspan = (0.0, 0.1))
  end
end

end # module


text:
### LP and MILP routines
function add_soc_constraints!(m::GlobalOptimizer, opt::T) where T
    for (func, set) in m._input_problem._conic_second_order
         MOI.add_constraint(opt, func, set)
    end
    return nothing
end

function optimize!(::SOCP, m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType}

    relaxed_optimizer = m.relaxed_optimizer
    MOI.empty!(relaxed_optimizer)

    m._relaxed_variable_index = add_variables(m, relaxed_optimizer, m._input_problem._variable_count)
    add_linear_constraints!(m, relaxed_optimizer)
    add_soc_constraints!(m, relaxed_optimizer)
    add_sv_or_aff_obj!(m, relaxed_optimizer,  m._input_problem._objective)
    MOI.set(relaxed_optimizer, MOI.ObjectiveSense(), m._input_problem._optimization_sense)

    if m._parameters.verbosity < 5
        MOI.set(relaxed_optimizer, MOI.Silent(), true)
    end
    m._parse_time = time() - m._start_time

    MOI.optimize!(relaxed_optimizer)

    #unpack_local_solve!(m, relaxed_optimizer)

    return
end

optimize!(::MISOCP, m::GlobalOptimizer) = optimize!(SOCP(), m)

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: MatthewStuber/EAGO


text:
export evaluate

"""
__Univariate__
```
evaluate(n::Int,x::Array{<:Real},a::AbstractVector{<:Real},b::AbstractVector{<:Real})
evaluate(n::Int,x::Real,a::AbstractVector{<:Real},b::AbstractVector{<:Real})
evaluate(n::Int,x::AbstractVector{<:Real},op::AbstractOrthoPoly)
evaluate(n::Int,x::Real,op::AbstractOrthoPoly)
```
Evaluate the `n`-th univariate basis polynomial at point(s) `x`
The function is multiply dispatched to facilitate its use with the composite type `AbstractOrthoPoly`

If several basis polynomials (stored in `ns`) are to be evaluated at points `x`, then call
```
evaluate(ns::AbstractVector{<:Int},x::AbstractVector{<:Real},op::AbstractOrthoPoly) = evaluate(ns,x,op.α,op.β)
evaluate(ns::AbstractVector{<:Int},x::Real,op::AbstractOrthoPoly) = evaluate(ns,[x],op)
```

If *all* basis polynomials are to be evaluated at points `x`, then call
```
evaluate(x::AbstractVector{<:Real},op::AbstractOrthoPoly) = evaluate(collect(0:op.deg),x,op)
evaluate(x::Real,op::AbstractOrthoPoly) = evaluate([x],op)
```
which returns an Array of dimensions `(length(x),op.deg+1)`.

!!! note
    - `n` is the degree of the univariate basis polynomial
    - `length(x) = N`, where `N` is the number of points
    - `(a,b)` are the recursion coefficients

__Multivariate__
```
evaluate(n::AbstractVector{<:Int},x::AbstractMatrix{<:Real},a::Vector{<:AbstractVector{<:Real}},b::Vector{<:AbstractVector{<:Real}})
evaluate(n::AbstractVector{<:Int},x::AbstractVector{<:Real},a::Vector{<:AbstractVector{<:Real}},b::Vector{<:AbstractVector{<:Real}})
evaluate(n::AbstractVector{<:Int},x::AbstractMatrix{<:Real},op::MultiOrthoPoly)
evaluate(n::AbstractVector{<:Int},x::AbstractVector{<:Real},op::MultiOrthoPoly)
```
Evaluate the n-th p-variate basis polynomial at point(s) x
The function is multiply dispatched to facilitate its use with the composite type `MultiOrthoPoly`

If several basis polynomials are to be evaluated at points `x`, then call
```
evaluate(ind::AbstractMatrix{<:Int},x::AbstractMatrix{<:Real},a::Vector{<:AbstractVector{<:Real}},b::Vector{<:AbstractVector{<:Real}})
evaluate(ind::AbstractMatrix{<:Int},x::AbstractMatrix{<:Real},op::MultiOrthoPoly)
```
where `ind` is a matrix of multi-indices.

If *all* basis polynomials are to be evaluated at points `x`, then call
```
evaluate(x::AbstractMatrix{<:Real},mop::MultiOrthoPoly) = evaluate(mop.ind,x,mop)
```
which returns an array of dimensions `(mop.dim,size(x,1))`.

!!! note
    - `n` is a multi-index
    - `length(n) == p`, i.e. a p-variate basis polynomial
    - `size(x) = (N,p)`, where `N` is the number of points
    - `size(a)==size(b)=p`.
"""
function evaluate(n::Int,x::AbstractArray{<:Real},a::AbstractVector{<:Real},b::AbstractVector{<:Real})
    @assert n >= 0 "Degree n has to be non-negative (currently n=$n)."
    # if length(a)==0 warn("Length of a is 0.") end
    @assert length(a) == length(b) "Inconsistent number of recurrence coefficients."
    @assert n <= length(a) "Specified degree is $n, but you only provided $(length(a)) coefficients."
    # recurrence relation for orthogonal polynomials
    nx = length(x)
    pminus, p = zeros(nx), ones(nx)
    if n==0
        if nx==1
            return first(p)
        else
            return p
        end
    end
    pplus = (x .- first(a)).*p .- first(b)*pminus
    for k in 2:n
        pminus = p
        p = pplus
        @inbounds pplus = (x .- a[k]).*p .- b[k]*pminus
    end
    nx == 1 ? first(pplus) : pplus
end
evaluate(n::Int,x::Real,a::AbstractVector{<:Real},b::AbstractVector{<:Real}) = evaluate(n,[x],a,b)
evaluate(n::Int,x::AbstractVector{<:Real},op::AbstractOrthoPoly) = evaluate(n,x,op.α,op.β)
evaluate(n::Int,x::Real,op::AbstractOrthoPoly) = evaluate(n,[x],op)

# univariate + several bases
function evaluate(ns,x::AbstractArray{<:Real},a::AbstractVector{<:Real},b::AbstractVector{<:Real})
    hcat(map(i->evaluate(i,x,a,b),ns)...)
end
evaluate(ns,x::Real,a::AbstractVector{<:Real},b::AbstractVector{<:Real}) = evaluate(ns,[x],a,b)

evaluate(ns,x::AbstractVector{<:Real},op::AbstractOrthoPoly) = evaluate(ns,x,op.α,op.β)
evaluate(ns,x::Real,op::AbstractOrthoPoly) = evaluate(ns,[x],op)
evaluate(x::AbstractVector{<:Real},op::AbstractOrthoPoly) = evaluate(collect(0:op.deg),x,op)
evaluate(x::Real,op::AbstractOrthoPoly) = evaluate([x],op)

# multivariate
function evaluate(n::AbstractVector{<:Int},x::AbstractMatrix{<:Real},a::AbstractVector{<:AbstractVector{<:Real}},b::AbstractVector{<:AbstractVector{<:Real}})
    @assert length(n) == size(x,2) "number of univariate bases (= $(length(n))) inconsistent with columns points x (= $(size(x,2)))"
    val = ones(Float64,size(x,1))
    for i in 1:length(n)
        @inbounds val = val.*evaluate(n[i],x[:,i],a[i],b[i])
    end
    return val
end
evaluate(n::AbstractVector{<:Int},x::AbstractVector{<:Real},a::AbstractVector{<:AbstractVector{<:Real}},b::AbstractVector{<:AbstractVector{<:Real}}) = evaluate(n,reshape(x,1,length(x)),a,b)
evaluate(n::AbstractVector{<:Int},x::AbstractMatrix{<:Real},op::MultiOrthoPoly) = evaluate(n,x,coeffs(op)...)
evaluate(n::AbstractVector{<:Int},x::AbstractVector{<:Real},op::MultiOrthoPoly) = evaluate(n,reshape(x,1,length(x)),op)

# using multi-index + multivariate
function evaluate(ind::AbstractMatrix{<:Int},x::AbstractMatrix{<:Real},a::AbstractVector{<:AbstractVector{<:Real}},b::AbstractVector{<:AbstractVector{<:Real}})
    vals = map(i->evaluate(ind[i,:],x,a,b),Base.OneTo(size(ind,1)))
    hcat(vals...) |> transpose |> Matrix
end

evaluate(ind::AbstractMatrix{<:Int},x::AbstractMatrix{<:Real},op::MultiOrthoPoly) = evaluate(ind,x,coeffs(op)...)
evaluate(x::AbstractMatrix{<:Real},mop::MultiOrthoPoly) = evaluate(mop.ind,x,mop)

evaluate(ind::AbstractMatrix{<:Int},x::AbstractVector{<:Real},a::AbstractVector{<:AbstractVector{<:Real}},b::AbstractVector{<:AbstractVector{<:Real}}) = evaluate(ind,reshape(x,1,length(x)),a,b)
evaluate(ind::AbstractMatrix{<:Int},x::AbstractVector{<:Real},op::MultiOrthoPoly) = evaluate(ind,reshape(x,1,length(x)),coeffs(op)...)
evaluate(x::AbstractVector{<:Real},mop::MultiOrthoPoly) = evaluate(mop.ind,reshape(x,1,length(x)),mop)