text:
## Boolean queries ##

# check the presence of vertex/vertices and hyperedge/hyperedges; unique is used to only check parallel hyperedges or repeated vertices of loops once
has_vertex(x::T, v) where {T<:AbstractHyperGraph} = v in vertices(x)
has_vertex(x::T, vs::AbstractVector) where {T<:AbstractHyperGraph} = [has_vertex(x, v) for v in unique(vs)]
has_vertices(x::T, vs::AbstractVector) where {T<:AbstractHyperGraph} = all(has_vertex(x, vs))
has_vertex(e::T, v) where {T<:AbstractHyperEdge} = v in vertices(e)
has_vertex(e::T, vs::AbstractVector) where {T<:AbstractHyperEdge} = [has_vertex(e, v) for v in unique(vs)]
has_vertices(e::T, vs::AbstractVector) where {T<:AbstractHyperEdge} = all(has_vertex(e, vs))
has_hyperedge(x::T, e::U) where {T<:AbstractHyperGraph, U<:AbstractHyperEdge} = e in hyperedges(x)
has_hyperedge(x::T, es::AbstractVector{U}) where {T<:AbstractHyperGraph, U<:AbstractHyperEdge} = [has_hyperedge(x, e) for e in unique(es)]
Base.in(e::T, x::U) where {T<:AbstractHyperEdge, U<:AbstractHyperGraph} = has_hyperedge(x, e)
has_hyperedges(x::T, es::AbstractVector{U}) where {T<:AbstractHyperGraph, U<:AbstractHyperEdge} = all(has_hyperedge(x, es))
in_src(v, e::T) where {T<:AbstractHyperEdge} = v in src(e)
in_src(vs::AbstractVector, e::T) where {T<:AbstractHyperEdge} = [in_src(v, e) for v in unique(vs)]
in_tgt(v, e::T) where {T<:AbstractHyperEdge} = v in tgt(e)
in_tgt(vs::AbstractVector, e::T) where {T<:AbstractHyperEdge} = [in_tgt(v, e) for v in unique(vs)]

# check incidence -- hyperedge e is incident on vertex v; either order is correct
isincident(e::AbstractHyperEdge, v) = has_vertex(e, v)
isincident(v, e::AbstractHyperEdge) = isincident(e, v)

# Base.isempty
Base.isempty(e::T) where {T<:AbstractHyperEdge} = length(vertices(e)) == 0
"""
    Base.isempty

A hyergraph is empty if it has no hyperedges.
"""
Base.isempty(x::T) where {T<:AbstractHyperGraph} = nhe(x) == 0

# a simple hypergraph has no loops or multiple hyperedges
issimple(x::T) where {T<:AbstractHyperGraph} = !has_loops(x) && !has_multi_hyperedges(x)

# check if a hypergraph may be represented as a graph, that is if the cardinality of all hyperedges is 2
isgraph(x::T) where {T<:AbstractHyperGraph} = cardinalities(x) == fill(2, nhe(x))

# check that a hyperedge is a subhyperedge of another hyperedge
@traitfn issubhyperedge(e_sub::T, e_super::T) where {T<:AbstractHyperEdge; !IsOriented{T}} = all([isincident(e_super, v) for v in vertices(e_sub)]) # unoriented
@traitfn issubhyperedge(e_sub::T, e_super::T) where {T<:AbstractHyperEdge; IsOriented{T}} = all(in_src(src(e_sub), e_super)) && all(in_tgt(tgt(e_sub), e_super)) # oriented

# check that a hyperedge is a subhyperedge of at most one hyperedge of a hypergraph (any hyperedge is a subhyperedge of itself)
isuniquesubhyperedge(e_sub::T, x::U) where {T<:AbstractHyperEdge, U<:AbstractHyperGraph} = sum([issubhyperedge(e_sub, e_super) for e_super in hyperedges(x)]) == 1

# loops
"""
    isloop

A hyperedge is a loop if it has repeated vertices.
"""
isloop(e::T) where {T<:AbstractHyperEdge} = length(vertices(e)) != length(unique(vertices(e)))
"""
    is_positive_loop

An oriented hyperedge is a positive loop wrt. a vertex if that vertex is both in the source set and in the target set.
"""
function is_positive_loop end
@traitfn is_positive_loop(e::AbstractHyperEdge::IsOriented, v) = has_vertex(e, v) ? v in intersect(src(e), tgt(e)) : error("vertex $v not in hyperedge")
@traitfn is_positive_loop(e::AbstractHyperEdge::IsOriented) = any([v in intersect(src(e), tgt(e)) for v in vertices(e)])
@traitfn is_positive_loop(e::AbstractHyperEdge::IsOriented, v, k) = cardinality(e) == k && is_positive_loop(e, v)
has_loops(x::AbstractHyperGraph) = any(isloop.(hyperedges(x)))

# uniformity, regularity
"""
    iskuniform

A hypergraph is k-uniform if all hyperedges have cardinality k. Equivalently, a hypergraph is k-uniform if its rank is equal to its co-rank [Berge1989, Bretto2013]; checking this is however slower than the method currently implemented.
"""
iskuniform(x::T, k::Int) where {T<:AbstractHyperGraph} = sum(cardinalities(hyperedges(x)) .== k) == nhe(x)
"""
    iskregular

A hypergraph is k-regular all its vertices have degree k.
"""
iskregular(x::T, k::Int) where {T<:AbstractHyperGraph} = sum(degrees(x) .== k) == nv(x)

# chemical hypergraph & hyperedge queries
function is_netstoich_null(e::ChemicalHyperEdge, v)
    s = stoichiometries(e, v)
    (length(first(s)) == length(last(s))) && all(first(s) .- last(s) .== 0)
end
iscatalyst(x::ChemicalHyperGraph, v) = v in catalysts(x)
iscatalyst(e::ChemicalHyperEdge, v) = v in catalysts(e)

# parallel and multi-hyperedges
has_parallel_hyperedges(x::T) where {T<:AbstractHyperGraph} = num_parallel_hyperedges(x) > 0
has_multi_hyperedges(x::T) where {T<:AbstractHyperGraph} = num_multi_hyperedges(x) > 0

## content queries ##

# incident hyperedges
# returns all hyperedges incident on a given vertex
function incident_hyperedges(x::T, v; check_presence=true) where {T<:AbstractHyperGraph}
    check_presence && !has_vertex(x, v) && error("vertex $v not found in hypergraph vertices")
    filter(e -> isincident(e, v), hyperedges(x))
end
incident_hyperedges(x::T, vs::AbstractVector) where {T<:AbstractHyperGraph} = union(vcat([incident_hyperedges(x, v) for v in vs])...)

# empty hyperedges, hypergraphs
empty_hyperedges(x::T) where {T<:AbstractHyperGraph} = filter(e -> isempty(e), hyperedges(x))
num_empty_hyperedges(x::T) where {T<:AbstractHyperGraph} = length(empty_hyperedges(x))
has_empty_hyperedges(x::T) where {T<:AbstractHyperGraph} = num_empty_hyperedges(x) > 0

# neighbors
function neighbors(x::T, v) where {T<:AbstractHyperGraph} # unoriented
    union(vcat(symdiff.(vertices(incident_hyperedges(x, v)), v))...)
end
@traitfn function inneighbors(x::T::IsOriented, v) where {T<:AbstractHyperGraph}

    # incident hyperedges where v is in the target set
    v_in_tgt = filter(e -> in_tgt(v, e), incident_hyperedges(x, v))

    union(vcat(src(v_in_tgt)...))
end
@traitfn function outneighbors(x::T::IsOriented, v) where {T<:AbstractHyperGraph}

    # incident hyperedges where v is in the source set
    v_in_src = filter(e -> in_src(v, e), incident_hyperedges(x, v))

    union(vcat(tgt(v_in_src)...))
end
@traitfn all_neighbors(x::T::(!IsOriented), v) where {T<:AbstractHyperGraph} = neighbors(x, v) # unoriented
@traitfn all_neighbors(x::T::IsOriented, v) where {T<:AbstractHyperGraph} = union(inneighbors(x, v), outneighbors(x, v)) # oriented

# loops
loops(es::AbstractVector{T}) where {T<:AbstractHyperEdge} = filter(e -> isloop(e), es)
loops(x::AbstractHyperGraph) = loops(hyperedges(x))
@traitfn positive_loops(es::AbstractVector{T}) where {T<:AbstractHyperEdge; IsOriented{T}} = filter(e -> is_positive_loop(e), es)
@traitfn positive_loops(x::AbstractHyperGraph::IsOriented) = positive_loops(hyperedges(x))
num_loops(x::AbstractHyperGraph) = length(loops(x))

# catalysts (chemical hypergraphs)
catalysts(e::ChemicalHyperEdge) = filter(v -> is_positive_loop(e, v) && is_netstoich_null(e, v), unique(vertices(e)))
catalysts(es::AbstractVector{ChemicalHyperEdge{T}}) where {T} = vcat(catalysts.(es)...)
catalysts(x::ChemicalHyperGraph) = vcat(catalysts(hyperedges(x))...)

# parallel and multi-hyperedges
"""
    parallel_hyperedges

There are parallel hyperedges in a hypergraph if any of its hyperedges is present more than once in the set of hyperedges."""
parallel_hyperedges(x::T) where {T<:AbstractHyperGraph} = filter(e -> sum(e == _e for _e in hyperedges(x)) > 1, hyperedges(x))
num_parallel_hyperedges(x::T) where {T<:AbstractHyperGraph} = length(parallel_hyperedges(x))
"""
    multi_hyperedges

There are multi-hyperedges in a hypergraph if any of its hyperedges is not a unique subedge i.e. there is at least another hyperedge of the same or greater length."""
multi_hyperedges(x::T) where {T<:AbstractHyperGraph} = filter(e -> !isuniquesubhyperedge(e, x), hyperedges(x))
num_multi_hyperedges(x::T) where {T<:AbstractHyperGraph} = length(multi_hyperedges(x))


text:
using Documenter, ReservoirComputing

makedocs(
    modules=[ReservoirComputing],
    clean=true,doctest=false,
    sitename = "ReservoirComputing.jl",
    pages = [
        "ReservoirComputing.jl" => "index.md",
        "Tutorials" => Any[
            "Basics"=>"examples/esn.md",
            "Using different layers"=>"examples/layers.md",
            "Using different linear methods"=>"examples/linear.md",
            "Double Activation Function ESN"=>"examples/dafesn.md",
             #"SVESM"=>"examples/svesm.md",
            "ESGP"=>"examples/esgp.md"

            ],
        "User Guide"=>Any[
            "Models"=>"user/models.md",
            "Layers"=>"user/layers.md",
            "Linear Models"=>"user/linear.md",
            "Special ESNs"=>"user/spec.md",
            "Non Linear Algorithm"=>"user/nla.md"
            ]
        ]
    )

deploydocs(
   repo = "github.com/SciML/ReservoirComputing.jl.git";
   push_preview = true
)

ext: jl


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

"""
    vandermonde_matrix(polynomial, coordinates)

Given some polynomial and coordinates points (1-3 dimensions), create a Vandermonde
matrix.

# Example

To genererate a Vandermonde matrix for a reference quadrangle `[-1.0, 1.0]^2` for
polynomial `p(u,v) = 1 + u + v + u*v`, one writes:

```julia
polynomial = :(1 + u + v + u*v)
coordinates = ((-1.0,-1.0), (1.0,-1.0), (1.0,1.0), (-1.0,1.0))
V = vandermonde_matrix(polynomial, coordinates)

# output

[
1.0 -1.0 -1.0  1.0
1.0  1.0 -1.0 -1.0
1.0  1.0  1.0  1.0
1.0 -1.0  1.0 -1.0
]

```

# References
- Wikipedia contributors. (2018, August 1). Vandermonde matrix. In Wikipedia, The Free Encyclopedia. Retrieved 10:00, August 20, 2018, from https://en.wikipedia.org/w/index.php?title=Vandermonde_matrix&oldid=852930962
"""
function vandermonde_matrix(polynomial::Expr, coordinates::NTuple{N, NTuple{D, T}}) where {N, D, T<:Number}
    A = zeros(N, N)
    first(polynomial.args) == :+ || error("Use only summation between terms of polynomial")
    args = polynomial.args[2:end]
    for i in 1:N
        X = coordinates[i]
        if D == 1
            data = (:u => X[1],)
        elseif D == 2
            data = (:u => X[1], :v => X[2])
        elseif D == 3
            data = (:u => X[1], :v => X[2], :w => X[3])
        end
        for (j, term) in enumerate(args)
            A[i,j] = subs(term, data)
        end
    end
    return A
end


text:
# load module from local directory

using SurfaceHopping
using Base.Test

# -----------------------------------------------------------------------------
H_d = H_diabatic(2.6,2.5,0.1)
U, U_nk = diagonalise(H_d,true)

ε = 1e-2
@test 1.0 ≈ U[1] atol = ε

# -----------------------------------------------------------------------------
H_ad = H_adiabatic(H_d)
H = *(U_nk',*(H_d,U_nk))

@test H[1] ≈ H_ad[1] atol = ε

# -----------------------------------------------------------------------------
phi_l,phi_r = diabatic_states(0.7,0.1)
phi_l2(R) = phi_l(R)^2; phi_r2(R) = phi_r(R)^2;

ε = 1e-4
integral = riemann(phi_l2,-10,10,1000)
@test 1.0 ≈ integral atol = ε

integral = riemann(phi_r2,-10,10,1000)
@test 1.0 ≈ integral atol = ε

# -----------------------------------------------------------------------------
psi_g, psi_e = adiabatic_states(phi_l, phi_r, U_nk)
psi_g2(R) = psi_g(R)^2; psi_e2(R) = psi_e(R)^2

integral = riemann(psi_g2,-10,10,1000)
@test 1.0 ≈ integral atol = ε

integral = riemann(psi_e2,-10,10,1000)
@test 1.0 ≈ integral atol = ε

# -----------------------------------------------------------------------------
dphi = diabatic_derivative()

d_ge, d_eg, d_gg, d_ee = NACV(phi_l, phi_r, dphi, U_nk, 2.6, -0.1, 0.1, 0.1)
#println(d_ge, "  ", d_eg, "  ", d_gg, "  ", d_ee)

# -----------------------------------------------------------------------------
cl = 1.0; cr = 0.0
ag, ae = a_mn(U_nk, cl, cr)
coeffs_diabatic = abs(cl)^2 + abs(cr)^2
@test 1.0 ≈ coeffs_diabatic atol = ε

coeffs_adiabatic = abs(ag)^2 + abs(ae)^2
@test 1.0 ≈ coeffs_adiabatic atol = ε

# -----------------------------------------------------------------------------
dpl, dpr = diabatic_potential(1.0,2.6,0.0)
apg, ape = adiabatic_potential(dpl, dpr, 2.0)
agn, aen = new_a_mn(ag, ae, apg, ape, d_ge, d_eg, d_gg, d_ee, 2.6, 0.1)

coeffs_adiabatic = abs(agn)^2 + abs(aen)^2
@test 1.0 ≈ coeffs_adiabatic atol = ε

# -----------------------------------------------------------------------------
println("Tests finished succesfully.")


text:
#This file provides a pseudocode for algorithm one that provides the architecture for the complete algorithm. Some of these objects can change.

include("PseudoDemand.jl")
include("PseudoClearingError.jl")
include("PseudoNeighbor.jl")
include("loops.jl")

function coursematch(M, k, βmax, besterror, t, Np)
    ptild = Array{Int}(M) #new search price
    pstar = Array[] #price that gives best error
    searcherror = Array[] #best error in a searchstart
    DoubleN = Array[] #neighbors of the search
    foundnextstep = Array{Bool}(1) #boolean. true if a new search price has been found (one that does not generate the same demand as one that is in the tabu list)
    dem = zeros(M) #the demand for courses at search price
    currenterror = Array[] #error in the current search.
    τ = Array[]
    starttime = Dates.Time(now())
    while (Dates.Time(now()) - starttime).value < t *1000000000
    #repeat from l.2 to 35. This is a do until runtime > t.
        p = Array[[0,0,0,0,0]] #Initial guess for p. Will be replaced also
        searcherror =  α.(d(p)) ##= 50 function from issue #1. For this person, will need to make sure it returns the absolute value
        c = 0
        # while loop from l.7 to 34
        while c < 5
            DoubleN = copy(Np)  # going to be neighb_fun
            #DoubleN = N(p) #Don't forget this needs to be sorted by clearing error (Issue #6)
            foundnextstep[1] = false
            # repeat from l. 10 to 16
            findnextprices!(ptild,  DoubleN, dem, foundnextstep, τ)
            # if from line 17 to 18
            if isempty(DoubleN) == true
                c = 5
            end
            #l. 19 to 33
            if foundnextstep[1] == true
                push!(τ,copy(dem))
                p = copy(ptild)
                currenterror =  α(dem)
                c = resetcounter(currenterror, searcherror, c)
                searcherror = replacesearcherror(currenterror, searcherror)
                pstar = replacepstar(currenterror, besterror, p, pstar)
                besterror = replacebesterror(currenterror, besterror)
            end
            println("we are at step $c")
        end
    end
    return pstar
end

M = 5 #Number of courses offered
k = 3 #Number of courses students take
βmax = 20 #Maximum budget
t = 1 #time in seconds
besterror = 100 #for now
Np = Array[[0,0,0,0,0],[9,40,5,2,0], [1,1,1,1,1], [1,2,3,0,0], [9,40,5,2,0], [32,40,100,1,2], [10,20,30,3,0]]
pstar = coursematch(M, k, βmax, besterror, t, Np)

max_stars_repo_name: ScPo-CompEcon/CourseMatch

max_issues_repo_name: ScPo-CompEcon/CourseMatch

max_forks_repo_name: ScPo-CompEcon/CourseMatch


text:
COMPRESSIONS = ["bz2", "bzip2", "gz", "gzip", "none", "unknown"]
KNOWN_COMPRESSIONS = ["bz2", "bzip2", "gz", "gzip", "none"]
ACTIONS = ["compress", "decompress", "unknown"]
KNOWN_ACTIONS = ["compress", "decompress"]

using DispatcherCache: get_compressor

@testset "Compression" begin
    mktempdir() do dir
        value = "some value"
        for (compression, action) in Iterators.product(COMPRESSIONS, ACTIONS)
            file = joinpath(abspath(dir), join("tmp", ".", compression))

            compression in ["bz2", "bzip2"] && action=="compress" && begin
                compressor = get_compressor(compression, action)
                @test compressor == Bzip2CompressorStream
                @test open(compressor, file, "w") do fid
                    write(fid, value)
                end > 0
            end

            compression in ["bz2", "bzip2"] && action=="decompress" && begin
                compressor = get_compressor(compression, action)
                @test compressor == Bzip2DecompressorStream
                @test open(compressor, file, "r") do fid
                    read(fid, typeof(value))
                end == value
            end

            compression in ["gz", "gzip"] && action=="compress" && begin
                compressor = get_compressor(compression, action)
                @test compressor == GzipCompressorStream
                @test open(compressor, file, "w") do fid
                    write(fid, value)
                end > 0
            end
            compression in ["gz", "gzip"] && action=="decompress" && begin
                compressor = get_compressor(compression, action)
                @test compressor == GzipDecompressorStream
                @test open(compressor, file, "r") do fid
                    read(fid, typeof(value))
                end == value
            end

            compression == "none"  && action in ["compress", "decompress"] &&
                @test get_compressor(compression, action) == NoopStream
                # No need to test reading/writing

            !(compression in KNOWN_COMPRESSIONS) || !(action in KNOWN_ACTIONS) &&
                @test_throws ErrorException get_compressor(compression, action)
        end
    end
end


text:
using Test, Random
using Quadmath

@testset "fp decomp" begin
    y = Float128(2.0)
    x,n = frexp(y)
    @test x == Float128(0.5)
    @test n == 2
    z = ldexp(Float128(0.5), 2)
    @test z == y
end

@testset "conversions" begin
@testset "conversion $T" for T in (Float16, Float32, Float64, Int32, Int64, Int128, UInt32, UInt64, UInt128, BigFloat, BigInt)
    @test Float128(T(1)) + Float128(T(2)) == Float128(T(3))
    @test Float128(T(1)) + Float128(T(2)) <= Float128(T(3))
    @test Float128(T(1)) + Float128(T(2)) != Float128(T(4))
    @test Float128(T(1)) + Float128(T(2)) < Float128(T(4))
    if isbitstype(T)
        @test T(Float128(T(1)) + Float128(T(2))) === T(3)
    else
        @test T(Float128(T(1)) + Float128(T(2))) == T(3)
    end
    if isbitstype(T) && T <: Integer
        nb = 8*sizeof(T) - 5
        x = T(9) << nb
        xf = Float128(9) * 2.0^nb
        @test Float128(x) == xf
        @test T(xf) == x
    end
end

@testset "conversion $T exceptions" for T in (Int32, Int64, UInt32, UInt64)
    x = Float128(typemax(T))
    @test_throws InexactError T(x+Float128(1))
    x = Float128(typemin(T))
    @test_throws InexactError T(x-Float128(1))
end

@testset "conversion $T exceptions" for T in (Float32, Float64)
    x = Float128(typemax(T))
    @test isinf(T(x+Float128(1)))
    x = Float128(typemin(T))
    @test isinf(T(x-Float128(1)))
end

@testset "BigFloat" begin
    x = parse(Float128, "0.1")
    y = parse(Float128, "0.2")
    @test Float64(x+y) == Float64(BigFloat(x) + BigFloat(y))
    @test x+y == Float128(BigFloat(x) + BigFloat(y))
end

@testset "BigInt" begin
    x = parse(Float128, "100.0")
    y = parse(Float128, "25.0")
    @test Float64(x+y) == Float64(BigInt(x) + BigInt(y))
    @test x+y == Float128(BigInt(x) + BigInt(y))
end
end

@test Base.exponent_one(Float128) == reinterpret(UInt128, Float128(1.0))

@testset "flipsign" begin
    x = Float128( 2.0)
    y = Float128(-2.0)
    @test x == flipsign(y, -one(Float128))
    @test y == flipsign(y,  1)
end

@testset "arithmetic" begin
    fpi = Float128(pi)
    finvpi = inv(fpi)
    @test (fpi + 3) - fpi == 3
    @test fpi * finvpi === one(Float128)
    @test finvpi / fpi == finvpi^2
end

@testset "modf" begin
    x = Float128(pi)
    fpart, ipart = modf(x)
    @test x == ipart + fpart
    @test signbit(fpart) == signbit(ipart) == false

    y = Float128(-pi)
    fpart, ipart = modf(y)
    @test y == ipart + fpart
    @test signbit(fpart) == signbit(ipart) == true

    z = x^3
    fpart, ipart = modf(x) .+ modf(y)
    @test x+y == ipart+fpart
end

isnan128(x) = isa(x, Float128) && isnan(x)
isinf128(x) = isa(x, Float128) && isinf(x)

@testset "nonfinite" begin
    Zero = Float128(0)
    One = Float128(1)
    huge = floatmax(Float128)
    myinf = huge + huge
    myminf = -myinf
    @test isinf128(myinf)
    @test isnan128(Zero / Zero)
    @test isinf128(One / Zero)
    @test isnan128(myinf - myinf)
    @test isnan128(myinf + myminf)
    @test Inf128 === myinf
    @test typemax(Float128) === myinf
    @test typemin(Float128) === myminf
end

@testset "transcendental etc. calls" begin
    # at least enough to cover all the wrapping code
    @testset "real" begin
        x = sqrt(Float128(2.0))
        xd = Float64(x)
        @test (x^Float128(4.0)) ≈ Float128(4.0)
        @test exp(x) ≈ exp(xd)
        @test abs(x) == x
        @test hypot(Float128(3),Float128(4)) == Float128(5)
        @test atan(x,x) ≈ Float128(pi) / 4
        if !Sys.iswindows()
            @test fma(x,x,Float128(-1.0)) ≈ Float128(1)
        end
    end
    @testset "complex" begin
        x = sqrt(ComplexF128(1.0 + 1.0im))
        xd = ComplexF64(x)
        @test x^x ≈ xd^xd
        @test exp(x) ≈ exp(xd)
        @test abs(x) ≈ abs(xd)
        @test sin(x) ≈ sin(xd)
        @test cos(x) ≈ cos(xd)
        @test tan(x) ≈ tan(xd)
        @test log(x) ≈ log(xd)
    end
end

@testset "misc. math" begin
    x = sqrt(Float128(2.0))
    @test abs(x^(-2) - Float128(0.5)) < 1.0e-32
    m = maxintfloat(Float128)
    @test m+one(Float128) == m
    @test m-one(Float128) != m
end

function hist(X, n)
    v = zeros(Int, n)
    for x in X
        v[floor(Int, x*n) + 1] += 1
    end
    v
end

@testset "random" begin
    # test for sanity and coarse uniformity
    @test typeof(rand(Float128)) == Float128
    for rng in [MersenneTwister(), RandomDevice()]
        counts = hist(rand(rng, Float128, 2000), 4)
        @test minimum(counts) > 300
        counts = hist([rand(rng, Float128) for i in 1:2000], 4)
        @test minimum(counts) > 300
    end
end

@testset "string conversion" begin
    s = string(Float128(3.0))
    p = r"3\.0+e\+0+"
    m = match(p, s)
    @test (m != nothing) && (m.match == s)
    @test parse(Float128,"3.0") == Float128(3.0)
end

@testset "irrationals" begin
    tiny = 2eps(Float128(1))
    @test abs(cos(Float128(pi)) + 1) < tiny
    @test abs(log(Float128(ℯ)) - 1) < tiny
    @test abs((2*Float128(MathConstants.golden) - 1)^2 - 5) < 5 * tiny
end

include("specfun.jl")

include("printf.jl")


text:
getrun() = loadtable(run_file_name, indexcols = :run)

function getdate()
    @label askdatetime
    println("Specify a date for this run\n[Enter: $(defaults[:date])]")
    _dt = strip(readline(stdin))
    dt = isempty(_dt) ? defaults[:datetime] : tryparsedatetime(_dt)
    if dt ≡ nothing
        @warn "the format of the date is wrong, try something like:" Date(now())
        @goto askdatetime
    end
    defaults[:date] = dt
    dt
end

function register_run()
    cdata = OrderedDict(k => "" for k in COMMON_FACTORS)
    run = getrun()
    experiment = getexperiment()
    cdata[:run] = uuid4()
    options = select(experiment, :experiment)
    defaults[:experiment] = first(select(experiment, :experiment))
    pushfirst!(options, defaults[:experiment])
    menu = RadioMenu(options)
    i = request("To which experiment does this run belong to?", menu)
    if i == 1
        i = findlast(isequal(menu.options[1]), menu.options)
    end
    cdata[:experiment] = experiment[i - 1].experiment
    cdata[:date] = getdate()
    println("Any comments about this run:")
    cdata[:comment] = strip(readline(stdin))
    println("Who was the experimenter?\n[Enter: $(defaults[:person])]")
    _person = strip(readline(stdin))
    cdata[:person] = isempty(_person) ? defaults[:person] : _person
    defaults[:person] = person
    data = OrderedDict(k => "" for k in FACTORS)
    for factor in FACTORS
        println("Specify $factor")
    _level = strip(readline(stdin))
    data[:factor] = isempty(_person) ? defaults[:person] : _person
    defaults[:person] = person
        data[factor] = strip(readline(stdin))
    end
    [(; cdata..., data...)] |> CSV.write(run_file_name, append = true)
end


text:
"""
    ImpactContact{T,N} <: Contact{T,N}

    contact object for impact (i.e., no friction)

    surface_normal_projector: mapping from world frame to surface tangent frame
    contact_point: position of contact on Body relative to center of mass
    offset: position of contact relative to contact_point
"""
mutable struct ImpactContact{T,N} <: Contact{T,N}
    surface_normal_projector::Adjoint{T,SVector{3,T}}
    contact_point::SVector{3,T}
    offset::SVector{3,T}

    function ImpactContact(body::Body{T}, normal::AbstractVector;
        contact_point=szeros(T, 3),
        offset::AbstractVector=szeros(T, 3)) where T
        V1, V2, V3 = orthogonal_columns(normal) #
        A = [V1 V2 V3]
        Ainv = inv(A)
        surface_normal_projector = Ainv[3,SA[1; 2; 3]]'
        new{Float64,2}(surface_normal_projector, contact_point, offset)
    end
end

function constraint(mechanism, contact::ContactConstraint{T,N,Nc,Cs}) where {T,N,Nc,Cs<:ImpactContact{T,N}}
    model = contact.model
    body = get_body(mechanism, contact.parent_id)
    x3, q3 = next_configuration(body.state, mechanism.timestep)
    SVector{1,T}(model.surface_normal_projector * (x3 + vector_rotate(model.contact_point,q3) - model.offset) - contact.impulses_dual[2][1])
end

function constraint_jacobian_configuration(model::ImpactContact, x3::AbstractVector, q3::Quaternion,
    x2::AbstractVector, v25::AbstractVector, q2::Quaternion, ϕ25::AbstractVector,
    λ, timestep)
    X = model.surface_normal_projector
    Q = model.surface_normal_projector * ∂vector_rotate∂q(model.contact_point, q3)
    return [X Q]
end

function constraint_jacobian_velocity(model::ImpactContact, x3::AbstractVector, q3::Quaternion,
    x2::AbstractVector, v25::AbstractVector, q2::Quaternion, ϕ25::AbstractVector,
    λ, timestep)
    V = model.surface_normal_projector * timestep
    Ω = model.surface_normal_projector * ∂vector_rotate∂q(model.contact_point, q3) * rotational_integrator_jacobian_velocity(q2, ϕ25, timestep)
    return [V Ω]
end

function force_mapping(model::ImpactContact, x::AbstractVector, q::Quaternion)
    X = model.surface_normal_projector
    return X
end

function set_matrix_vector_entries!(mechanism::Mechanism, matrix_entry::Entry, vector_entry::Entry,
    contact::ContactConstraint{T,N,Nc,Cs,N½}) where {T,N,Nc,Cs<:ImpactContact{T,N},N½}
    # ∇impulses[dual .* impulses - μ; g - s] = [diag(dual); -diag(0,1,1)]
    # ∇dual[dual .* impulses - μ; g - s] = [diag(impulses); -diag(1,0,0)]
    γ = contact.impulses[2] +  REG * neutral_vector(contact.model)
    s = contact.impulses_dual[2] + REG * neutral_vector(contact.model)

    ∇s = hcat(γ, -Diagonal(sones(N½)))
    ∇γ = hcat(s, -Diagonal(szeros(N½)))
    matrix_entry.value = hcat(∇s, ∇γ)

    # [-dual .* impulses + μ; -g + s]
    vector_entry.value = vcat(-complementarityμ(mechanism, contact), -constraint(mechanism, contact))
    return
end


text:
@testset "LinearCombination" begin
    @testset "scalar 1d" begin
        @testset for TC in rnc((Int,Float64))
            C = rand(TC,3)

            p = @inferred(LinearCombination(C, Monomials(3)))
            @test ndims(p) == ndims(typeof(p)) == 1

            @test @inferred(p(1)) ≈ sum(C)
            @test @inferred(p([1,2])) ≈ [sum(C), sum(C.*2 .^ (0:2))]
        end
    end

    @testset "scalar 2d" begin
        @testset for TC in rnc((Int,Float64))
            C = rand(TC,2,2)

            p = @inferred(LinearCombination(C, Monomials.((2,2))))
            @test @inferred(LinearCombination(C, Monomials(2))) == p
            @test ndims(p) == ndims(typeof(p)) == 2

            @test @inferred(p(1,1)) ≈ sum(C)
            @test @inferred(p(1,1:2)) ≈ [1,1]'*C*[1 1; 1 2]
            @test @inferred(p(1:2,1:2)) ≈ [1 1; 1 2]'*C*[1 1; 1 2]
        end
    end

    @testset "matrix 1d" begin
        C = rand(5)
        p = LinearCombination(C,Monomials(length(C)))
        x = LinRange(-1,1,11)
        M = Diagonal(collect(x))
        @test diag(@inferred(p(M))) == p.(x)
    end

    @testset "matrix times vector 1d" begin
        C = rand(5)
        p = LinearCombination(C,Monomials(length(C)))
        x = LinRange(-1,1,11)
        M = Diagonal(collect(x))
        v = one.(x)
        @test @inferred(((M,v)->@evaluate p(M) * v)(M,v)) == p.(x)
    end
end


text:
"""
    get_author_title(pdf::String) -> Tuple{String, String}

Extract the author and title from the PDF.

# Arguments
- `pdf::String`: $(PDF_ARGUMENT)

# Returns
- `Tuple{String, String}`: the author and title

# Throws
- [`FileDoesNotExist`](@ref): $(FILE_DOES_NOT_EXIST_EXCEPTION)
- [`NotPDF`](@ref): $(NOT_PDF_EXCEPTION)

# Example
```jldoctest; output = false
using PDFHighlights

path_to_pdf = joinpath(
    pathof(PDFHighlights) |> dirname |> dirname,
    "test",
    "pdf",
    "TestPDF.pdf",
)

get_author_title(path_to_pdf) == ("Pavel Sobolev", "A dummy PDF for tests")

# output

true
```
"""
function get_author_title(pdf::String)::Tuple{String, String}

    !isfile(pdf) && throw(FileDoesNotExist(pdf))
    !endswith(pdf, ".pdf") && throw(NotPDF(pdf))

    # Compute the check sum
    checksum = ""
    open(pdf, "r") do io
        checksum = bytes2hex(sha1(io))
    end

    get!(GET_AUTHOR_TITLE_OUTPUTS, checksum) do

        author_ref = Ref{Cstring}()
        title_ref = Ref{Cstring}()

        ccall(
            (:get_author_title, PDFHighlightsWrapper),
            Cvoid,
            (Cstring, Ref{Cstring}, Ref{Cstring}),
            pdf, author_ref, title_ref,
        )

        return unsafe_string(author_ref[]), unsafe_string(title_ref[])

    end

end


text:
module Utils

export is_integer, is_alpha

# include("helper.jl")

"""
    is_integer(s)
Return Whether or not the string s is an integer.
"""
function is_integer(s)
    return isa(tryparse(Int, s), Int)
end

"""
    is_alpha(s::String)
    Faster version of isletter function. Only works for ASCI, not unicode.
Return True or False.
"""
function is_alpha(text)
    all(c -> 'a' <= c <= 'z' || 'A' <= c <= 'Z', text)
end

end

max_stars_repo_name: ChristianMichelsen/reference-free

max_issues_repo_name: ChristianMichelsen/reference-free

max_forks_repo_name: ChristianMichelsen/reference-free


text:
using MRMP: dist

@testset "dist" begin
    p1_1 = [-1.0, -2.0]
    p1_2 = [1.0, 2.0]
    p2_1 = [1.0, -2.0]
    p2_2 = [-1.0, 2.0]

    @test dist(p1_1, p1_2, p2_1, p2_2) == 0
end


text:
using Crystallography:
    PrimitiveCubic,
    FaceCenteredCubic,
    BodyCenteredCubic,
    PrimitiveHexagonal,
    RCenteredHexagonal,
    PrimitiveTetragonal,
    BodyCenteredTetragonal,
    PrimitiveOrthorhombic,
    BCenteredOrthorhombic,
    ACenteredOrthorhombic,
    FaceCenteredOrthorhombic,
    BodyCenteredOrthorhombic,
    PrimitiveMonoclinic,
    CCenteredMonoclinic,
    BCenteredMonoclinic,
    PrimitiveTriclinic

import Crystallography: Bravais, Lattice

export Ibrav

struct Ibrav
    i::Int
    function Ibrav(i)
        @assert i in union(0:1:14, (-3, -5, -9, 91, -12, -13)) "Bravais lattice undefined for `ibrav = $i`!"
        return new(i)
    end
end

function Bravais(ibrav::Ibrav)
    i = ibrav.i
    if i == 1
        return PrimitiveCubic(true)
    elseif i == 2
        return FaceCenteredCubic(true)
    elseif i == 3
        return BodyCenteredCubic(true)
    elseif i == -3
        return BodyCenteredCubic(false)
    elseif i == 4
        return PrimitiveHexagonal(true)
    elseif i == 5
        return RCenteredHexagonal(true)
    elseif i == -5
        return RCenteredHexagonal(false)
    elseif i == 6
        return PrimitiveTetragonal(true)
    elseif i == 7
        return BodyCenteredTetragonal(true)
    elseif i == 8
        return PrimitiveOrthorhombic(true)
    elseif i == 9
        return BCenteredOrthorhombic(true)
    elseif i == -9
        return BCenteredOrthorhombic(false)
    elseif i == 91
        return ACenteredOrthorhombic(true)  # New in QE 6.5
    elseif i == 10
        return FaceCenteredOrthorhombic(true)
    elseif i == 11
        return BodyCenteredOrthorhombic(true)
    elseif i == 12
        return PrimitiveMonoclinic(true)
    elseif i == -12
        return PrimitiveMonoclinic(false)
    elseif i == 13
        return CCenteredMonoclinic(true)
    elseif i == -13
        return BCenteredMonoclinic(true)  # New in QE 6.5
    elseif i == 14
        return PrimitiveTriclinic(true)
    else
        throw(ArgumentError("ambiguous `ibrav` number $i is given!"))
    end
end

"""
    Lattice(::Bravais, p)

Create a Bravais lattice from the exact lattice type and cell parameters `p` (not `celldm`!).

The first elements of `p` are `a`, `b`, `c`; the last 3 are `α`, `β`, `γ` (in radians).
"""
Lattice(::PrimitiveCubic, p) = Lattice(p[1] * [[1, 0, 0], [0, 1, 0], [0, 0, 1]]...)
Lattice(::FaceCenteredCubic, p) = Lattice(p[1] / 2 * [[-1, 0, 1], [0, 1, 1], [-1, 1, 0]]...)
function Lattice(bravais::BodyCenteredCubic, p)
    if bravais.obverse
        return Lattice(p[1] / 2 * [[1, 1, 1], [-1, 1, 1], [-1, -1, 1]]...)
    else  # -3
        return Lattice(p[1] / 2 * [[-1, 1, 1], [1, -1, 1], [1, 1, -1]]...)
    end
end
Lattice(::PrimitiveHexagonal, p) =
    Lattice(p[1] * [[1, 0, 0], [-1 / 2, √3 / 2, 0], [0, 0, p[3]]]...)
function Lattice(bravais::RCenteredHexagonal, p)
    cosγ = p[4]
    ty = sqrt((1 - cosγ) / 6)
    tz = sqrt((1 + 2cosγ) / 3)
    if bravais.obverse
        tx = sqrt((1 - cosγ) / 2)
        return Lattice(p[1] * [[tx, -ty, tz], [0, 2ty, tz], [-tx, -ty, tz]]...)
    else  # -5
        a′ = p[1] / √3
        u = tz - 2 * √2 * ty
        v = tz + √2 * ty
        return Lattice(a′ * [[u, v, v], [v, u, v], [v, v, u]]...)
    end
end
Lattice(::PrimitiveTetragonal, p) = Lattice(p[1] * [[1, 0, 0], [0, 1, 0], [0, 0, p[3]]]...)
function Lattice(::BodyCenteredTetragonal, celldm)
    r = celldm[3]
    return Lattice(celldm[1] / 2 * [[1, -1, r], [1, 1, r], [-1, -1, r]]...)
end
Lattice(::PrimitiveOrthorhombic, p) =
    Lattice(p[1] * [[1, 0, 0], [0, p[2], 0], [0, 0, p[3]]]...)
function Lattice(bravais::BCenteredOrthorhombic, p)
    a, b, c = p[1] .* (1, p[2], p[3])
    if bravais.obverse
        return Lattice([[a / 2, b / 2, 0], [-a / 2, b / 2, 0], [0, 0, c]]...)
    else
        return Lattice([[a / 2, -b / 2, 0], [a / 2, b / 2, 0], [0, 0, c]]...)
    end
end
function Lattice(::ACenteredOrthorhombic, p)
    a, r1, r2 = p[1:3]
    return Lattice(a * [[1, 0, 0], [0, r1 / 2, -r2 / 2], [0, r1 / 2, r2 / 2]]...)
end  # New in QE 6.4
function Lattice(::FaceCenteredOrthorhombic, p)
    a, b, c = p[1], p[1] * p[2], p[1] * p[3]
    return Lattice(1 / 2 * [[a, 0, c], [a, b, 0], [0, b, c]]...)
end
function Lattice(::BodyCenteredOrthorhombic, p)
    a, b, c = p[1] .* (1, p[2], p[3])
    return Lattice(1 / 2 * [[a, b, c], [-a, b, c], [-a, -b, c]]...)
end
function Lattice(bravais::PrimitiveMonoclinic, p)
    if bravais.obverse
        a, r1, r2, cosγ = p[1:4]
        return Lattice(a * [[1, 0, 0], [r1 * cosγ, r1 * sin(acos(cosγ)), 0], [0, 0, r2]]...)
    else
        a, r1, r2, _, cosβ = p[1:5]
        return Lattice(a * [[1, 0, 0], [0, r1, 0], [r2 * cosβ, 0, r2 * sin(acos(cosβ))]]...)
    end
end
function Lattice(::CCenteredMonoclinic, p)
    a, r1, r2, cosγ = p[1:4]
    return Lattice(
        a *
        [[1 / 2, 0, -r2 / 2], [r1 * cosγ, r1 * sin(acos(cosγ)), 0], [1 / 2, 0, r2 / 2]]...,
    )
end
function Lattice(::BCenteredMonoclinic, p)
    a, r1, r2, _, cosβ = p[1:3]
    return Lattice(
        a *
        [[1 / 2, r1 / 2, 0], [-1 / 2, r1 / 2, 0], [r2 * cosβ, 0, r2 * sin(acos(cosβ))]]...,
    )
end
function Lattice(::PrimitiveTriclinic, p)
    a, r1, r2, cosα, cosβ, cosγ = p[1:6]  # Every `p` that is an iterable can be used
    sinγ = sin(acos(cosγ))
    δ = r2 * sqrt(1 + 2 * cosα * cosβ * cosγ - cosα^2 - cosβ^2 - cosγ^2) / sinγ
    return Lattice(
        a * [
            [1, 0, 0],
            [r1 * cosγ, r1 * sinγ, 0],
            [r2 * cosβ, r2 * (cosα - cosβ * cosγ) / sinγ, δ],
        ]...,
    )
end


text:
function isshaped(samples::BAT.DensitySampleVector)
    isa(varshape(samples), NamedTupleShape) ? (return true) : (return false)
end

function asindex(samples::BAT.DensitySampleVector, name::Union{Expr, Symbol})
    if isshaped(samples)
        return asindex(varshape(samples), name)
    else
        throw(ArgumentError("Samples are unshaped. Key :$name cannot be matched. Use index instead."))
    end
end

function asindex(
    x::Union{DensitySampleVector, NamedTupleDist, MarginalDist},
    key::Integer
)
    return key
end

#for MarginalDist
function asindex(marg::MarginalDist, name::Union{Expr, Symbol})
    idx = asindex(marg.origvalshape, name)
    if idx in marg.dims
        return idx
    else
        throw(ArgumentError("Key :$name not in MarginalDist"))
    end
end

function getstring(samples::BAT.DensitySampleVector, idx::Integer)
    if isshaped(samples)
        vs = varshape(samples)
        names = all_active_names(vs)
        return names[idx]
    else
        throw(ArgumentError("Samples are unshaped. Key :$name cannot be matched. Use index instead."))
    end
end

function getstring(marg::MarginalDist, idx::Integer)
    vs = marg.dist.shape
    names = all_active_names(vs)
    if checkbounds(Bool, names, idx)
        return names[idx]
    else
        throw(ArgumentError("Index $idx not in $marg"))
    end
end


text:
"""
    air_density(Tair,pressure; ...)

Air density of moist air from air temperature and pressure_

# Arguments
- Tair:      Air temperature (deg C)
- pressure:  Atmospheric pressure (kPa)
optional
- `constants=`[`BigleafConstants`](@ref)`()`: Dictionary with entries
   Kelvin, kPa2Pa

# Details
Air density ``\\rho`` is calculated as:

  ``\\rho = {pressure \\over Rd * Tair}``

# Value
air density (kg m-3)

# Examples
```@example doc
# air density at 25degC and standard pressure (101.325kPa)
air_density(25,101.325)
```

# References
Foken, T, 2008: Micrometeorology. Springer, Berlin, Germany.
"""
function air_density(Tair,pressure; constants=BigleafConstants())
  Tair_K     = Tair + oftype(Tair,constants.Kelvin)
  pressure_Pa = pressure * oftype(pressure,constants.kPa2Pa)
  rho = pressure_Pa / (oftype(Tair_K,constants.Rd) * Tair_K)
end

"""
    pressure_from_elevation(elev,Tair,VPD=missing;...)

An estimate of mean pressure at a given elevation as predicted by the
             hypsometric equation.

# Arguments
- elev:      Elevation a_s_l_ (m)
- Tair:      Air temperature (deg C)
- VPD:       Vapor pressure deficit (kPa); optional
optional
- `constants=`[`BigleafConstants`](@ref)`()`: Dictionary with entries
    Kelvin, pressure0, Rd, g, Pa2kPa

# Details
Atmospheric pressure is approximated by the hypsometric equation:

``pressure = pressure_0 / (exp(g * elevation / (Rd Temp)))``

The hypsometric equation gives an estimate of the standard pressure
      at a given altitude.
      If VPD is provided, humidity correction is applied and the
      virtual temperature instead of air temperature is used_ VPD is
      internally converted to specific humidity.

# Value
Atmospheric pressure (kPa)

# References
Stull B_, 1988: An Introduction to Boundary Layer Meteorology.
            Kluwer Academic Publishers, Dordrecht, Netherlands.

# Examples
```@example doc
# mean pressure at 500m altitude at 25 deg C and VPD of 1 kPa
pressure_from_elevation(500,25,1)
```
"""
function pressure_from_elevation(elev,Tair,VPD=missing; constants=BigleafConstants())
  Tair_K     = Tair + oftype(Tair,constants.Kelvin)
  if ismissing(VPD)
    pressure = oftype(Tair,constants.pressure0) / exp(oftype(Tair,constants.g) * elev / (oftype(Tair,constants.Rd)*Tair_K))
  else
    pressure1   = oftype(Tair,constants.pressure0) / exp(oftype(Tair,constants.g) * elev / (oftype(Tair,constants.Rd)*Tair_K))
    Tv          = virtual_temp(
      Tair_K - oftype(Tair,constants.Kelvin), pressure1 * constants.Pa2kPa, VPD;
      Esat_formula=Sonntag1990(),constants) + oftype(Tair,constants.Kelvin)
    pressure    = oftype(Tair,constants.pressure0) / exp(oftype(Tair,constants.g) * elev / (oftype(Tair,constants.Rd)*Tv))
  end
  pressure = pressure * constants.Pa2kPa
end

"""
    psychrometric_constant(Tair,pressure; ...)

Computes the psychrometric 'constant'.

# Argumens
- Tair:      Air temperature (deg C)
- pressure:  Atmospheric pressure (kPa)
optional
- `constants=`[`BigleafConstants`](@ref)`()`: Dictionary with entries
    cp, eps

# Details
The psychrometric constant (``\\gamma``) is given as:

``\\gamma = cp * pressure / (eps * \\lambda)``

 where ``\\lambda`` is the latent heat of vaporization (J kg-1),
 as calculated from [`latent_heat_vaporization`](@ref).

# Value
the psychrometric constant (kPa K-1)

# References
Monteith J.L., Unsworth M.H., 2008: Principles of Environmental Physics.
            3rd Edition. Academic Press, London.

# Examples
```@example doc
psychrometric_constant.(5.0:5.0:25.0, 100)
```
"""
function psychrometric_constant(Tair,pressure; constants=BigleafConstants())
  lambda = latent_heat_vaporization(Tair)
  gamma  = (oftype(pressure,constants.cp) * pressure) / (oftype(lambda,constants.eps) * lambda)
end

"""
    latent_heat_vaporization(Tair)

Latent heat of vaporization as a function of air temperature
using

``\\lambda = (2.501 - 0.00237 \\, Tair) 10^6``.

# Arguments:
- Tair:  Air temperature (deg C)

# Value
``\\lambda``: Latent heat of vaporization (J kg-1)

# References
- Stull, B_, 1988: An Introduction to Boundary Layer Meteorology (p_641)
            Kluwer Academic Publishers, Dordrecht, Netherlands
- Foken, T, 2008: Micrometeorology_ Springer, Berlin, Germany

```@example
latent_heat_vaporization.(5:5:45)
```
"""
function latent_heat_vaporization(Tair)
  k1 = 2.501
  k2 = 0.00237
  lambda = ( k1 - k2 * Tair ) * 1e+06
end

"""
    wetbulb_temp(Tair, pressure, VPD; ...)

Calculate the wet bulb temperature, i_e_ the temperature
             that the air would have if it was saturated

# Arguments
- Tair:      Air temperature (deg C)
- pressure:  Atmospheric pressure (kPa)
- VPD:       Vapor pressure deficit (kPa)
optional
- accuracy:  Accuracy of the result (deg C)
- `Esat_formula`: formula used in [`Esat_from_Tair`](@ref)
- `constants=`[`BigleafConstants`](@ref)`()`: Dictionary with entries
     cp, eps,Pa2kPa

# Details
Wet-bulb temperature (Tw) is calculated from the following expression:

``e = Esat(Tw) - gamma* (Tair - Tw)``

The equation is optimized for Tw.
Actual vapor pressure e (kPa) is calculated from VPD using [`VPD_to_e`](@ref).
The psychrometric constant gamma (kPa K-1) is calculated using
[`psychrometric_constant`](@ref).

# Value
wet-bulb temperature (degC)

# References
Monteith J.L., Unsworth M.H., 2008: Principles of Environmental Physics.
            3rd edition. Academic Press, London.

# Examples
```@example doc
wetbulb_temp.([20,25.0], 100, [1,1.6])
```
"""
function wetbulb_temp(Tair, pressure, VPD; accuracy=1e-03,
  Esat_formula=Sonntag1990(), constants=BigleafConstants())
  gamma  = psychrometric_constant(Tair,pressure)
  ea     = VPD_to_e(VPD,Tair;Esat_formula)
  wetbulb_temp_from_e_Tair_gamma(ea,Tair,gamma; accuracy,Esat_formula,constants)
end,
function wetbulb_temp_from_e_Tair_gamma(ea, Tair, gamma; accuracy=1e-03,
  Esat_formula=Sonntag1990(), constants=BigleafConstants())
  if accuracy > one(accuracy)
    @warn ("'accuracy' is set to 1 degC")
    accuracy = one(accuracy)
  end
  fopt(Tw) = abs(ea - (Esat_from_Tair(Tw; Esat_formula = Esat_formula,constants) -
    0.93*gamma*(Tair - Tw)))
  resopt = optimize(fopt, -100, 100, abs_tol = accuracy)
  roundmult(resopt.minimizer, accuracy)
end

"""
    dew_point(Tair,VPD; ...)
    dew_point_from_e(ea; ...)

Calculate the dew point, the temperature to which air must be
             cooled to become saturated (i_e_ e = Esat(Td))

# Arguments
- Tair:     Air temperature (degC)
- VPD:      Vapor pressure deficit (kPa)
- ea:       actual water vapor pressure (kPa)
optional
- accuracy = 1e-03: Accuracy of the result (deg C)
- `Esat_formula`: formula used in [`Esat_from_Tair`](@ref)
- `constants=`[`BigleafConstants`](@ref)`()`: Dictionary with entries
    Pa2kPa

# Details
Dew point temperature (Td) is defined by the temperature for which
saturaed vapour pressure equals current vapour pressure, i.e. below which
water would start to condensate.

``e = Esat(Td)``

where e is vapor pressure of the air and Esat is the vapor pressure deficit.
         This equation is solved for Td by optimization.

# Value
dew point temperature (degC)

# References
Monteith J.L., Unsworth M.H., 2008: Principles of Environmental Physics.
            3rd edition. Academic Press, London.

# Examples
```@example doc
Tair = 20.0
VPD = 1.5
Td = dew_point(Tair, VPD; accuracy = 1e-2)
(e = VPD_to_e(VPD,Tair), esat_Td = Esat_from_Tair(Td))
```
"""
function dew_point(Tair, VPD; accuracy=1e-03,Esat_formula=Sonntag1990(),
                      constants=BigleafConstants())
  ea = VPD_to_e(VPD,Tair;Esat_formula)
  dew_point_from_e(ea; accuracy, Esat_formula, constants)
end,
function dew_point_from_e(ea;accuracy=1e-03,Esat_formula=Sonntag1990(), constants=BigleafConstants())
  if accuracy > one(accuracy)
    @warn ("'accuracy' is set to 1 degC")
    accuracy = one(accuracy)
  end
  fopt(Td) = abs(ea - Esat_from_Tair(Td;Esat_formula = Esat_formula,constants))
  resopt = optimize(fopt, -100, 100, abs_tol = accuracy)
  roundmult(resopt.minimizer, accuracy)
end
# https://discourse.julialang.org/t/rounding-to-a-multiple-of-float/21295/21
roundmult(val, prec) = (inv_prec = 1 / prec; round(val * inv_prec) / inv_prec)

"""
    virtual_temp(Tair,pressure,VPD; ...)

Virtual temperature, defined as the temperature at which dry air would have the same
             density as moist air at its actual temperature.

# Arguments
- `Tair`:      Air temperature (deg C)
- `pressure`:  Atmospheric pressure (kPa)
- `VPD`:       Vapor pressure deficit (kPa)
- `Esat_formula`: formula used in [`Esat_from_Tair`](@ref)
optional
- `constants=`[`BigleafConstants`](@ref)`()`: Dictionary with entries
  - :Kelvin - conversion degree Celsius to Kelvin
  - :eps - ratio of the molecular weight of water vapor to dry air (-)

# Details
The virtual temperature is given by:

```Tv = Tair / (1 - (1 - eps) e/pressure)```

 where Tair is in Kelvin (converted internally)_ Likewise, VPD is converted
 to actual vapor pressure (e in kPa) with [`VPD_to_e`](@ref) internally.

# Value
virtual temperature (deg C)

# References
- Monteith J.L., Unsworth M.H., 2008: Principles of Environmental Physics.
            3rd edition. Academic Press, London.

# Examples
```jldoctest; output = false
Tair,pressure,VPD = 25.0,100.0,1.5
vt = virtual_temp(Tair,pressure,VPD)
≈(vt, 26.9, atol =0.1)
# output
true
```
"""
function virtual_temp(Tair,pressure,VPD;Esat_formula=Sonntag1990(),
                         constants=BigleafConstants())
  e    = VPD_to_e(VPD,Tair;Esat_formula)
  Tair_Kelvin = Tair + oftype(Tair,constants.Kelvin)
  Tv_Kelvin = Tair_Kelvin / (1 - (1 - oftype(pressure,constants.eps)) * e/pressure)
  Tv = Tv_Kelvin - oftype(Tair,constants.Kelvin)
end

"""
    kinematic_viscosity(Tair,pressure; ...)

Calculate the kinematic viscosity of air.

# Parameters
- Tair      Air temperature (deg C)
- pressure  Atmospheric pressure (kPa)
optional
- `constants=`[`BigleafConstants`](@ref)`()`: Dictionary with entries
    Kelvin, pressure0, Tair0, kPa2Pa

# Details
Eq where v is the kinematic viscosity of the air (m2 s-1),
         given by (Massman 1999b):

``v = 1.327 * 10^-5(pressure0/pressure)(Tair/Tair0)^{1.81}``

# Value
kinematic viscosity of air (m2 s-1)

# References
Massman, W.J., 1999b: Molecular diffusivities of Hg vapor in air,
            O2 and N2 near STP and the kinematic viscosity and thermal diffusivity
            of air near STP. Atmospheric Environment 33, 453-457.

# Examples
```jldoctest; output = false
Tair,pressure = 25.0,100.0
vis = kinematic_viscosity(Tair,pressure)
≈(vis, 1.58e-5, atol =1e-7)
# output
true
```
"""
function kinematic_viscosity(Tair,pressure; constants=BigleafConstants())
  (ismissing(Tair) || ismissing(pressure)) && return(missing)
  Tair_Kelvin     = Tair + oftype(Tair, constants.Kelvin)
  pressure_kPa = pressure * oftype(pressure, constants.kPa2Pa)
  v  = oftype(Tair,1.327e-05)*(oftype(pressure,constants.pressure0)/pressure_kPa) *
    (Tair_Kelvin/oftype(Tair,constants.Tair0))^oftype(Tair,1.81)
end


text:
export generateCanonicalFG_Honeycomb!

global _honeycombRecipe = Dict{Symbol,Symbol}(
  :l6   =>:l0,
  :l18  =>:l0,
  :l42  =>:l0,
  :l57  =>:l0,
  :l63  =>:l0,
  :l70  =>:l0,
  :l10  =>:l1,
  :l50  =>:l1,
  :l56  =>:l1,
  :l69  =>:l1,
  :l43  =>:l2,
  :l49  =>:l2,
  :l61  =>:l2,
  :l36  =>:l3,
  :l54  =>:l3,
  :l28  =>:l4,
  :l34  =>:l4,
  :l47  =>:l4,
  :l12  =>:l5,
  :l19  =>:l5,
  :l26  =>:l5,
  :l33  =>:l5,
  :l40  =>:l5,
  :l13  =>:l7,
  :l25  =>:l7,
  :l17  =>:l8,
  :l58  =>:l9,
  :l64  =>:l9,
  :l35  =>:l11,
  :l41  =>:l11,
  :l48  =>:l11,
  :l55  =>:l11,
  :l62  =>:l11,
  :l20  =>:l14,
  :l32  =>:l14,
  :l24  =>:l15,
  :l27  =>:l21,
  :l39  =>:l21,
  :l31  =>:l22,
  :l38  =>:l29,
  :l46  =>:l37,
  :l53  =>:l44,
  :l60  =>:l51,

  # pose offset legs
  :x41  =>:left,
  :x63  =>:left,
  :x78  =>:left,
)

# piecewise posecount construction
# 7,14,21,28,35,43,50

function _addLandmarkBeehive!(fg,
                              lastPose::Symbol;
                              refKey::Symbol=:simulated,
                              solvable::Int=1,
                              graphinit::Bool=true,
                              landmarkRegex::Regex=r"l\d+",
                              srcNumber::Integer = match(r"\d+", string(lastPose)).match |> x->parse(Int,x),
                              atol::Real=1,
                              _doHoneycomb::Bool=false )
  #
  newFactor = RoME.Pose2Point2BearingRange(Normal(0,0.03), Normal(20,0.5))
  isAlready, simPPE, genLabel = IIF._checkVariableByReference(fg, lastPose, landmarkRegex, RoME.Point2, newFactor; atol=atol)

  # FIXME, oddball option that should be refactored (breaks Beehive usecase if true)
  if _doHoneycomb
    # force isAlready until fixed _checkVariableByReference parametric solution at -pi Optim issue
    global _honeycombRecipe
    isAlready = false
    genLabel = Symbol(:l, srcNumber)
    if haskey(_honeycombRecipe, genLabel)
      isAlready = true
      genLabel = _honeycombRecipe[genLabel]
      simPPE = getPPE(fg, genLabel, refKey)
    end
  end

  # maybe add new variable
  if !isAlready
    @info "New variable with simPPE" genLabel round.(simPPE.suggested,digits=2)
    newVar = addVariable!(fg, genLabel, RoME.Point2, solvable=solvable, tags=[:LANDMARK;])
    addFactor!(fg, [lastPose; genLabel], newFactor, solvable=solvable, graphinit=graphinit)

    # also set :simulated PPE for similar future usage
    setPPE!(newVar, :simulated, typeof(simPPE), simPPE)   # TODO this API can be improved
  else
    @info "Adding simulated loop closure with perfect data association" lastPose genLabel
    addFactor!(fg, [lastPose; genLabel], newFactor, solvable=solvable, graphinit=graphinit)
  end

  #
  return genLabel
end

function _driveHex!(fgl::AbstractDFG,
                    posecount::Int;
                    poseCountTarget::Real=Inf,
                    graphinit::Bool=false,
                    gaugePrior::Symbol=:x0f1,
                    refKey::Symbol=:simulated,
                    addLandmarks::Bool=true,
                    landmarkSolvable::Int=1,
                    _doHoneycomb::Bool=false,
                    atol::Real=1,
                    postpose_cb::Function=(fg_,latestpose)->()  )
  #

  # Drive around in a hexagon
  for i in (posecount):(posecount+5)
    if poseCountTarget <= posecount
      break
    end
    psym = Symbol("x$i")
    # nsym = Symbol("x$(i+1)")
    pp = Pose2Pose2(MvNormal([10.0;0;pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
    posecount += 1
    v_n = _addPoseCanonical!(fgl, psym, posecount, pp, refKey=refKey, graphinit=graphinit, postpose_cb=postpose_cb)
    nsym = getLabel(v_n)
    # add a new landmark (if not yet present)
    !addLandmarks ? nothing : _addLandmarkBeehive!(fgl, nsym; refKey=refKey, solvable=landmarkSolvable, atol=atol, graphinit=false, _doHoneycomb=_doHoneycomb)
  end

  return posecount
end

function _offsetHexLeg( fgl::AbstractDFG,
                        posecount::Int;
                        poseCountTarget::Real=Inf,
                        direction::Symbol=:right,
                        graphinit::Bool=false,
                        psym = Symbol("x$(posecount)"),
                        # nsym = Symbol("x$(posecount+1)"),
                        refKey::Symbol=:simulated,
                        guagePrior::Symbol=:x0f1,
                        addLandmarks::Bool=true,
                        landmarkSolvable::Int=1,
                        atol::Real=1,
                        _doHoneycomb::Bool=false,
                        postpose_cb::Function=(fg_,latestpose)->()   )
  #
  #skip out early if pose count target has been reached
  if poseCountTarget <= posecount
    return posecount
  end
  dirsign = if direction == :right
    -1
  elseif direction == :left
    +1
  else
    error("unknown direction symbol $direction")
  end
  pp = Pose2Pose2(MvNormal([10.0;0;dirsign*pi/3], diagm([0.1;0.1;0.1].^2)))

  posecount += 1
  v_n = _addPoseCanonical!(fgl, psym, posecount, pp, refKey=refKey, graphinit=graphinit, postpose_cb=postpose_cb)
  nsym = getLabel(v_n)

  # add a new landmark (if not yet present)
  !addLandmarks ? nothing : _addLandmarkBeehive!(fgl, nsym; refKey=refKey, solvable=landmarkSolvable, atol=atol, graphinit=false, _doHoneycomb=_doHoneycomb )

  return posecount
end

function generateCanonicalFG_Honeycomb!(poseCountTarget::Int=36;
                                        graphinit::Bool = false,
                                        dfg::AbstractDFG = LightDFG{SolverParams}(solverParams=SolverParams(graphinit=graphinit)),
                                        useMsgLikelihoods::Bool=getSolverParams(dfg).useMsgLikelihoods,
                                        direction::Symbol = :right,
                                        solvable::Int=1,
                                        refKey::Symbol=:simulated,
                                        addLandmarks::Bool=true,
                                        landmarkSolvable::Int=0,
                                        atol::Real=1,
                                        postpose_cb::Function=(fg_,latestpose)->()     )
  #
  global _honeycombRecipe

  # does anything exist in the graph yet
  posecount = if :x0 in ls(dfg)
    # what is the last pose
    lastPose = (ls(dfg, r"x\d+") |> sortDFG)[end]
    # get latest posecount number
    match(r"\d+", string(lastPose)).match |> x->parse(Int,x)
  else
    # initial zero pose
    generateCanonicalFG_ZeroPose(dfg=dfg, varType=RoME.Pose2, graphinit=graphinit, postpose_cb=postpose_cb) # , μ0=[0;0;1e-5] # tried for fix NLsolve on wrap issue

    # # reference ppe on :x0
    # refVal = zeros(3)
    # ppe = DFG.MeanMaxPPE(refKey, refVal, refVal, refVal)
    # setPPE!(dfg[:x0], refKey, DFG.MeanMaxPPE, ppe)

    # add a new landmark (if not yet present)
    !addLandmarks ? nothing : _addLandmarkBeehive!(dfg, :x0; refKey=refKey, solvable=landmarkSolvable, atol=atol, _doHoneycomb=true, graphinit=false)

    # staring posecount (i.e. :x0)
    0
  end

  # keep adding poses until the target number is reached
  while posecount < poseCountTarget
    posecount = _driveHex!(dfg, posecount; graphinit=graphinit, landmarkSolvable=landmarkSolvable, atol=atol, poseCountTarget=poseCountTarget, postpose_cb=postpose_cb, _doHoneycomb=true)
    # drive the offset legs
    lastPose = Symbol(:x, posecount)
    if haskey(_honeycombRecipe, lastPose)
      posecount = _offsetHexLeg(dfg, posecount; direction=_honeycombRecipe[lastPose], graphinit=graphinit, atol=atol, landmarkSolvable=landmarkSolvable, poseCountTarget=poseCountTarget, postpose_cb=postpose_cb, _doHoneycomb=true)
    end
    posecount = _offsetHexLeg(dfg, posecount; direction=direction, graphinit=graphinit, atol=atol, landmarkSolvable=landmarkSolvable, poseCountTarget=poseCountTarget, postpose_cb=postpose_cb, _doHoneycomb=true)
  end

  # NOTE solvable forced for everything at this time
  setSolvable!.(dfg, ls(dfg),  solvable)
  setSolvable!.(dfg, lsf(dfg), solvable)

  return dfg
end


text:
module NeuralStyle

using BSON: @load, @save
using Dates
using Flux
using Images
using Plots
using Random

# Exports from BSON
export @load, @save

# Basic utilities
include("utils.jl")
export nstrides1d, nstrides2d, pad1d, pad2d, tensor2image, image2tensor

# Provides locations for included contents and styles
include("datasets.jl")
export list_contents, list_styles

# Homegrown models
include("models/abstractmodel.jl")
export AutoEncoderModel
include("models/aem5.jl")
export AEM5
include("models/pvgg.jl")
export pVGG

end # module


text:
module LazySequences

export seq, first, rest, cons, Seqable, @lazyseq, @lazycat, cat, map, zip, take

import Base.start, Base.next, Base.done, Base.map, Base.zip, Base.cat,
       Base.first

abstract Seqable

seq(s::Seqable) = s
seq(::Nothing) = nothing

# Cons
# ----

type Cons <: Seqable
    first

    # This must be something with a seq function.
    rest
end

first(s::Cons) = s.first

function rest(s::Cons)
    # This effectively caches lazy-lists.
    s.rest = seq(s.rest)
end

cons(first, rest) = Cons(first, rest)

# Common operations
# -----------------

type Cat <: Seqable
    a::Seqable
    b # Something with a seq function.
end

first(s::Cat) = first(s.a)

function rest(s::Cat)
    a = rest(s.a)
    if a === nothing
        s.b = seq(s.b) # caching lazy seqs
    else
        Cat(a, s.b)
    end
end

function getindex(s::Cat, n::Int64)
    sequence = Any[]
    for i in take(n, s)
        push!(sequence, i)
    end
    sequence[end]
end

cat(a) = seq(a)
cat(a, b) = Cat(seq(a), b)
can(a, b, cs...) = Cat(seq(a), cat(b, cs...))

immutable Take <: Seqable
    n::Int
    a::Seqable
end

first(s::Take) = first(s.a)

function rest(s::Take)
    a = rest(s.a)
    a === nothing || s.n <= 1 ? nothing : Take(s.n - 1, a)
end

function take(n::Int, a)
    if n <= 0
        nothing
    else
        Take(n, seq(a))
    end
end

# Lazy-seq
# --------

immutable LazySeq
    realize::Function
end

seq(s::LazySeq) = s.realize()

macro lazyseq(body)
    quote
        LazySeq(() -> $(esc(body)))
    end
end

# Lazy cat
# --------

macro lazycat(ss...)
    quote
        cat($([:(@lazyseq $(s)) for s in ss]...))
    end
end

# Backward compatibility
# ----------------------

start(s::LazySeq) = s.realize()
start(s::Seqable) = s

function next(::Seqable, s::Seqable)
    first(s), rest(s)
end

done(::Seqable, s::Nothing) = true
done(::Seqable, s::Any) = false

# Realizations
# ------------

immutable AbstractArraySeq <: Seqable
    xs::AbstractArray
    i::Int
end

function seq(xs::AbstractArray)
    AbstractArraySeq(xs, 1)
end

first(s::AbstractArraySeq) = s.xs[s.i]
rest(s::AbstractArraySeq) =
    s.i >= length(s.xs) ? nothing :  AbstractArraySeq(s.xs, s.i + 1)

immutable Zip <: Seqable
    ts::Vector
end

function first(s::Zip)
    tuple({first(t) for t in s.ts}...)
end

function rest(s::Zip)
    ts = {rest(t) for t in s.ts}
    if any(t -> t === nothing, ts)
        nothing
    else
        Zip(ts)
    end
end

function zip(ss::Union(LazySeq, Seqable)...)
    ss = Seqable[seq(s) for s in ss]
    if any(s -> s === nothing, ss)
        nothing
    else
        Zip(ss)
    end
end

map(f::Function, ::Nothing) = nothing
map(f::Function, s::Seqable) = cons(f(first(s)), @lazyseq map(f, rest(s)))

function map(f::Function, s0::Seqable, s1::Seqable, ss::Seqable...)
    map(v -> apply(f, v), zip(s0, s1, ss...))
end

# Lazy character streams from IO objects.
# ---------------------------------------

function seq(io::IO)
    eof(io) ? nothing : cons(read(io, Char), @lazyseq seq(io))
end

function seqlines(io::IO)
    eof(io) ? nothing : cons(readline(io), @lazyseq seqlines(io))
end

end # end module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct CreateDataLakeStoreAccountParameters <: SwaggerModel
    location::Any # spec type: Union{ Nothing, String } # spec name: location
    tags::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: tags
    identity::Any # spec type: Union{ Nothing, EncryptionIdentity } # spec name: identity
    properties::Any # spec type: Union{ Nothing, CreateDataLakeStoreAccountProperties } # spec name: properties

    function CreateDataLakeStoreAccountParameters(;location=nothing, tags=nothing, identity=nothing, properties=nothing)
        o = new()
        validate_property(CreateDataLakeStoreAccountParameters, Symbol("location"), location)
        setfield!(o, Symbol("location"), location)
        validate_property(CreateDataLakeStoreAccountParameters, Symbol("tags"), tags)
        setfield!(o, Symbol("tags"), tags)
        validate_property(CreateDataLakeStoreAccountParameters, Symbol("identity"), identity)
        setfield!(o, Symbol("identity"), identity)
        validate_property(CreateDataLakeStoreAccountParameters, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type CreateDataLakeStoreAccountParameters

const _property_map_CreateDataLakeStoreAccountParameters = Dict{Symbol,Symbol}(Symbol("location")=>Symbol("location"), Symbol("tags")=>Symbol("tags"), Symbol("identity")=>Symbol("identity"), Symbol("properties")=>Symbol("properties"))
const _property_types_CreateDataLakeStoreAccountParameters = Dict{Symbol,String}(Symbol("location")=>"String", Symbol("tags")=>"Dict{String, String}", Symbol("identity")=>"EncryptionIdentity", Symbol("properties")=>"CreateDataLakeStoreAccountProperties")
Base.propertynames(::Type{ CreateDataLakeStoreAccountParameters }) = collect(keys(_property_map_CreateDataLakeStoreAccountParameters))
Swagger.property_type(::Type{ CreateDataLakeStoreAccountParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_CreateDataLakeStoreAccountParameters[name]))}
Swagger.field_name(::Type{ CreateDataLakeStoreAccountParameters }, property_name::Symbol) =  _property_map_CreateDataLakeStoreAccountParameters[property_name]

function check_required(o::CreateDataLakeStoreAccountParameters)
    (getproperty(o, Symbol("location")) === nothing) && (return false)
    true
end

function validate_property(::Type{ CreateDataLakeStoreAccountParameters }, name::Symbol, val)
end


text:
struct ConsistencyTest{E<:CalibrationErrorEstimator,P,T,V} <: HypothesisTests.HypothesisTest
    """Calibration estimator."""
    estimator::E
    """Predictions."""
    predictions::P
    """Targets."""
    targets::T
    """Calibration estimate."""
    estimate::V
end

function ConsistencyTest(
    estimator::CalibrationErrorEstimator,
    predictions::AbstractVector,
    targets::AbstractVector,
)
    estimate = estimator(predictions, targets)
    return ConsistencyTest(estimator, predictions, targets, estimate)
end

# HypothesisTests interface

HypothesisTests.default_tail(::ConsistencyTest) = :right

HypothesisTests.testname(::ConsistencyTest) = "Consistency resampling test"

function HypothesisTests.pvalue(test::ConsistencyTest; kwargs...)
    return pvalue(Random.GLOBAL_RNG, test; kwargs...)
end
function HypothesisTests.pvalue(
    rng::Random.AbstractRNG, test::ConsistencyTest; bootstrap_iters::Int=1_000
)
    bootstrap_iters > 0 || error("number of bootstrap samples must be positive")

    predictions = test.predictions
    sampledpredictions = similar(predictions)
    sampledtargets = similar(test.targets)
    samples = StructArrays.StructArray((sampledpredictions, sampledtargets))

    estimate = test.estimate
    estimator = test.estimator
    sampler = Random.Sampler(rng, ConsistencyResampling.Consistent(predictions))
    n = 0
    for _ in 1:bootstrap_iters
        # perform consistency resampling
        Random.rand!(rng, samples, sampler)

        # evaluate the calibration error
        sampledestimate = estimator(sampledpredictions, sampledtargets)

        # check if the estimate for the resampled data is ≥ the original estimate
        if sampledestimate ≥ estimate
            n += 1
        end
    end

    return n / bootstrap_iters
end

# parameter of interest: name, value under H0, point estimate
function HypothesisTests.population_param_of_interest(test::ConsistencyTest)
    return nameof(typeof(test.estimator)), zero(test.estimate), test.estimate
end

HypothesisTests.show_params(io::IO, test::ConsistencyTest, ident="") = nothing


text:
# direct implementations
e(x::AbstractMatrix) = x
tr(x::AbstractMatrix) = collect(x')
fx(x::AbstractMatrix) = reverse(x; dims=1)
fy(x::AbstractMatrix) = reverse(x; dims=2)
r180(x::AbstractMatrix) = reverse(x)

# compound implementations
r90(x::AbstractMatrix) = fy(x') # fy(x') is slightly more efficient than fy(tr(x))
r270(x::AbstractMatrix) = fx(x') # fx(x') is slightly more efficient than fx(tr(x))
trr(x::AbstractMatrix) = r180(x') # r180(x') is slightly more efficient than r180(tr(x))

# group operation for elements of the symmetry group
Base.:∘(::Val{e}, ::Val{e}) = e
Base.:∘(::Val{e}, ::Val{fx}) = fx
Base.:∘(::Val{e}, ::Val{fy}) = fy
Base.:∘(::Val{e}, ::Val{tr}) = tr
Base.:∘(::Val{e}, ::Val{trr}) = trr
Base.:∘(::Val{e}, ::Val{r90}) = r90
Base.:∘(::Val{e}, ::Val{r180}) = r180
Base.:∘(::Val{e}, ::Val{r270}) = r270

Base.:∘(::Val{fx}, ::Val{e}) = fx
Base.:∘(::Val{fx}, ::Val{fx}) = e
Base.:∘(::Val{fx}, ::Val{fy}) = r180
Base.:∘(::Val{fx}, ::Val{tr}) = r270
Base.:∘(::Val{fx}, ::Val{trr}) = r90
Base.:∘(::Val{fx}, ::Val{r90}) = trr
Base.:∘(::Val{fx}, ::Val{r180}) = fy
Base.:∘(::Val{fx}, ::Val{r270}) = tr

Base.:∘(::Val{fy}, ::Val{e}) = fy
Base.:∘(::Val{fy}, ::Val{fx}) = r180
Base.:∘(::Val{fy}, ::Val{fy}) = e
Base.:∘(::Val{fy}, ::Val{tr}) = r90
Base.:∘(::Val{fy}, ::Val{trr}) = r270
Base.:∘(::Val{fy}, ::Val{r90}) = tr
Base.:∘(::Val{fy}, ::Val{r180}) = fx
Base.:∘(::Val{fy}, ::Val{r270}) = trr

Base.:∘(::Val{tr}, ::Val{e}) = tr
Base.:∘(::Val{tr}, ::Val{fx}) = r90
Base.:∘(::Val{tr}, ::Val{fy}) = r270
Base.:∘(::Val{tr}, ::Val{tr}) = e
Base.:∘(::Val{tr}, ::Val{trr}) = r180
Base.:∘(::Val{tr}, ::Val{r90}) = fx
Base.:∘(::Val{tr}, ::Val{r180}) = trr
Base.:∘(::Val{tr}, ::Val{r270}) = fy

Base.:∘(::Val{trr}, ::Val{e}) = trr
Base.:∘(::Val{trr}, ::Val{fx}) = r270
Base.:∘(::Val{trr}, ::Val{fy}) = r90
Base.:∘(::Val{trr}, ::Val{tr}) = r180
Base.:∘(::Val{trr}, ::Val{trr}) = e
Base.:∘(::Val{trr}, ::Val{r90}) = fy
Base.:∘(::Val{trr}, ::Val{r180}) = tr
Base.:∘(::Val{trr}, ::Val{r270}) = fx

Base.:∘(::Val{r90}, ::Val{e}) = r90
Base.:∘(::Val{r90}, ::Val{fx}) = tr
Base.:∘(::Val{r90}, ::Val{fy}) = trr
Base.:∘(::Val{r90}, ::Val{tr}) = fy
Base.:∘(::Val{r90}, ::Val{trr}) = fx
Base.:∘(::Val{r90}, ::Val{r90}) = r180
Base.:∘(::Val{r90}, ::Val{r180}) = r270
Base.:∘(::Val{r90}, ::Val{r270}) = e

Base.:∘(::Val{r180}, ::Val{e}) = r180
Base.:∘(::Val{r180}, ::Val{fx}) = fy
Base.:∘(::Val{r180}, ::Val{fy}) = fx
Base.:∘(::Val{r180}, ::Val{tr}) = trr
Base.:∘(::Val{r180}, ::Val{trr}) = tr
Base.:∘(::Val{r180}, ::Val{r90}) = r270
Base.:∘(::Val{r180}, ::Val{r180}) = e
Base.:∘(::Val{r180}, ::Val{r270}) = r90

Base.:∘(::Val{r270}, ::Val{e}) = r270
Base.:∘(::Val{r270}, ::Val{fx}) = trr
Base.:∘(::Val{r270}, ::Val{fy}) = tr
Base.:∘(::Val{r270}, ::Val{tr}) = fx
Base.:∘(::Val{r270}, ::Val{trr}) = fy
Base.:∘(::Val{r270}, ::Val{r90}) = e
Base.:∘(::Val{r270}, ::Val{r180}) = r90
Base.:∘(::Val{r270}, ::Val{r270}) = r180

Base.:∘(f::Function, g::Function) = ∘(Val(f), Val(g))


text:
using RefractionCorrection
using Test
# using LazIO

@testset "RefractionCorrection.jl" begin

    @testset "Correct file" begin
        #MAIN#
        # ds = LazIO.open(filename)
        #
        # #subset
        # pp = points[1500000:1520000]
        #
        # #Array with coordinates
        # coords = coordinates(pp,h)
        #
        # #Separate laser pulses
        # laser_pulses = laserpulses(coords)
        #
        # #Return corrected points
        # corrected_laser_pulses, indexes_corr, indexes_un_corr  = corrected_points(laser_pulses,coords)
        #
        # #write pointcloud
        # write_pointcloud(ds,h,pp,indexes_corr,indexes_un_corr)
        @test 1 == 1
    end

    @testset "Correct" begin
        minz, maxz = -6, -5
        p1 = 0, 0, maxz
        p2 = 0, 0, minz
        _, _, z = correct(p1, p2)
        @test minz < z < maxz
    end

    @testset "Polar" begin
        p = 1., 1., 1.
        x = cart2pol(p...)
        pp = pol2cart(x...)
        @test all(p .≈ pp)
    end
end


text:
using Logging
try
    using PackageCompiler
catch
    @warn "PackageCompiler is not installed, attempting to install"
    using Pkg
    Pkg.add("PackageCompiler")
    using PackageCompiler
end

create_app("..", "fierapp", filter_stdlibs = true)

max_stars_repo_name: SERVIR/fier-cli

max_issues_repo_name: SERVIR/fier-cli

max_forks_repo_name: SERVIR/fier-cli


text:
"""
    mutable struct EdwardsInteraction <: AbstractInteraction

An interaction representing an Edward's excluded volume interaction between like-species.

The interaction energy is of the form:
``H = (1/2) ∑_i u_i ∫dr ϕ_i^2(r)``
"""
mutable struct EdwardsInteraction <: AbstractInteraction
    mids   :: Set{Int}
    uvals  :: Dict{Int,Float64}
    system :: Option{FieldSystem}
    EdwardsInteraction() = new(Set(), Dict(), nothing)
end

#==============================================================================#

Base.show(io::IO, itx::EdwardsInteraction) = @printf(io, "EdwardsInteraction(%d monomers)", length(itx.mids))

function setup!(itx::EdwardsInteraction, sys::FieldSystem)
    itx.system = sys
    return nothing
end

function set_interaction!(itx::EdwardsInteraction, mid::Integer, u::Real)
    if mid in itx.mids
        @warn "Replacing interaction for monomer with id = $(mid)."
    end
    push!(itx.mids, mid)
    itx.uvals[mid] = u
    return nothing
end

function energy(itx::EdwardsInteraction)
    @assert !isnothing(itx.system)
    sys = itx.system

    energy = 0.0
    for mid in itx.mids
        if hasmonomer(sys, mid)
            u = itx.uvals[mid]
            v = sys.monomers[mid].vol
            rho = sys.density[mid]

            @simd for i in eachindex(rho)
                @inbounds energy += (u/2)*(rho[i]/v)^2
            end
        end
    end

    return energy/ngrid(sys)
end

function energy_bulk(itx::EdwardsInteraction)
    @assert !isnothing(itx.system)
    sys = itx.system

    energy = 0.0
    for mid in itx.mids
        if hasmonomer(sys, mid)
            u = itx.uvals[mid]
            v = sys.monomers[mid].vol
            rho = sys.monomer_fracs[mid]

            energy += (u/2)*(rho/v)^2
        end
    end

    return energy
end

function potential!(itx::EdwardsInteraction, mid::Integer, pot::FieldGrid)
    @assert !isnothing(itx.system)
    sys = itx.system

    if mid in itx.mids
        u = itx.uvals[mid]
        v = sys.monomers[mid].vol
        rho = sys.density[mid]

        @simd for i in eachindex(pot)
            @inbounds pot[i] += u*(rho[i]/v)
        end
    end

    return nothing
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: mediastore_data

using AWS.Compat
using AWS.UUIDs
"""
    DeleteObject()

Deletes an object at the specified path.

# Required Parameters
- `Path`: The path (including the file name) where the object is stored in the container. Format: <folder name>/<folder name>/<file name>

"""

delete_object(Path; aws_config::AWSConfig=global_aws_config()) = mediastore_data("DELETE", "/$(Path)"; aws_config=aws_config)
delete_object(Path, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore_data("DELETE", "/$(Path)", args; aws_config=aws_config)

"""
    DescribeObject()

Gets the headers for an object at the specified path.

# Required Parameters
- `Path`: The path (including the file name) where the object is stored in the container. Format: <folder name>/<folder name>/<file name>

"""

describe_object(Path; aws_config::AWSConfig=global_aws_config()) = mediastore_data("HEAD", "/$(Path)"; aws_config=aws_config)
describe_object(Path, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore_data("HEAD", "/$(Path)", args; aws_config=aws_config)

"""
    GetObject()

Downloads the object at the specified path. If the object’s upload availability is set to streaming, AWS Elemental MediaStore downloads the object even if it’s still uploading the object.

# Required Parameters
- `Path`: The path (including the file name) where the object is stored in the container. Format: <folder name>/<folder name>/<file name> For example, to upload the file mlaw.avi to the folder path premiumcanada in the container movies, enter the path premium/canada/mlaw.avi. Do not include the container name in this path. If the path includes any folders that don't exist yet, the service creates them. For example, suppose you have an existing premium/usa subfolder. If you specify premium/canada, the service creates a canada subfolder in the premium folder. You then have two subfolders, usa and canada, in the premium folder.  There is no correlation between the path to the source and the path (folders) in the container in AWS Elemental MediaStore. For more information about folders and how they exist in a container, see the AWS Elemental MediaStore User Guide. The file name is the name that is assigned to the file that you upload. The file can have the same name inside and outside of AWS Elemental MediaStore, or it can have the same name. The file name can include or omit an extension.

# Optional Parameters
- `Range`: The range bytes of an object to retrieve. For more information about the Range header, see http://www.w3.org/Protocols/rfc2616/rfc2616-sec14.html#sec14.35. AWS Elemental MediaStore ignores this header for partially uploaded objects that have streaming upload availability.
"""

get_object(Path; aws_config::AWSConfig=global_aws_config()) = mediastore_data("GET", "/$(Path)"; aws_config=aws_config)
get_object(Path, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore_data("GET", "/$(Path)", args; aws_config=aws_config)

"""
    ListItems()

Provides a list of metadata entries about folders and objects in the specified folder.

# Optional Parameters
- `MaxResults`: The maximum number of results to return per API request. For example, you submit a ListItems request with MaxResults set at 500. Although 2,000 items match your request, the service returns no more than the first 500 items. (The service also returns a NextToken value that you can use to fetch the next batch of results.) The service might return fewer results than the MaxResults value. If MaxResults is not included in the request, the service defaults to pagination with a maximum of 1,000 results per page.
- `NextToken`: The token that identifies which batch of results that you want to see. For example, you submit a ListItems request with MaxResults set at 500. The service returns the first batch of results (up to 500) and a NextToken value. To see the next batch of results, you can submit the ListItems request a second time and specify the NextToken value. Tokens expire after 15 minutes.
- `Path`: The path in the container from which to retrieve items. Format: <folder name>/<folder name>/<file name>
"""

list_items(; aws_config::AWSConfig=global_aws_config()) = mediastore_data("GET", "/"; aws_config=aws_config)
list_items(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = mediastore_data("GET", "/", args; aws_config=aws_config)

"""
    PutObject()

Uploads an object to the specified path. Object sizes are limited to 25 MB for standard upload availability and 10 MB for streaming upload availability.

# Required Parameters
- `Body`: The bytes to be stored.
- `Path`: The path (including the file name) where the object is stored in the container. Format: <folder name>/<folder name>/<file name> For example, to upload the file mlaw.avi to the folder path premiumcanada in the container movies, enter the path premium/canada/mlaw.avi. Do not include the container name in this path. If the path includes any folders that don't exist yet, the service creates them. For example, suppose you have an existing premium/usa subfolder. If you specify premium/canada, the service creates a canada subfolder in the premium folder. You then have two subfolders, usa and canada, in the premium folder.  There is no correlation between the path to the source and the path (folders) in the container in AWS Elemental MediaStore. For more information about folders and how they exist in a container, see the AWS Elemental MediaStore User Guide. The file name is the name that is assigned to the file that you upload. The file can have the same name inside and outside of AWS Elemental MediaStore, or it can have the same name. The file name can include or omit an extension.

# Optional Parameters
- `Cache-Control`: An optional CacheControl header that allows the caller to control the object's cache behavior. Headers can be passed in as specified in the HTTP at https://www.w3.org/Protocols/rfc2616/rfc2616-sec14.html#sec14.9. Headers with a custom user-defined value are also accepted.
- `Content-Type`: The content type of the object.
- `x-amz-storage-class`: Indicates the storage class of a Put request. Defaults to high-performance temporal storage class, and objects are persisted into durable storage shortly after being received.
- `x-amz-upload-availability`: Indicates the availability of an object while it is still uploading. If the value is set to streaming, the object is available for downloading after some initial buffering but before the object is uploaded completely. If the value is set to standard, the object is available for downloading only when it is uploaded completely. The default value for this header is standard. To use this header, you must also set the HTTP Transfer-Encoding header to chunked.
"""

put_object(Body, Path; aws_config::AWSConfig=global_aws_config()) = mediastore_data("PUT", "/$(Path)", Dict{String, Any}("Body"=>Body); aws_config=aws_config)
put_object(Body, Path, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediastore_data("PUT", "/$(Path)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Body"=>Body), args)); aws_config=aws_config)


text:
module RecipesBase

export
    @recipe,
    @series,
    @userplot,
    @shorthands,
    RecipeData,
    AbstractBackend,
    AbstractPlot,
    AbstractLayout

# Common abstract types for the Plots ecosystem
abstract type AbstractBackend end
abstract type AbstractPlot{T<:AbstractBackend} end
abstract type AbstractLayout end

# a placeholder to establish the name so that other packages (Plots.jl for example)
# can add their own definition of RecipesBase.plot since RecipesBase is the common
# dependency of the Plots ecosystem
function plot end
function plot! end

# a placeholder to establish the name so that other packages (Plots.jl for example)
# can add their own definition of RecipesBase.animate since RecipesBase is the common
# dependency of the Plots ecosystem. Plots.jl will handle the basic cases, while
# other packages can now extend for their types
function animate end

# a placeholder to establish the name so that other packages (Plots.jl for example)
# can add their own definition of RecipesBase.is_key_supported(k::Symbol)
function is_key_supported end

# a placeholder to establish the name so that other packages (Plots.jl for example)
# can add their own definition of RecipesBase.group_as_matrix(t)
group_as_matrix(t) = false

# This holds the recipe definitions to be dispatched on
# the function takes in an attribute dict `d` and a list of args.
# This default definition specifies the "no-arg" case.
apply_recipe(plotattributes::AbstractDict{Symbol,Any}) = ()

# Is a key explicitly provided by the user?
# Should be overridden for subtypes representing plot attributes.
is_explicit(d::AbstractDict{Symbol,Any}, k) = haskey(d, k)

const _debug_recipes = Bool[false]
function debug(v::Bool = true)
    _debug_recipes[1] = v
end

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

# this holds the data and attributes of one series, and is returned from apply_recipe
struct RecipeData
    plotattributes::AbstractDict{Symbol,Any}
    args::Tuple
end

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

@inline to_symbol(s::Symbol) = s
@inline to_symbol(qn::QuoteNode) = qn.value

@inline wrap_tuple(tup::Tuple) = tup
@inline wrap_tuple(v) = (v,)

# check for flags as part of the `-->` expression
function _is_arrow_tuple(expr::Expr)
    expr.head == :tuple && !isempty(expr.args) &&
        isa(expr.args[1], Expr) &&
        expr.args[1].head == :(-->)
end

function _equals_symbol(arg::Symbol, sym::Symbol)
    arg == sym
end
function _equals_symbol(arg::Expr, sym::Symbol) #not sure this method is necessary anymore on 0.7
    arg.head == :quote && arg.args[1] == sym
end
function _equals_symbol(arg::QuoteNode, sym::Symbol)
    arg.value == sym
end
_equals_symbol(x, sym::Symbol) = false

# build an apply_recipe function header from the recipe function header
function get_function_def(func_signature::Expr, args::Vector)
    front = func_signature.args[1]
    if func_signature.head == :where
        Expr(:where, get_function_def(front, args), esc.(func_signature.args[2:end])...)
    elseif func_signature.head == :call
        func = Expr(:call, :(RecipesBase.apply_recipe), esc.([:(plotattributes::AbstractDict{Symbol, Any}); args])...)
        if isa(front, Expr) && front.head == :curly
            Expr(:where, func, esc.(front.args[2:end])...)
        else
            func
        end
    else
        error("Expected `func_signature = ...` with func_signature as a call or where Expr... got: $func_signature")
    end
end

function create_kw_body(func_signature::Expr)
    # get the arg list, stripping out any keyword parameters into a
    # bunch of get!(kw, key, value) lines
    func_signature.head == :where && return create_kw_body(func_signature.args[1])
    args = func_signature.args[2:end]
    kw_body = Expr(:block)
    cleanup_body = Expr(:block)
    if isa(args[1], Expr) && args[1].head == :parameters
        for kwpair in args[1].args
            k, v = kwpair.args
            if isa(k, Expr) && k.head == :(::)
                k = k.args[1]
                @warn("Type annotations on keyword arguments not currently supported in recipes. Type information has been discarded")
            end
            push!(kw_body.args, :($k = get!(plotattributes, $(QuoteNode(k)), $v)))
            push!(cleanup_body.args, :(RecipesBase.is_key_supported($(QuoteNode(k))) || delete!(plotattributes, $(QuoteNode(k)))))
        end
        args = args[2:end]
    end
    args, kw_body, cleanup_body
end

# process the body of the recipe recursively.
# when we see the series macro, we split that block off:
    # let
    #   d2 = copy(d)
    #   <process_recipe_body on d2>
    #   RecipeData(d2, block_return)
    # end
# and we push this block onto the series_blocks list.
# then at the end we push the main body onto the series list
function process_recipe_body!(expr::Expr)
    for (i,e) in enumerate(expr.args)
        if isa(e,Expr)

            # process trailing flags, like:
            #   a --> b, :quiet, :force
            quiet, require, force = false, false, false
            if _is_arrow_tuple(e)
                for flag in e.args
                    if _equals_symbol(flag, :quiet)
                        quiet = true
                    elseif _equals_symbol(flag, :require)
                        require = true
                    elseif _equals_symbol(flag, :force)
                        force = true
                    end
                end
                e = e.args[1]
            end

            # the unused operator `:=` will mean force: `x := 5` is equivalent to `x --> 5, force`
            # note: this means "x is defined as 5"
            if e.head == :(:=)
                force = true
                e.head = :(-->)
            end

            # we are going to recursively swap out `a --> b, flags...` commands
            # note: this means "x may become 5"
            if e.head == :(-->)
                k, v = e.args
                if isa(k, Symbol)
                    k = QuoteNode(k)
                end

                set_expr = if force
                    # forced override user settings
                    :(plotattributes[$k] = $v)
                else
                    # if the user has set this keyword, use theirs
                    :(RecipesBase.is_explicit(plotattributes, $k) || (plotattributes[$k] = $v))
                end

                expr.args[i] = if quiet
                    # quietly ignore keywords which are not supported
                    :(RecipesBase.is_key_supported($k) ? $set_expr : nothing)
                elseif require
                    # error when not supported by the backend
                    :(RecipesBase.is_key_supported($k) ? $set_expr : error("In recipe: required keyword ", $k, " is not supported by backend $(backend_name())"))
                else
                    set_expr
                end

            elseif e.head == :return
                # To allow `return` in recipes just extract the returned arguments.
                expr.args[i] = first(e.args)

            elseif e.head != :call
                # we want to recursively replace the arrows, but not inside function calls
                # as this might include things like Dict(1=>2)
                process_recipe_body!(e)
            end
        end
    end
end

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

"""
This handy macro will process a function definition, replace `-->` commands, and
then add a new version of `RecipesBase.apply_recipe` for dispatching on the arguments.

This functionality is primarily geared to turning user types and settings into the
data and attributes that describe a Plots.jl visualization.

Set attributes using the `-->` command, and return a comma separated list of arguments that
should replace the current arguments.

An example:

```
using RecipesBase

# Our custom type that we want to display
type T end

@recipe function plot{N<:Integer}(t::T, n::N = 1; customcolor = :green)
    markershape --> :auto, :require
    markercolor --> customcolor, :force
    xrotation --> 5
    zrotation --> 6, :quiet
    rand(10,n)
end

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

# Plots will be the ultimate consumer of our recipe in this example
using Plots; gr()

# This call will implicitly call `RecipesBase.apply_recipe` as part of the Plots
# processing pipeline (see the Pipeline section of the Plots documentation).
# It will plot 5 line plots, all with black circles for markers.
# The markershape argument must be supported, and the zrotation argument's warning
# will be suppressed.  The user can override all arguments except markercolor.
plot(T(), 5; customcolor = :black, shape=:c)
```

In this example, we see lots of the machinery in action.  We create a new type `T` which
we will use for dispatch, and an optional argument `n`, which will be used to determine the
number of series to display.  User-defined keyword arguments are passed through, and the
`-->` command can be trailed by flags:

- quiet:   Suppress unsupported keyword warnings
- require: Error if keyword is unsupported
- force:   Don't allow user override for this keyword
"""
macro recipe(funcexpr::Expr)
    func_signature, func_body = funcexpr.args

    if !(funcexpr.head in (:(=), :function))
        error("Must wrap a valid function call!")
    end
    if !(isa(func_signature, Expr) && func_signature.head in (:call, :where))
        error("Expected `func_signature = ...` with func_signature as a call or where Expr... got: $func_signature")
    end
    if length(func_signature.args) < 2
        error("Missing function arguments... need something to dispatch on!")
    end

    args, kw_body, cleanup_body = create_kw_body(func_signature)
    func = get_function_def(func_signature, args)

    # this is where the receipe func_body is processed
    # replace all the key => value lines with argument setting logic
    # and break up by series.
    process_recipe_body!(func_body)

    # now build a function definition for apply_recipe, wrapping the return value in a tuple if needed.
    # we are creating a vector of RecipeData objects, one per series.
    funcdef = Expr(:function, func, esc(quote
        @nospecialize
        if RecipesBase._debug_recipes[1]
            println("apply_recipe args: ", $args)
        end
        $kw_body
        $cleanup_body
        series_list = RecipesBase.RecipeData[]
        func_return = $func_body
        if func_return != nothing
            push!(series_list, RecipesBase.RecipeData(plotattributes, RecipesBase.wrap_tuple(func_return)))
        end
        series_list
    end))
    funcdef
end

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

"""
Meant to be used inside a recipe to add additional RecipeData objects to the list:

```
@recipe function f(::T)
    # everything get this setting
    linecolor --> :red

    @series begin
        # this setting is only for this series
        fillcolor := :green

        # return the args, just like in recipes
        rand(10)
    end

    # this is the main series... though it can be skipped by returning nothing.
    # note: a @series block returns nothing
    rand(100)
end
```
"""
macro series(expr::Expr)
    esc(quote
        let plotattributes = copy(plotattributes)
            args = $expr
            push!(series_list, RecipesBase.RecipeData(plotattributes, RecipesBase.wrap_tuple(args)))
            nothing
        end
    end)
end

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

"""
You can easily define your own plotting recipes with convenience methods:
```
@userplot GroupHist

@recipe function f(gh::GroupHist)
    # set some attributes, add some series, using gh.args as input
end
# now you can plot like:
grouphist(rand(1000,4))
```
"""
macro userplot(expr)
    _userplot(expr)
end

function _userplot(expr::Expr)
    if expr.head != :struct
        error("Must call userplot on a [mutable] struct expression.  Got: $expr")
    end

    typename = gettypename(expr.args[2])
    funcname = Symbol(lowercase(string(typename)))
    funcname2 = Symbol(funcname, "!")

    # return a code block with the type definition and convenience plotting methods
    esc(quote
        $expr
        export $funcname, $funcname2
        Core.@__doc__ $funcname(args...; kw...) = RecipesBase.plot($typename(args); kw...)
        Core.@__doc__ $funcname2(args...; kw...) = RecipesBase.plot!($typename(args); kw...)
        Core.@__doc__ $funcname2(plt::RecipesBase.AbstractPlot, args...; kw...) = RecipesBase.plot!(plt, $typename(args); kw...)
    end)
end

function _userplot(sym::Symbol)
    _userplot(:(mutable struct $sym
            args
    end))
end

gettypename(sym::Symbol) = sym

function gettypename(expr::Expr)
    expr.head == :curly || @error "Unexpected struct name: $expr"
    expr.args[1]
end

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

"""
    @shorthands(funcname::Symbol)

Defines and exports shorthand plotting method definitions (`\$funcname` and `\$funcname!`).
Pass the series type (as a symbol) to the macro.

## Examples

```julia
# define some series type
@recipe function f(::Type{Val{:myseriestype}}, x, y)
    # some implementation here
end
# docstrings are forwarded
\"\"\"
    myseriestype(x, y)
Plot my series type!
\"\"\"
@shorthands myseriestype
```
"""
macro shorthands(funcname::Symbol)
    funcname2 = Symbol(funcname, "!")
    esc(quote
        export $funcname, $funcname2
        Core.@__doc__ $funcname(args...; kw...) = RecipesBase.plot(args...; kw..., seriestype = $(Meta.quot(funcname)))
        Core.@__doc__ $funcname2(args...; kw...) = RecipesBase.plot!(args...; kw..., seriestype = $(Meta.quot(funcname)))
    end)
end

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

# allow usage of type recipes without depending on StatsPlots

"""
`recipetype(s, args...)`

Use this function to refer to type recipes by their symbol, without taking a dependency.

# Example

```julia
import RecipesBase: recipetype
recipetype(:groupedbar, 1:10, rand(10, 2))
```

instead of

```julia
import StatsPlots: GroupedBar
GroupedBar((1:10, rand(10, 2)))
```
"""
recipetype(s, args...) = recipetype(Val(s), args...)

function recipetype(s::Val{T}, args...) where T
    error("No type recipe defined for type $T. You may need to load StatsPlots")
end

include("precompile.jl")

end # module