text:
"""
    module Shapes

Definitions of various shapes.

All shapes have the signature: `shape(parameters::Tuple, x0, y0, θ=0)`, and some
have keyword argument `reference` which determines which point `(x0,y0)` is referring to.

The length of `parameters` depends on the shape. For example, `Circle` has just one
parameter, the radius, while `Ellipse` has two, the semi-major and semi-minor axes.
"""
module Shapes
abstract type AbstractDisk <: AbstractShape{2,1} end
abstract type AbstractQuadrilateral <: AbstractShape{2,4} end

export Circle,
Ellipse,
DeformedDisk,
Annulus,
Square,
Rectangle,
Parallelogram,
normal_distance,
get_number_of_sides,
rotate,
unrotate

export AbstractShape,
AbstractDisk,
AbstractQuadrilateral

using ...Defaults
using Formatting
using LinearAlgebra
using NLopt
using RecipesBase
using StaticArrays

abstract type AbstractShape{N} end # N is the number of sides
abstract type AbstractDisk <: AbstractShape{1} end
abstract type AbstractQuadrilateral <: AbstractShape{4} end

get_number_of_sides(s::AbstractShape{N}) where N = N

"""
    struct Circle <: AbstractDisk <: AbstractShape{1}

    Circle((R,), x0, y0, θ=0; reference=:center) -> circle
    Circle(R, x0, y0, θ=0; reference=:center) -> circle

`R` radius
`x0, y0` is location of `reference`

-------------------

    (::Circle)(x,y) --> is_in_circle::Bool
"""
struct Circle <: AbstractDisk
    R::Float64
    x0::Float64
    y0::Float64
    θ::Float64
    sinθ::Float64
    cosθ::Float64
    corners::Tuple{}
    models::Array{Opt,1}

    Circle((R,), args...; kwargs...) = Circle(R[1], args...; kwargs...)
    function Circle(R::Real,x0::Real,y0::Real,θ::Real=0; reference::Symbol=:center)
        if reference ∈ [:bottom, :Bottom, :b, :B]
            x0,y0 = (x0,y0) .+ rotate( 0, R,cos(θ),sin(θ))
        elseif reference ∈ [:top, :Top, :T, :t]
            x0,y0 = (x0,y0) .+ rotate( 0,-R,cos(θ),sin(θ))
        elseif reference ∈ [:left, :Left, :L, :l]
            x0,y0 = (x0,y0) .+ rotate( R, 0,cos(θ),sin(θ))
        elseif reference ∈ [:right, :Right, :R, :r]
            x0,y0 = (x0,y0) .+ rotate(-R, 0,cos(θ),sin(θ))
        end
        c = new(R,x0,y0,θ,sin(θ),cos(θ))
        models = build_models(c)
        return new(R,x0,y0,θ,sin(θ),cos(θ),(),models)
    end

    (c::Circle)(x,y) = hypot((x-c.x0), (y-c.y0)) < c.R-eps(c.R)

    function Base.show(io::IO, circle::Circle)
        print(io, "Circle(R=", fmt("2.2f",circle.R), ", x0=", fmt("2.2f",circle.x0), ", y0=", fmt("2.2f",circle.y0), ")")
    end
end
"""
    perimeter(t,shape,side) -> x,y,dx/dt,dy/dt
"""
function perimeter(t::Number,c::Circle,side::Int)
    sn,cs = sincos(t)
    x = c.x0 + c.R*cs
    y = c.y0 + c.R*sn
    ẋ = -c.R*sn
    ẏ =  c.R*cs
    return x,y,ẋ,ẏ
end

"""
    struct Ellipse <: AbstractDisk <: AbstractShape{1}

    Ellipse((a,b), x0, y0, θ=0; reference=:center) -> ellipse
    Ellipse(a, b, x0, y0, θ=0; reference=:center) -> ellipse

`a, b` axis lengths
`x0, y0` is location of `reference`
`θ` angle of `a` axis

----------------

    (::Ellipse)(x,y) -> is_in_ellipse::Bool
"""
struct Ellipse <: AbstractDisk
    a::Float64
    b::Float64
    x0::Float64
    y0::Float64
    θ::Float64
    sinθ::Float64
    cosθ::Float64
    corners::Tuple{}
    models::Array{Opt,1}

    Ellipse((a,b),args...;kwargs...) = Ellipse(a,b,args...;kwargs...)
    function Ellipse(a::Number,b::Number,x0::Number,y0::Number,θ::Number=0; reference::Symbol=:center)
        if reference ∈ [:bottom, :Bottom, :b, :B]
            x0,y0 = (x0,y0) .+ rotate( 0, b,cos(θ),sin(θ))
        elseif reference ∈ [:top, :Top, :T, :t]
            x0,y0 = (x0,y0) .+ rotate( 0,-b,cos(θ),sin(θ))
        elseif reference ∈ [:left, :Left, :L, :l]
            x0,y0 = (x0,y0) .+ rotate( a, 0,cos(θ),sin(θ))
        elseif reference ∈ [:right, :Right, :R, :r]
            x0,y0 = (x0,y0) .+ rotate(-a, 0,cos(θ),sin(θ))
        end
         e = new(a,b,x0,y0,θ,sin(θ),cos(θ))
         models = build_models(e)
         return new(a,b,x0,y0,θ,sin(θ),cos(θ),(),models)
    end

    Ellipse(c::Circle) = Ellipse(c.R,c.R,c.x0,c.y0,c.θ)

    function (e::Ellipse)(x,y)
        xurot, yurot = unrotate(x, y, e)
        return hypot((xurot-e.x0)/e.a, (yurot-e.y0)/e.b) < 1-eps()
    end

    function Base.show(io::IO, ellipse::Ellipse)
        print(io, "Ellipse(a=", fmt("2.2f",ellipse.a), ", b=", fmt("2.2f",ellipse.b), ", x0=", fmt("2.2f",ellipse.x0), ", y0=", fmt("2.2f",ellipse.y0), ", θ=∠", fmt("3.2f",(mod2pi(ellipse.θ))*180/π), "°)")
    end
end
function perimeter(t::Number,e::Ellipse,side::Int)
    sn,cs = sincos(t)
    xur,yur =    e.a*cs, e.b*sn
    xpur,ypur = -e.a*sn, e.b*cs
    x,y = (e.x0,e.y0) .+ rotate(xur,yur,e.cosθ,e.sinθ)
    ẋ,ẏ = rotate(xpur,ypur,e.cosθ,e.sinθ)
    return x,y,ẋ,ẏ
end

"""
    struct DeformedDisk{N} <: AbstractDisk <: AbstractShape{1}

    DeformedDisk{N}(R, x0, y0, M, a, φ, θ=0) -> deformeddisk
    DeformedDisk{N}((R,M,a,φ), x0, y0, θ=0) -> deformeddisk

`R` is radius,
`x0` and `y0` is center of circle of radius `R`,
`M` is array of length `N` of multipole integers,
`a` is array of length `N` amplitudes,
`φ` is array of length `N` of angles
`θ` is overall rotation angle

-----------------

    (::DeformedDisk)(x,y) -> is_in_disk::Bool
"""
struct DeformedDisk{N} <: AbstractDisk
    R::Float64
    x0::Float64
    y0::Float64
    θ::Float64
    M::SArray{Tuple{N},Int,1,N}
    a::SArray{Tuple{N},Float64,1,N}
    φ::SArray{Tuple{N},Float64,1,N}
    sinθ::Float64
    cosθ::Float64
    corners::Tuple{}
    models::Array{Opt,1}

    DeformedDisk{n}((R,M,a,φ),x0,y0,args...;kwargs...) where n= DeformedDisk{n}(R,x0,y0,M,a,φ,args...;kwargs...)
    function DeformedDisk{n}(R::Number,x0::Number,y0::Number,M,a,φ,θ=0) where n
        @assert n==length(M)==length(a)==length(φ) "parameter N in DeformedDisk{N}(...) must be equal to length(M)"
        d = new{n}(R,x0,y0,θ,SVector{n}(M),SVector{n}(a),SVector{n}(φ.+θ),sin(θ),cos(θ))
        models = build_models(d)
        return new{n}(R,x0,y0,d.θ,d.M,d.a,d.φ,d.sinθ,d.cosθ,(),models)
    end

    DeformedDisk(c::Circle) = DeformedDisk{0}(c.R,c.x0,c.y0,[],[],[],c.θ)

    function (d::DeformedDisk)(x,y)
        θ = atan(y-d.y0,x-d.x0)
        r = hypot(x-d.x0,y-d.y0) < d.R + sum(d.a.*(map((m,φ)->cos(m*(θ-φ)),d.M,d.φ))) - eps(d.R)
        return r
    end

    function Base.show(io::IO, d::DeformedDisk{N}) where N
         print(io,"DeformedDisk with $N multipoles")
    end
end
function perimeter(t::Number,d::DeformedDisk,side::Int)
    r = d.R + sum(d.a.*(map((m,φ)->cos(m*(t-φ)),d.M,d.φ)))
    ṙ = sum(d.a.*(map((m,φ)->-m*sin(m*(t-φ)),d.M,d.φ)))
    sn,cs = sincos(t)
    x,y = d.x0 + r*cs, d.y0 + r*sn
    ẋ = ṙ*cs - r*sn
    ẏ = ṙ*sn + r*cs
    return x,y,ẋ,ẏ
end

# plotting for all disks
@recipe function f(d::AbstractDisk)
    θ = LinRange(0,2π,201)
    x = Array{Float64,1}(undef,length(θ))
    y = Array{Float64,1}(undef,length(θ))
    for i ∈ eachindex(θ)
        x[i],y[i] = perimeter(θ[i],d,1)
    end
    alpha --> 0
    seriestype --> :path
    fillcolor --> SHAPE_COLOR
    fillrange --> d.y0
    fillalpha --> SHAPE_FILL_ALPHA
    aspect_ratio --> 1
    legend-->false
    (x,y)
end

"""
    struct Annulus <: AbstractShape{2}

    Annulus((R1,R2),x0,y0,θ=0) -> annulus
    Annulus(R1, R2, x0, y0, θ=0) -> annulus

`R1` is inner radius, `R2` outer
`x0, y0` is center

-------------
    (::Annulus)(x,y) -> is_in_annulus::Bool
"""
struct Annulus <: AbstractShape{2}
    R1::Float64
    R2::Float64
    x0::Float64
    y0::Float64
    θ::Float64
    sinθ::Float64
    cosθ::Float64
    corners::Tuple{}
    models::Array{Opt,1}

    Annulus((R1,R2),args...;kwargs...) = Annulus(R1,R2,args...;kwargs...)
    function Annulus(R1::Number,R2::Number,x0::Number,y0::Number,θ::Number=0)
        @assert R1≤R2 "R1=$R1, R2=$R2 must satisfy R1≤R2"
        a = new(R1,R2,x0,y0,0,0,1)
        models = build_models(a)
        return new(R1,R2,x0,y0,0,sin(θ),cos(θ),(),models)
    end

    (c::Annulus)(x,y) = c.R1+eps(c.R1) < hypot( (x-c.x0),(y-c.y0) ) < c.R2-eps(c.R2)

    function Base.show(io::IO, a::Annulus)
        print(io, "Annulus(R1=$(fmt("2.2f",a.R1)), R2=$(fmt("2.2f",a.R2)), x0=$(fmt("2.2f",a.x0)), y0=$(fmt("2.2f",a.y0)))")
    end
end
function perimeter(t::Number,a::Annulus,side::Int)
    sn,cs = sincos(t)
    R = side==1 ? a.R1 : a.R2
    x = a.x0 + R*cs
    y = a.y0 + R*sn
    ẋ = -R*sn
    ẏ =  R*cs
    return x,y,ẋ,ẏ
end
@recipe function f(cr::Annulus)
    @series Circle(cr.R1,cr.x0,cr.y0)
    @series Circle(cr.R2,cr.x0,cr.y0)
end

"""
    struct Square <: AbstractQuadrilateral <: AbstractShape{4}

    Square((a,),x0,y0,θ=0; reference=:center) -> square
    Square(a,x0,y0,θ=0); reference=:center -> square

`a` length of sides
`x0, y0` is location of `reference`
`θ` angle of `a` side

------------

    (::Square)(x,y) -> is_in_square::Bool
"""
struct Square <: AbstractQuadrilateral
    a::Float64
    x0::Float64
    y0::Float64
    θ::Float64
    sinθ::Float64
    cosθ::Float64
    corners::NTuple{4,NTuple{2,Float64}}
    models::Array{Opt,1}

    Square((a,),args...;kwargs...) = Square(a,args...;kwargs...)
    function Square(a::Number,x0::Number,y0::Number,θ::Number=0; reference::Symbol=:center)
        if reference ∈ [:bottom, :Bottom, :b, :B]
            x0,y0 = (x0,y0) .+ rotate( 0, a/2,cos(θ),sin(θ))
        elseif reference ∈ [:top, :Top, :T, :t]
            x0,y0 = (x0,y0) .+ rotate( 0,-a/2,cos(θ),sin(θ))
        elseif reference ∈ [:left, :Left, :L, :l]
            x0,y0 = (x0,y0) .+ rotate( a/2, 0,cos(θ),sin(θ))
        elseif reference ∈ [:right, :Right, :R, :r]
            x0,y0 = (x0,y0) .+ rotate(-a/2, 0,cos(θ),sin(θ))
        elseif reference ∈ [:topright, :TopRight, :TR, :tr, :ne, :NE]
            x0,y0 = (x0,y0) .+ rotate(-a/2,-a/2,cos(θ),sin(θ))
        elseif reference ∈ [:topleft, :TopLeft, :TL, :tl, :nw, :NW]
            x0,y0 = (x0,y0) .+ rotate( a/2,-a/2,cos(θ),sin(θ))
        elseif reference ∈ [:bottomleft, :BottomLeft, :BL, :bl, :sw, :SW]
            x0,y0 = (x0,y0) .+ rotate( a/2, a/2,cos(θ),sin(θ))
        elseif reference ∈ [:bottomright, :BottomRight, :BR, :br, :se, :SE]
            x0,y0 = (x0,y0) .+ rotate(-a/2, a/2,cos(θ),sin(θ))
        end
        s = new(a,x0,y0,θ,sin(θ),cos(θ))
        c1 = rotate(s.x0+s.a/2,s.y0+s.a/2, s)
        c2 = rotate(s.x0-s.a/2,s.y0+s.a/2, s)
        c3 = rotate(s.x0-s.a/2,s.y0-s.a/2, s)
        c4 = rotate(s.x0+s.a/2,s.y0-s.a/2, s)
        models = build_models(s)
        return new(a,x0,y0,θ,sin(θ),cos(θ),(c1,c2,c3,c4),models)
    end

    function (s::Square)(x,y)
        xurot, yurot = unrotate(x, y, s)
        return -s.a/2+10eps(s.a/2) < xurot-s.x0 < s.a/2-10eps(s.a/2)  &&  -s.a/2+10eps(s.a/2) < yurot-s.y0 < s.a/2-10eps(s.a/2)
    end

    function Base.show(io::IO, square::Square)
        print(io, "Square(a=$(fmt("2.2f",square.a)), x0=$(fmt("2.2f",square.x0)), y0=$(fmt("2.2f",square.y0)), θ=∠", fmt("3.2f",(mod2pi(square.θ))*180/π), "°)")
    end
end
function perimeter(t::Number,s::Square,side::Int)
    if side==1
        x = -s.a/2
        y = -s.a/2 + t*s.a
        ẋ = 0
        ẏ = s.a
    elseif side==2
        x =  s.a/2
        y = -s.a/2 + t*s.a
        ẋ = 0
        ẏ = s.a
    elseif side==3
        x = -s.a/2 + t*s.a
        y = -s.a/2
        ẋ = s.a
        ẏ = 0
    elseif side==4
        x = -s.a/2 + t*s.a
        y =  s.a/2
        ẋ = s.a
        ẏ = 0
    end
    x += s.x0
    y += s.y0
    xrot,yrot = rotate(x,y,s)
    ẋrot,ẏrot = rotate(ẋ,ẏ,s.cosθ,s.sinθ)
    return xrot,yrot,ẋrot,ẏrot
end

"""
    struct Rectangle <: AbstractQuadrilateral <: AbstractShape{4}

    Rectangle((a,b), x0, y0, θ=0; reference=:center) -> rectangle
    Rectangle(a, b, x0, y0, θ=0; reference=:center) -> rectangle

`a, b` length of sides
`x0, y0` is location of `reference`
`θ` angle of `a` side

-----------------------

    (::Rectangle)(x,y) -> is_in_rectangle::Bool
"""
struct Rectangle <: AbstractQuadrilateral
    a::Float64
    b::Float64
    x0::Float64
    y0::Float64
    θ::Float64
    sinθ::Float64
    cosθ::Float64
    corners::NTuple{4,NTuple{2,Float64}}
    models::Array{Opt,1}

    Rectangle((a,b),args...;kwargs...) = Rectangle(a,b,args...;kwargs...)
    function Rectangle(a::Number,b::Number,x0::Number,y0::Number,θ::Number=0; reference::Symbol=:center)
        if reference ∈ [:bottom, :Bottom, :b, :B]
            x0,y0 = (x0,y0) .+ rotate( 0, b/2,cos(θ),sin(θ))
        elseif reference ∈ [:top, :Top, :T, :t]
            x0,y0 = (x0,y0) .+ rotate( 0,-b/2,cos(θ),sin(θ))
        elseif reference ∈ [:left, :Left, :L, :l]
            x0,y0 = (x0,y0) .+ rotate( a/2, 0,cos(θ),sin(θ))
        elseif reference ∈ [:right, :Right, :R, :r]
            x0,y0 = (x0,y0) .+ rotate(-a/2, 0,cos(θ),sin(θ))
        elseif reference ∈ [:topright, :TopRight, :TR, :tr, :ne, :NE]
            x0,y0 = (x0,y0) .+ rotate(-a/2,-b/2,cos(θ),sin(θ))
        elseif reference ∈ [:topleft, :TopLeft, :TL, :tl, :nw, :NW]
            x0,y0 = (x0,y0) .+ rotate( a/2,-b/2,cos(θ),sin(θ))
        elseif reference ∈ [:bottomleft, :BottomLeft, :BL, :bl, :sw, :SW]
            x0,y0 = (x0,y0) .+ rotate( a/2, b/2,cos(θ),sin(θ))
        elseif reference ∈ [:bottomright, :BottomRight, :BR, :br, :se, :SE]
            x0,y0 = (x0,y0) .+ rotate(-a/2, b/2,cos(θ),sin(θ))
        end
        r = new(a,b,x0,y0,θ,sin(θ),cos(θ))
        models = build_models(r)
        corners_x, corners_y = Array{Float64}(undef,4), Array{Float64}(undef,4)
        c1 = (r.x0,r.y0) .+ rotate( r.a/2, r.b/2,r.cosθ,r.sinθ)
        c2 = (r.x0,r.y0) .+ rotate(-r.a/2, r.b/2,r.cosθ,r.sinθ)
        c3 = (r.x0,r.y0) .+ rotate(-r.a/2,-r.b/2,r.cosθ,r.sinθ)
        c4 = (r.x0,r.y0) .+ rotate( r.a/2,-r.b/2,r.cosθ,r.sinθ)
        return new(a,b,x0,y0,θ,sin(θ),cos(θ),(c1,c2,c3,c4),models)
    end

    Rectangle(s::Square) = Rectangle(s.a,s.a,s.x0,s.y0,s.θ)

    function (r::Rectangle)(x,y)
        xurot, yurot = unrotate(x, y, r)
        return -r.a/2+20eps(r.a/2) < xurot-r.x0 < r.a/2-20eps(r.a/2)  &&  -r.b/2+20eps(r.b/2) < yurot-r.y0 < r.b/2-20eps(r.b/2)
    end

    function Base.show(io::IO, rect::Rectangle)
        print(io, "Rectangle(a=$(fmt("2.2f",rect.a)), b=$(fmt("2.2f",rect.b)), x0=$(fmt("2.2f",rect.x0)), y0=$(fmt("2.2f",rect.y0)), θ=∠", fmt("3.2f",(mod2pi(rect.θ))*180/π), "°)")
    end
end
function perimeter(t::Number,s::Rectangle,side::Int)
    if side==1
        x = -s.a/2
        y = -s.b/2 + t*s.b
        ẋ = 0
        ẏ = s.b
    elseif side==2
        x =  s.a/2
        y = -s.b/2 + t*s.b
        ẋ = 0
        ẏ = s.b
    elseif side==3
        x = -s.a/2 + t*s.a
        y = -s.b/2
        ẋ = s.a
        ẏ = 0
    elseif side==4
        x = -s.a/2 + t*s.a
        y =  s.b/2
        ẋ = s.a
        ẏ = 0
    end
    x += s.x0
    y += s.y0
    xrot,yrot = rotate(x,y,s)
    ẋrot,ẏrot = rotate(ẋ,ẏ,s.cosθ,s.sinθ)
    return xrot,yrot,ẋrot,ẏrot
end

# """
#     Parallelogram(a, b, α, x0, y0, θ)
# """
# struct Parallelogram <: AbstractParallelogram
#     a::Float64
#     b::Float64
#     α::Float64
#     x0::Float64
#     y0::Float64
#     θ::Float64
#     sinθ::Float64
#     cosθ::Float64
#     tanα::Float64
#     cosα::Float64
#     models::Array{Opt,1}
#
#     function Parallelogram(a::Number, b::Number, α::Number, x0::Number, y0::Number, θ::Number)
#         new(a,b,α,x0,y0,θ,cos(θ),sin(θ),tan(α),cos(α))
#     end
#
#     function (p::Parallelogram)(x,y)
#         xrot, yrot = rotate(x-p.x0, y-p.y0, p.cosθ, p.sinθ)
#         return p.tanα*(xrot-p.a) < yrot < p.tanα*xrot  &&  0 < yrot < p.cosα*p.b
#     end
#
#     Base.show(io::IO, par::Parallelogram) = print(io, "Parallelogram(a=$(fmt("2.2f",par.a)), b=$(fmt("2.2f",par.b)), α=$(fmt("2.2f",par.α)), x0=$(fmt("2.2f",par.x0)), y0=$(fmt("2.2f",par.y0)), θ=$(fmt("2.2f",par.θ)))")
#
#     @recipe function f(pg::Parallelogram)
#         x = cumsum([0, pg.a, +pg.b*cos(pg.α), -pg.a, -pg.b*cos(pg.α)])
#         y = cumsum([0, 0, +pg.b*sin(pg.α), 0, -pg.b*sin(pg.α)])
#         alpha --> 0
#         seriestype --> :path
#         fillcolor --> SHAPE_COLOR
#         fillrange --> 0
#         fillalpha --> SHAPE_FILL_ALPHA
#         aspect_ratio --> 1
#         legend --> false
#         xrot, yrot = rotate(x,y,pg.cosθ,-pg.sinθ)
#         pg.x0 .+ xrot, pg.y0 .+ yrot
#     end
# end

# plot all quadrilaterals
@recipe function f(q::AbstractQuadrilateral)
    t = LinRange(0,1,201)
    x1,y1 = Array{Float64}(undef,length(t)), Array{Float64}(undef,length(t))
    x2,y2 = Array{Float64}(undef,length(t)), Array{Float64}(undef,length(t))
    x3,y3 = Array{Float64}(undef,length(t)), Array{Float64}(undef,length(t))
    x4, y4 = Array{Float64}(undef,length(t)), Array{Float64}(undef,length(t))
    for i ∈ eachindex(t)
        x1[i],y1[i] = perimeter(t[i],q,1)
        x2[i],y2[i] = perimeter(t[i],q,2)
        x3[i],y3[i] = perimeter(t[i],q,3)
        x4[i],y4[i] = perimeter(t[i],q,4)
    end
    alpha --> 0
    seriestype --> :path
    fillcolor --> SHAPE_COLOR
    fillrange --> q.y0
    fillalpha --> SHAPE_FILL_ALPHA
    aspect_ratio --> 1
    legend --> false
    vcat(x1,x4,reverse(x3)),vcat(y1,y4,reverse(y3))
end

"""
    rotate(x,y,cosθ,sinθ) -> xrot, yrot

rotate by `θ` about origin

----------------
    rotate(x,y,shape) -> xrot, yrot

rotate about `shape`'s center by `shape`'s angle
"""
rotate(x::Real,y::Real,cosθ::Real,sinθ::Real) = iszero(sinθ) ? (x,y) : (cosθ*x-sinθ*y, sinθ*x+cosθ*y)
rotate(x::Real,y::Real,s::AbstractShape) = (s.x0,s.y0) .+ rotate(x-s.x0,y-s.y0,s.cosθ,s.sinθ)
function rotate(x::AbstractArray,y::AbstractArray,cosθ::Real,sinθ::Real)
    xrot = Array{Float64}(undef,size(xyrot))
    yrot = Array{Float64}(undef,size(xyrot))
    rotate!(xrot,yrot,x,y,cosθ,sinθ)
    return xrot, yrot
end
function rotate(x::AbstractArray,y::AbstractArray,s::AbstractShape)
    xrot,yrot = rotate(x.-s.x0,y.-y.s0,s.cosθ,s.sinθ)
    for i ∈ eachindex(xrot)
        xrot[i] += s.x0
        yrot[i] += s.y0
    end
    return xrot,yrot
end
function rotate!(xrot::AbstractArray,yrot::AbstractArray,x::AbstractArray,y::AbstractArray,cosθ::Real,sinθ::Real)
    for i ∈ eachindex(x)
        xrot[i],yrot[i] = rotate(x[i],y[i],cosθ,sinθ)
    end
    return nothing
end

"""
    unrotate(x,y,cosθ,sinθ) -> xrot, yrot

rotate by `-θ` about origin

----------------
    unrotate(x,y,shape) -> xrot, yrot

rotate about `shape`'s center by the negative of `shape`'s angle
"""
unrotate(x,y,cosθ::Real,sinθ::Real) = rotate(x,y,cosθ,-sinθ)
unrotate(x::Real,y::Real,s::AbstractShape) = (s.x0,s.y0) .+ unrotate(x-s.x0,y-s.y0,s.cosθ,s.sinθ)
function unrotate(x::AbstractArray,y::AbstractArray,s::AbstractShape)
    xrot,yrot = unrotate(x.-s.x0,y.-y.s0,s.cosθ,s.sinθ)
    for i ∈ eachindex(xrot)
        xrot[i] += s.x0
        yrot[i] += s.y0
    end
    return xrot,yrot
end

function build_models(s::AbstractShape{N}) where N
    models = Array{Opt}(undef,N)
    for i ∈ eachindex(models)
        models[i] = Opt(NL_NORMAL_ALGORITHM, 1)
        models[i].xtol_rel = NL_NORMAL_XTOL
        models[i].ftol_rel = NL_NORMAL_FTOL
        models[i].maxeval = NL_NORMAL_MAXEVAL
        if typeof(s)<:AbstractQuadrilateral
            models[i].lower_bounds = 0
            models[i].upper_bounds = 1
        end
    end
    return models
end

function model_objective(v::Vector,grad::Vector,X::Number,Y::Number,s::AbstractShape,side::Int)
    x,y,ẋ,ẏ = perimeter(v[1],s,side)
    if length(grad) > 0
        grad[1] = 2(x-X)*ẋ + 2(y-Y)*ẏ
    end
    return (x-X)^2 + (y-Y)^2
end

"""
    normal_distance(shape,x,y) -> normal, tangent, distance, side

Get the `normal` vector from the surface of `shape` to a point `x,y`,
and the associated `tangent` vector and distance from surface `side`.

The returned arguments are sorted in ascending order by `distance`, so that, e.g.,
`normal[1]` is the normal vector associated with the closest side.
"""
function normal_distance(s::TS,x::Number,y::Number) where TS<:AbstractShape
    n = Array{Array{Float64,1}}(undef,length(s.models))
    d = Array{Float64}(undef,length(s.models))
    sides = 1:length(s.models)
    for i ∈ eachindex(s.models)
        model = s.models[i]
        model.min_objective = (v,g)->model_objective(v,g,x,y,s,i)
        minf,mint,ret = optimize(model,[t0_shape(x,y,s)])
        X,Y = perimeter(mint[1],s,i)
        n[i] = [x-X, y-Y]
        n[i] = n[i]/norm(n[i])
        d[i] = sqrt(minf)
    end
    perm = sortperm(d)
    n, d = n[perm], d[perm]
    NORMAL = Array{Array{Float64,1}}(undef,length(s.models))
    TANGENT = Array{Array{Float64,1}}(undef,length(s.models))
    for i ∈ eachindex(NORMAL)
        nx, ny = n[i][1],n[i][2]
        NORMAL[i] = [nx, ny]
        TANGENT[i] = [ny,-nx]
    end
    return NORMAL, TANGENT, d, sides[perm]
end

t0_shape(x::Number,y::Number,s::AbstractQuadrilateral) = .5
t0_shape(x::Number,y::Number,s::AbstractDisk) = atan(y-s.y0,x-s.x0)

end #module


text:
# This file contains a function
# to write input data to a file

"""
    write(input::InputStruct, file::AbstractString)

Write input data from an instance of [`InputStruct`](@ref) to a file.

# Usage

```jldoctest; output = false
using Scats
s = Scats.API()
file, _ = mktemp()
s.Input.write(file)

# output

```
"""
function write(input::InputStruct, file::AbstractString)

    open(file, "w") do f
        # Print
        println(f, "Sample size")
        println(f, input.N)
        println(f, "\nSample step")
        println(f, input.Δt)
        println(f, "\nSignificance level")
        println(f, input.q)
        println(f, "\nTime array")
        println(f, input.t)
        println(f, "\nValues array")
        println(f, input.x)
    end

end

ext: jl


text:
# test filtering

T = Float32
fs = 100.  # sampling rate
N = 2^14   # number of points
Nwin = 10
freqmin = 10. # minimum frequency
freqmax = 20. # maximum frequency
corners = 4   # number of corners
zerophase = true # do zerophase filtering
starttime = d2u(DateTime(Date(now()))) # starttime for SeisChanel / SeisData

# test bandpass
@testset "bandpass" begin

# create white noise
A = rand(T,N,Nwin)
B = bandpass(A,freqmin,freqmax,fs,corners=corners,zerophase=zerophase)
# check size and type stability
@test size(B) == size(A)
@test eltype(B) == eltype(A)
bandpass!(A,freqmin,freqmax,fs,corners=corners,zerophase=zerophase)
@test A == B

# test power in pass bands
freqs = FFTW.rfftfreq(N,fs)
lowband = findall(0 + eps(T) .< freqs .< freqmin)
passband = findall(freqmin .<= freqs .<= freqmax)
highband = findall(freqs .> freqmax)
F = rfft(A,1)
Fpow = abs.(F)
# test lowband
@test all(median(Fpow[lowband,:],dims=1) .< median(Fpow[passband,:],dims=1))

# test highband
@test all(median(Fpow[highband,:],dims=1) .< median(Fpow[passband,:],dims=1))

# test zerophase
A = rand(T,N,Nwin)
@test bandpass(A,freqmin,freqmax,fs,corners=corners,zerophase=true) != bandpass(A,freqmin,
               freqmax,fs,corners=corners,zerophase=false)

# test corners
@test bandpass(A,freqmin,freqmax,fs,corners=4,zerophase=true) != bandpass(A,freqmin,
               freqmax,fs,corners=3,zerophase=true)

# test high frequency
highfreq = fs
B = highpass(A,freqmin,fs,corners=corners,zerophase=zerophase)
bandpass!(A,freqmin,highfreq,fs,corners=4,zerophase=true)
@test A == B

# test error with freqmin > fs / 2
@test_throws ArgumentError bandpass!(A,highfreq,freqmax,fs,corners=4,zerophase=true)

# test RawData filtering
A = rand(T,N,Nwin)
R = RawData()
R.x = A
R.fs = fs
Rnew = bandpass(R,freqmin,freqmax)
RInt = bandpass(R,Int(freqmin),Int(freqmax))
@test Rnew.fs == R.fs
@test Rnew.freqmin == freqmin
@test Rnew.freqmax == freqmax
@test RInt.freqmin == freqmin
@test RInt.freqmax == freqmax
bandpass!(R,freqmin,freqmax)
@test R.fs == R.fs
@test R.freqmin == freqmin
@test R.freqmax == freqmax
@test R.x == Rnew.x
@test R.x == RInt.x

# test zerophase
A = rand(T,N,Nwin)
R.x = A
@test bandpass(R,freqmin,freqmax,corners=corners,zerophase=true).x != bandpass(R,freqmin,
               freqmax,corners=corners,zerophase=false).x

# test corners
@test bandpass(R,freqmin,freqmax,corners=4,zerophase=true).x != bandpass(R,freqmin,
              freqmax,corners=3,zerophase=true).x

# test CorrData filtering
C = CorrData()
C.corr = A
C.fs = fs
Cnew = bandpass(C,freqmin,freqmax)
CInt = bandpass(C,Int(freqmin),Int(freqmax))
@test Cnew.fs == C.fs
@test Cnew.freqmin == freqmin
@test Cnew.freqmax == freqmax
@test CInt.freqmin == freqmin
@test CInt.freqmax == freqmax
bandpass!(C,freqmin,freqmax)
@test C.fs == C.fs
@test C.freqmin == freqmin
@test C.freqmax == freqmax
@test C.corr == Cnew.corr
@test C.corr == CInt.corr

# test zerophase
A = rand(T,N,Nwin)
C.corr = A
@test bandpass(C,freqmin,freqmax,corners=corners,zerophase=true).corr != bandpass(C,freqmin,
               freqmax,corners=corners,zerophase=false).corr

# test corners
@test bandpass(C,freqmin,freqmax,corners=4,zerophase=true).corr != bandpass(C,freqmin,
              freqmax,corners=3,zerophase=true).corr

# test filtering SeisChannel
C = SeisChannel()
C.x = rand(T,N)
C.fs = fs
C.t = [1 starttime * 1e6;length(C.x) 0]
Cnew = bandpass(C,freqmin,freqmax)
CInt = bandpass(C,Int(freqmin),Int(freqmax))
@test Cnew.fs == fs
bandpass!(C,freqmin,freqmax)
@test C.fs == C.fs
@test C.x == Cnew.x
@test C.x == CInt.x

# test zerophase -> there is no option for zerophase here
A = rand(T,N)
C.x = A
@test bandpass(C,freqmin,freqmax,corners=corners,zerophase=true).x == bandpass(C,freqmin,
               freqmax,corners=corners,zerophase=false).x

# test corners
@test bandpass(C,freqmin,freqmax,corners=4,zerophase=true).x != bandpass(C,freqmin,
              freqmax,corners=3,zerophase=true).x

# test filtering SeisData
S = SeisData(2)
S.x[1] = rand(T,N)
S.fs[1] = fs
S.t[1] = [1 starttime * 1e6;length(C.x) 0]
S.x[2] = rand(T,N)
S.fs[2] = fs
S.t[2] = [1 starttime * 1e6;length(C.x) 0]
Snew = bandpass(S,freqmin,freqmax)
SInt = bandpass(S,Int(freqmin),Int(freqmax))
@test Snew.fs == S.fs
bandpass!(S,freqmin,freqmax)
@test all(S.fs .== fs)
@test S.x[1] == Snew.x[1]
@test S.x[2] == Snew.x[2]
@test S.x[1] == SInt.x[1]
@test S.x[2] == SInt.x[2]

# test zerophase -> there is no option for zerophase here
A = rand(T,N,Nwin)
S.x[1] = rand(T,N)
S.x[2] = rand(T,N)
@test bandpass(S,freqmin,freqmax,corners=corners,zerophase=true).x[1] == bandpass(S,freqmin,
               freqmax,corners=corners,zerophase=false).x[1]
@test bandpass(S,freqmin,freqmax,corners=corners,zerophase=true).x[2] == bandpass(S,freqmin,
              freqmax,corners=corners,zerophase=false).x[2]

# test corners
@test bandpass(S,freqmin,freqmax,corners=4,zerophase=true).x[1] != bandpass(S,freqmin,
              freqmax,corners=3,zerophase=true).x[1]
@test bandpass(S,freqmin,freqmax,corners=4,zerophase=true).x[2] != bandpass(S,freqmin,
            freqmax,corners=3,zerophase=true).x[2]

# test bandpass with Int inputs
A = rand(T,N,Nwin)
B = bandpass(A,Int(freqmin),Int(freqmax),Int(fs),corners=corners,zerophase=zerophase)
# check size and type stability
@test size(B) == size(A)
@test eltype(B) == eltype(A)
bandpass!(A,freqmin,freqmax,fs,corners=corners,zerophase=zerophase)
@test A == B

# NEEDED: test GPUArray
end

@testset "bandstop" begin

# create white noise
A = rand(T,N,Nwin)
B = bandstop(A,freqmin,freqmax,fs,corners=corners,zerophase=zerophase)
# check size and type stability
@test size(B) == size(A)
@test eltype(B) == eltype(A)
bandstop!(A,freqmin,freqmax,fs,corners=corners,zerophase=zerophase)
@test A == B

# test power in pass bands
freqs = FFTW.rfftfreq(N,fs)
lowband = findall(0 + eps(T) .< freqs .< freqmin)
passband = findall(freqmin .<= freqs .<= freqmax)
highband = findall(freqs .> freqmax)
F = rfft(A,1)
Fpow = abs.(F)
# test lowband
@test all(median(Fpow[lowband,:],dims=1) .> median(Fpow[passband,:],dims=1))

# test highband
@test all(median(Fpow[highband,:],dims=1) .> median(Fpow[passband,:],dims=1))

# test zerophase
A = rand(T,N,Nwin)
@test bandstop(A,freqmin,freqmax,fs,corners=corners,zerophase=true) != bandstop(A,freqmin,
               freqmax,fs,corners=corners,zerophase=false)

# test corners
@test bandstop(A,freqmin,freqmax,fs,corners=4,zerophase=true) != bandstop(A,freqmin,
               freqmax,fs,corners=3,zerophase=true)

# test high freq
highfreq = fs / 2
@test_throws ArgumentError bandstop!(A,freqmin,highfreq,fs,corners=4,zerophase=true)

# test low freq
highfreq = fs / 2
@test_throws ArgumentError bandstop!(A,highfreq,freqmax,fs,corners=4,zerophase=true)

# test RawData filtering
A = rand(T,N,Nwin)
R = RawData()
R.x = A
R.fs = fs
Rnew = bandstop(R,freqmin,freqmax)
RInt = bandstop(R,Int(freqmin),Int(freqmax))
@test Rnew.fs == R.fs
bandstop!(R,freqmin,freqmax)
@test R.fs == R.fs
@test R.x == Rnew.x
@test R.x == RInt.x

# test zerophase
A = rand(T,N,Nwin)
R.x = A
@test bandstop(R,freqmin,freqmax,corners=corners,zerophase=true).x != bandstop(R,freqmin,
               freqmax,corners=corners,zerophase=false).x

# test corners
@test bandstop(R,freqmin,freqmax,corners=4,zerophase=true).x != bandstop(R,freqmin,
              freqmax,corners=3,zerophase=true).x

# test CorrData filtering
C = CorrData()
C.corr = A
C.fs = fs
Cnew = bandstop(C,freqmin,freqmax)
CInt = bandstop(C,Int(freqmin),Int(freqmax))
@test Cnew.fs == C.fs
bandstop!(C,freqmin,freqmax)
@test C.fs == C.fs
@test C.corr == Cnew.corr
@test C.corr == CInt.corr

# test zerophase
A = rand(T,N,Nwin)
C.corr = A
@test bandstop(C,freqmin,freqmax,corners=corners,zerophase=true).corr != bandstop(C,freqmin,
               freqmax,corners=corners,zerophase=false).corr

# test corners
@test bandstop(C,freqmin,freqmax,corners=4,zerophase=true).corr != bandstop(C,freqmin,
              freqmax,corners=3,zerophase=true).corr

# test filtering SeisChannel
C = SeisChannel()
C.x = rand(T,N)
C.fs = fs
C.t = [1 starttime * 1e6;length(C.x) 0]
Cnew = bandstop(C,freqmin,freqmax)
CInt = bandstop(C,Int(freqmin),Int(freqmax))
@test Cnew.fs == fs
bandstop!(C,freqmin,freqmax)
@test C.fs == C.fs
@test C.x == Cnew.x
@test C.x == CInt.x

# test zerophase -> there is no option for zerophase here
A = rand(T,N)
C.x = A
@test bandstop(C,freqmin,freqmax,corners=corners,zerophase=true).x == bandstop(C,freqmin,
               freqmax,corners=corners,zerophase=false).x

# test corners
@test bandstop(C,freqmin,freqmax,corners=4,zerophase=true).x != bandpass(C,freqmin,
              freqmax,corners=3,zerophase=true).x

# test filtering SeisData
S = SeisData(2)
S.x[1] = rand(T,N)
S.fs[1] = fs
S.t[1] = [1 starttime * 1e6;length(C.x) 0]
S.x[2] = rand(T,N)
S.fs[2] = fs
S.t[2] = [1 starttime * 1e6;length(C.x) 0]
Snew = bandstop(S,freqmin,freqmax)
SInt = bandstop(S,Int(freqmin),Int(freqmax))
@test Snew.fs == S.fs
bandstop!(S,freqmin,freqmax)
@test all(S.fs .== fs)
@test S.x[1] == Snew.x[1]
@test S.x[2] == Snew.x[2]
@test S.x[1] == SInt.x[1]
@test S.x[2] == SInt.x[2]

# test zerophase -> there is no option for zerophase here
A = rand(T,N,Nwin)
S.x[1] = rand(T,N)
S.x[2] = rand(T,N)
@test bandstop(S,freqmin,freqmax,corners=corners,zerophase=true).x[1] == bandstop(S,freqmin,
               freqmax,corners=corners,zerophase=false).x[1]
@test bandstop(S,freqmin,freqmax,corners=corners,zerophase=true).x[2] == bandstop(S,freqmin,
              freqmax,corners=corners,zerophase=false).x[2]

# test corners
@test bandstop(S,freqmin,freqmax,corners=4,zerophase=true).x[1] != bandstop(S,freqmin,
              freqmax,corners=3,zerophase=true).x[1]
@test bandstop(S,freqmin,freqmax,corners=4,zerophase=true).x[2] != bandstop(S,freqmin,
            freqmax,corners=3,zerophase=true).x[2]

# test bandstop with Int inputs
A = rand(T,N,Nwin)
B = bandstop(A,Int(freqmin),Int(freqmax),Int(fs),corners=corners,zerophase=zerophase)
# check size and type stability
@test size(B) == size(A)
@test eltype(B) == eltype(A)
bandstop!(A,freqmin,freqmax,fs,corners=corners,zerophase=zerophase)
@test A == B

# NEEDED: test GPUArray
end

# test bandpass
@testset "lowpass" begin

# create white noise
A = rand(T,N,Nwin)
B = lowpass(A,freqmax,fs,corners=corners,zerophase=zerophase)
# check size and type stability
@test size(B) == size(A)
@test eltype(B) == eltype(A)
lowpass!(A,freqmax,fs,corners=corners,zerophase=zerophase)
@test A == B

# test power in pass bands
freqs = FFTW.rfftfreq(N,fs)
lowband = findall(0 + eps(T) .< freqs .< freqmax)
highband = findall(freqs .>= freqmax)
F = rfft(A,1)
Fpow = abs.(F)
# test filter
@test all(median(Fpow[lowband,:],dims=1) .> median(Fpow[highband,:],dims=1))

# test zerophase
A = rand(T,N,Nwin)
@test lowpass(A,freqmax,fs,corners=corners,zerophase=true) != lowpass(A,
               freqmax,fs,corners=corners,zerophase=false)

# test corners
@test lowpass(A,freqmax,fs,corners=4,zerophase=true) != lowpass(A,
               freqmax,fs,corners=3,zerophase=true)

# test low freq
highfreq = fs / 2
B = lowpass(A,highfreq,fs,corners=4,zerophase=true)
@test B == A
lowpass!(A,highfreq,fs,corners=4,zerophase=true)
@test B == A

# test RawData filtering
R = RawData()
R.x = A
R.fs = fs
Rnew = lowpass(R,freqmax)
@test Rnew.fs == R.fs
@test Rnew.freqmax == freqmax
lowpass!(R,freqmax)
@test R.fs == R.fs
@test R.freqmax == freqmax
@test R.x == Rnew.x

# test zerophase
A = rand(T,N,Nwin)
R.x = A
@test lowpass(R,freqmax,corners=corners,zerophase=true).x != lowpass(R,
               freqmax,corners=corners,zerophase=false).x

# test corners
@test lowpass(R,freqmax,corners=4,zerophase=true).x != lowpass(R,
              freqmax,corners=3,zerophase=true).x

# test CorrData filtering
C = CorrData()
C.corr = A
C.fs = fs
Cnew = lowpass(C,freqmax)
CInt = lowpass(C,Int(freqmax))
@test Cnew.fs == C.fs
@test Cnew.freqmax == freqmax
@test CInt.freqmax == freqmax
lowpass!(C,freqmax)
@test C.fs == C.fs
@test C.freqmax == freqmax
@test C.corr == Cnew.corr
@test C.corr == CInt.corr

# test zerophase
A = rand(T,N,Nwin)
C.corr = A
@test lowpass(C,freqmax,corners=corners,zerophase=true).corr != lowpass(C,
               freqmax,corners=corners,zerophase=false).corr

# test corners
@test lowpass(C,freqmax,corners=4,zerophase=true).corr != lowpass(C,
              freqmax,corners=3,zerophase=true).corr

# test filtering SeisChannel
C = SeisChannel()
C.x = rand(T,N)
C.fs = fs
C.t = [1 starttime * 1e6;length(C.x) 0]
Cnew = lowpass(C,freqmax)
CInt = lowpass(C,Int(freqmax))
@test Cnew.fs == fs
lowpass!(C,freqmax)
@test C.fs == C.fs
@test C.x == Cnew.x
@test C.x == CInt.x

# test zerophase -> there is no option for zerophase here
A = rand(T,N)
C.x = A
@test lowpass(C,freqmax,corners=corners,zerophase=true).x == lowpass(C,
               freqmax,corners=corners,zerophase=false).x

# test corners
@test lowpass(C,freqmax,corners=4,zerophase=true).x != lowpass(C,
              freqmax,corners=3,zerophase=true).x

# test filtering SeisData
S = SeisData(2)
S.x[1] = rand(T,N)
S.fs[1] = fs
S.t[1] = [1 starttime * 1e6;length(C.x) 0]
S.x[2] = rand(T,N)
S.fs[2] = fs
S.t[2] = [1 starttime * 1e6;length(C.x) 0]
Snew = lowpass(S,freqmax)
SInt = lowpass(S,Int(freqmax))
@test Snew.fs == S.fs
lowpass!(S,freqmax)
@test all(S.fs .== fs)
@test S.x[1] == Snew.x[1]
@test S.x[2] == Snew.x[2]
@test S.x[1] == SInt.x[1]
@test S.x[2] == SInt.x[2]

# test zerophase -> there is no option for zerophase here
A = rand(T,N,Nwin)
S.x[1] = rand(T,N)
S.x[2] = rand(T,N)
@test lowpass(S,freqmax,corners=corners,zerophase=true).x[1] == lowpass(S,
               freqmax,corners=corners,zerophase=false).x[1]
@test lowpass(S,freqmax,corners=corners,zerophase=true).x[2] == lowpass(S,
              freqmax,corners=corners,zerophase=false).x[2]

# test corners
@test lowpass(S,freqmax,corners=4,zerophase=true).x[1] != lowpass(S,
              freqmax,corners=3,zerophase=true).x[1]
@test lowpass(S,freqmax,corners=4,zerophase=true).x[2] != lowpass(S,
            freqmax,corners=3,zerophase=true).x[2]

# test lowpass with Int inputs
A = rand(T,N,Nwin)
B = lowpass(A,Int(freqmax),Int(fs),corners=corners,zerophase=zerophase)
# check size and type stability
@test size(B) == size(A)
@test eltype(B) == eltype(A)
lowpass!(A,freqmax,fs,corners=corners,zerophase=zerophase)
@test A == B

# NEEDED: test GPUArray
end

@testset "highpass" begin

# create white noise
A = rand(T,N,Nwin)
B = highpass(A,freqmin,fs,corners=corners,zerophase=zerophase)
# check size and type stability
@test size(B) == size(A)
@test eltype(B) == eltype(A)
highpass!(A,freqmin,fs,corners=corners,zerophase=zerophase)
@test A == B

# test power in pass bands
freqs = FFTW.rfftfreq(N,fs)
lowband = findall(0 + eps(T) .< freqs .< freqmax)
highband = findall(freqs .>= freqmax)
F = rfft(A,1)
Fpow = abs.(F)
# test filter
@test all(median(Fpow[lowband,:],dims=1) .< median(Fpow[highband,:],dims=1))

# test zerophase
A = rand(T,N,Nwin)
@test highpass(A,freqmin,fs,corners=corners,zerophase=true) != highpass(A,
               freqmin,fs,corners=corners,zerophase=false)

# test corners
@test highpass(A,freqmin,fs,corners=4,zerophase=true) != highpass(A,
               freqmin,fs,corners=3,zerophase=true)

# test low freq
highfreq = fs / 2
@test_throws ArgumentError highpass!(A,highfreq,fs,corners=4,zerophase=true)

# test RawData filtering
R = RawData()
R.x = A
R.fs = fs
Rnew = highpass(R,freqmin)
RInt = highpass(R,Int(freqmin))
@test Rnew.fs == R.fs
@test Rnew.freqmin == freqmin
highpass!(R,freqmin)
@test R.fs == R.fs
@test R.freqmin == freqmin
@test RInt.freqmin == freqmin
@test R.x == Rnew.x
@test R.x == RInt.x

# test zerophase
A = rand(T,N,Nwin)
R.x = A
@test highpass(R,freqmin,corners=corners,zerophase=true).x != highpass(R,
               freqmin,corners=corners,zerophase=false).x

# test corners
@test highpass(R,freqmin,corners=4,zerophase=true).x != highpass(R,
              freqmin,corners=3,zerophase=true).x

# test CorrData filtering
C = CorrData()
C.corr = A
C.fs = fs
Cnew = highpass(C,freqmin)
CInt = highpass(C,Int(freqmin))
@test Cnew.fs == C.fs
@test Cnew.freqmin == freqmin
@test CInt.freqmin == freqmin
highpass!(C,freqmin)
@test C.fs == C.fs
@test C.freqmin == freqmin
@test C.corr == Cnew.corr
@test C.corr == CInt.corr

# test zerophase
A = rand(T,N,Nwin)
C.corr = A
@test highpass(C,freqmin,corners=corners,zerophase=true).corr != highpass(C,
               freqmin,corners=corners,zerophase=false).corr

# test corners
@test highpass(C,freqmin,corners=4,zerophase=true).corr != highpass(C,
              freqmin,corners=3,zerophase=true).corr

# test filtering SeisChannel
C = SeisChannel()
C.x = rand(T,N)
C.fs = fs
C.t = [1 starttime * 1e6;length(C.x) 0]
Cnew = highpass(C,freqmin)
CInt = highpass(C,Int(freqmin))
@test Cnew.fs == fs
highpass!(C,freqmin)
@test C.fs == C.fs
@test C.x == Cnew.x
@test C.x == CInt.x

# test zerophase -> there is no option for zerophase here
A = rand(T,N)
C.x = A
@test highpass(C,freqmin,corners=corners,zerophase=true).x == highpass(C,
               freqmin,corners=corners,zerophase=false).x

# test corners
@test highpass(C,freqmin,corners=4,zerophase=true).x != highpass(C,
              freqmin,corners=3,zerophase=true).x

# test filtering SeisData
S = SeisData(2)
S.x[1] = rand(T,N)
S.fs[1] = fs
S.t[1] = [1 starttime * 1e6;length(C.x) 0]
S.x[2] = rand(T,N)
S.fs[2] = fs
S.t[2] = [1 starttime * 1e6;length(C.x) 0]
Snew = highpass(S,freqmin)
SInt = highpass(S,Int(freqmin))
@test Snew.fs == S.fs
highpass!(S,freqmin)
@test all(S.fs .== fs)
@test S.x[1] == Snew.x[1]
@test S.x[2] == Snew.x[2]
@test S.x[1] == SInt.x[1]
@test S.x[2] == SInt.x[2]

# test zerophase -> there is no option for zerophase here
A = rand(T,N,Nwin)
S.x[1] = rand(T,N)
S.x[2] = rand(T,N)
@test highpass(S,freqmin,corners=corners,zerophase=true).x[1] == highpass(S,
               freqmin,corners=corners,zerophase=false).x[1]
@test highpass(S,freqmin,corners=corners,zerophase=true).x[2] == highpass(S,
              freqmin,corners=corners,zerophase=false).x[2]

# test corners
@test highpass(S,freqmin,corners=4,zerophase=true).x[1] != highpass(S,
              freqmin,corners=3,zerophase=true).x[1]
@test highpass(S,freqmin,corners=4,zerophase=true).x[2] != highpass(S,
            freqmin,corners=3,zerophase=true).x[2]

# test highpass with Int inputs
A = rand(T,N,Nwin)
B = highpass(A,Int(freqmin),Int(fs),corners=corners,zerophase=zerophase)
# check size and type stability
@test size(B) == size(A)
@test eltype(B) == eltype(A)
highpass!(A,freqmin,fs,corners=corners,zerophase=zerophase)
@test A == B

# NEEDED: test GPUArray
end

# test envelope
@testset "envelope" begin
t = collect(0.:π/100:10π)
t = T.(t)
A = sin.(t) .* cos.(5 .* t)
up, down = envelope(A)
@test size(up) == size(A)
@test size(down) == size(A)
@test eltype(up) == eltype(A)
@test eltype(down) == eltype(A)
@test all(up .- A .> -1e-6)
@test all(A .- down .> -1e-6)

end


text:
using LPSjl
using Test

@testset "LPSjl.jl" begin

end


text:
using DSP: conv

@doc raw"""
    synthesis(wt::WaveletTransformed)

Implements the continuous wavelet synthesis on some `WaveletTransformed` object. That is, it
reconstructs the signal, `wt` was created from. The optional keyword argument `hilbert` specifies
whether the complex Hilbert transformed of the original signal or the real signal is returned. By
default the real reconstruction is returned.
"""
function synthesis(wt::ContinuousWavelet.ContinuousWaveletTransformed; hilbert=false)
    N, M = size(wt)
    res = zeros(Complex{Float64}, N)
    # First convolution with smallest scale
    K = Int32(exp2(ceil(log2(2*cutoff_time(wt.wavelet)*wt.scales[1])))); off = K÷2
    kern = map(t -> eval_synthesis(wt.wavelet,
        (t-K/2)/wt.scales[1])/wt.scales[1]/wt.scales[1], 1:K)
    last = conv(kern, wt[:,1])[(off):(off+N-1)]

    current = similar(last)
    for i in 2:M
        # eval kernel for current scale
        K = Int32(exp2(ceil(log2(2*cutoff_time(wt.wavelet)*wt.scales[i])))); off = K÷2
        kern = map(t -> eval_synthesis(wt.wavelet,
            (t-K/2)/wt.scales[i])/wt.scales[i]/wt.scales[i], 1:K)
        # convolution of i-th voice with kernel at i-th scale
        current[:] = conv(kern, wt[:,i])[(off):(off+N-1)]
        # integrate by mid-point evaluation
        res[:] += ( (current+last) .* ((wt.scales[i]-wt.scales[i-1]))/2 )
        # swap current <-> last w/o copy
        tmp = last;
        last = current;
        current = tmp;
    end
    # If Hilbert transformed result is requested -> complex otherwise -> real
    (hilbert) ? res : 2*real.(res)
end


text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at https://mozilla.org/MPL/2.0/.

"""
    default_container(indices)

If `indices` is a [`NestedIterator`](@ref), return a
[`SparseAxisArray`](@ref). Otherwise, `indices` should be
a `VectorizedProductIterator` and the function returns
`Array` if all iterators of the product are `Base.OneTo` and retunrs
[`DenseAxisArray`](@ref) otherwise.
"""
function default_container end

"""
    container(f::Function, indices, ::Type{C})

Create a container of type `C` with indices `indices` and values at given
indices given by `f`.

    container(f::Function, indices)

Create a container with indices `indices` and values at given indices given by
`f`. The type of container used is determined by [`default_container`](@ref).

## Examples

```@jldoctest; setup=:(using JuMP)
julia> Containers.container((i, j) -> i + j, Containers.vectorized_product(Base.OneTo(3), Base.OneTo(3)))
3×3 Array{Int64,2}:
 2  3  4
 3  4  5
 4  5  6

julia> Containers.container((i, j) -> i + j, Containers.vectorized_product(1:3, 1:3))
2-dimensional DenseAxisArray{Int64,2,...} with index sets:
    Dimension 1, 1:3
    Dimension 2, 1:3
And data, a 3×3 Array{Int64,2}:
 2  3  4
 3  4  5
 4  5  6

julia> Containers.container((i, j) -> i + j, Containers.vectorized_product(2:3, Base.OneTo(3)))
2-dimensional DenseAxisArray{Int64,2,...} with index sets:
    Dimension 1, 2:3
    Dimension 2, Base.OneTo(3)
And data, a 2×3 Array{Int64,2}:
 3  4  5
 4  5  6

julia> Containers.container((i, j) -> i + j, Containers.nested(() -> 1:3, i -> i:3, condition = (i, j) -> isodd(i) || isodd(j)))
SparseAxisArray{Int64,2,Tuple{Int64,Int64}} with 5 entries:
  [1, 2]  =  3
  [2, 3]  =  5
  [3, 3]  =  6
  [1, 1]  =  2
  [1, 3]  =  4
```
"""
function container end

function container(f::Function, indices)
    return container(f, indices, default_container(indices))
end

const ArrayIndices{N} = VectorizedProductIterator{NTuple{N,Base.OneTo{Int}}}
default_container(::ArrayIndices) = Array
function container(f::Function, indices::ArrayIndices, ::Type{Array})
    return map(I -> f(I...), indices)
end

function _oneto(indices::AbstractVector{<:Integer})
    if indices == 1:length(indices)
        return Base.OneTo(length(indices))
    end
    return error("Index set for array is not one-based interval.")
end

function _oneto(::Any)
    return error("Index set for array is not one-based interval.")
end

_oneto(indices::Base.OneTo) = indices

function container(
    f::Function,
    indices::VectorizedProductIterator,
    ::Type{Array},
)
    return container(
        f,
        vectorized_product(_oneto.(indices.prod.iterators)...),
        Array,
    )
end
default_container(::VectorizedProductIterator) = DenseAxisArray
function container(
    f::Function,
    indices::VectorizedProductIterator,
    ::Type{DenseAxisArray},
)
    return DenseAxisArray(map(I -> f(I...), indices), indices.prod.iterators...)
end
default_container(::NestedIterator) = SparseAxisArray
# Returns the element type. If it is unknown but it is known to be `N`-tuples,
# returns `NTuple{N, Any}`.
_eltype_or_any(indices::Array) = eltype(indices)
function container(f::Function, indices, ::Type{SparseAxisArray})
    # Same as `map` but does not allocate the resulting vector.
    mappings = Base.Generator(I -> I => f(I...), indices)
    # Same as `Dict(mapping)` but it will error if two indices are the same.
    data = NoDuplicateDict(mappings)
    dict = data.dict
    if isempty(dict)
        # If `dict` is empty, it was not able to determine the type
        # of the key hence the type of `dict` is `Dict{Any, Any}`.
        # This is an issue since `SparseAxisArray` needs the key
        # type to be a tuple.
        dict = Dict{_eltype_or_any(indices),Any}()
    end
    return SparseAxisArray(dict)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@snn_kw struct IF2Parameter{FT=Float32} <: AbstractIFParameter
    Ee::FT = 0mV
    Ei::FT = 0mV
end

@snn_kw mutable struct IF2 <: AbstractIF
    param::IF2Parameter = IF2Parameter()
    N::Int32 = 100
end

"""
[Integrate-And-Fire Neuron](https://neuronaldynamics.epfl.ch/online/Ch1.S3.html)
"""
IF2

function integrate!(p::IF2, param::IF2Parameter, dt::Float32)
    @unpack N = p
    @unpack Ee, Ei = param
    @inbounds for i = 1:N
        v[i] += dt * (ge[i] * (Ee - v[i]) + gi[i] * (Ei - v[i]) - (v[i] - El)) / τm
        ge[i] += dt * -ge[i] / τe
        gi[i] += dt * -gi[i] / τi
    end
    @inbounds for i = 1:N
        fire[i] = v[i] > Vt
        v[i] = ifelse(fire[i], Vr, v[i])
    end
end


text:
# ---
# title: 1366. Rank Teams by Votes
# id: problem1366
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Sort
# link: <https://leetcode.com/problems/rank-teams-by-votes/description/>
# hidden: true
# ---
#
# In a special ranking system, each voter gives a rank from highest to lowest to
# all teams participated in the competition.
#
# The ordering of teams is decided by who received the most position-one votes.
# If two or more teams tie in the first position, we consider the second
# position to resolve the conflict, if they tie again, we continue this process
# until the ties are resolved. If two or more teams are still tied after
# considering all positions, we rank them alphabetically based on their team
# letter.
#
# Given an array of strings `votes` which is the votes of all voters in the
# ranking systems. Sort all teams according to the ranking system described
# above.
#
# Return _a string of all teams_ **sorted** by the ranking system.
#
#
#
# **Example 1:**
#
#
#
#     Input: votes = ["ABC","ACB","ABC","ACB","ACB"]
#     Output: "ACB"
#     Explanation: Team A was ranked first place by 5 voters. No other team was voted as first place so team A is the first team.
#     Team B was ranked second by 2 voters and was ranked third by 3 voters.
#     Team C was ranked second by 3 voters and was ranked third by 2 voters.
#     As most of the voters ranked C second, team C is the second team and team B is the third.
#
#
# **Example 2:**
#
#
#
#     Input: votes = ["WXYZ","XYZW"]
#     Output: "XWYZ"
#     Explanation: X is the winner due to tie-breaking rule. X has same votes as W for the first position but X has one vote as second position while W doesn't have any votes as second position.
#
#
# **Example 3:**
#
#
#
#     Input: votes = ["ZMNAGUEDSJYLBOPHRQICWFXTVK"]
#     Output: "ZMNAGUEDSJYLBOPHRQICWFXTVK"
#     Explanation: Only one voter so his votes are used for the ranking.
#
#
# **Example 4:**
#
#
#
#     Input: votes = ["BCA","CAB","CBA","ABC","ACB","BAC"]
#     Output: "ABC"
#     Explanation:
#     Team A was ranked first by 2 voters, second by 2 voters and third by 2 voters.
#     Team B was ranked first by 2 voters, second by 2 voters and third by 2 voters.
#     Team C was ranked first by 2 voters, second by 2 voters and third by 2 voters.
#     There is a tie and we rank teams ascending by their IDs.
#
#
# **Example 5:**
#
#
#
#     Input: votes = ["M","M","M","M"]
#     Output: "M"
#     Explanation: Only team M in the competition so it has the first rank.
#
#
#
#
# **Constraints:**
#
#   * `1 <= votes.length <= 1000`
#   * `1 <= votes[i].length <= 26`
#   * `votes[i].length == votes[j].length` for `0 <= i, j < votes.length`.
#   * `votes[i][j]` is an English **upper-case** letter.
#   * All characters of `votes[i]` are unique.
#   * All the characters that occur in `votes[0]` **also  occur** in `votes[j]` where `1 <= j < votes.length`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# #=##############################################################################################
# Filename: couple_vlm_bem.jl
# Author: Ryan Anderson
# Contact: rymanderson@gmail.com
# README: `<: Action` function calculates aerodynamic forces on a blade-element momentum-vortex lattice system
# =###############################################################################################

# """
# couple_vlm_bem <: Action

# Solves for the aerodynamic force distribution at each step.

# NOTE: THIS ACTION DOES NOT SOLVE THE VORTEX LATTICE. Call `solve_wing_CF_CM` prior to calling this action.

# Inputs:

# * `aircraft::Aircraft` : `Aircraft` system object
# * `parameters<:Parameters` `Parameters` struct
# * `freestream::Freestream` : `Freestream` object
# * `environment::Environment` `Environment` object
# * `step_range::AbstractArray` : array of steps for which the simulation is run
# * `stepi::Int` : index of the current step
# * `step_symbol::String` : defines the step, e.g. `alpha` or `time`

# `parameters <: Parameters` requires the following elements:

# * `cfs::Vector{Array{Float64,2}}` : vector with length equal to the number of lifting surfaces, each member containing an array of size (3,nspanwisepanels) of force coefficients cd, cy, cl
# * `cms::Vector{Array{Float64,2}}` : vector with length equal to the number of lifting surfaces, each member containing an array of size (3,nspanwisepanels) of moment coefficients cmx, cmy, cmz
# * `cls::Vector{Array{Float64,2}}` : each element is an array of size (nspanwisepanels, nsteps) containing local lift coefficients at each lifting line section, corresponding to each lifting surface
# * `cds::Vector{Array{Float64,2}}` : each element is an array of size (nspanwisepanels, nsteps) containing local drag coefficients at each lifting line section, corresponding to each lifting surface
# * `cys::Vector{Array{Float64,2}}` : each element is an array of size (nspanwisepanels, nsteps) containing local side force coefficients at each lifting line section, corresponding to each lifting surface
# * `cmxs::Vector{Array{Float64,2}}` : each element is an array of size (nspanwisepanels, nsteps) containing local x-axis (roll) moment coefficients at each lifting line section, corresponding to each lifting surface
# * `cmys::Vector{Array{Float64,2}}` : each element is an array of size (nspanwisepanels, nsteps) containing local y-axis (pitch) moment coefficients at each lifting line section, corresponding to each lifting surface
# * `cmzs::Vector{Array{Float64,2}}` : each element is an array of size (nspanwisepanels, nsteps) containing local z-axis (yaw) moment force coefficients at each lifting line section, corresponding to each lifting surface
# * `omegas::Vector{Float64}` : a vector of length `length(step_range)` containing a vector of rotational velocities for each rotor
# * `Js::Vector{Vector{Float64}}` : a vector of length `length(step_range)` containing a vector of advance ratios for each rotor
# * `CTs::Vector{Vector{Float64}}` : a vector of length `length(step_range)` containing a vector of thrust coefficients for each rotor
# * `CQs::Vector{Vector{Float64}}` : a vector of length `length(step_range)` containing a vector of torque coefficients for each rotor
# * `ηs::Vector{Vector{Float64}}` : a vector of length `length(step_range)` containing a vector of propulsive efficiencies for each rotor

# """
# function couple_vlm_bem(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol)
#     # extract lift and moment distribution
#     cfs, cms = VL.lifting_line_coefficients(aircraft.wingsystem.system, aircraft.wingsystem.lifting_line_rs, aircraft.wingsystem.lifting_line_chords)
#     # store to `parameters`
#     nwings = length(aircraft.wingsystem.system.surfaces)
#     for iwing in 1:nwings
#         parameters.cls[iwing][:,stepi] = cfs[iwing][3,:]
#         parameters.cds[iwing][:,stepi] = cfs[iwing][1,:]
#         parameters.cys[iwing][:,stepi] = cfs[iwing][2,:]
#         parameters.cmxs[iwing][:,stepi] = cms[iwing][1,:]
#         parameters.cmys[iwing][:,stepi] = cms[iwing][2,:]
#         parameters.cmzs[iwing][:,stepi] = cms[iwing][3,:]
#     end
#     return false
# end

# """
# lift_moment_distribution(system, step_range)

# Method returns initialized elements required for the `parameters <: Parameters` struct during simulation.

# Inputs:

# * `aircraft::Aircraft` : system to be simulated
# * `step_range::AbstractArray` : defines each step of the simulation

# Outputs:

# * `cls::Vector{Array{Float64,2}}` : each element is an array of size (nspanwisepanels, nsteps) containing local lift coefficients at each lifting line section, corresponding to each lifting surface
# * `cds::Vector{Array{Float64,2}}` : each element is an array of size (nspanwisepanels, nsteps) containing local drag coefficients at each lifting line section, corresponding to each lifting surface
# * `cys::Vector{Array{Float64,2}}` : each element is an array of size (nspanwisepanels, nsteps) containing local side force coefficients at each lifting line section, corresponding to each lifting surface
# * `cmxs::Vector{Array{Float64,2}}` : each element is an array of size (nspanwisepanels, nsteps) containing local x-axis (roll) moment coefficients at each lifting line section, corresponding to each lifting surface
# * `cmys::Vector{Array{Float64,2}}` : each element is an array of size (nspanwisepanels, nsteps) containing local y-axis (pitch) moment coefficients at each lifting line section, corresponding to each lifting surface
# * `cmzs::Vector{Array{Float64,2}}` : each element is an array of size (nspanwisepanels, nsteps) containing local z-axis (yaw) moment force coefficients at each lifting line section, corresponding to each lifting surface

# """
# function lift_moment_distribution(aircraft, step_range)

#     nwings = length(aircraft.wingsystem.system.surfaces)
#     nspanwisepanels = [size(surface)[2] for surface in aircraft.wingsystem.system.surfaces]
#     cls = [zeros(nspanwisepanels[i], length(step_range)) for i in 1:nwings]
#     cds = deepcopy(cls)
#     cys = deepcopy(cls)
#     cmxs = deepcopy(cls)
#     cmys = deepcopy(cls)
#     cmzs = deepcopy(cls)

#     return cls, cds, cys, cmxs, cmys, cmzs
# end

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
module TestMLJParameters

using MLJBase
using Test

import LightGBM

parameters = Dict(
    :boosting => "goss",
    :drop_rate => 0.6,
    :max_drop => 100,
    :skip_drop => 0.9,
    :xgboost_dart_mode => true,
    :uniform_drop => true,
    :drop_seed => 10,
    :top_rate => 0.7,
    :other_rate => 0.2,
)

# The tests for regressor and classifier is to check
# if the MLJ parameters are present in the LGBMEstimators

@testset "MLJ parameters -- regressor" begin

    X, y = @MLJBase.load_boston;

    # Passed the dictionary of parameters to the regressor model
    model = LightGBM.MLJInterface.LGBMRegressor(;parameters...)

    # Extract the estimator
    fit_result, _, _ = MLJBase.fit(model, 0, X, y)
    estimator, _, _ = fit_result

    for (k, v) in parameters
        @test getproperty(estimator, k) == v
    end

end

@testset "MLJ parameters -- classifier" begin

    X, y = @MLJBase.load_iris;
    model = LightGBM.MLJInterface.LGBMClassifier(;parameters...)
    fit_result, _, _ = MLJBase.fit(model, 0, X, y)
    estimator, _, _ = fit_result

    for (k, v) in parameters
        @test getproperty(estimator, k) == v
    end

end

end #end module


text:
using Documenter
using SeisTest
using Pkg
Pkg.add("SeisProcessing")
Pkg.add("DSP")
Pkg.add("Plots")
makedocs(
    sitename = "SeisTest",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing)== "true"),
    modules = [SeisTest],
	pages = [
          "Home" => "index.md",
		  "Library" => Any[
 	 			"Public" => "public.md",
 	 ],
	      "Example"=>"example.md",
		  "Publication"=>"publication.md"
		],
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(repo = "github.com/Ji-seismic/SeisTest.git")

max_stars_repo_name: Ji-seismic/SeisTest

max_issues_repo_name: Ji-seismic/SeisTest

max_forks_repo_name: Ji-seismic/SeisTest


text:
using Choices
using Test
using Base.Iterators: take, countfrom

# test iterating over the first few emirp numbers
# https://en.wikipedia.org/wiki/Emirp

# first some common definitions:
divides(x::Integer, y::Integer) = isinteger(x / y)
isprime(x::Integer) = !any(y -> divides(x, y), 2:trunc(Integer, √x))
reverse_digits(x::Int) = parse(Int, reverse(string(x)))

# the conventional way to do it with a Generator:
function isemirp(x::Int)
  isprime(x) || return false
  rev_x = reverse_digits(x)
  return (x != rev_x && isprime(rev_x))
end
generated_emirps = [take((x for x ∈ countfrom() if isemirp(x)), 10)...]

# a way to do it with Choices:
emirps_chooser =
  choose(countfrom()) |> x ->
  (isprime(x) ? reverse_digits(x) : fail()) |> rev_x ->
  ((x != rev_x && isprime(rev_x)) ? x : fail())
chosen_emirps = [take(emirps_chooser, 10)...]
cut()

expected_emirps = [13, 17, 31, 37, 71, 73, 79, 97, 107, 113]
@test generated_emirps == expected_emirps
@test chosen_emirps == expected_emirps


text:
#=
Visualisation of a the central limit theorem convergence, (https://en.wikipedia.org/wiki/Central_limit_theorem).

The mean of n samples from "any" Independent Identically Distributed random variable converges towards a gaussian for n -> Inf
=#

using Javis
using Colors

# Following modules don't belong to Javis dependencies so to run them
# one needs to add them to the enviroment used to run this script
using StatsBase
using Distributions

"""
    ground(args...)

Function defining the background of the animation.
"""
function ground(args...)
    background("black")
    sethue("white")
end

"""
    fixed_gaussian(loc_hist; color)

Returns an anonymous function to be given to barchart as labelfunction.
It takes care of normalizing the values of the limit gaussian distribution according
to the ```loc_hist``` histogram and render these values as circles. Also draws the
bottom tickline.
"""
function fixed_gaussian(loc_hist; color)
    (values, i, lowpos, highpos, barwidth, scaledvalue) -> begin

        minvalue, maxvalue = extrema(loc_hist.weights)
        barchartheight = boxheight(boundingbox) - 2margin
        minbarrange = minvalue - abs(minvalue)
        maxbarrange = maxvalue + abs(maxvalue)
        @layer begin

            sethue(color)
            if i <= length(loc_hist.edges[1])
                scaledgaussvalue =
                    rescale(
                        pdf(
                            gauss,
                            loc_hist.edges[1][i + 1] - loc_hist.edges[1].step.hi / 2,
                        ),
                        minbarrange,
                        maxbarrange,
                    ) * barchartheight
                circle(lowpos - (0, scaledgaussvalue), 3, :fill)
            end
        end

        tickline(
            boxbottomleft(boundingbox) - Point(0, margin),
            boxbottomright(boundingbox) - Point(0, margin),
            startnumber = loc_hist.edges[1][1],
            finishnumber = loc_hist.edges[1][end],
        )
    end
end

"""
    hist_bar

Utility function only used to change the histogram bars color.
"""
function hist_bar(; color)
    (values, i, lowpos, highpos, barwidth, scaledvalue) -> begin
        @layer begin
            sethue(color)
            Luxor.setline(barwidth)
            line(lowpos, highpos, :stroke)
        end
    end
end

# This can be changed to any distribution from Distribtuions.jl
dist = Bernoulli(0.3)

# Parameters for the barchart function
boundingbox = BoundingBox(O + (-250, -120), O + (250, 120))
margin = 5

# number of frames
n_frames = 700

# Number of histograms shown. To make the animation
# slower set this to a lower value, less histograms
# will be shown but they will reach the same n for
# convergence as set by n_frames
n_hists = 700

# number of samples at each n
n_samples = 10000

# Change to adjust bar colors and gauss colors
barcolor = HSV(colorant"goldenrod1")
gausscolor = HSV(240, barcolor.s, barcolor.v)

# Sample the mean r.v. for increasing values of n (1 to n_hists)
samples = map(1:n_hists) do n
    map(1:n_samples) do _
        (mean(rand(dist, n)) - mean(dist)) * sqrt(n)
    end
end

finalmin, finalmax = extrema(samples[end])

# Turn the samples into histograms
steps = map(samples) do sampling
    hist = fit(
        Histogram,
        sampling,
        weights(ones(length(sampling))),
        range(finalmin, finalmax, length = 100),
    )
    hist = StatsBase.normalize(hist, mode = :pdf)
    hist
end

# Define the gaussian distribution where they should converge
gauss = Normal(0, StatsBase.std(dist))

my_video = Video(600, 500)
Background(1:n_frames, ground)

step_size = n_frames ÷ n_hists
frame_brakes = 1:step_size:n_frames

# Fix some point where writings will be shown
titlepoint = Point(0, -100)
distpoint = Point(-200, 0)
counterpoint = Point(200, 0)

# The last and thus hopefully most closely converged histogram in our sequence
final_hist = steps[end]

for (frame_n, hist) in zip(frame_brakes, steps)

    Object(
        frame_n:(frame_n + step_size - 1),
        @JShape begin
            barchart(
                hist.weights,
                boundingbox = boundingbox,
                labels = true,

                # Provide the hist_bar we defined as the barfunction
                barfunction = hist_bar(color = barcolor),

                # Provide the fixed_gaussian we defined as the labelfunction
                # it will plot the gaussian dots and the ticks on the bottom
                labelfunction = fixed_gaussian(final_hist, color = gausscolor),
            )
        end
    )

    # The counter digits
    Object(
        frame_n:(frame_n + step_size - 1),
        @JShape begin
            sethue(barcolor)
            fontsize(15)
            text(string(frame_n ÷ step_size), counterpoint, halign = :center)
        end
    )
end

# All the writings except the changing digit.
Object(
    1:n_frames,
    @JShape begin
        fontsize(40)
        text("Central Limit Theorem", titlepoint, halign = :center)

        fontsize(20)
        label("N", :N, counterpoint, offset = 20)
        label("Distribution", :N, distpoint, offset = 20)

        sethue(barcolor)
        fontsize(15)
        text(
            # Gather the distribution name and parameters
            # only works if dist is from Distributions.jl
            join([string(typeof(dist).name.name); string(params(dist))], ""),
            distpoint,
            halign = :center,
        )
    end
)

render(my_video, pathname = "central_limit_theorem.gif")


text:
using DifferentialEquations

@testset "DifferentialEquations solve API using an initial-value problem"
    prob = @ivp(x' = 1.0x, x(0) ∈ [1/2])

    sol = DifferentialEquations.solve(prob, tspan=(0.0, 1.0))
    @test sol isa OrdinaryDiffEq.ODECompositeSolution
    @test first(sol.t) == 0.0 && last(sol.t) == 1.0

    sol = DifferentialEquations.solve(prob, T=1.0)
    @test sol isa OrdinaryDiffEq.ODECompositeSolution
    @test first(sol.t) == 0.0 && last(sol.t) == 1.0

    # inplace (default)
    sol = DifferentialEquations.solve(prob, T=1.0, inplace=true)

    # out-of-place
    sol = DifferentialEquations.solve(prob, T=1.0, inplace=false)
end


text:
using DynamicPPL
using AbstractMCMC
using AbstractPPL
using Bijectors
using Distributions
using DistributionsAD
using Documenter
using ForwardDiff
using MacroTools
using MCMCChains
using Tracker
using Zygote
using Setfield

using Distributed
using LinearAlgebra
using Pkg
using Random
using Serialization
using Test

using DynamicPPL: getargs_dottilde, getargs_tilde, Selector

const DIRECTORY_DynamicPPL = dirname(dirname(pathof(DynamicPPL)))
const DIRECTORY_Turing_tests = joinpath(DIRECTORY_DynamicPPL, "test", "turing")
const GROUP = get(ENV, "GROUP", "All")

Random.seed!(100)

include("test_util.jl")

@testset "DynamicPPL.jl" begin
    if GROUP == "All" || GROUP == "DynamicPPL"
        @testset "interface" begin
            include("utils.jl")
            include("compiler.jl")
            include("varinfo.jl")
            include("simple_varinfo.jl")
            include("model.jl")
            include("sampler.jl")
            include("prob_macro.jl")
            include("independence.jl")
            include("distribution_wrappers.jl")
            include("contexts.jl")
            include("context_implementations.jl")

            include("threadsafe.jl")

            include("serialization.jl")

            include("loglikelihoods.jl")
        end

        @testset "compat" begin
            include(joinpath("compat", "ad.jl"))
        end

        @testset "doctests" begin
            DocMeta.setdocmeta!(
                DynamicPPL, :DocTestSetup, :(using DynamicPPL); recursive=true
            )
            doctestfilters = [
                # Older versions will show "0 element Array" instead of "Type[]".
                r"(Any\[\]|0-element Array{.+,[0-9]+})",
                # Older versions will show "Array{...,1}" instead of "Vector{...}".
                r"(Array{.+,\s?1}|Vector{.+})",
                # Older versions will show "Array{...,2}" instead of "Matrix{...}".
                r"(Array{.+,\s?2}|Matrix{.+})",
                # Errors from macros sometimes result in `LoadError: LoadError:`
                # rather than `LoadError:`, depending on Julia version.
                r"ERROR: (LoadError:\s)+",
            ]
            doctest(DynamicPPL; manual=false, doctestfilters=doctestfilters)
        end
    end

    if GROUP == "All" || GROUP == "Downstream"
        @testset "turing" begin
            try
                # activate separate test environment
                Pkg.activate(DIRECTORY_Turing_tests)
                Pkg.develop(PackageSpec(; path=DIRECTORY_DynamicPPL))
                Pkg.instantiate()

                # make sure that the new environment is considered `using` and `import` statements
                # (not added automatically on Julia 1.3, see e.g. PR #209)
                if !(joinpath(DIRECTORY_Turing_tests, "Project.toml") in Base.load_path())
                    pushfirst!(LOAD_PATH, DIRECTORY_Turing_tests)
                end

                include(joinpath("turing", "runtests.jl"))
            catch err
                err isa Pkg.Resolve.ResolverError || rethrow()
                # If we can't resolve that means this is incompatible by SemVer and this is fine
                # It means we marked this as a breaking change, so we don't need to worry about
                # Mistakenly introducing a breaking change, as we have intentionally made one
                @info "Not compatible with this release. No problem." exception = err
            end
        end
    end
end

max_stars_repo_name: TuringLang/MicroPPL

max_issues_repo_name: TuringLang/MicroPPL

max_forks_repo_name: TuringLang/MicroPPL


text:
using Flux
using CSV
using DataFrames
using Random
using MLDataUtils
using Statistics
using Printf
using BSON
include("dataloader.jl")

using CUDAapi
if has_cuda()
    @info "CUDA is on"
    import CuArrays
    CuArrays.allowscalar(false)
else
    @error "Couldn't find CUDA"
    exit()
end

# Parameter setup
csv_path = "" # ------------------------ TO BE SET
dataset_size = 10_000
batch_size = 64
epochs = 50
patience = 10
train_test_split = 0.8
random_seed = 2
Random.seed!(random_seed)

@info("Preparing the dataset")
dataset_csv = CSV.File(csv_path; header=["col_path", "col_depth", "col_hmo", "col_freqpeak", "col_thetapeak"])
df = DataFrame(dataset_csv)
df = df[shuffle(1:size(df, 1)),:]  # Shuffling
df = df[1:dataset_size, [:col_path, :col_depth]]  # Reducing

for i in 1:size(df, 1)
    # CSV doesn't have the file extension
    df.col_path[i] = string(df.col_path[i], ".h5")
end

# Make data generator
function load_dataset(dataset)
    x = dataset.col_path
    y = dataset.col_depth / 10
    return x,  convert(Array{Float32}, y)  # Y
end

x_data, y_data = load_dataset(df);
(x_train, y_train), (x_test, y_test) = splitobs((x_data, y_data), at = train_test_split)
trainloader = DataLoader((x_train, y_train), batchsize=batch_size, shuffle=true)
testloader = DataLoader((x_test, y_test), batchsize=batch_size, shuffle=true)

@info("Building the model")
model = Chain(
    Conv((3, 3), 4=>16, pad=(1, 1), relu),
    Conv((3, 3), 16=>32, pad=(1, 1), relu),
    x -> reshape(x, :, size(x, 4)),
    Dense(51200, 256, relu),
    Dense(256, 256, relu),
    Dense(256, 1, relu)
)

# Load onto GPU if available
@info("Loading onto GPU")
model = gpu(model)
# trainloader = gpu(trainloader)
# testloader = gpu(testloader)

opt = ADAM(1e-05, (0.99, 0.999))

function loss(x, y)
    return Flux.mse(model(x), y)
end

function evaluate(testloader)
    errors = []
    for (x, y) in testloader
        mse = Flux.mse(cpu(model(x)), cpu(y))
        append!(errors, mse)
    end
    return mean(errors)
end

@info("Training")
best_val = 999
last_improvement = 0
for epoch_idx in 1:epochs
    global best_val, last_improvement
    Flux.train!(loss, Flux.params(model), trainloader, opt)

    new_val = evaluate(testloader)
    @info(@sprintf("[%d]: Validation MSE: %.4f", epoch_idx, new_val))
    if new_val < best_val
        best_val = new_val
        last_improvement = epoch_idx
        filename = string(epoch_idx, "__", new_val, ".bson")
        @info(" -> New best! Saving model out to ", filename)
        BSON.@save joinpath(dirname(@__FILE__), filename) model epoch_idx new_val
    end

    # Learning rate scheduler
    if epoch_idx == 40 || epoch_idx == 60 || epoch_idx == 80
        opt.eta /= 10.0
        @info(@sprintf("Epoch [%d]: Dropping learning rate to %.1e", epoch_idx, opt.eta))
    end

    # Early stopping
    if epoch_idx - last_improvement > patience
        @info("Early stopping...")
        exit()
    end

end


text:
#
# AnnotatedArrays.jl -
#
# Objects that combine values stored in an array and properties stored in a
# dictionary or a named tuple.
#

module AnnotatedArrays

export
    AnnotatedArray,
    nkeys,
    properties

using ..ArrayTools.PseudoArrays

# Method `hasproperty` was introduced in Julia 1.2.
if isdefined(Base, :hasproperty)
    import Base: hasproperty
end

"""

`AnnotatedArrays.Properties` is the union of allowed types to store the
properties of `AnnotatedArray` instances.  It can be a `NamedTuple`, a
dictionary with symbolic keys or a dictionary indexed by strings.  The former
case provides the fastest access to properties but they are immutable, the
latter case provides the slowest access to properties.

If other types that symbols or strings are allowed, they must be restricted to
avoid ambiguities with array indices (and ranges, etc.).

For symbolic properties names (stored by a `NamedTuple` or an
`AbstractDict{Symbol}`), the properties can be accessed like ordinary fields.

"""
const Properties = Union{AbstractDict,NamedTuple}
const SymbolicProperties = Union{AbstractDict{Symbol},NamedTuple}

"""

`AbstractAnnotatedArray{T,N,P,S}` is the super-type of annotated arrays which
store values of type `T` as in a `N`-dimensional array with indexing style `S`
and properties in an object of type `P` (a dictionary or a named tuple).

Two useful aliases are defined:

- `DynamicallyAnnotatedArray{T,N,K,V,S}` for annotated arrays whose properties
  can be modified because they are stored in a dictionary with key type `K` and
  value type `V`.

- `StaticallyAnnotatedArray{T,N,S}` for annotated arrays which have immutable
  properties (because they are stored in a named tuple).  To build such
  annotated arrays, provide a named tuple as the last argument of the
  `AnnotatedArray` constructor:

  ```julia
  A = AnnotatedArray(arr, (key=val,))
  B = AnnotatedArray{T}(undef, (dim1, dim2, ...), (key1=val1, key2=val2, ...))
  ```

  Do not forget to have a trailing comma if the named tuple has a single
  element (as for `A` in the above example).

The abstract type `AbstractAnnotatedArray{T,N,P,S}` is defined (but not
exported by `AnnotatedArrays`) so that other types than `AnnotatedArray` can be
derived with similar behavior.  For derived types, say `CustomType`, two
methods should be specialized:

```julia
using ArrayTools.AnnotatedArrays
Base.parent(obj::CustomType) -> the array storing the values of obj
AnnotatedArrays.properties(obj::CustomType) -> the object storing the properties of obj
```

"""
abstract type AbstractAnnotatedArray{T,N,P<:Properties,S} <: PseudoArray{T,N,S} end

const DynamicallyAnnotatedArray{T,N,K,V,S} = AbstractAnnotatedArray{T,N,<:AbstractDict{K,V},S}
@doc @doc(AbstractAnnotatedArray) DynamicallyAnnotatedArray

const StaticallyAnnotatedArray{T,N,S} = AbstractAnnotatedArray{T,N,<:NamedTuple,S}
@doc @doc(AbstractAnnotatedArray) StaticallyAnnotatedArray

struct AnnotatedArray{T,N,P<:Properties,
                      A<:AbstractArray{T,N},
                      S} <: AbstractAnnotatedArray{T,N,P,S}
    data::A
    prop::P
    function AnnotatedArray{T,N,P,A,S}(data::A,
                                       prop::P) where {T,N,P<:NamedTuple,
                                                       A<:AbstractArray{T,N},S}
        @assert IndexStyle(data) === S()
        return new{T,N,P,A,S}(data, prop)
    end
    function AnnotatedArray{T,N,P,A,S}(data::A,
                                       prop::P) where {T,N,K,V,
                                                       P<:AbstractDict{K,V},
                                                       A<:AbstractArray{T,N},S}
        @assert IndexStyle(data) === S()
        K === Any && error("key type must be more specialized than `Any`")
        K <: Integer && error("key type must not be an integer")
        K <: CartesianIndex && error("key type must not be a Cartesian index")
        K <: AbstractRange{<:Integer} && error("key type must not be an integer range")
        K <: Colon && error("key type must not be a colon")
        return new{T,N,P,A,S}(data, prop)
    end
end

# As expected by the PseudoArray interface, extend Base.parent() to return the
# array backing the storage of values.
@inline Base.parent(A::AnnotatedArray) = Base.getfield(A, :data)

# Generic outer constructor with given data and properties.
function AnnotatedArray(data::A, prop::P) where {T,N,P<:Properties,
                                                 A<:AbstractArray{T,N}}
    S = typeof(IndexStyle(data))
    return AnnotatedArray{T,N,P,A,S}(data, prop)
end

# Get rid of matching type-parameter when initial array is specified.  The
# pairs of definitions are needed to disentangle ambiguities.
AnnotatedArray{T,N}(data::AbstractArray{T,N}, args...; kwds...) where {T,N} =
    AnnotatedArray(data, args...; kwds...)
AnnotatedArray{T,N}(data::AbstractArray{T,N}; kwds...) where {T,N} =
    AnnotatedArray(data; kwds...)
AnnotatedArray{T}(data::AbstractArray{T}, args...; kwds...) where {T} =
    AnnotatedArray(data, args...; kwds...)
AnnotatedArray{T}(data::AbstractArray{T}; kwds...) where {T} =
    AnnotatedArray(data; kwds...)

# Constructor for properties provided by keywords.
AnnotatedArray(data::AbstractArray; kwds...) =
    AnnotatedArray(data, Dict{Symbol,Any}(kwds...))

# Constructors with initial array and initial properties specified as key-value
# pairs.  When initial properties are specified as key-value pairs, we want to
# have a dictionary whose key type is specialized if possible (for efficiency)
# but avoid having value type specialized (for flexibility).  If one wants
# specific properties key and value types, it is always possible to explicitly
# specify a dictionary.
AnnotatedArray(data::AbstractArray, args::Pair{K,<:Any}...) where {K} =
    AnnotatedArray(data, Dict{K,Any}(args...))
AnnotatedArray(data::AbstractArray, args::Pair{<:AbstractString,<:Any}...) =
    AnnotatedArray(data, Dict{String,Any}(args...))
AnnotatedArray(data::AbstractArray, args::Pair...) =
    AnnotatedArray(data, _initialproperties(Dict(args...)))

_initialproperties(prop::Dict{K,Any}) where {K} = prop
_initialproperties(prop::Dict{K}) where {K} = convert(Dict{K,Any}, prop)

# Constructors that allocate the initial array.  If only keywords are
# specified, the dimensions of the array may be provided by integer arguments.
AnnotatedArray{T,N}(init, dims::Integer...; kwds...) where {T,N} =
    AnnotatedArray{T,N}(init, dims; kwds...)
AnnotatedArray{T}(init, dims::Integer...; kwds...) where {T} =
    AnnotatedArray{T}(init, dims; kwds...)
AnnotatedArray{T,N}(init, dims::NTuple{N,Integer}, args...; kwds...) where {T,N} =
    AnnotatedArray(Array{T,N}(init, dims), args...; kwds...)
AnnotatedArray{T}(init, dims::Tuple{Vararg{Integer}}, args...; kwds...) where {T} =
    AnnotatedArray(Array{T}(init, dims), args...; kwds...)

# The following ones are to avoid common errors.
AnnotatedArray(::UndefInitializer, dims::Integer...; kwds...) =
    throw_missing_type_parameter()
AnnotatedArray(::UndefInitializer, dims::Tuple{Vararg{Integer}}, args...; kwds...) =
    throw_missing_type_parameter()

"""

```julia
properties(A)
```

yields the properties associated with `A`, an instance of `AnnotatedArray`.

"""
@inline properties(A::AnnotatedArray) = Base.getfield(A, :prop)

"""

```julia
propertyname(T, sym)
```

converts symbol `sym` to a suitable key for an instance of an annotated array
of type `T` (a sub-type of `AbstractAnnotatedArray`), throwing an error if this
conversion is not supported.

"""
propertyname(::Type{<:DynamicallyAnnotatedArray{<:Any,<:Any,Symbol}}, sym::Symbol) = sym
propertyname(::Type{<:StaticallyAnnotatedArray}, sym::Symbol) = sym
propertyname(::Type{<:DynamicallyAnnotatedArray{<:Any,<:Any,String}}, sym::Symbol) =
    String(sym)
@noinline propertyname(::Type{T}, sym::Symbol) where {T<:AbstractAnnotatedArray} =
    error(string("converting symbolic key to ", keytype(T), " is not supported"))

# FIXME: The following specialized method extensions are omitted because
#        `obj.key` is much slower than `obj[key]`
#
#propertyname(::Type{<:DynamicallyAnnotatedArray{<:Any,<:Any,String}}, sym::Symbol) =
#    String(sym)
#propertyname(::Type{<:DynamicallyAnnotatedArray{<:Any,<:Any,K}}, sym::Symbol) where {K} =
#    convert(K, sym)

# Extend methods so that syntax `obj.field` can be used.

@inline Base.getproperty(A::StaticallyAnnotatedArray, sym::Symbol) =
    getproperty(properties(A), sym)

@inline Base.getproperty(A::T, sym::Symbol) where {T<:DynamicallyAnnotatedArray} =
    getindex(properties(A), propertyname(T, sym))

Base.setproperty!(A::StaticallyAnnotatedArray, sym::Symbol, val) =
    throw_immutable_properties()

@inline Base.setproperty!(A::T, sym::Symbol, val) where {T<:DynamicallyAnnotatedArray} =
    setindex!(properties(A), val, propertyname(T, sym))

Base.propertynames(A::StaticallyAnnotatedArray, private::Bool=false) =
    propertynames(properties(A), private)

# FIXME: Should return `Tuple(keys(properties(A)))` to conform to the doc. of
#        `propertynames` but this is slower and for most purposes, an iterable
#        is usually needed.
Base.propertynames(A::DynamicallyAnnotatedArray, private::Bool=false) where {T,N,P} =
    keys(properties(A))

#
# Notes:
#
# * keytype() and valtype() do not take NamedTuple as argument.  So to avoid
#   type-piracy, we define our own keytype() and keytype() methods and just
#   override Base.keytype() and Base.valtype() for AbstractAnnotatedArray.
#
# * The following methods to determine the result of valtype() for a named
#   tuple should cover most of the needs. Maybe we should use the same
#   algorithm as the one used by dictionary constructors to determine the value
#   type.
#
keytype(x) = Base.keytype(x) # use Base definition by default

keytype(::NamedTuple) = Symbol
keytype(::Type{<:NamedTuple}) = Symbol

keytype(::T) where {T<:AbstractAnnotatedArray} = keytype(T)
keytype(::Type{<:AbstractAnnotatedArray{<:Any,<:Any,P}}) where {P} = keytype(P)

Base.keytype(::T) where {T<:AbstractAnnotatedArray} = keytype(T)
Base.keytype(::Type{T}) where {T<:AbstractAnnotatedArray} = keytype(T)

valtype(x) = Base.valtype(x) # use Base definition by default

valtype(::T) where {T<:NamedTuple} = valtype(T)
valtype(::Type{<:NamedTuple{<:Any,T}}) where {T} = _eltype(T)
_eltype(::Type{<:Tuple{Vararg{T}}}) where {T} = T
_eltype(::Type{<:Tuple}) = Any

valtype(::T) where {T<:AbstractAnnotatedArray} = valtype(T)
valtype(::Type{<:AbstractAnnotatedArray{<:Any,<:Any,P}}) where {P} = valtype(P)

Base.valtype(::T) where {T<:AbstractAnnotatedArray} = valtype(T)
Base.valtype(::Type{T}) where {T<:AbstractAnnotatedArray} = valtype(T)

"""

```julia
nkeys(A)
```

yields the number of keys in the dictionary `A` or in the dictionary associated
with `A` if it is an instance of `AttributeArray`.

"""
nkeys(A::AbstractAnnotatedArray) = nkeys(properties(A))
nkeys(A::NamedTuple) = length(A)
nkeys(A::AbstractDict) = length(A)

# Extend methods to access an instance of `AttributeArray` like a dictionary.
Base.haskey(A::DynamicallyAnnotatedArray, key) = haskey(properties(A), key)
Base.haskey(A::StaticallyAnnotatedArray, key::Symbol) = haskey(properties(A), key)
Base.haskey(A::StaticallyAnnotatedArray, key) = false

Base.keys(A::AbstractAnnotatedArray) = keys(properties(A))
Base.values(A::AbstractAnnotatedArray) = values(properties(A))

Base.getkey(A::DynamicallyAnnotatedArray, key, def) = getkey(properties(A), key, def)
Base.getkey(A::StaticallyAnnotatedArray, key::Symbol, def) = (haskey(A, key) ? key : def)
Base.getkey(A::StaticallyAnnotatedArray, key, def) = def

Base.get(A::DynamicallyAnnotatedArray, key, def) = get(properties(A), key, def)
Base.get(A::StaticallyAnnotatedArray, key::Symbol, def) = get(properties(A), key, def)
Base.get(A::StaticallyAnnotatedArray, key, def) = def

Base.get!(A::DynamicallyAnnotatedArray, key, def) = get!(properties(A), key, def)
Base.get!(A::StaticallyAnnotatedArray, key, def) = throw_immutable_properties()

# FIXME: Check whether key type must be narrowed in getindex() and setindex!()
#        to avoid conflict with array indexing.
Base.getindex(A::DynamicallyAnnotatedArray{<:Any,<:Any,K}, key::K) where {K} =
    getindex(properties(A), key)
Base.getindex(A::DynamicallyAnnotatedArray{<:Any,<:Any,String}, key::AbstractString) =
    getindex(properties(A), key)
Base.getindex(A::StaticallyAnnotatedArray, key::Symbol) =
    getindex(properties(A), key)

@inline function Base.setindex!(A::DynamicallyAnnotatedArray{<:Any,<:Any,K},
                                val, key::K) where {K}
    properties(A)[key] = val
    A
end

Base.setindex!(A::StaticallyAnnotatedArray, val, key::Symbol) =
    throw_immutable_properties()

Base.delete!(A::DynamicallyAnnotatedArray, key) = begin
    delete!(properties(A), key)
    return A
end
Base.delete!(A::StaticallyAnnotatedArray, key) = throw_immutable_properties()

Base.pop!(A::DynamicallyAnnotatedArray, key) = pop!(properties(A), key)
Base.pop!(A::DynamicallyAnnotatedArray, key, def) = pop!(properties(A), key, def)
Base.pop!(A::StaticallyAnnotatedArray, key) = throw_immutable_properties()
Base.pop!(A::StaticallyAnnotatedArray, key, def) = throw_immutable_properties()

Base.pairs(A::AbstractAnnotatedArray) = pairs(properties(A))

Base.merge(A::AbstractAnnotatedArray, others...) =
    merge(properties(A), others...)
Base.merge(A::AbstractDict, B::AbstractAnnotatedArray, others...) =
    merge(A, properties(B), others...)
Base.merge(combine::Function, A::AbstractAnnotatedArray, others...) =
    merge(combine, properties(A), others...)
Base.merge(combine::Function, A::AbstractDict, B::AbstractAnnotatedArray, others...) =
    merge(combine, A, properties(B), others...)

Base.merge!(A::DynamicallyAnnotatedArray, others...) = begin
    merge!(properties(A), others...)
    return A
end
Base.merge!(A::StaticallyAnnotatedArray, others...) = throw_immutable_properties()
Base.merge!(A::AbstractDict, B::AbstractAnnotatedArray, others...) =
    merge!(A, properties(B), others...)

Base.merge!(combine::Function, A::DynamicallyAnnotatedArray, others...) = begin
    merge!(combine, properties(A), others...)
    return A
end
Base.merge!(combine::Function, A::StaticallyAnnotatedArray, others...) =
    throw_immutable_properties()
Base.merge!(combine::Function, A::AbstractDict, B::AbstractAnnotatedArray, others...) =
    merge!(combine, A, properties(B), others...)

# Errors.
throw_missing_type_parameter() = error("type parameter T must be specified")
throw_immutable_properties() = error("properties are immutable")

end # module


text:
export  convert2affinePCE,
        calculateAffinePCE,
        assign2multi,
        evaluatePCE,
        sampleMeasure,
        samplePCE,
        mean,
        var,
        std

# auxiliary functions
function _createMethodVector(n::Int, word::String="adaptiverejection")
    [word for i in 1:n]
end

function _createMethodVector(m::ProductMeasure, word::String="adaptiverejection")
    _createMethodVector(length(m.measures), word)
end

function _createMethodVector(mop::MultiOrthoPoly, word::String="adaptiverejection")
    _createMethodVector(mop.measure, word)
end

function _checkStandardDevation(σ::Real)
    σ < 0 && throw(DomainError(σ, "σ has to be non-negative"))
    σ < 1e-4 && @warn "σ is close to zero (σ = $σ)"
end

function _checkKind(kind::String)
    lowercase(kind) ∉ ["lbub", "μσ"] && throw(DomainError(kind, "this kind is not supported"))
    lowercase(kind)
end

function _checkBounds(lb::Real, ub::Real)
    lb >= ub && throw(DomainError((lb, ub), "inconsistent bounds"))
end

function _checkNumberOfSamples(n::Int)
    n < 1 && throw(DomainError(n, "invalid number of samples"))
end
"""
    calculateAffinePCE(α::AbstractVector{<:Real})
Computes the affine PCE coefficients ``x_0`` and ``x_1`` from recurrence coefficients ``\alpha``.
"""
function calculateAffinePCE(α::AbstractVector{<:Real})
    length(α) < 1 && throw(DomainError(length(α), "not enough recursion coefficients"))
    return [α[1], 1.]
end

calculateAffinePCE(op::AbstractOrthoPoly) = calculateAffinePCE(op.α)

function calculateAffinePCE(i::Int, mop::MultiOrthoPoly)
    p = length(mop.name)
    i > p && throw(InconsistencyError("mop is $p-variate, but the PC for the $i-variate entry was requested"))
    calculateAffinePCE(mop.uni[i])
end

"""
```
convert2affinePCE(mu::Real, sigma::Real, op::AbstractCanonicalOrthoPoly; kind::String)
```
Computes the affine PCE coefficients ``x_0`` and ``x_1`` from
```math
X = a_1 + a_2 \\Xi = x_0 + x_1 \\phi_1(\\Xi),
```
where ``\\phi_1(t) = t-\\alpha_0`` is the first-order monic basis polynomial.

Works for subtypes of AbstractCanonicalOrthoPoly. The keyword `kind in ["lbub", "μσ"]`
specifies whether `p1` and `p2` have the meaning of lower/upper bounds or mean/standard deviation.
"""
function convert2affinePCE(a1::Real,a2::Real,α0::Real)
    [ a1 + α0 * a2; a2 ]
end

function convert2affinePCE(mu::Real, sigma::Real, op::GaussOrthoPoly)
    _checkStandardDevation(sigma)
    convert2affinePCE(mu, sigma, first(op.α))
end

function convert2affinePCE(par1::Real, par2::Real, op::Uniform01OrthoPoly; kind::String="lbub")
    kind = _checkKind(kind)
    a1, a2 =
        if kind == "lbub"
            _checkBounds(par1, par2)
            par1, par2 - par1
        elseif kind == "μσ"
            _checkStandardDevation(par2)
            par1 - sqrt(3)*par2 , 2*sqrt(3)*par2
        end
    convert2affinePCE(a1, a2, first(op.α))
end

function convert2affinePCE(par1::Real, par2::Real, op::Uniform_11OrthoPoly; kind::String="lbub")
    kind = _checkKind(kind)
    a1, a2 =
        if kind == "lbub"
            _checkBounds(par1, par2)
            0.5 * (par1 + par2), 0.5 * (par2 - par1)
        elseif kind == "μσ"
            _checkStandardDevation(par2)
            par1, sqrt(3) * par2
        end
    convert2affinePCE(a1, a2, first(op.α))
end

function convert2affinePCE(p1::Real, p2::Real, op::Beta01OrthoPoly; kind::String="lbub")
    kind = _checkKind(kind)
    α, β = op.measure.ashapeParameter, op.measure.bshapeParameter
    a1, a2 =
        if kind == "lbub"
            _checkBounds(p1, p2)
            a1, a2 = p1, p2-p1
        elseif kind == "μσ"
            _checkStandardDevation(p2)
            a1, a2 = p1 - sqrt(α / β) * sqrt(1 + α + β) * p2, (α + β) * sqrt((α + β + 1) / (α * β)) * p2
        end
    convert2affinePCE(a1, a2, first(op.α))
end

function convert2affinePCE(p1::Real, p2::Real, op::GammaOrthoPoly)
    throw(error("convert2affine not yet implemented for $(typeof(op))"))
end

function convert2affinePCE(p1::Real, p2::Real, op::LogisticOrthoPoly)
    _checkStandardDevation(p2)
    convert2affinePCE(p1, p2, first(op.α))
end

function assign2multi(x::AbstractVector{<:Real},i::Int,ind::AbstractMatrix{<:Int})
    l, p = size(ind)
    nx, deg = length(x), ind[end,end]
    nx > deg + 1 && throw(DomainError(nx, "inconsistent number of coefficients ($nx vs $(deg+1))"))
    i > p && throw(DomainError((i,p), "basis is $p-variate, you requested $i-variate"))
    myind = findUnivariateIndices(i,ind)[1:nx]
    y = spzeros(Float64,l)
    y[myind] = x
    return y
end

"""
__Univariate__
```
sampleMeasure(n::Int,name::String,w::Function,dom::Tuple{<:Real,<:Real},symm::Bool,d::Dict;method::String="adaptiverejection")
sampleMeasure(n::Int,m::Measure;method::String="adaptiverejection")
sampleMeasure(n::Int,op::AbstractOrthoPoly;method::String="adaptiverejection")
```
Draw `n` samples from the measure `m` described by its
- `name`
- weight function `w`,
- domain `dom`,
- symmetry property `symm`,
- and, if applicable, parameters stored in the dictionary `d`.
By default an adaptive rejection sampling method is used (from [AdaptiveRejectionSampling.jl](https://github.com/mauriciogtec/AdaptiveRejectionSampling.jl)),
unless it is a common random variable for which [Distributions.jl](https://github.com/JuliaStats/Distributions.jl) is used.

The function is dispatched to accept `AbstractOrthoPoly`.

__Multivariate__
```
sampleMeasure(n::Int,m::ProductMeasure;method::Vector{String}=["adaptiverejection" for i=1:length(m.name)])
sampleMeasure(n::Int,mop::MultiOrthoPoly;method::Vector{String}=["adaptiverejection" for i=1:length(mop.meas.name)])
```
Multivariate extension which provides array of samples with `n` rows and
as many columns as the multimeasure has univariate measures.
"""
function sampleMeasure(n::Int, w::Function, dom::Tuple{<:Real,<:Real}; method::String="adaptiverejection")
    _checkNumberOfSamples(n)
    method = lowercase(method)

    if method == "adaptiverejection"
        # works only if w is log-concave
        sampler = RejectionSampler(w, dom)
        return run_sampler!(sampler, n)
    elseif method == "rejection"
        # all purpose method but needs a solid envelope PDF
        throw(error("method $method not yet implemented"))
    elseif method == "inversecdf"
        throw(error("method $method not yet implemented"))
        # requires CDF and reliable root-finding method
    else
        throw(error("method $method not implemented"))
    end
end

function sampleMeasure(n::Int, meas::AbstractMeasure; method::String="adaptiverejection")
    sampleMeasure(n, meas.w, meas.dom; method=method)
end

sampleMeasure(n::Int, op::AbstractOrthoPoly; method::String="adaptiverejection") = sampleMeasure(n,op.measure; method=method)
sampleMeasure(n::Int, op::AbstractCanonicalOrthoPoly) = sampleMeasure(n, op.measure::AbstractCanonicalMeasure)

function sampleMeasure(n::Int, dist::Distribution{Univariate,Continuous})
    _checkNumberOfSamples(n)
    rand(dist, n)
end

sampleMeasure(n::Int, meas::GaussMeasure) = sampleMeasure(n, Normal())
sampleMeasure(n::Int, meas::Uniform01Measure) = sampleMeasure(n, Uniform())
sampleMeasure(n::Int, meas::Beta01Measure) = sampleMeasure(n, Beta(meas.ashapeParameter, meas.bshapeParameter))
sampleMeasure(n::Int, meas::GammaMeasure) = sampleMeasure(n, Gamma(meas.shapeParameter, 1/meas.rateParameter))
sampleMeasure(n::Int, meas::LogisticMeasure) = sampleMeasure(n, Logistic())

function sampleMeasure(n::Int, meas::AbstractCanonicalMeasure; method::String="adaptiverejection")
    @warn "ignoring keyword method; sampling from Distributions.jl instead"
    sampleMeasure(n, meas)
end

function sampleMeasure(n::Int, measure::ProductMeasure; method::Vector{String}=_createMethodVector(measure))
    samples = Matrix{Float64}(undef, n, 0)
    for (k, unimeasure) in enumerate(measure.measures)
        samples = hcat(samples, sampleMeasure(n, unimeasure; method = method[k]))
    end
    samples
end

sampleMeasure(n::Int, mop::MultiOrthoPoly; method::Vector{String}=_createMethodVector(mop)) = sampleMeasure(n,mop.measure; method=method)

"""
    evaluatePCE(x::AbstractVector{<:Real},ξ::AbstractVector{<:Real},α::AbstractVector{<:Real},β::AbstractVector{<:Real})
Evaluation of polynomial chaos expansion
```math
\\mathsf{x} = \\sum_{i=0}^{L} x_i \\phi_i{\\xi_j},
```
where `L+1 = length(x)` and ``x_j`` is the ``j``th sample where ``j=1,\\dots,m``
with `m = length(ξ)`.
"""
function evaluatePCE(x::AbstractVector{<:Real},ξ::AbstractVector{<:Real},α::AbstractVector{<:Real},β::AbstractVector{<:Real})
    length(α) != length(β) && throw(InconsistencyError("inconsistent number of recurrence coefficients"))
    Nsmpl = length(ξ)
    _checkNumberOfSamples(Nsmpl)
    Nx, Nrec = length(x), length(α)
    Nx > Nrec && throw(InconsistencyError("not enough recursion coefficients"))
    if Nrec > Nx α,β = α[1:Nx], β[1:Nx] end
    ϕ = zeros(Float64, Nsmpl, Nx)
    for n in 1:Nx
        ϕ[:,n] = evaluate(n-1,ξ,α,β)
    end
    ϕ*x
end
evaluatePCE(x::AbstractVector{<:Real},ξ::Real,α::AbstractVector{<:Real},β::AbstractVector{<:Real}) = evaluatePCE(x,[ξ],α,β)
evaluatePCE(x::AbstractVector{<:Real},ξ::AbstractVector{<:Real},op::AbstractOrthoPoly) = evaluatePCE(x,ξ,op.α,op.β)

function evaluatePCE(x::AbstractVector{<:Real},ξ::Matrix{<:Real},α::AbstractVector{<:AbstractVector{<:Real}},β::AbstractVector{<:AbstractVector{<:Real}},ind::AbstractMatrix{Int})
    Nsmpl = size(ξ,1)
    _checkNumberOfSamples(Nsmpl)
    !(length(α) == length(β) == size(ξ,2) == size(ind,2)) && throw(InconsistencyError("inconsistent number of coefficients"))
    Nx = length(x)
    Nx > size(ind,1) && throw(InconsistencyError("too few pc coefficients (resp: too small basis)"))
    ϕ = zeros(Float64,Nsmpl,Nx)
    for n in 1:Nx
        ϕ[:,n] = evaluate(ind[n,:],ξ,α,β)
    end
    ϕ*x
end
function evaluatePCE(x::AbstractVector{<:Real},ξ::Matrix{<:Real},mOP::MultiOrthoPoly)
    a,b = coeffs(mOP)
    evaluatePCE(x,ξ,a,b,mOP.ind)
end

"""
__Univariate__
```
samplePCE(n::Int,x::AbstractVector{<:Real},op::AbstractOrthoPoly;method::String="adaptiverejection")
```
Combines [`sampleMeasure`](@ref) and [`evaluatePCE`](@ref), i.e. it first draws `n` samples
from the measure, then evaluates the PCE for those samples.

__Multivariate__
```
samplePCE(n::Int,x::AbstractVector{<:Real},mop::MultiOrthoPoly;method::Vector{String}=["adaptiverejection" for i=1:length(mop.meas.name)])
```
"""
function samplePCE(n::Int,x::AbstractVector{<:Real},op::AbstractOrthoPoly;method::String="adaptiverejection")
    ξ = sampleMeasure(n,op;method=method)
    evaluatePCE(x,ξ,op)
end

function samplePCE(n::Int,x::AbstractVector{<:Real},op::AbstractCanonicalOrthoPoly;method::String="adaptiverejection")
    ξ = sampleMeasure(n,op)
    evaluatePCE(x,ξ,op)
end

function samplePCE(n::Int,x::AbstractVector{<:Real},mop::MultiOrthoPoly;method::Vector{String} = _createMethodVector(mop))
    ξ = sampleMeasure(n, mop; method=method)
    evaluatePCE(x, ξ, mop)
end

"""
__Univariate__
```
mean(x::AbstractVector,op::AbstractOrthoPoly)
```
__Multivariate__
```
mean(x::AbstractVector,mop::MultiOrthoPoly)
```
compute mean of random variable with PCE `x`
"""
mean(x::AbstractVector,op::AbstractOrthoPoly) = x[1]*computeSP2(0,op.β)

function mean(x::AbstractVector,mop::MultiOrthoPoly)
    nunc = length(mop.uni)
    x[1]*computeSP(zeros(Int64,nunc),mop)
end

"""
__Univariate__
```
var(x::AbstractVector,op::AbstractOrthoPoly)
var(x::AbstractVector,t2::Tensor)
```
__Multivariate__
```
var(x::AbstractVector,mop::MultiOrthoPoly)
var(x::AbstractVector,t2::Tensor)
```
compute variance of random variable with PCE `x`
"""
function var(x::AbstractVector,op::AbstractOrthoPoly)
    t = computeSP2(op)
    # length(t2) > length(x) && throw(InconsistencyError("cannot compute variance; too many PCE coefficients"))
    sum( x[i]^2*t[i] for i in 2:length(x) )
end

function var(x::AbstractVector,t2::AbstractTensor)
    sum( x[i]^2*t2.get([i-1,i-1]) for i in 2:length(x) )
end

function var(x::AbstractVector,mop::MultiOrthoPoly)
    length(x) > size(mop.ind,1) && throw(InconsistencyError("cannot compute variance; too many PCE coefficients"))
    var(x,Tensor(2,mop))
end

"""
__Univariate__
```
std(x::AbstractVector,op::AbstractOrthoPoly)
```
__Multivariate__
```
std(x::AbstractVector,mop::MultiOrthoPoly)
```
compute standard deviation of random variable with PCE `x`
"""
std(x::AbstractVector,op::AbstractOrthoPoly) = sqrt(var(x,op))
std(x::AbstractVector,mop::MultiOrthoPoly) = sqrt(var(x,mop))