text:
"""
Base abstract class for implementing regions in `N`-dimensional space with
dimensions of type `T`.
"""
abstract type Region{T,N} end

dimtype(::Type{<:Region{T}}) where T = T
dimtype(r::Region) = dimtype(typeof(r))
Base.ndims(::Type{<:Region{<:Any, N}}) where N = N
Base.ndims(r::Region) = ndims(typeof(r))

"""
`N`-dimensional point.
"""
struct Point{T, N} <: Region{T, N}
    coord::NTuple{N, T}
end

"""
    empty(::Type{T}) where T<:Region

Generate empty (uninitialized) region of type `T`.
"""
empty(::Type{Point{T,N}}) where {T,N} =
    Point{T,N}(ntuple(_ -> convert(T, NaN), N))
empty(::Type{Point{T,N}}) where {T<:Integer,N} =
    Point{T,N}(ntuple(_ -> zero(T), N))

"""
    area(a::Region)

`N`-dimensional "area" (volume etc) of `a`.
"""
area(a::Point{T}) where T = zero(typeof(zero(T)*zero(T)))

"""
    perimeter(a::Region)

The sum of the `a` sides.
"""
perimeter(a::Point{T}) where T = zero(T)

"""
    isvalid(a::Region)

Check that the parameters of `a` are valid and it defines a proper region.
"""
isvalid(a::Point) = all(!isnan, a.coord)

"""
Rectangular region constrained by `low[i]`...`high[i]` in each of `N` dimensions.
"""
struct Rect{T, N} <: Region{T, N}
    low::NTuple{N, T}
    high::NTuple{N, T}
end

Rect(low::NTuple{N, T1}, high::NTuple{N, T2}) where {N, T1<:Number, T2<:Number} =
    Rect{Base.promote_type(T1, T2), N}(low, high)

Rect{T,N}(pt::Point{T,N}) where {T,N} = Rect{T,N}(pt.coord, pt.coord)
Rect(pt::Point) = Rect(pt.coord, pt.coord)

# "empty" (uninitialized) rectangles
empty(::Type{Rect{T,N}}) where {T,N} =
    Rect{T,N}(ntuple(_ -> convert(T, NaN), N), ntuple(_ -> convert(T, NaN), N))
empty(::Type{Rect{T,N}}) where {T<:Integer,N} =
    Rect{T,N}(ntuple(_ -> typemax(T), N), ntuple(_ -> typemax(T), N))

"""
    combine(a::Region, b::Region)

MBR that contains both `a` and `b` regions.
"""
combine(a::Rect{T,N}, b::Rect{T,N}) where {T,N} =
    @inbounds Rect{T,N}(ntuple(i -> min(a.low[i], b.low[i]), N),
                        ntuple(i -> max(a.high[i], b.high[i]), N))

# check that coords are not NaN and side lengths not negative.
@generated isvalid(a::Rect{T,N}) where {T,N} =
    quote
        Base.Cartesian.@nall $N i -> !isnan(a.low[i]) && !isnan(a.high[i]) && a.low[i] <= a.high[i]
    end

"""
    intersect(a::Region, b::Region)

MBR that is an intersection of `a` and `b` regions.
"""
function intersect(a::Rect{T,N}, b::Rect{T,N}) where {T,N}
    res = @inbounds Rect{T,N}(ntuple(i -> max(a.low[i], b.low[i]), N),
                              ntuple(i -> min(a.high[i], b.high[i]), N))
    return isvalid(res) ? res : empty(typeof(a))
end

@generated area(a::Rect{T,N}) where {T,N} =
    quote
        @inbounds res = a.high[1] - a.low[1]
        @inbounds Base.Cartesian.@nexprs $(N-1) i -> res *= (a.high[i+1] - a.low[i+1])
        return res
    end

@generated perimeter(a::Rect{T,N}) where {T,N} =
    quote
        @inbounds res = a.high[1] - a.low[1]
        @inbounds Base.Cartesian.@nexprs $(N-1) i -> res += (a.high[i+1] - a.low[i+1])
        return res*ndims(a)
    end

"""
    overlap_area(a::Region, b::Region)

The area of MBR for `a` and `b` intersection.
"""
@generated overlap_area(a::Rect{T,N}, b::Rect{T,N}) where {T,N} =
    quote
        @inbounds res = min(a.high[1], b.high[1]) - max(a.low[1], b.low[1])
        res <= 0.0 && return zero(typeof(res))
        @inbounds Base.Cartesian.@nexprs $(N-1) i -> begin
            width_i = min(a.high[i+1], b.high[i+1]) - max(a.low[i+1], b.low[i+1])
            width_i <= zero(typeof(width_i)) && return zero(typeof(res))
            res *= width_i
        end
        return res
    end

"""
    combined_area(a::Region, b::Region)

The area of MBR for the union of `a` and `b`.
"""
@generated combined_area(a::Rect{T,N}, b::Rect{T,N}) where {T,N} =
    quote
        @inbounds res = max(a.high[1], b.high[1]) - min(a.low[1], b.low[1])
        @inbounds Base.Cartesian.@nexprs $(N-1) i -> res *= max(a.high[i+1], b.high[i+1]) - min(a.low[i+1], b.low[i+1])
        return res
    end

"""
    enlargement(a::Region, b::Region)

How much `a` grows when combined with `b`.

The difference between the MBR area of the `a` and `b` union and the MBR area of `a`.
"""
enlargement(a::Rect{T,N}, b::Rect{T,N}) where {T,N} = combined_area(a, b) - area(a)

"""
    center(a::Region)

The center point of `a`.
"""
center(a::Rect{T,N}) where {T,N} =
    Point(ntuple(i -> 0.5*(a.low[i] + a.high[i]), N))

@generated sqrdistance(a::Point{T,N}, b::Point{T,N}) where {T,N} =
    quote
        @inbounds res = abs2(a.coord[1] - b.coord[1])
        @inbounds Base.Cartesian.@nexprs $(N-1) i -> res += abs2(a.coord[i+1] - b.coord[i+1])
        return res
    end

"""
    intersects(a::Region, b::Region)

Check whether `a` intersects with `b`.
"""
@generated intersects(a::Rect{T,N}, b::Rect{T,N}) where {T,N} =
    quote
        @inbounds Base.Cartesian.@nall $N i -> (a.low[i] <= b.low[i] <= a.high[i]) || (b.low[i] <= a.low[i] <= b.high[i])
    end

"""
    contains(a::Region, b::Region)

Check whether `a` contains `b`.
"""
@generated contains(a::Rect{T,N}, b::Rect{T,N}) where {T,N} =
    quote
        @inbounds Base.Cartesian.@nall $N i -> (a.low[i] <= b.low[i]) && (a.high[i] >= b.high[i])
    end

@generated contains(a::Rect{T,N}, b::Point{T,N}) where {T,N} =
    quote
        @inbounds Base.Cartesian.@nall $N i -> (a.low[i] <= b.coord[i] <= a.high[i])
    end

"""
    in(a::Region, b::Region)

Check whether `a` is contained inside `b`.
"""
Base.in(a::Rect, b::Rect) = contains(b, a)
Base.in(a::Point, b::Rect) = contains(b, a)

# point is equal to the MBR if its low and high corners coincides
Base.:(==)(a::Point, b::Rect) = a.coord == b.low == b.high
Base.:(==)(a::Rect, b::Point) = b == a

"""
    touches(a::Rect, b::Rect)

Check whether `a` and `b` rectangles touches from the inside
(i.e. any `low` side touches `low` or `high` touches `high`).
"""
@generated touches(a::Rect{T,N}, b::Rect{T,N}) where {T,N} =
    quote
        @inbounds Base.Cartesian.@nany $N i -> (a.low[i] == b.low[i]) || (a.high[i] == b.high[i])
    end


text:
module adapters_postgresql_functions_test

using Test # @test
using Octo.Adapters.PostgreSQL # Repo Schema SELECT AVG FROM

Repo.debug_sql()

Repo.connect(
    adapter = Octo.Adapters.PostgreSQL,
    user = "postgres",
)

struct User
end
Schema.model(User, table_name="users")

u = from(User)

Repo.execute([DROP TABLE IF EXISTS User])
Repo.execute(Raw("""CREATE TABLE IF NOT EXISTS users (
                     ID SERIAL,
                     name VARCHAR(255),
                     salary FLOAT(8),
                     PRIMARY KEY (ID) )
                 """))

Repo.insert!(User, [
    (Name="Jeremy",  Salary=10000.50),
    (Name="Cloris",  Salary=20000.50),
    (Name="John",    Salary=30000.50),
    (Name="Hyunden", Salary=40000.50),
    (Name="Justin",  Salary=50000.50),
    (Name="Tom",     Salary=60000.50),
])

df = Repo.query([SELECT (AVG(u.salary), MAX(u.salary), MIN(u.salary)) FROM u])
@test Pretty.table(df) == """
|       avg |       max |       min |
| --------- | --------- | --------- |
|   35000.5 |   60000.5 |   10000.5 |
1 row."""

df = Repo.query([SELECT (u.name, u.salary) FROM u ORDER BY u.salary DESC])
@test Pretty.table(df) == """
| name      |    salary |
| --------- | --------- |
| Tom       |   60000.5 |
| Justin    |   50000.5 |
| Hyunden   |   40000.5 |
| John      |   30000.5 |
| Cloris    |   20000.5 |
| Jeremy    |   10000.5 |
6 rows."""

end # module adapters_postgresql_functions_test

ext: jl


text:
module DustExtinction

using Unitful
using UnitfulAstro
using DataDeps
using Parameters

export redden,
       deredden,
       CCM89,
       CAL00,
       OD94,
       GCC09,
       VCG04,
       SFD98Map,
       ebv_galactic,
       FM90

"""
    DustExtinction.ExtinctionLaw

The abstract super-type for dust extinction laws. See the extended help (`??DustExtinction.ExtinctionLaw` from the REPL) for more information about the interface.

# Extended Help

## Interface

Here's how to make a new extinction law, called `MyLaw`
1. Create your subtype
       struct MyLaw <: DustExtinction.ExtinctionLaw end
2. Define the limits. This will default to `(0, Inf)`. Currently this doesn't do anything internally.
       DustExtinction.bounds(::Type{<:MyLaw}) = (min, max)
3. Define the law. You only need to provide one function which takes wavelength as angstrom. If your law is naturally written for inverse-micron, there is a helper function `aa_to_invum`.
       (::MyLaw)(wavelength::Real)
4. (Optional) enable `Unitful.jl` support by adding this function. If you are building a new law within `DustExtinction.jl` you can add your law to the code-gen list inside `DustExtinction.jl/src/DustExtinction.jl`.
       (l::MyLaw)(wavelength::Unitful.Quantity) = l(ustrip(u"angstrom", wavelength)) * u"mag"
"""
abstract type ExtinctionLaw end

Base.broadcastable(law::ExtinctionLaw) = Ref(law)

"""
    DustExtinction.bounds(::ExtinctionLaw)::Tuple
    DustExtinction.bounds(::Type{<:ExtinctionLaw})::Tuple

Get the natural wavelengths bounds for the extinction law, in angstrom
"""
bounds(::E) where {E <: ExtinctionLaw} = bounds(E)
bounds(::Type{<:ExtinctionLaw}) = (0, Inf)

checkbounds(::E, wave) where {E <: ExtinctionLaw} = checkbounds(E, wave)
function checkbounds(E::Type{<:ExtinctionLaw}, wave)
    b = bounds(E)
    return b[1] ≤ wave ≤ b[2]
end

"""
    redden(::ExtinctionLaw, wave, flux; Av=1)
    redden(::Type{ExtinctionLaw}, wave, flux; Av=1, law_kwargs...)

Redden the given `flux` using the given extinction law at the given wavelength.

If `wave` is `<:Real` then it is expected to be in angstrom and if it is `<:Unitful.Quantity` it will be automatically converted. `Av` is the total extinction value. The extinction law can be a constructed struct or a `Type`. If it is a `Type`, `law_kwargs` will be passed to the constructor.

# Examples

```jldoctest
julia> wave = 3000; flux = 1000;

julia> redden(CCM89, wave, flux; Rv=3.1)
187.38607779757183

julia> redden(CCM89(Rv=3.1), wave, flux; Av=2)
35.11354215235764
```

# See Also
[`deredden`](@ref)
"""
redden(L::Type{<:ExtinctionLaw}, wave, flux; Av = 1, kwargs...) = redden(L(values(kwargs)...), wave, flux; Av = Av)
redden(law::ExtinctionLaw, wave::Real, flux; Av = 1) = flux * 10^(-0.4 * Av * law(wave))
redden(law::ExtinctionLaw, wave::Quantity, flux::Real; Av = 1) = redden(law, ustrip(u"Å", wave), flux; Av = Av)
redden(law::ExtinctionLaw, wave::Quantity, flux::Quantity; Av = 1) = flux * (Av * law(wave))

"""
    deredden(::ExtinctionLaw, wave, flux; Av=1)
    deredden(::Type{ExtinctionLaw}, wave, flux; Av=1, law_kwargs...)

Deredden the given `flux` using the given extinction law at the given wavelength.

If `wave` is `<:Real` then it is expected to be in angstrom and if it is `<:Unitful.Quantity` it will be automatically converted. `Av` is the total extinction value. The extinction law can be a constructed struct or a `Type`. If it is a `Type`, `law_kwargs` will be passed to the constructor.

# Examples

```jldoctest
julia> wave = 3000; flux = 187.386;

julia> deredden(CCM89, wave, flux; Rv=3.1)
999.9995848273642

julia> deredden(CCM89(Rv=3.1), wave, flux; Av=2)
5336.573541539394
```

# See Also
[`redden`](@ref)
"""
deredden(L::Type{<:ExtinctionLaw}, wave, flux; Av = 1, kwargs...) = deredden(L(values(kwargs)...), wave, flux; Av = Av)
deredden(law::ExtinctionLaw, wave::Real, flux; Av = 1) = flux / 10^(-0.4 * Av * law(wave))
deredden(law::ExtinctionLaw, wave::Quantity, flux::Real; Av = 1) = deredden(law, ustrip(u"Å", wave), flux; Av = Av)
deredden(law::ExtinctionLaw, wave::Quantity, flux::Quantity; Av = 1) = flux / (Av * law(wave))

# --------------------------------------------------------------------------------
# bring in the support

include("deprecate.jl")
include("color_laws.jl")
include("dust_maps.jl")
include("fittable_laws.jl")

# --------------------------------------------------------------------------------
# Here be codegen!

# generate unitful support for the following laws
# this can be removed when julia support is pinned to 1.3 or higher,
# at which point adding `(l::ExtinctionLaw)(wave::Quantity)` is possible, until then
# using this code-gen does the trick but requires manually editing
# instead of providing support for all <: ExtinctionLaw
for law in [CCM89, OD94, CAL00, GCC09, VCG04]
    (l::law)(wavelength::Quantity) = l(ustrip(u"Å", wavelength)) * u"mag"
end

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

function __init__()
    # register our data dependencies
    register(DataDep("sfd98_map",
    """
    SFD98 Galactic Dust Maps
    Website: https://sncosmo.github.io
    """,
    ["https://sncosmo.github.io/data/dust/SFD_dust_4096_ngp.fits",
    "https://sncosmo.github.io/data/dust/SFD_dust_4096_sgp.fits"],
    ["50b6aaad0b880762d0fd081177802dcc17c39d7044a410dd5649e2dfd0503e97",
    "84891a59054adab44a7be54051e4dcf0e66e3f13eee0d845ce3739242f553b83"]))
end

end # module


text:
# canonical data version: 1.2.0

using Test

include("binary-search.jl")

@testset "default binary search" begin
    @testset "value in array" begin
        @test binarysearch([6], 6) == 1:1
        @test binarysearch([1, 3, 4, 6, 8, 9, 11], 6) == 4:4
        @test binarysearch([1, 3, 4, 6, 8, 9, 11], 1) == 1:1
        @test binarysearch([1, 3, 4, 6, 8, 9, 11], 11) == 7:7
        @test binarysearch([1, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 634], 144) == 10:10
        @test binarysearch([1, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377], 21) == 6:6
    end
    @testset "value not in array" begin
        @test binarysearch([1, 3, 4, 6, 8, 9, 11], 7) == 5:4
        @test binarysearch([1, 3, 4, 6, 8, 9, 11], 0) == 1:0
        @test binarysearch([1, 3, 4, 6, 8, 9, 11], 13) == 8:7
        @test binarysearch([], 1) == 1:0
    end
end

@testset "bonus tasks" begin
    @testset "reverse search" begin
        @testset "value in array" begin
            @test binarysearch([6], 6, rev = true) == 1:1
            @test binarysearch([11, 9, 8, 6, 4, 3, 1], 6, rev = true) == 4:4
            @test binarysearch([11, 9, 8, 6, 4, 3, 1], 1, rev = true) == 7:7
            @test binarysearch([11, 9, 8, 6, 4, 3, 1], 11, rev = true) == 1:1
            @test binarysearch([634, 377, 233, 144, 89, 55, 34, 21, 13, 8, 5, 3, 1], 144, rev = true) == 4:4
            @test binarysearch([377, 233, 144, 89, 55, 34, 21, 13, 8, 5, 3, 1], 21, rev = true) == 7:7
        end
        @testset "value not in array" begin
            @test binarysearch([11, 9, 8, 6, 4, 3, 1], 7, rev = true) == 4:3
            @test binarysearch([11, 9, 8, 6, 4, 3, 1], 0, rev = true) == 8:7
            @test binarysearch([11, 9, 8, 6, 4, 3, 1], 13, rev = true) == 1:0
            @test binarysearch([], 1, rev = true) == 1:0
        end
    end

    @testset "apply transformation" begin
        @testset "value in array" begin
            @test binarysearch([5.5], 6, by = round) == 1:1
            @test binarysearch([1.1, 2.9, 4.4, 5.5, 8.1, 9.0, 10.8], 6, by = round) == 4:4
            @test binarysearch([1.1, 2.9, 4.4, 5.5, 8.1, 9.0, 10.8], 1, by = round) == 1:1
            @test binarysearch([1.1, 2.9, 4.4, 5.5, 8.1, 9.0, 10.8], 11, by = round) == 7:7
            @test binarysearch([1.1, 2.9, 4.4, 5.5, 8.1, 9.0, 10.8], 11, by =  abs2∘round) == 7:7
            @test binarysearch([1, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 634], 144.4, by = round) == 10:10
            @test binarysearch([1, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377], 20.6, by = round) == 6:6
        end
        @testset "value not in array" begin
            @test binarysearch([1.1, 2.9, 4.4, 5.5, 8.1, 9.0, 10.8], 7, by = round) == 5:4
            @test binarysearch([1.1, 2.9, 4.4, 5.5, 8.1, 9.0, 10.8], 0, by = round) == 1:0
            @test binarysearch([1.1, 2.9, 4.4, 5.5, 8.1, 9.0, 10.8], 13, by = round) == 8:7
            @test binarysearch([], 1, by = round) == 1:0
        end
    end

    @testset "compare with > instead of <" begin
        # this is equivalent to searching in reverse order
        @testset "value in array" begin
            @test binarysearch([6], 6, lt = >) == 1:1
            @test binarysearch([11, 9, 8, 6, 4, 3, 1], 6, lt = >) == 4:4
            @test binarysearch([11, 9, 8, 6, 4, 3, 1], 1, lt = >) == 7:7
            @test binarysearch([11, 9, 8, 6, 4, 3, 1], 11, lt = >) == 1:1
            @test binarysearch([634, 377, 233, 144, 89, 55, 34, 21, 13, 8, 5, 3, 1], 144, lt = >) == 4:4
            @test binarysearch([377, 233, 144, 89, 55, 34, 21, 13, 8, 5, 3, 1], 21, lt = >) == 7:7
        end
        @testset "value not in array" begin
            @test binarysearch([11, 9, 8, 6, 4, 3, 1], 7, lt = >) == 4:3
            @test binarysearch([11, 9, 8, 6, 4, 3, 1], 0, lt = >) == 8:7
            @test binarysearch([11, 9, 8, 6, 4, 3, 1], 13, lt = >) == 1:0
            @test binarysearch([], 1, lt = >) == 1:0
        end
    end
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
for N in [Float64, Rational{Int}, Float32]
    I1 = Interval(N(-1), N(1))
    I2 = Interval(N(2), N(3))
    I3 = Interval(N(-2), N(-1))
    B1 = BallInf(N[1, 1], N(2))
    B2 = BallInf(N[3, 1], N(2))
    B3 = BallInf(N[3, 3], N(2))

    # constructor
    RI1 = Rectification(I1)
    RI2 = Rectification(I2)
    RI3 = Rectification(I3)
    RB1 = Rectification(B1)
    RB2 = Rectification(B2)
    RB3 = Rectification(B3)

    # set
    @test set(RI1) == I1
    @test set(RI2) == I2
    @test set(RI3) == I3
    @test set(RB1) == B1
    @test set(RB2) == B2
    @test set(RB3) == B3

    # dimension
    @test dim(RI1) == dim(RI2) == dim(RI3) == 1
    @test dim(RB1) == dim(RB2) == dim(RB3) == 2

    # support vector
    # hyperrectangles
    @test σ(N[1], RI1) == N[1] && σ(N[-1], RI1) == N[0]
    @test σ(N[1], RI2) == N[3] && σ(N[-1], RI2) == N[2]
    @test σ(N[1], RI3) == N[0] && σ(N[-1], RI3) == N[0]
    @test σ(N[1, 1], RB1) == N[3, 3] && σ(N[-1, 1], RB1) == N[0, 3]
    @test σ(N[1, 1], RB2) == N[5, 3] && σ(N[-1, 1], RB2) == N[1, 3]
    @test σ(N[1, 1], RB3) == N[5, 5] && σ(N[-1, 1], RB3) == N[1, 5]
    # Cartesian products
    RC1 = Rectification(I1 × I2)
    @test σ(N[1, 1], RC1) == N[1, 3] && σ(N[-1, -1], RC1) == N[0, 2]
    RC2 = Rectification(CartesianProductArray([I1, I2, I3]))
    @test σ(N[1, 1, 1], RC2) == N[1, 3, 0] &&
          σ(N[-1, -1, -1], RC2) == N[0, 2, 0]
    # other sets in 1D fall back to interval conversion
    @test σ(N[1], Rectification(Ball1(N[0], N(1)))) == N[1]
    # other sets in higher dimensions compute lazy linear maps and intersections
    @test σ(N[-1, -1], Rectification(Ball1(N[-2, -2], N(1)))) == N[0, 0]

    # an_element
    x = an_element(RI1)
    @test x isa Vector{N} && length(x) == 1 && N(0) <= x[1] <= N(1)
    x = an_element(RB2)
    @test x isa Vector{N} && length(x) == 2 && N(1) <= x[1] <= N(3) &&
          N(0) <= x[2] <= N(3)

    # membership
    @test N[-1, 1] ∉ RB1
    @test N[1, 1] ∈ RB1
    @test N[0, 1] ∈ RB1
    @test_throws ErrorException N[0, 4] ∈ RB1

    # emptiness
    @test !isempty(RI1)
    @test isempty(Rectification(EmptySet{N}(1)))

    # boundedness
    @test isbounded(RI1) && isboundedtype(typeof(RI1))
    @test isbounded(Rectification(HalfSpace(N[1], N(0))))
    RI4 = Rectification(Universe{N}(2))
    @test !isbounded(RI4) && !isboundedtype(typeof(RI4))
    @test !isbounded(Rectification(HalfSpace(N[-1], N(0))))

    # conversion
    @test convert(Interval, RI1) == Interval(N(0), N(1))
    @test convert(Hyperrectangle, RB1) ==
                  Hyperrectangle(low=N[0, 0], high=N[3, 3])

    # support function
    # upper right quadrant
    B = Ball1(N[2, 2], N(1))
    RB = Rectification(B)
    for d in [N[1, 0], N[-1, 0], N[0, 1], N[0, -1]]
        @test ρ(d, RB) == ρ(d, B)
    end
    # lower right quadrant
    B = Ball1(N[2, -2], N(1))
    RB = Rectification(B)
    for d in [N[1, 0], N[-1, 0]]
        @test ρ(d, RB) == ρ(d, B)
    end
    @test ρ(N[0, 1], RB) == ρ(N[0, -1], RB) == N(0)
    # upper left quadrant
    B = Ball1(N[-2, 2], N(1))
    RB = Rectification(B)
    for d in [N[0, 1], N[0, -1]]
        @test ρ(d, RB) == ρ(d, B)
    end
    @test ρ(N[1, 0], RB) == ρ(N[-1, 0], RB) == N(0)
    # lower left quadrant
    B = Ball1(N[-2, -2], N(1))
    RB = Rectification(B)
    @test ρ(N[1, 0], RB) == ρ(N[-1, 0], RB) == ρ(N[0, 1], RB) ==
          ρ(N[0, -1], RB) == N(0)

    # concretize
    @test concretize(RI1) == rectify(I1)
    @test concretize(RI2) == rectify(I2)
    @test concretize(RI3) == rectify(I3)
end

# tests that only work with Float64
for N in [Float64]
    # support function
    # some of the tests do not work because of insufficient precision in the
    # intersection; if the precision changes, these tests can be replaced by
    # their true (commented out) expected results

    # two right quadrants
    B = Ball1(N[2, 0], N(1))
    RB = Rectification(B)
    for d in [N[1, 0], N[-1, 0], N[0, 1]]
        @test ρ(d, RB) ≈ ρ(d, B)
    end
#     @test ρ(N[0, -1], RB) ≈ N(0)
    @test N(0) ≤ ρ(N[0, -1], RB) ≤ N(1e-9)
    # all four quadrants
    P = VPolygon([N[-1, 1], N[-1.5, 0.5], N[1.5, 0.5], N[1, -0.5]])
    RP = Rectification(P)
    @test ρ(N[1, 0], RP) ≈ N(1.5)
    @test ρ(N[0, 1], RP) ≈ N(1)
    @test ρ(N[1, 1], RP) ≈ ρ(N[1, 1], P) == N(2)
#     @test ρ(N[-1, 0], RP) ≈ N(0)
    @test N(0) ≤ ρ(N[-1, 0], RP) ≤ N(1e-8)
#     @test ρ(N[0, -1], RP) ≈ N(0)
    @test N(0) ≤ ρ(N[0, -1], RP) ≤ N(1e-8)
#     @test ρ(N[-1, -1], RP) ≈ N(0)
    @test N(0) ≤ ρ(N[-1, -1], RP) ≤ N(1.1e-1)

    @test σ(N[1, 1], Rectification(Ball1(N[0, 0], N(1)))) == N[0, 1]
end


text:
using Test, Defer

mutable struct E end
mutable struct F end

const closables = Set{Any}()
mutable struct Closable
  closed::Bool
  name::String
  Closable(name::String) = let a = new(false, name); push!(closables, a); a; end
end

Base.close(c::Closable) = (c.closed = true; nothing)

mutable struct Unclosable
  closed::Bool
  name::String
  Unclosable(name::String) = let a = new(false, name); push!(closables, a); a; end
end

Base.close(c::Unclosable) = (c.closed = true; throw(E()))

function verify_closed()
  setdiff!(closables, filter(x->x.closed, closables))
  if !isempty(closables)
    error(string("Unclosed objects: ", join(map(x->x.name, collect(closables)), ", ")))
  end
end

Closable("A")
@test_throws ErrorException verify_closed()
empty!(closables)

Unclosable("B")
@test_throws ErrorException verify_closed()
empty!(closables)

@scope begin
  @defer close(Closable("1"))
end
verify_closed()

@test_throws F @scope begin
  throw(F())
  @defer close(Closable("2"))
end
verify_closed()

@test_throws F @scope begin
  @defer close(Closable("3"))
  throw(F())
end
verify_closed()

@test_throws E @scope begin
  @defer close(Unclosable("4"))
end
verify_closed()

@test_throws F @scope begin
  throw(F())
  @defer close(Unclosable("5"))
end
verify_closed()

@test_throws CompositeException @scope begin
  @defer close(Unclosable("6"))
  throw(F())
end
verify_closed()

@test_throws CompositeException @scope begin
  @defer close(Unclosable("7"))
  @defer close(Unclosable("8"))
end
verify_closed()

@scope @! Closable("9")
verify_closed()

scope_nogc() do
  @! Closable("10")
end
verify_closed()

@scope let a = Closable("11")
  @! Closable("12")
end
verify_closed()

@scope let a = Closable("11a"),  b = Closable("11b")
  @! Closable("12a")
end
verify_closed()

@scope f1(s::String) = @! Closable(s)
@scope function f2(t::String) @! Closable(t) end
f1("13")
f2("14")
verify_closed()


text:
using Test, Suppressor
import HMCResearchRandomWalks.Arlett
using HMCResearchRandomWalks.Arlett:
    approach_wall!, boundary_check,
    north!, eastwest!, evaluate_proposed,
    one_step!, insafebounds, inescapebounds,
    randseed

# see arlett_frontend_tests.jl for run_sim.jl
# see arlett_backend_tests.jl for all other location ("in...") functions

const smallest_stepsize = min(Arlett.WATER_STEP_SIZE,
                            Arlett.ENZ_STEP_SIZE,
                            Arlett.PPD_STEP_SIZE)
const significant_err = smallest_stepsize / Arlett.MAX_STEPS_PER_WALK
const diagonal = sqrt(2) / 2

function test_xy(testxy::Array{Float64,1},
                            correctxy::Array{Float64,1})
    @test isapprox(testxy[1], correctxy[1], atol=significant_err)
    @test isapprox(testxy[2], correctxy[2], atol=significant_err)
end

# Vertical Proposed Point via North!
# enz_water
@testset "vertical displacement at water-enz" begin
    Arlett.set_FLOW_BIAS(false)
    correct_below = [0.0, 1.6475]
    correct_above = [0.0, 1.7]
    wrong_below = [0.0, 1.6]
    wrong_above = [0.0, 1.6525]
    test_below = north!(wrong_below, correct_above, 0.0, -1.0, "water", "enz")
    test_above = north!(wrong_above, correct_below, 0.0, 1.0, "enz", "water")
    test_below_ep = evaluate_proposed(correct_above, 0.0, -1.0)[1]
    test_above_ep = evaluate_proposed(correct_below, 0.0, 1.0)[1]
    test_below_bc = boundary_check(correct_above, 0.0, -1.0)[1]
    test_above_bc = boundary_check(correct_below, 0.0, 1.0)[1]
    test_xy(test_below, correct_below)
    test_xy(test_above, correct_above)
    test_xy(test_below_ep, correct_below)
    test_xy(test_above_ep, correct_above)
    test_xy(test_below_bc, correct_below)
    test_xy(test_above_bc, correct_above)
end

@testset "positive slope displacement at enz-ppd" begin
    Arlett.set_FLOW_BIAS(false)
    correct_below = [0.0, 0.15 - 0.001 * diagonal / 2]
    correct_above = [0.0 + (0.001 + 0.005) * diagonal / 2,
                        0.15 + 0.005 * diagonal / 2]
    wrong_above = [0 + 0.001 * diagonal,
                        0.15 + 0.001 * diagonal / 2]
    wrong_below = [0.0 + (0.001 - 0.005) * diagonal / 2,
                        0.15 - 0.005 * diagonal / 2]
    test_above = north!(wrong_above, correct_below,
                         diagonal, diagonal, "ppd", "enz")
    test_below = north!(wrong_below, correct_above,
                         -1 * diagonal, -1 * diagonal, "enz", "ppd")
    test_above_ep = evaluate_proposed(correct_below, diagonal, diagonal)[1]
    test_below_ep = evaluate_proposed(correct_above, -1 * diagonal, -1 * diagonal)[1]
    test_above_bc = boundary_check(correct_below, diagonal, diagonal)[1]
    test_below_bc = boundary_check(correct_above, -1 * diagonal, -1 * diagonal)[1]
    test_xy(test_below, correct_below)
    test_xy(test_above, correct_above)
    test_xy(test_below_ep, correct_below)
    test_xy(test_above_ep, correct_above)
    test_xy(test_below_bc, correct_below)
    test_xy(test_above_bc, correct_above)
end

@testset "negative slope displacement at water-ppd" begin
    Arlett.set_FLOW_BIAS(false)
    for x in [-90, -45, 45, 90]
        correct_above = [x, 0.15 + 0.1 * diagonal / 2]
        correct_below = [x + (0.1 + 0.001) * diagonal / 2,
                            0.15 - 0.001 * diagonal / 2]
        wrong_below = [x + 0.1 * diagonal,
                            0.15 - 0.1 * diagonal / 2]
        wrong_above = [x + (0.1 - 0.001) * diagonal / 2,
                            0.15 + 0.001 * diagonal / 2]
        test_above = north!(wrong_above, correct_below,
                            -1 * diagonal, diagonal, "ppd", "water")
        test_below = north!(wrong_below, correct_above,
                            diagonal, -1 * diagonal, "water", "ppd")
        test_above_ep = evaluate_proposed(correct_below, -1 * diagonal, diagonal)[1]
        test_below_ep = evaluate_proposed(correct_above, diagonal, -1 * diagonal)[1]
        test_above_bc = boundary_check(correct_below, -1 * diagonal, diagonal)[1]
        test_below_bc = boundary_check(correct_above, diagonal, -1 * diagonal)[1]
        test_xy(test_below, correct_below)
        test_xy(test_above, correct_above)
        test_xy(test_below_ep, correct_below)
        test_xy(test_above_ep, correct_above)
        test_xy(test_below_bc, correct_below)
        test_xy(test_above_bc, correct_above)
    end
end

@testset "displacement at enz-water east" begin
    correct_enz = [13.4975, 1.6]
    correct_water = [13.55, 1.6]
    test_enz = boundary_check(correct_water, -1.0, 0.0)[1]
    test_water = boundary_check(correct_enz, 1.0, 0.0)[1]
    test_xy(test_enz, correct_enz)
    test_xy(test_water, correct_water)
end

@testset "displacement at enz-water west" begin
    correct_enz = [-13.4975, 1.6]
    correct_water = [-13.55, 1.6]
    test_enz = boundary_check(correct_water, 1.0, 0.0)[1]
    test_water = boundary_check(correct_enz, -1.0, 0.0)[1]
    test_xy(test_enz, correct_enz)
    test_xy(test_water, correct_water)
end

@testset "universal location functions in one_step!.jl" begin
    for safexy in [[0.0, 2.0],[122.2, 122.2],[-122.2,112.2]]
        @test Arlett.insafebounds(safexy)
        @test Arlett.inescapebounds(safexy)
    end
    for escapexy in [[0.0, 0.0],[122.5, 122.5],[-122.5, 122.5]]
        @test Arlett.inescapebounds(escapexy)
        @test !Arlett.insafebounds(escapexy)
    end
    for notescapexy in [[122.6, 122.5],[122.6, 122.6],[-122.5, 122.6],
                            [-122.6, 122.5],[0.0, -2.0]]
        @test !Arlett.inescapebounds(notescapexy)
        @test !Arlett.insafebounds(notescapexy)
    end
end

@testset "approch_wall!.jl and wallcases()" begin
    correct_wall_side = [57.42 + 0.1 * 0.5, 1.0]
    test_wall_side_aw = approach_wall!([57.42, 1.0], 1.0, 0.0, 0.1)
    test_wall_side_wc = Arlett.wallcases([57.42, 1.0], 1.0, 0.0, 0.1)[1]
    test_wall_side_bc = boundary_check([57.42, 1.0], 1.0, 0.0)[1]
    test_xy(test_wall_side_aw, correct_wall_side)
    test_xy(test_wall_side_wc, correct_wall_side)
    test_xy(test_wall_side_bc, correct_wall_side)

    Arlett.set_CATALASE_ON_WALLS(false)
    Arlett.set_FLOW_BIAS(false)
    correct_wall_top = [15, 1.51 - 0.1 * 0.5 * 0.5 * 0.5 * 0.5]
    test_wall_top_aw = approach_wall!([15, 1.51], 0.0, -1.0, 0.1)
    test_wall_top_wc = Arlett.wallcases([15, 1.51], 0.0, -1.0, 0.1)[1]
    test_wall_top_bc = boundary_check([15, 1.51], 0.0, -1.0)[1]
    test_xy(test_wall_top_aw, correct_wall_top)
    test_xy(test_wall_top_wc, correct_wall_top)
    test_xy(test_wall_top_bc, correct_wall_top)

    Arlett.set_CATALASE_ON_WALLS(true)
    for wallx in [-115, -60, -30, 30, 60, 115]
        correct_catalase_top = undef
        test_catalase_top, test_report = boundary_check([wallx, 1.51], 0.0, -1.0)
        @test test_catalase_top == undef
        @test test_report == "top wall"
    end

    side_corner_coinflip = false
    top_corner_coinflip = false
    while (!side_corner_coinflip || !top_corner_coinflip)
        corner_case, corner_report = boundary_check([32.51, 1.51], -1 * diagonal, -1 * diagonal)
        if corner_report == "side wall" side_corner_coinflip = true end
        if corner_report == "top wall" top_corner_coinflip = true end
    end
    @test side_corner_coinflip && top_corner_coinflip
end

@testset "onestep!" begin
    init_xy = [0.0, 50.0]
    copy_xy = copy(init_xy)

    next_xy = one_step!(copy_xy)[1]

    @test init_xy != copy_xy
    @test copy_xy === next_xy

    @test -0.1 <= next_xy[1] <= 0.1
    @test 49.9 <= next_xy[2] <= 50.1
    d = sqrt(
        (next_xy[1] - init_xy[1])^2 +
        (next_xy[2] - init_xy[2])^2
    )
    @test isapprox(d, 0.1, atol=significant_err)
end

@testset "randseed" begin
    @test typeof(Arlett.randseed()) == Int
    @test 0 <= Arlett.randseed() <= 10000
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
begin
# Initializing the system

# Hamiltonian of the system
hamiltonian = sigmax
# Initial State
init_state = normalize!(QuArray([0.5+0.1im, 0.2+0.2im]))
# Time step
tlist = 0.:0.1:2*pi
# Initial State density matrix
init_state_dm = init_state*init_state'

const solver = @compat Dict{Any, Any}(:qu_euler => QuEuler, :qu_cn => QuCrankNicolson, :qu_krylov => QuKrylov,
                                      :qu_ode45 => QuODE45, :qu_ode78 => QuODE78, :qu_ode23s => QuODE23s,
                                      :quexpm_expo => QuExpokit, :quexpm_expmv => QuExpmV)

for (key,value) in solver
    @eval begin
        $key = QuStateEvolution(hamiltonian, init_state, tlist, $value())
    end
end

const next_state = @compat Dict{Any, Any}(:next_state_euler => qu_euler, :next_state_cn => qu_cn, :next_state_krylov => qu_krylov,
                                          :next_state_ode45 => qu_ode45, :next_state_ode78 => qu_ode78, :next_state_ode23s => qu_ode23s,
                                          :next_state_quexpm_expo => quexpm_expo, :next_state_quexpm_expmv => quexpm_expmv)

for (key,value) in next_state
    @eval begin
        $key = next($value, start($value))
    end
end
next_state_actual = expm(-im*sigmax*0.1)*init_state

@test_approx_eq_eps coeffs(next_state_euler[1][2]) coeffs(next_state_actual) 1e-2
@test_approx_eq_eps coeffs(next_state_cn[1][2]) coeffs(next_state_actual) 1e-4
@test_approx_eq coeffs(next_state_krylov[1][2]) coeffs(next_state_actual)
@test_approx_eq coeffs(next_state_ode45[1][2]) coeffs(next_state_actual)
@test_approx_eq coeffs(next_state_ode78[1][2]) coeffs(next_state_actual)
@test_approx_eq_eps coeffs(next_state_ode23s[1][2]) coeffs(next_state_actual) 1.0e-5
@test_approx_eq coeffs(next_state_quexpm_expo[1][2]) coeffs(next_state_actual)
@test_approx_eq coeffs(next_state_quexpm_expmv[1][2]) coeffs(next_state_actual)

# Tests for the different constructs for Liouville von Neumann Equation

lvn = QuLiouvillevonNeumannEq(QuDynamics.liouvillian_op(sigmax))

qexpokit_lvn = QuStateEvolution(lvn, init_state*init_state', tlist, QuExpokit())
qexpokit_dm = QuStateEvolution(sigmax, init_state*init_state', tlist, QuExpokit())
next_state_qexpokit_lvn = next(qexpokit_lvn, start(qexpokit_lvn))
next_state_qexpokit_dm = next(qexpokit_dm, start(qexpokit_dm))
@assert next_state_qexpokit_lvn[1][2] == next_state_qexpokit_dm[1][2]
qode45_lvn = QuStateEvolution(lvn, init_state*init_state', tlist, QuODE45())
qode45_dm = QuStateEvolution(sigmax, init_state*init_state', tlist, QuODE45())
next_state_qode45_lvn = next(qode45_lvn, start(qode45_lvn))
next_state_qode45_dm = next(qode45_dm, start(qode45_dm))
@assert next_state_qode45_lvn[1][2] == next_state_qode45_dm[1][2]
qeuler_lvn = QuStateEvolution(lvn, init_state*init_state', tlist, QuEuler())
qeuler_dm = QuStateEvolution(sigmax, init_state*init_state', tlist, QuEuler())
next_state_qeuler_lvn = next(qeuler_lvn, start(qeuler_lvn))
next_state_qeuler_dm = next(qeuler_dm, start(qeuler_dm))
@assert next_state_qeuler_lvn[1][2] == next_state_qeuler_dm[1][2]

# Tests comparing results for Liouville von Neumann Equation

next_state_actual = expm(-im*QuDynamics.liouvillian_op(sigmax)*0.1)*vec(init_state_dm)
@test_approx_eq coeffs(vec(next_state_qexpokit_lvn[1][2])) coeffs(next_state_actual)
@test_approx_eq_eps coeffs(vec(next_state_qode45_dm[1][2])) coeffs(next_state_actual) 1e-10
@test_approx_eq_eps coeffs(vec(next_state_qeuler_lvn[1][2])) coeffs(next_state_actual) 1e-2

# Tests for the different constructs for Lindblad Master Equation

c_ops = [lowerop(2)]
lmeq = QuLindbladMasterEq(sigmax, c_ops)

qexpokit_lmeq = QuStateEvolution(lmeq, init_state*init_state', tlist, QuExpokit())
qexpokit_dm = QuStateEvolution(sigmax, c_ops, init_state*init_state', tlist, QuExpokit())
next_state_qexpokit_lmeq = next(qexpokit_lmeq, start(qexpokit_lmeq))
next_state_qexpokit_dm = next(qexpokit_dm, start(qexpokit_dm))
@assert next_state_qexpokit_lmeq[1][2] == next_state_qexpokit_dm[1][2]
qode78_lmeq = QuStateEvolution(lmeq, init_state*init_state', tlist, QuODE78())
qode78_dm = QuStateEvolution(sigmax, c_ops, init_state*init_state', tlist, QuODE78())
next_state_qode78_lmeq = next(qode78_lmeq, start(qode78_lmeq))
next_state_qode78_dm = next(qode78_dm, start(qode78_dm))
@assert next_state_qode78_lmeq[1][2] == next_state_qode78_dm[1][2]
qcn_lmeq = QuStateEvolution(lmeq, init_state*init_state', tlist, QuCrankNicolson())
qcn_dm = QuStateEvolution(sigmax, c_ops, init_state*init_state', tlist, QuCrankNicolson())
next_state_qcn_lmeq = next(qcn_lmeq, start(qcn_lmeq))
next_state_qcn_dm = next(qcn_dm, start(qcn_dm))
@assert next_state_qcn_lmeq[1][2] == next_state_qcn_dm[1][2]

# Tests comparing results for Lindblad Master Equation

next_state_actual = expm(-im*QuDynamics.lindblad_op(sigmax, [lowerop(2)])*0.1)*vec(init_state_dm)
@test_approx_eq coeffs(vec(next_state_qexpokit_lmeq[1][2])) coeffs(next_state_actual)
@test_approx_eq coeffs(vec(next_state_qexpokit_dm[1][2])) coeffs(next_state_actual)
@test_approx_eq coeffs(vec(next_state_qexpokit_lmeq[1][2])) coeffs(next_state_actual)

# TODO : QuBase.jl complex WIP. Issue 38. https://github.com/JuliaQuantum/QuBase.jl/pull/38
Base.complex{B<:QuBase.OrthonormalBasis}(qarr::QuBase.AbstractQuArray{B}) = QuBase.similar_type(qarr)(complex(coeffs(qarr)), bases(qarr))

# Monte Carlo Wave Function Method

rho = init_state_dm
rhos = Array(typeof(rho), length(tlist)-1)
ps = Array(typeof(rho), length(tlist)-1)
qumcwfen = QuMCWFEnsemble(complex(init_state), 1000)

for i=1:length(tlist)-1
    rhos[i] = complex(zeros(rho))
    ps[i] = complex(zeros(rho))
end

for psi0 in qumcwfen
    i = 1
    mcwf = QuMCWF()
    for (t,psi) in QuStateEvolution(hamiltonian, c_ops, psi0, tlist, mcwf)
        rhos[i] = rhos[i] + (psi*psi')/length(qumcwfen)/norm(psi)^2
        i = i + 1
    end
end

prop_expm = QuStateEvolution(hamiltonian, c_ops, rho, tlist, QuExpmV())
j = 1
for (t, rhot) in prop_expm
    ps[j] = rhot
    j = j + 1
end

ind=rand(1:(length(tlist)-1))
@test_approx_eq_eps coeffs(vec(rhos[ind])) coeffs(vec(ps[ind])) 1e-1

# tests for  evolution opeator.
evolved_state = coeffs(propagate(QuExpokit(), QuSchrodingerEq(hamiltonian), tlist[2], tlist[1], init_state))
@test_approx_eq coeffs(QuEvolutionOp(sigmax, 0.1)*init_state) evolved_state
@test_approx_eq coeffs(QuEvolutionOp(sigmax, tlist[5], tlist[4])*init_state) evolved_state
@test_approx_eq coeffs(QuEvolutionOp(QuSchrodingerEq(sigmax), tlist[5], tlist[4])*init_state) evolved_state

# caching and uncached version testing
@assert QuLindbladMasterEqUncached(sigmax, [sigmax]).lindblad == nothing
@assert QuLindbladMasterEq(sigmax, [sigmax]).lindblad == QuDynamics.operator(QuLindbladMasterEqUncached(sigmax,  [sigmax]))
end


text:
function SparseRealVector()
    return SparseRealVector(())
end


text:
using JLD2, FileIO
using Plots
gr()

# load primegapsdicts from disk
gapsdictarray = []
for i = 2:13
    append!(gapsdictarray, [load_object("gapsdict10x$(i).jld2")])
end

# fit linear polynomial to log(N_gaps)
using Polynomials

function linfittodict(gapsdict)
    x = Int64.(collect(keys(gapsdict)))
    y = log.(collect(values(gapsdict)))
    fit(x, y, 1)
end

function polygradiant(poly)
    poly[1]
end

linpolyarray = map(linfittodict, gapsdictarray)
gradiantarray = map(polygradiant, linpolyarray)

function polyplotrange(poly)
    [2, floor(roots(poly)[1])]
end

# plot primegaps and linpolyfits
begin
    fig = plot(title="Primegaps for all Primes p < N", yaxis=:log)

    xlabel!(fig, "size of gap")
    ylabel!(fig, "number of occurences")

    for i = 1:12
        scatter!(fig, gapsdictarray[i], label="N=10^$(i+1)")

        xs = polyplotrange(linpolyarray[i])
        ys = exp.(linpolyarray[i].(xs))

        plot!(fig, xs, ys, label=nothing, line=(:dash, :red))
    end

    vline!(fig, [6, 30, 210], label=nothing, line=(:dot, :black))

    plot!(fig, dpi=300)
end

# plot gradient of polyfits
begin
    fig_gradient = plot(exp10.(2:13), gradiantarray, marker=:dot)

    yaxis!(fig_gradient, "gradiant of linear fit")
    xaxis!(fig_gradient, "N")

    xticks!(fig_gradient, [10^i for i in 2:13])

    plot!(fig_gradient, legend=false, xaxis=:log10, dpi=300)
end

# compare ratio of gaps length 6 and 30
function sixtothirty(gapsdict)
    get(gapsdict, 6, 1)/get(gapsdict, 30, 1)
end

sixtothirtyarray = map(sixtothirty, gapsdictarray)

begin
    fig_sixtothirty = plot([10^i for i in 4:13], sixtothirtyarray[3:end], marker=:dot)

    yaxis!(fig_sixtothirty, "ratio of occurence 6/30")
    xaxis!(fig_sixtothirty, "N")

    xticks!(fig_sixtothirty, [10^i for i in 4:13])
    # yticks!(fig_sixtothirty, vcat([1], collect(0:5:25)))
    # yticks!(fig_sixtothirty, 0:5:25)

    ylims!(fig_sixtothirty, (0, ylims(fig_sixtothirty)[2]))

    # hline!(fig_sixtothirty, [1], line=(:dot, :black))

    plot!(fig_sixtothirty, legend=false, xaxis=:log10, dpi=300)
end

# save plots to disk
savefig(fig, "Primegaps.png")
savefig(fig_sixtothirty, "sixtothirty.png")
savefig(fig_gradient, "gradient.png")

max_stars_repo_name: chortbauer/primegaps

max_issues_repo_name: chortbauer/primegaps

max_forks_repo_name: chortbauer/primegaps


text:
# NOTE: This script is only supposed to be executed by our own build script.
# see: https://docs.travis-ci.com/user/languages/julia/#dependency-management

import Pkg

@static if VERSION < v"1.4"
    Pkg.add([Pkg.PackageSpec(name="Documenter", version=v"0.26"),
             Pkg.PackageSpec(name="Symbolics", version=v"0.1.2"),
             Pkg.PackageSpec(name="TaylorModels", version=v"0.0.1")])
end

@static if VERSION > v"1.4"
    Pkg.add([Pkg.PackageSpec(name="Javis", version=v"0.5"),
             Pkg.PackageSpec(name="Symbolics", version=v"0.1.32")])
end


text:
# NOTE: This file contains machine-generated code by translate_exported_functions.jl.

# hFILE *hopen(const char *filename, const char *mode, ...) HTS_RESULT_USED;
# ErrorException("unknown token: '.'")

# hFILE *hdopen(int fd, const char *mode) HTS_RESULT_USED;
@defun hdopen(fd::Cint, mode::Cstring)::Ptr{hFILE}

# int hisremote(const char *filename) HTS_RESULT_USED;
@defun hisremote(filename::Cstring)::Cint

# char *haddextension(struct kstring_t *buffer, const char *filename, int replace, const char *extension) HTS_RESULT_USED;
@defun haddextension(buffer::Ptr{kstring_t}, filename::Cstring, replace::Cint, extension::Cstring)::Ptr{Cchar}

# int hclose(hFILE *fp) HTS_RESULT_USED;
@defun hclose(fp::Ptr{hFILE})::Cint

# void hclose_abruptly(hFILE *fp);
@defun hclose_abruptly(fp::Ptr{hFILE})::Cvoid

# off_t hseek(hFILE *fp, off_t offset, int whence) HTS_RESULT_USED;
@defun hseek(fp::Ptr{hFILE}, offset::Coff_t, whence::Cint)::Coff_t

# ssize_t hgetdelim(char *buffer, size_t size, int delim, hFILE *fp) HTS_RESULT_USED;
@defun hgetdelim(buffer::Ptr{Cchar}, size::Csize_t, delim::Cint, fp::Ptr{hFILE})::Cssize_t

# char *hgets(char *buffer, int size, hFILE *fp) HTS_RESULT_USED;
@defun hgets(buffer::Ptr{Cchar}, size::Cint, fp::Ptr{hFILE})::Ptr{Cchar}

# ssize_t hpeek(hFILE *fp, void *buffer, size_t nbytes) HTS_RESULT_USED;
@defun hpeek(fp::Ptr{hFILE}, buffer::Ptr{Cvoid}, nbytes::Csize_t)::Cssize_t

# int hflush(hFILE *fp) HTS_RESULT_USED;
@defun hflush(fp::Ptr{hFILE})::Cint

# char *hfile_mem_get_buffer(hFILE *file, size_t *length);
@defun hfile_mem_get_buffer(file::Ptr{hFILE}, length::Ptr{Csize_t})::Ptr{Cchar}

# char *hfile_mem_steal_buffer(hFILE *file, size_t *length);
@defun hfile_mem_steal_buffer(file::Ptr{hFILE}, length::Ptr{Csize_t})::Ptr{Cchar}


text:
#=
getindex.jl
Indexing support for LinearMapAX
2018-01-19, Jeff Fessler, University of Michigan
=#

# [end]
function Base.lastindex(A::LinearMapAX)
    return prod(size(A._lmap))
end

# [?,end] and [end,?]
function Base.lastindex(A::LinearMapAX, d::Int)
    return size(A._lmap, d)
end

# A[i,j]
function Base.getindex(A::LinearMapAX, i::Int, j::Int)
    T = eltype(A)
    e = zeros(T, size(A._lmap,2)); e[j] = one(T)
    tmp = A._lmap * e
    return tmp[i]
end

# A[:,j]
# it is crucial to provide this function rather than to inherit from
# Base.getindex(A::AbstractArray, ::Colon, ::Int)
# because Base.getindex does this by iterating (I think).
function Base.getindex(A::LinearMapAX, ::Colon, j::Int)
    T = eltype(A)
    e = zeros(T, size(A,2)); e[j] = one(T)
    return A._lmap * e
end

# A[ii,j]
Base.getindex(A::LinearMapAX, ii::Indexer, j::Int) = A[:,j][ii]

# A[i,jj]
Base.getindex(A::LinearMapAX, i::Int, jj::Indexer) = A[i,:][jj]

# A[:,jj]
# this one is also important for efficiency
Base.getindex(A::LinearMapAX, ::Colon, jj::AbstractVector{Bool}) =
    A[:,findall(jj)]
Base.getindex(A::LinearMapAX, ::Colon, jj::Indexer) =
    hcat([A[:,j] for j in jj]...)

# A[ii,:]
# trick: cannot use A' for a FunctionMap with no fc
function Base.getindex(A::LinearMapAX, ii::Indexer, ::Colon)
    if (:fc in propertynames(A._lmap)) && isnothing(A._lmap.fc)
        return hcat([A[ii,j] for j in 1:size(A,2)]...)
    else
        return A'[:,ii]'
    end
end

# A[ii,jj]
Base.getindex(A::LinearMapAX, ii::Indexer, jj::Indexer) = A[:,jj][ii,:]

# A[k]
function Base.getindex(A::LinearMapAX, k::Int)
    c = CartesianIndices(size(A._lmap))[k] # is there a more elegant way?
    return A[c[1], c[2]]
end

# A[kk]
Base.getindex(A::LinearMapAX, kk::AbstractVector{Bool}) = A[findall(kk)]
Base.getindex(A::LinearMapAX, kk::Indexer) = [A[k] for k in kk]

# A[i,:]
# trick: one row slice returns a 1D ("column") vector
Base.getindex(A::LinearMapAX, i::Int, ::Colon) = A[[i],:][:]

# A[:,:] = Matrix(A)
Base.getindex(A::LinearMapAX, ::Colon, ::Colon) = Matrix(A._lmap)

# A[:]
Base.getindex(A::LinearMapAX, ::Colon) = A[:,:][:]


text:
# # Predefined parameter sets

## load packages
using Photosynthesis
using PlotPlants
FT = Float32;
#------------------------------------------------------------------------------

# ## Jmax
_td_1 = JmaxTDBernacchi(FT);
_td_2 = JmaxTDCLM(FT);
_td_3 = JmaxTDLeuning(FT);
_ts   = collect(FT, 273:1:323);
_jm_1 = photo_TD_from_val.([_td_1], FT(100), _ts);
_jm_2 = photo_TD_from_val.([_td_2], FT(100), _ts);
_jm_3 = photo_TD_from_val.([_td_3], FT(100), _ts);

_fig,_axes = create_canvas("Jmax");
_ax1 = _axes[1];
_ax1.plot(_ts .- 273.15, _jm_1, "k-", label="Bernacchi");
_ax1.plot(_ts .- 273.15, _jm_2, "k:", label="CLM");
_ax1.plot(_ts .- 273.15, _jm_3, "k--", label="Leuning");
_ax1.set_xlabel("Leaf temperature (°C)");
_ax1.set_ylabel("Jcmax (μmol m⁻² s⁻¹)");
_ax1.legend(loc="upper left");
_fig
#------------------------------------------------------------------------------

# ## Kc
_td_1 = KcTDBernacchi(FT);
_td_2 = KcTDCLM(FT);
_ts   = collect(FT, 273:1:323);
_kc_1 = photo_TD_from_set.([_td_1], _ts);
_kc_2 = photo_TD_from_set.([_td_2], _ts);

_fig,_axes = create_canvas("Kc");
_ax1 = _axes[1];
_ax1.plot(_ts .- 273.15, _kc_1, "k-", label="Bernacchi");
_ax1.plot(_ts .- 273.15, _kc_2, "k:", label="CLM");
_ax1.set_xlabel("Leaf temperature (°C)");
_ax1.set_ylabel("Kc (Pa)");
_ax1.legend(loc="upper left");
_fig
#------------------------------------------------------------------------------

# ## Ko
_td_1 = KoTDBernacchi(FT);
_td_2 = KoTDCLM(FT);
_ts   = collect(FT, 273:1:323);
_ko_1 = photo_TD_from_set.([_td_1], _ts);
_ko_2 = photo_TD_from_set.([_td_2], _ts);

_fig,_axes = create_canvas("Ko");
_ax1 = _axes[1];
_ax1.plot(_ts .- 273.15, _ko_1, "k-", label="Bernacchi");
_ax1.plot(_ts .- 273.15, _ko_2, "k:", label="CLM");
_ax1.set_xlabel("Leaf temperature (°C)");
_ax1.set_ylabel("Ko (Pa)");
_ax1.legend(loc="upper left");
_fig
#------------------------------------------------------------------------------

# ## Kpep
_td_1 = KpepTDBoyd(FT);
_td_2 = KpepTDCLM(FT);
_ts   = collect(FT, 273:1:323);
_kp_1 = photo_TD_from_set.([_td_1], _ts);
_kp_2 = photo_TD_from_set.([_td_2], _ts);

_fig,_axes = create_canvas("Kpep");
_ax1 = _axes[1];
_ax1.plot(_ts .- 273.15, _kp_1, "k-", label="Boyd");
_ax1.plot(_ts .- 273.15, _kp_2, "k:", label="CLM");
_ax1.set_xlabel("Leaf temperature (°C)");
_ax1.set_ylabel("Kpep (Pa)");
_ax1.legend(loc="upper left");
_fig
#------------------------------------------------------------------------------

# ## Respiration
_td_1 = RespirationTDBernacchi(FT);
_td_2 = RespirationTDCLM(FT);
_td_3 = Q10TDAngiosperm(FT);
_td_4 = Q10TDGymnosperm(FT);
_ts   = collect(FT, 273:1:323);
_rd_1 = photo_TD_from_val.([_td_1], FT(1), _ts);
_rd_2 = photo_TD_from_val.([_td_2], FT(1), _ts);
_rd_3 = photo_TD_from_val.([_td_3], FT(1), _ts);
_rd_4 = photo_TD_from_val.([_td_4], FT(1), _ts);

_fig,_axes = create_canvas("Respiration");
_ax1 = _axes[1];
_ax1.plot(_ts .- 273.15, _rd_1, "k-", label="Bernacchi");
_ax1.plot(_ts .- 273.15, _rd_2, "k:", label="CLM");
_ax1.plot(_ts .- 273.15, _rd_3, "r-", label="Q10 Angiosperm");
_ax1.plot(_ts .- 273.15, _rd_4, "r:", label="Q10 Gymnosperm");
_ax1.set_xlabel("Leaf temperature (°C)");
_ax1.set_ylabel("Respiration (μmol m⁻² s⁻¹)");
_ax1.legend(loc="upper left");
_fig
#------------------------------------------------------------------------------

# ## Vcmax
_td_1 = VcmaxTDBernacchi(FT);
_td_2 = VcmaxTDCLM(FT);
_td_3 = VcmaxTDLeuning(FT);
_ts   = collect(FT, 273:1:323);
_vc_1 = photo_TD_from_val.([_td_1], FT(100), _ts);
_vc_2 = photo_TD_from_val.([_td_2], FT(100), _ts);
_vc_3 = photo_TD_from_val.([_td_3], FT(100), _ts);

_fig,_axes = create_canvas("Vcmax");
_ax1 = _axes[1];
_ax1.plot(_ts .- 273.15, _vc_1, "k-", label="Bernacchi");
_ax1.plot(_ts .- 273.15, _vc_2, "k:", label="CLM");
_ax1.plot(_ts .- 273.15, _vc_3, "k--", label="Leuning");
_ax1.set_xlabel("Leaf temperature (°C)");
_ax1.set_ylabel("Vcmax (μmol m⁻² s⁻¹)");
_ax1.legend(loc="upper left");
_fig
#------------------------------------------------------------------------------

# ## Vomax
_td_1 = VomaxTDBernacchi(FT);
_ts   = collect(FT, 273:1:323);
_vo_1 = photo_TD_from_val.([_td_1], FT(100), _ts);

_fig,_axes = create_canvas("Vomax");
_ax1 = _axes[1];
_ax1.plot(_ts .- 273.15, _vo_1, "k-", label="Bernacchi");
_ax1.set_xlabel("Leaf temperature (°C)");
_ax1.set_ylabel("Vomax (μmol m⁻² s⁻¹)");
_ax1.legend(loc="upper left");
_fig
#------------------------------------------------------------------------------

# ## Vpmax
_td_1 = VpmaxTDBoyd(FT);
_ts   = collect(FT, 273:1:323);
_vp_1 = photo_TD_from_val.([_td_1], FT(100), _ts);

_fig,_axes = create_canvas("Vpmax");
_ax1 = _axes[1];
_ax1.plot(_ts .- 273.15, _vp_1, "k-", label="Boyd");
_ax1.set_xlabel("Leaf temperature (°C)");
_ax1.set_ylabel("Vpmax (μmol m⁻² s⁻¹)");
_ax1.legend(loc="upper left");
_fig
#------------------------------------------------------------------------------

# ## Γ*
_td_1 = ΓStarTDBernacchi(FT);
_td_2 = ΓStarTDCLM(FT);
_ts   = collect(FT, 273:1:323);
_Γs_1 = photo_TD_from_set.([_td_1], _ts);
_Γs_2 = photo_TD_from_set.([_td_2], _ts);

_fig,_axes = create_canvas("Γ*");
_ax1 = _axes[1];
_ax1.plot(_ts .- 273.15, _Γs_1, "k-", label="Bernacchi");
_ax1.plot(_ts .- 273.15, _Γs_2, "k:", label="CLM");
_ax1.set_xlabel("Leaf temperature (°C)");
_ax1.set_ylabel("Γ* (Pa)");
_ax1.legend(loc="upper left");
_fig
#------------------------------------------------------------------------------


text:
struct LinearLS{T<:Real}
    x :: Vector{T}  # x data
    y :: Vector{T}  # y data
    f :: Function   # y = m*f(x) + q
    w :: Vector{T}  # weights (inverse a priori error on the `y`'s)
    m :: Float64    # slope
    q :: Float64    # intercept
    function LinearLS(x::AbstractVector{<:Real}, y::AbstractVector{<:Real},
        f::Function, w::AbstractVector{<:Real}, m::Real, q::Real)
        T = promote_type(eltype(x), eltype(y), eltype(w))
        new{T}(Vector{T}(x), Vector{T}(y), f, Vector{T}(w), m, q)
    end
end

function Base.show(io::IO, ls::LinearLS)
    println(io, "Linear fit Y ~ m*f(X) + q")
    println(io, "\tf = ", ls.f)
    println(io, "\tm = ", ls.m)
    println(io, "\tq = ", ls.q)
end

"""
    linearls(x::AbstractVector, y::AbstractVector, f::Function=identity;
        weights=ones(length(x)))

Linear fit for `y = m*f(x) + q`.
Return a `LinearLS` object storing the estimated `m, q`
Optionally provide `weights` for Weighted LS

"""
function linearls(x::AbstractVector, y::AbstractVector, f::Function=identity;
        weights=ones(length(x)))
    X = diagm(weights) * [f.(x) ones(length(x))]
    Y = weights .* y
    m, q = X \ Y
    LinearLS(x, y, f, weights, m, q)
end

"""
    predict(ls::LinearLS, x::AbstractVector)

Predict new values based on the linear fit object `ls`. See `linearls`

"""
function predict(ls::LinearLS, x::AbstractVector)
    ls.f.(x) .* ls.m .+ ls.q
end


text:
using Sundials
using Base.Test

@testset "CVODE" begin
    @testset "Roberts CVODE Simplified" begin include("cvode_Roberts_simplified.jl") end
    @testset "Roberts CVODE Direct" begin include("cvode_Roberts_dns.jl") end
end
@testset "ARK" begin
    @testset "Roberts ARKODE Direct" begin include("arkode_Roberts_dns.jl") end
end
@testset "IDA" begin
    @testset "Roberts IDA Simplified" begin include("ida_Roberts_simplified.jl") end
    @testset "Roberts IDA Direct" begin include("ida_Roberts_dns.jl") end
    @testset "Heat IDA Direct" begin include("ida_Heat2D.jl") end
    #@testset "Cable IDA Direct" begin include("ida_Cable.jl") end # requires Grid
end
@testset "Kinsol" begin
    @testset "Kinsol Simplified" begin include("kinsol_mkin_simplified.jl") end
    @testset "Kinsol MKin" begin include("kinsol_mkinTest.jl") end
    @testset "Kinsol Banded" begin include("kinsol_banded.jl") end
end
@testset "Handle Tests" begin include("handle_tests.jl") end

@testset "Common Interface" begin
    @testset "CVODE" begin include("common_interface/cvode.jl") end
    @testset "ARKODE" begin include("common_interface/arkode.jl") end
    @testset "IDA" begin include("common_interface/ida.jl") end
    @testset "Jacobians" begin include("common_interface/jacobians.jl") end
    @testset "Callbacks" begin include("common_interface/callbacks.jl") end
    @testset "Iterator" begin include("common_interface/iterators.jl") end
    @testset "Errors" begin include("common_interface/errors.jl") end
end


text:
using Pickle.Torch
using Strided

const torch_tests = [
  "load",
  "load_zip",
  "save",
  "prebuild"
]

@testset "Torch" begin

  for t in torch_tests
    fp = joinpath(dirname(@__FILE__), "torch_$t.jl")
    @info "Test Torch $(uppercase(t))"
    include(fp)
  end

end


text:
nargin(func) = [method.nargs-1 for method in methods(func).ms]

# To treat lists and nonlists the same
struct List end
struct NonList end

islist(::Type{<:AbstractVector}) = List()
islist(::Type{<:Tuple}) = List()
islist(::Type{<:Number}) = NonList()
islist(::Type{<:Nothing}) = NonList()

aslist(x::T) where T = aslist(islist(T), x)
aslist(::List, x) = x
aslist(::NonList, x) = [x]

max_stars_repo_name: jonniedie/ComposableDiffEqs

max_issues_repo_name: jonniedie/ComposableDiffEqs

max_forks_repo_name: jonniedie/ComposableDiffEqs


text:
"""
    label_generators(generators::Array{T}; alphabet::Vector{Char}) -> labelled_generators::Bijection{String, T}

Labels the given `generators` using the given `alphabet`, a vector of characters and returns the labelled generators.
The `alphabet` can be bigger than the number of generators, only the first `Ω = length(generators)` characters are used.
"""
function label_generators(generators::Array{T}; alphabet = default_alphabet())::Bijection{String, T} where T
    Ω = length(generators)
    labelled_generators = Bijection(Dict(map(=>, alphabet[1:Ω] .* "", generators)))
    return labelled_generators
end

"""
    default_alphabet()

Defines a default list of characters that are used to label the generators. Maximum is currently 103 symbols.
The default alphabet consists of `'a':'z'`, `'A':'Z'`, `'𝒜':'𝒵'` and `'α':'ω'`.
"""
function default_alphabet()::Vector{Char}
    latin_lc = collect('a':'z')   # 26 symbols
    latin_uc = collect('A':'Z')  # 26 symbols
    scr_uc = collect('𝒜':'𝒵')   # 26 symbols
    greek_lc = collect('α':'ω')   # 25 symbols
    return union(latin_lc, latin_uc, scr_uc, greek_lc)
end

"""
    labelled_group_generator_simple(labelled_generators::Bijection{String, T}, num_max_elements::Int64; prnt = false) -> labelled_group::Bijection{String, T}

Naive Method: Generates all group elements and labels them with the first occuring label. It stops after generating `num_max_elements` elements. The resulting label may not be the shortest possible label.
The generators can be given to the function either as a labelled `Bijection{String, T}` or as an unlabelled `Array{T}` that will automatically labelled.

"""
function labelled_group_generator_simple(labelled_generators::Bijection{String, T}, num_max_elements::Int64; prnt = false)::Bijection{String, T} where T
    labelled_group = Bijection(copy(labelled_generators))
    labelled_generators_keys = sort(String.(keys(labelled_generators)))
    while length(labelled_group) < num_max_elements
        for one in labelled_generators_keys, two in sort(String.(keys(labelled_group)))
            # calculate element
            el = labelled_generators[one] * labelled_group[two]

            # create the new label
            el_label = one * two

            # create the Pair
            el = el_label => el

            # try to store the Pair in the Bijection, fails if key or value is already in Bijection
            try push!(labelled_group, el) catch end

            # break if all elements are labelled
            if length(labelled_group) == num_max_elements
                break
            end
        end
    end
    return labelled_group
end

function labelled_group_generator_simple(generators::Array{T}, num_max_elements::Int64; prnt = false)::Bijection{String, T} where T
    labelled_generators = label_generators(generators)
    return labelled_group_generator_simple(labelled_generators, num_max_elements; prnt)
end

"""
    labelled_group_generator_shortest(labelled_generators::Bijection{String, T}; prnt = false)::Tuple{Bijection{String, T}, Int64} -> labelled_group::Bijection{String, T}, num_all_mult::Int64

For the general method and description of `prnt` see `group_generator_basic()`.
This method generates all other elements of the closed group and labels them using the lexicographically first appearing (and shortest) label.
As Bijections are slower than Sets, this method can be really slow compared to `group_generator_basic()`.
"""
function labelled_group_generator_shortest(labelled_generators::Bijection{String, T}; prnt = false)::Tuple{Bijection{String, T}, Int64} where T

    # arrays to store generated elements in different levels of generation
    group_prev_level = Bijection{String, T}()
    group_current_level = Bijection(copy(labelled_generators))
    group_next_level = [1]

    level = 0
    size_after = 0
    new_elements = 0

    num_all_mult = 0

    labelled_generators_keys = sort(String.(keys(labelled_generators)))

    while length(group_current_level) > 0
        level += 1
        num_multiplications = 0

        # calculate the number of multiplications in this level
        num_max_mult = length(group_current_level) * length(labelled_generators)

        if prnt
            prog = Progress(num_max_mult)
        end

        # empty next for new iteration
        group_next_level = Bijection{String, T}()

        # multiplicate current with the generators and save new elements to next
        for one in labelled_generators_keys
            for two in sort(String.(keys(group_current_level)))
                num_multiplications += 1

                # calculate element
                el = labelled_generators[one] * group_current_level[two]
                # create the new label
                el_label = one * two

                # create the Pair
                el_pair = el_label => el

                # try to store the Pair in the Bijection, fails if key or value is already in Bijection
                try push!(group_next_level, el_pair) catch end

                if prnt
                    size_next_level = length(group_next_level)
                    update!(prog, num_multiplications; showvalues = [(:i,"$num_multiplications/$num_max_mult"),(:level, level),(:size,size_after),(:new_elements, new_elements),(:size_next_level, "$size_next_level ($(size_next_level/num_multiplications * 100)%)")])
                end
            end
        end

        if prnt print("current → prev...") end

        # current is now done, move to previous
        # union!() doesn't work for Bijections
        for (el_label, el) in group_current_level
            try push!(group_prev_level, el_label => el) catch end
        end

        if prnt println(" done.") end

        if prnt print("next → current...") end

        # move all from next to current that are not in previous
        # group_current_level = setdiff(group_next_level, group_prev_level)
        group_current_level = Bijection{String, T}()
        for (el_label, el) in group_next_level
            vals = values(group_prev_level)
            if !(el in vals)
                # element not in previous, store in current
                push!(group_current_level, el_label => el)
            end # element is in previous with shorter label
        end

        if prnt println(" done.") end

        # some stats to display
        size_after = length(group_prev_level)
        new_elements = length(group_current_level)
        num_all_mult += num_multiplications
        println("level #",level," size ", size_after, " and ", new_elements," new elements (",num_multiplications," multiplications)")
        if prnt
            println("level #",level," size ", size_after, " and ", new_elements," new elements (",num_multiplications," multiplications)")
            # countdown to prevent directly overwriting with progressbar
            print("5..."); sleep(1); print("4..."); sleep(1); print("3..."); sleep(1); print("2..."); sleep(1); println("1..."); sleep(1);
        end
    end

    println("\nEnded after ", level, " iterations. \nThe resulting group has ", size_after, " elements.")

    # returns the final group and the number of multiplications needed
    return group_prev_level, num_all_mult
end

function labelled_group_generator_shortest(generators::Array{T}; prnt = false)::Tuple{Bijection{String, T}, Int64} where T
    labelled_generators = label_generators(generators)
    return labelled_group_generator_shortest(labelled_generators; prnt)
end


text:
using XGBoost

D = 10
N = 200
NT = 100

x = float32(randn(N, D))
function rastrigin(x)
  D = length(x)
  10 * D + sum( x.^2 ) - 10 * sum( cos( 2 * π * x ) )
end
ya = mapslices(rastrigin, x, 2)
y = float32(vec(ya))
xx = float32(randn(NT, D))
yya = mapslices(rastrigin, xx, 2)
yy = float32(vec(yya))

bt = xgboost(x, 100, label = y, eta=1, max_depth=4, objective="reg:linear")
ypred = predict(bt, xx)
medape(ypred, y) = 100.0 * median(abs((ypred .- y) ./ y))
medape(ypred, yy)


text:
using Test
using ErrorFlow

@testset "@catch" begin
    err = @catch error("oh no")
    @test err.value == ErrorException("oh no")
    @test err.hints == []

    err = @catch error("oh no") "bad"
    @test err.value == ErrorException("oh no")
    @test err.hints == ["bad"]
end

@testset "@chain example" begin
    msg_fail_load_config = "Failed to load config"
    msg_fail_parse_config = "Failed to parse config"

    app(succeed) = @chain load_config(succeed) msg_fail_load_config
    load_config(succeed) = @chain parse_line157(succeed) msg_fail_parse_config
    function parse_line157(succeed)
        if succeed
            return 42
        else
            error("Cannot parse line 157")
        end
    end

    @test app(true) == 42
    @inferred Union{Int, Error} app(true)
    err = app(false)
    @test err.value == ErrorException("Cannot parse line 157")
    s = sprint(showerror, err)
    @test occursin(msg_fail_load_config, s)
    @test occursin(msg_fail_parse_config, s)
    @test occursin("parse_line157", s)
end


text:
## sampler ###################################################################

# allows to call `Sampler` only when the the arg isn't a Sampler itself

sampler(::Type{RNG}, X, n::Repetition) where {RNG<:AbstractRNG} =
    Sampler(RNG, X, n)

sampler(::Type{RNG}, ::Type{X}, n::Repetition) where {RNG<:AbstractRNG,X} =
    Sampler(RNG, X, n)

sampler(::Type{RNG}, X::Sampler, n::Repetition) where {RNG<:AbstractRNG} = X

sampler(rng::AbstractRNG, X, n::Repetition) = sampler(typeof(rng), X, n)

## SamplerTrivial, SamplerSimple, SamplerTag #################################

# we redefine (copy-paste) SamplerTrivial, SamplerSimple, SamplerTag from
# Random so that they inherit from SamplerReset

struct SamplerTrivial{T,E} <: SamplerReset{E}
    self::T
end

"""
    SamplerTrivial(x)

Create a sampler that just wraps the given value `x`. This is the default fall-back for
values.
The `eltype` of this sampler is equal to `eltype(x)`.

The recommended use case is sampling from values without precomputed data.
"""
SamplerTrivial(x::T) where {T} = SamplerTrivial{T,gentype(T)}(x)

Sampler(::Type{<:AbstractRNG}, x::Distribution, ::Repetition) = SamplerTrivial(x)

Base.getindex(sp::SamplerTrivial) = sp.self

# simple sampler carrying data (which can be anything)
struct SamplerSimple{T,S,E} <: SamplerReset{E}
    self::T
    data::S
end

"""
    SamplerSimple(x, data)

Create a sampler that wraps the given value `x` and the `data`.
The `eltype` of this sampler is equal to `eltype(x)`.

The recommended use case is sampling from values with precomputed data.
"""
SamplerSimple(x::T, data::S) where {T,S} = SamplerSimple{T,S,gentype(T)}(x, data)

Base.getindex(sp::SamplerSimple) = sp.self

# simple sampler carrying a (type) tag T and data
struct SamplerTag{T,S,E} <: SamplerReset{E}
    data::S
    SamplerTag{T}(s::S) where {T,S} = new{T,S,gentype(T)}(s)
end