text:
# common facilities

# scalar functions
"""
    xlogx(x::Real)

Return `x * log(x)` for `x ≥ 0`, handling `x = 0` by taking the downward limit.

```jldoctest
julia> StatsFuns.xlogx(0)
0.0
```
"""
xlogx(x::Real) = x > zero(x) ? x * log(x) : zero(log(x))

"""
    xlogy(x::Real, y::Real)

Return `x * log(y)` for `y > 0` with correct limit at `x = 0`.
"""
xlogy(x::T, y::T) where {T<:Real} = x > zero(T) ? x * log(y) : zero(log(x))
xlogy(x::Real, y::Real) = xlogy(promote(x, y)...)

"""
    logistic(x::Real)

The [logistic](https://en.wikipedia.org/wiki/Logistic_function) sigmoid function mapping a real number to a value in the interval [0,1],
```math
\\sigma(x) = \\frac{1}{e^{-x} + 1} = \\frac{e^x}{1+e^x}.
```

Its inverse is the [`logit`](@ref) function.
"""
logistic(x::Real) = inv(exp(-x) + one(x))

"""
    logit(p::Real)

The [logit](https://en.wikipedia.org/wiki/Logit) or log-odds transformation,
```math
\\log\\left(\\frac{x}{1-x}\\right), \\text{where} 0 < x < 1
```
Its inverse is the [`logistic`](@ref) function.
"""
logit(x::Real) = log(x / (one(x) - x))

"""
    log1psq(x::Real)

Return `log(1+x^2)` evaluated carefully for `abs(x)` very small or very large.
"""
log1psq(x::Real) = log1p(abs2(x))
function log1psq(x::Union{Float32,Float64})
    ax = abs(x)
    ax < maxintfloat(x) ? log1p(abs2(ax)) : 2 * log(ax)
end

"""
    log1pexp(x::Real)

Return `log(1+exp(x))` evaluated carefully for largish `x`.

This is also called the ["softplus"](https://en.wikipedia.org/wiki/Rectifier_(neural_networks))
transformation, being a smooth approximation to `max(0,x)`. Its inverse is [`logexpm1`](@ref).
"""
log1pexp(x::Real) = x < 18.0 ? log1p(exp(x)) : x < 33.3 ? x + exp(-x) : oftype(exp(-x), x)
log1pexp(x::Float32) = x < 9.0f0 ? log1p(exp(x)) : x < 16.0f0 ? x + exp(-x) : oftype(exp(-x), x)

"""
    log1mexp(x::Real)

Return `log(1 - exp(x))`

See:
 * Martin Maechler (2012) "Accurately Computing log(1 − exp(− |a|))",
   http://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf

Note: different than Maechler (2012), no negation inside parentheses
"""
log1mexp(x::Real) = x < loghalf ? log1p(-exp(x)) : log(-expm1(x))

"""
    log2mexp(x::Real)

Return `log(2 - exp(x))` evaluated as `log1p(-expm1(x))`
"""
log2mexp(x::Real) = log1p(-expm1(x))

"""
    logexpm1(x::Real)

Return `log(exp(x) - 1)` or the "invsoftplus" function.
It is the inverse of [`log1pexp`](@ref) (aka "softplus").
"""
logexpm1(x::Real) = x <= 18.0 ? log(expm1(x)) : x <= 33.3 ? x - exp(-x) : oftype(exp(-x), x)
logexpm1(x::Float32) = x <= 9f0 ? log(expm1(x)) : x <= 16f0 ? x - exp(-x) : oftype(exp(-x), x)

const softplus = log1pexp
const invsoftplus = logexpm1

"""
    log1pmx(x::Float64)

Return `log(1 + x) - x`

Use naive calculation or range reduction outside kernel range.  Accurate ~2ulps for all `x`.
"""
function log1pmx(x::Float64)
    if !(-0.7 < x < 0.9)
        return log1p(x) - x
    elseif x > 0.315
        u = (x-0.5)/1.5
        return _log1pmx_ker(u) - 9.45348918918356180e-2 - 0.5*u
    elseif x > -0.227
        return _log1pmx_ker(x)
    elseif x > -0.4
        u = (x+0.25)/0.75
        return _log1pmx_ker(u) - 3.76820724517809274e-2 + 0.25*u
    elseif x > -0.6
        u = (x+0.5)*2.0
        return _log1pmx_ker(u) - 1.93147180559945309e-1 + 0.5*u
    else
        u = (x+0.625)/0.375
        return _log1pmx_ker(u) - 3.55829253011726237e-1 + 0.625*u
    end
end

"""
    logmxp1(x::Float64)

Return `log(x) - x + 1` carefully evaluated.
"""
function logmxp1(x::Float64)
    if x <= 0.3
        return (log(x) + 1.0) - x
    elseif x <= 0.4
        u = (x-0.375)/0.375
        return _log1pmx_ker(u) - 3.55829253011726237e-1 + 0.625*u
    elseif x <= 0.6
        u = 2.0*(x-0.5)
        return _log1pmx_ker(u) - 1.93147180559945309e-1 + 0.5*u
    else
        return log1pmx(x - 1.0)
    end
end

# The kernel of log1pmx
# Accuracy within ~2ulps for -0.227 < x < 0.315
function _log1pmx_ker(x::Float64)
    r = x/(x+2.0)
    t = r*r
    w = @horner(t,
                6.66666666666666667e-1, # 2/3
                4.00000000000000000e-1, # 2/5
                2.85714285714285714e-1, # 2/7
                2.22222222222222222e-1, # 2/9
                1.81818181818181818e-1, # 2/11
                1.53846153846153846e-1, # 2/13
                1.33333333333333333e-1, # 2/15
                1.17647058823529412e-1) # 2/17
    hxsq = 0.5*x*x
    r*(hxsq+w*t)-hxsq
end

"""
    logaddexp(x::Real, y::Real)

Return `log(exp(x) + exp(y))`, avoiding intermediate overflow/undeflow, and handling non-finite values.
"""
function logaddexp(x::T, y::T) where T<:Real
    # x or y is  NaN  =>  NaN
    # x or y is +Inf  => +Inf
    # x or y is -Inf  => other value
    isfinite(x) && isfinite(y) || return max(x,y)
    x > y ? x + log1p(exp(y - x)) : y + log1p(exp(x - y))
end
logaddexp(x::Real, y::Real) = logaddexp(promote(x, y)...)

Base.@deprecate logsumexp(x::Real, y::Real) logaddexp(x,y)

"""
    logsumexp(X)

Compute `log(sum(exp, X))`, evaluated avoiding intermediate overflow/undeflow.

`X` should be an iterator of real numbers.
"""
function logsumexp(X)
    isempty(X) && return log(sum(X))
    reduce(logaddexp, X)
end
function logsumexp(X::AbstractArray{T}; dims=:) where {T<:Real}
    u = reduce(max, X, dims=dims, init=log(zero(T)))
    u isa AbstractArray || isfinite(u) || return float(u)
    let u=u # avoid https://github.com/JuliaLang/julia/issues/15276
        # TODO: remove the branch when JuliaLang/julia#31020 is merged.
        if u isa AbstractArray
            u .+ log.(sum(exp.(X .- u); dims=dims))
        else
            u + log(sum(x -> exp(x-u), X))
        end
    end
end

"""
    softmax!(r::AbstractArray, x::AbstractArray)

Overwrite `r` with the `softmax` (or _normalized exponential_) transformation of `x`

That is, `r` is overwritten with `exp.(x)`, normalized to sum to 1.

See the [Wikipedia entry](https://en.wikipedia.org/wiki/Softmax_function)
"""
function softmax!(r::AbstractArray{R}, x::AbstractArray{T}) where {R<:AbstractFloat,T<:Real}
    n = length(x)
    length(r) == n || throw(DimensionMismatch("Inconsistent array lengths."))
    u = maximum(x)
    s = 0.
    @inbounds for i = 1:n
        s += (r[i] = exp(x[i] - u))
    end
    invs = convert(R, inv(s))
    @inbounds for i = 1:n
        r[i] *= invs
    end
    r
end

"""
    softmax(x::AbstractArray{<:Real})

Return the [`softmax transformation`](https://en.wikipedia.org/wiki/Softmax_function) applied to `x`
"""
softmax!(x::AbstractArray{<:AbstractFloat}) = softmax!(x, x)
softmax(x::AbstractArray{<:Real}) = softmax!(similar(x, Float64), x)


text:
using Plots
x = 1:10; y = rand(10); # These are the plotting data
gr()
display(plot(x,y))
savefig("myplot.png")

ext: jl

max_stars_repo_name: sergev/vak-opensource

max_issues_repo_name: sergev/vak-opensource

max_forks_repo_name: sergev/vak-opensource


text:
using Stheno: EmpiricalMean, EmpiricalKernel
using FillArrays: Zeros, Fill

@testset "conversion" begin

    # Mean function conversion.
    let
        rng, N = MersenneTwister(123456), 11
        x = randn(rng, N)

        # Check that base mean functions convert properly when made finite.
        for (μ, T) in [(CustomMean(sin), Vector{Float64}),
                       (ZeroMean{Float64}(), Zeros{Float64}),
                       (ConstantMean(randn(rng)), Fill)]
            μ_finite = FiniteMean(μ, x)
            @test AbstractVector(μ_finite) isa T
            @test AbstractVector(μ_finite) == map(μ, x)
        end

        # Test `EmpiricalMean` converts correctly.
        @test AbstractVector(EmpiricalMean(x)) === x

        # Check that `BlockMean`s behave well in both shallow finite scenarios.
        x1, x2 = randn(rng, 11), randn(rng, 7)
        μ1, μ2 = ConstantMean(randn(rng)), ZeroMean{Float64}()
        μ1f, μ2f = FiniteMean(μ1, x1), FiniteMean(μ2, x2)
        bf_μ = FiniteMean(BlockMean([μ1, μ2]), BlockData([x1, x2]))
        b_μf = BlockMean([μ1f, μ2f])

        @test AbstractVector(bf_μ) isa BlockVector
        @test AbstractVector(bf_μ) == vcat(AbstractVector(μ1f), AbstractVector(μ2f))

        @test AbstractVector(b_μf) isa BlockVector
        @test AbstractVector(bf_μ) == vcat(AbstractVector(μ1f), AbstractVector(μ2f))
    end

    # Kernel conversion.
    let
        rng, N = MersenneTwister(123456), 11
        x = randn(rng, N)

        # Test that base kernel functions convert properly when made finite.
        for (k, T) in [(ZeroKernel{Float64}(), Zeros),
                       (ConstantKernel(5.0), Fill),
                       (EQ(), Matrix{Float64}),
                       (Linear(0.0), Matrix{Float64}),
                       (Noise(1.0), Diagonal)]
            k_finite = FiniteKernel(k, x)
            @test AbstractMatrix(k_finite) isa LazyPDMat{Float64, <:T}
            @test AbstractMatrix(k_finite) == pairwise(k, x)
        end

        A = randn(5, 5)
        Σ = LazyPDMat(A' * A + eps() * I)
        @test AbstractMatrix(EmpiricalKernel(Σ)) === Σ

        # Check that `BlockKernel`s play nicely when made finite.
        x1, x2 = randn(rng, 11), randn(rng, 7)
        k1, k2, k12 = EQ(), Linear(1.0), ZeroKernel{Float64}()
        ks_off = Matrix{Stheno.CrossKernel}(undef, 2, 2)
        ks_off[1, 2] = k12
        bk = BlockKernel([k1, k2], ks_off)
        bf_k = FiniteKernel(bk, BlockData([x1, x2]))

        @test AbstractMatrix(bf_k) isa LazyPDMat{Float64, <:Symmetric{Float64, <:BlockMatrix}}
        @test AbstractMatrix(bf_k) == pairwise(bk, BlockData([x1, x2]))
    end

    # CrossKernel conversions.
    let
        # Check that CrossKernels play nicely when made finite.
        rng = MersenneTwister(123456)
        x1, x2, x3, x4 = randn(rng, 11), randn(rng, 7), randn(rng, 3), randn(rng, 13)
        k1, k2, k12, k21 = EQ(), EQ(), EQ(), EQ()
        bk = BlockCrossKernel([k1 k12; k21 k2])
        bf_k = FiniteCrossKernel(bk, BlockData([x1, x2]), BlockData([x3, x4]))

        @test AbstractMatrix(bf_k) isa BlockMatrix
        @test AbstractMatrix(bf_k) == pairwise(bk, BlockData([x1, x2]), BlockData([x3, x4]))
    end
end


text:
###################################################
## Save and load networks
###################################################

"""
    load a network previously saved as a template
"""
function loadTemplate(name::AbstractString)
    path = "$(Pkg.dir("RoutingNetworks"))/saved/$(name).jld"
    if isfile(path)
        return JLD.load(path, "n")
    else
        error("template does not exist")
    end
end

"""
    Save a network as a template
"""
function saveTemplate(n::Network, name::AbstractString)
    path = "$(Pkg.dir("RoutingNetworks"))/saved/$(name).jld"
    JLD.save(path, "n", n)
end

# """
#     Create some common templates
# """
# function createTemplates()
#     #Manhattan
#     if !isTemplate("manhattan")
#         MANHATTAN_BORDER = [(40.87851017592601,-73.927903175354);(40.87743937506905,-73.92168045043945);(40.87516792197097,-73.91871929168701);(40.87393481479329,-73.91193866729736);(40.87182549927143,-73.90970706939697);(40.866016421491814,-73.91168117523193);(40.859655059077824,-73.91764640808105);(40.84855366761019,-73.9266586303711);(40.84394377141786,-73.92949104309082);(40.83550228531863,-73.9339542388916);(40.82816381252365,-73.9339542388916);(40.81972031701224,-73.93318176269531);(40.808677200957106,-73.93335342407227);(40.801920506109774,-73.92803192138672);(40.7958778790764,-73.92820358276367);(40.78210123234386,-73.94107818603516);(40.77547182731999,-73.93953323364258);(40.7701418259051,-73.94588470458984);(40.75831029512206,-73.95669937133789);(40.744395800976775,-73.96888732910156);(40.71746884168838,-73.97077560424805);(40.70927151739564,-73.9764404296875);(40.7066689811733,-73.99497985839844);(40.697299008636755,-74.01008605957031);(40.69899090674371,-74.02381896972656);(40.74907763805906,-74.01351928710938);(40.77950154452169,-73.99532318115234);(40.81822635589172,-73.96751403808594);(40.853163243121564,-73.94845962524414)]
#         n = queryOsmPolygon(MANHATTAN_BORDER)
#         saveTemplate(n, "manhattan")
#     end
# end

"""
    Check if template exists
"""
function isTemplate(name::AbstractString)
    path = "$(Pkg.dir("RoutingNetworks"))/saved/$(name).jld"
    return isfile(path)
end

max_stars_repo_name: andrew-zhm/RoutingNetworks

max_issues_repo_name: andrew-zhm/RoutingNetworks

max_forks_repo_name: andrew-zhm/RoutingNetworks


text:
module ClearStacktrace

using Crayons

const MODULECRAYONS = Ref(
        [
            crayon"blue",
            crayon"yellow",
            crayon"red",
            crayon"green",
            crayon"cyan",
            crayon"magenta",
        ]
)
const TYPECOLORS = Ref(Any[0x706060, 0x607060, 0x606080])
const CRAYON_HEAD = Ref(Crayon(bold = true))
const CRAYON_HEADSEP = Ref(Crayon(bold = true))
const CRAYON_FUNCTION = Ref(Crayon())
const CRAYON_FUNC_EXT = Ref(Crayon(foreground = :dark_gray))
const CRAYON_LOCATION = Ref(Crayon(foreground = :dark_gray))
const CRAYON_NUMBER = Ref(Crayon(foreground = :blue))
const NUMPAD = Ref(1)
const FUNCPAD = Ref(2)
const MODULEPAD = Ref(2)

function expandbasepath(str)
    if !isnothing(match(r"^\./\w+\.jl$", str))
        sourcestring = Base.find_source_file(str[3:end]) # cut off ./
        replace(sourcestring, raw"\\" => '/')
    else
        str
    end
end

function replaceuserpath(str)
    str1 = replace(str, homedir() => "~")
    # seems to be necessary for some paths with small letter drive c:// etc
    replace(str1, lowercasefirst(homedir()) => "~")
end

getline(frame) = frame.line
getfile(frame) = string(frame.file) |> expandbasepath |> replaceuserpath
getfunc(frame) = string(frame.func)
getmodule(frame) = try; string(frame.linfo.def.module) catch; "" end

getsig(frame) = try;  join("::" .* repr.(frame.linfo.specTypes.parameters[2:end]), ", ") catch; "" end

function convert_trace(trace)
    files = getfile.(trace)
    lines = getline.(trace)
    funcs = getfunc.(trace)
    moduls = getmodule.(trace)
    signatures = getsig.(trace)
    inlineds = getfield.(trace, :inlined)

    # replace empty modules if there is another frame from the same file
    for (i_this, mo) in enumerate(moduls)
        if mo == ""
            for i_other in 1:length(trace)
                if files[i_this] == files[i_other] && moduls[i_other] != ""
                    moduls[i_this] = moduls[i_other]
                end
            end
        end
    end

    (files = files, lines = lines, funcs = funcs,
        moduls = moduls, signatures = signatures, inlineds = inlineds)
end

function printtrace(io::IO, stacktrace)

    files, lines, funcs, moduls, signatures, inlineds = convert_trace(stacktrace)

    numbers = ["[" * string(i) * "]" for i in 1:length(stacktrace)]
    numwidth = maximum(length, numbers)

    exts = [inl ? " [i]" : "" for inl in inlineds]
    funcs_w_ext = funcs .* exts
    funcwidth = maximum(length, funcs_w_ext)

    modulwidth = max(maximum(length, moduls), length("Module"))

    umoduls = setdiff(unique(moduls), [""])
    modcrayons = Dict(u => c for (u, c) in
        Iterators.zip(umoduls, Iterators.cycle(MODULECRAYONS[])))

    print(io, rpad("", numwidth + NUMPAD[]))
    print(io, CRAYON_HEAD[](rpad("Function", funcwidth + FUNCPAD[])))
    print(io, CRAYON_HEAD[](rpad("Module", modulwidth + MODULEPAD[])))
    print(io, CRAYON_HEAD[]("Signature"))
    println(io)

    print(io, rpad("", numwidth + NUMPAD[]))
    print(io, CRAYON_HEADSEP[](rpad("────────", funcwidth + FUNCPAD[])))
    print(io, CRAYON_HEADSEP[](rpad("──────", modulwidth + MODULEPAD[])))
    print(io, CRAYON_HEADSEP[]("─────────"))
    println(io)

    for (i, (num, func, ext, modul, file, line, signature)) in enumerate(
            zip(numbers, funcs, exts, moduls, files, lines, signatures))
        print(io, CRAYON_NUMBER[](rpad(num, numwidth + NUMPAD[])))

        print(io, CRAYON_FUNCTION[](func))

        print(io, CRAYON_FUNC_EXT[](ext * (" " ^ (FUNCPAD[] + funcwidth - length(funcs_w_ext[i])))))

        mcrayon = get(modcrayons, modul, crayon"white")
        print(io, mcrayon(rpad(modul, modulwidth + MODULEPAD[])))

        print_signature(io, signature, TYPECOLORS[])

        println(io)

        println(io, CRAYON_LOCATION[](string(file) * ":" * string(line)))
    end
end

function print_signature(io::IO, sig, colors)
    # split before and after curly braces and commas

    if length(colors) != 3
        error("""
        Three colors are needed to color a type signature without collisions.
        You supplied $(length(colors)): $colors
        """)
    end

    if isempty(sig)
        return
    end

    regex = r"(?<=[\{\}\,])|(?=[\{\}\,])"

    parts = split(sig, regex)

    # pretend that we used two colors already to avoid the edge case
    colorstack = colors[1:2]

    print(io, Crayon(foreground = :dark_gray)("("))
    for p in parts
        color = nothing

        if p == "{"
            # opening brace has same color as previous word
            color = colorstack[end]
            # now add the color of the word before that to the end of the stack
            # so the next chosen color will be neither of these two
            push!(colorstack, colorstack[end-1])
        elseif p == "}"
            # closing brace gets the color from the previous stack level
            # remove that one
            pop!(colorstack)
            color = colorstack[end]
        elseif p == ","
            color = :dark_gray
        else
            # change color for every element, choose not the last two ones
            color = setdiff(colors, colorstack[end-1:end])[1]
            # then change the last color to the current ones
            colorstack[end] = color
        end

        cray = Crayon(foreground = color)
        print(io, cray(p))
    end
    print(io, Crayon(foreground = :dark_gray)(")"))
end

@warn "Overloading Base.show_backtrace(io::IO, t::Vector) with custom version"
function Base.show_backtrace(io::IO, t::Vector)

    ### this part is copied from the original function
    resize!(Base.LAST_SHOWN_LINE_INFOS, 0)
    filtered = Base.process_backtrace(t)
    isempty(filtered) && return

    if length(filtered) == 1 && StackTraces.is_top_level_frame(filtered[1][1])
        f = filtered[1][1]
        if f.line == 0 && f.file == Symbol("")
            # don't show a single top-level frame with no location info
            return
        end
    end
    ###

    println(io, "\nStacktrace:")

    # process_backtrace returns a Tuple{Frame, Int}
    frames = first.(filtered)

    printtrace(io, frames)
end

end # module


text:
using STLFile
using Test

@testset "STLFile.jl" begin
    # Write your tests here.
end


text:
using Dates

mutable struct SwapRateHelper{PrT <: Dates.Period, PrS <: Dates.Period} <: RateHelper
  rate::Quote
  tenor::PrT
  fwdStart::PrS
  swap::VanillaSwap
end

function SwapRateHelper(rate::Float64, tenor::PrT, cal::C, fixedFrequency::F, fixedConvention::B, fixedDayCount::DC, iborIndex::IborIndex, spread::Float64, fwdStart::PrS,
                    pricingEngine::P = DiscountingSwapEngine(), settlementDays::Int = iborIndex.fixingDays, nominal::Float64 = 1.0, swapT::ST = Payer(),
                    fixedRate::Float64 = 0.0) where {PrT <: Dates.Period, C <: BusinessCalendar, F <: Frequency, B <: BusinessDayConvention, DC <: DayCount, PrS <: Dates.Period, P <: PricingEngine, ST <: SwapType}
  # do stuff
  fixedCal = cal
  floatingCal = cal
  floatTenor = iborIndex.tenor
  fixedTenor = QuantLib.Time.TenorPeriod(fixedFrequency)
  fixedTermConvention = fixedConvention
  floatConvention = iborIndex.convention
  floatTermConvention = iborIndex.convention
  fixedRule = DateGenerationBackwards()
  floatRule = DateGenerationBackwards()
  floatDayCount = iborIndex.dc
  # fixed_rate = 0.0

  ref_date = adjust(floatingCal, floatConvention, settings.evaluation_date)
  spot_date = advance(Dates.Day(settlementDays), floatingCal, ref_date, floatConvention)
  start_date = adjust(floatingCal, floatConvention, spot_date + fwdStart)
  ## TODO Float end of month (defaults to false)
  end_date = start_date + tenor

  # build schedules
  fixed_schedule = Schedule(start_date, end_date, fixedTenor, fixedConvention, fixedTermConvention, fixedRule, false, fixedCal)
  float_schedule = Schedule(start_date, end_date, floatTenor, floatConvention, floatTermConvention, floatRule, false, floatingCal)

  swap = VanillaSwap(swapT, nominal, fixed_schedule, fixedRate, fixedDayCount, iborIndex, spread, float_schedule, floatDayCount, pricingEngine, fixedConvention)

  return SwapRateHelper{PrT, PrS}(Quote(rate), tenor, fwdStart, swap)
end

maturity_date(sh::SwapRateHelper) = maturity_date(sh.swap)

struct DepositRateHelper{B <: BusinessCalendar, C <: BusinessDayConvention, DC <: DayCount} <: RateHelper
  rate::Quote
  tenor::TenorPeriod
  fixingDays::Int
  calendar::B
  convention::C
  endOfMonth::Bool
  dc::DC
  iborIndex::IborIndex
  evaluationDate::Date
  referenceDate::Date
  earliestDate::Date
  maturityDate::Date
  fixingDate::Date
end

function DepositRateHelper(rate::Quote, tenor::TenorPeriod, fixingDays::Integer, calendar::B, convention::C, endOfMonth::Bool, dc::DC) where {B <: BusinessCalendar, C <: BusinessDayConvention, DC <: DayCount}
  ibor_index = IborIndex("no-fix", tenor, fixingDays, NullCurrency(), calendar, convention, endOfMonth, dc)
  evaluation_date = settings.evaluation_date
  reference_date = adjust(calendar, convention, evaluation_date)
  earliest_date = advance(Dates.Day(fixingDays), calendar, reference_date, convention)
  maturity_d = maturity_date(ibor_index, earliest_date)
  fix_date = fixing_date(ibor_index, earliest_date)
  return DepositRateHelper{B, C, DC}(rate, tenor, fixingDays, calendar, convention, endOfMonth, dc, ibor_index, evaluation_date, reference_date, earliest_date, maturity_d, fix_date)
end

maturity_date(rate::RateHelper) = rate.maturityDate

value(rate::RateHelper) = rate.rate.value

function implied_quote(swap_helper::SwapRateHelper)
  swap = swap_helper.swap
  recalculate!(swap)
  #
  # println("Floating Leg NPV ", floating_leg_NPV(swap))
  # println("Floating Leg BPS ", floating_leg_BPS(swap))
  # println("Fixed Leg BPS ", fixed_leg_BPS(swap))
  # println("Swap spread ", swap.spread)

  floatingLegNPV = floating_leg_NPV(swap)
  spread = swap.spread
  spreadNPV = floating_leg_BPS(swap) / basisPoint * spread
  totNPV = -(floatingLegNPV + spreadNPV)

  return totNPV / (fixed_leg_BPS(swap) / basisPoint)
end

function implied_quote(rh::RateHelper)
  return fixing(rh.iborIndex, rh.iborIndex.ts, rh.fixingDate, true)
end

struct FraRateHelper{D <: Dates.Period, II <: IborIndex} <: RateHelper
  rate::Quote
  evaluationDate::Date
  periodToStart::D
  iborIndex::II
  fixingDate::Date
  earliestDate::Date
  latestDate::Date
end

function FraRateHelper(rate::Quote, monthsToStart::Int, monthsToEnd::Int, fixingDays::Int, calendar::BusinessCalendar, convention::BusinessDayConvention, endOfMonth::Bool, dc::DayCount)
  periodToStart = Dates.Month(monthsToStart)
  iborIndex = IborIndex("no-fix", TenorPeriod(Dates.Month(monthsToEnd - monthsToStart)), fixingDays, NullCurrency(), calendar, convention, endOfMonth, dc)
  evaluationDate = settings.evaluation_date

  # initialize dates
  refDate = adjust(iborIndex.fixingCalendar, evaluationDate)
  spotDate = advance(Dates.Day(iborIndex.fixingDays), calendar, refDate)
  earliestDate = advance(periodToStart, calendar, spotDate, convention)
  latestDate = maturity_date(iborIndex, earliestDate)
  fixingDate = fixing_date(iborIndex, earliestDate)

  return FraRateHelper{Dates.Month, typeof(iborIndex)}(rate, evaluationDate, periodToStart, iborIndex, fixingDate, earliestDate, latestDate)
end

maturity_date(fra::FraRateHelper) = fra.latestDate

# Clone functions #
function clone(depo::DepositRateHelper, ts::TermStructure = depo.iborIndex.ts)
  # first we have to clone a new index
  newIdx = clone(depo.iborIndex, ts)
  #now we build a new depo helper
  return DepositRateHelper(depo.rate, depo.tenor, depo.fixingDays, depo.calendar, depo.convention, depo.endOfMonth, depo.dc, newIdx, depo.evaluationDate, depo.referenceDate,
                          depo.earliestDate, depo.maturityDate, depo.fixingDate)
end

function clone(fra::FraRateHelper, ts::TermStructure = fra.iborIndex.ts)
  # clone new index
  newIdx = clone(fra.iborIndex, ts)

  # build new fra helper
  return FraRateHelper(fra.rate, fra.evaluationDate, fra.periodToStart, newIdx, fra.fixingDate, fra.earliestDate, fra.latestDate)
end

function clone(swapHelper::SwapRateHelper, ts::TermStructure = swapHelper.swap.iborIndex.ts)
  # first we need a new PE for the swap
  newPE = clone(swapHelper.swap.pricingEngine, ts)
  # then we have to clone the swap
  newSwap = clone(swapHelper.swap, newPE, ts)

  # now we can clone the helper
  return SwapRateHelper(swapHelper.rate, swapHelper.tenor, swapHelper.fwdStart, newSwap)
end

update_termstructure(rateHelper::RateHelper, ts::TermStructure) = clone(rateHelper, ts)


text:
# Unit tests for JUDI linear operators (without PDE solves)
# Philipp Witte (pwitte.slim@gmail.com)
# May 2018
#
# Mathias Louboutin, mlouboutin3@gatech.edu
# Updated July 2020

# Tests
function test_transpose(Op)
    @test isequal(size(Op), size(conj(Op)))
    @test isequal(reverse(size(Op)), size(transpose(Op)))
    @test isequal(reverse(size(Op)), size(adjoint(Op)))
    @test isequal(reverse(size(Op)), size(transpose(Op)))
    return true
end

function test_getindex(Op, nsrc)
    # requires: Op.info.nsrc == 2
    Op_sub = Op[1]
    @test isequal(Op_sub.info.nsrc, 1)
    @test isequal(Op_sub.model, Op.model)
    @test isequal(size(Op_sub), convert(Tuple{Int64, Int64}, size(Op) ./ nsrc))

    inds = nsrc > 1 ? (1:nsrc) : 1
    Op_sub = Op[inds]
    @test isequal(Op_sub.info.nsrc, nsrc)
    @test isequal(Op_sub.model, Op.model)
    @test isequal(size(Op_sub), size(Op))
    return true
end

########################################################### judiModeling ###############################################

@testset "judiModeling Unit Test with $(nsrc) sources" for nsrc=[1, 2]

    info = example_info(nsrc=nsrc)
    model = example_model()
    F_forward = judiModeling(info, model; options=Options())
    F_adjoint = judiModelingAdjoint(info, model; options=Options())

    @test isequal(typeof(F_forward), judiModeling{Float32, Float32})
    @test isequal(typeof(F_adjoint), judiModelingAdjoint{Float32, Float32})

    # conj, transpose, adjoint
    @test test_transpose(F_forward)
    @test test_transpose(F_adjoint)

    # get index
    @test test_getindex(F_forward, nsrc)
    @test test_getindex(F_adjoint, nsrc)

    if VERSION>v"1.2"
        a = randn(Float32, model.n...)
        F2 = F_forward(;m=a)
        @test isapprox(F2.model.m, a)
        F2 = F_forward(Model(model.n, model.d, model.o, a))
        @test isapprox(F2.model.m, a)
        @test F2.model.n == model.n
    end

end

############################################################# judiPDE ##################################################

@testset "judiPDE Unit Test with $(nsrc) sources" for nsrc=[1, 2]

    info = example_info(nsrc=nsrc)
    model = example_model()
    rec_geometry = example_rec_geometry(nsrc=nsrc)

    PDE_forward = judiPDE("PDE", info, model, rec_geometry; options=Options())
    PDE_adjoint = judiPDEadjoint("PDEadjoint", info, model, rec_geometry; options=Options())

    @test isequal(typeof(PDE_forward), judiPDE{Float32, Float32})
    @test isequal(typeof(PDE_adjoint), judiPDEadjoint{Float32, Float32})

    # conj, transpose, adjoint
    @test test_transpose(PDE_forward)
    @test test_transpose(PDE_adjoint)

    # get index
    @test test_getindex(PDE_forward, nsrc)
    @test test_getindex(PDE_adjoint, nsrc)

    # Multiplication w/ judiProjection
    src_geometry = example_src_geometry()
    Ps = judiProjection(info, src_geometry)

    PDE = PDE_forward*transpose(Ps)
    @test isequal(typeof(PDE), judiPDEfull{Float32, Float32})
    @test isequal(PDE.recGeometry, rec_geometry)
    @test isequal(PDE.srcGeometry, src_geometry)

    PDEad = PDE_adjoint*transpose(Ps)
    @test isequal(typeof(PDEad), judiPDEfull{Float32, Float32})
    @test isequal(PDEad.srcGeometry, rec_geometry)
    @test isequal(PDEad.recGeometry, src_geometry)

end

######################################################### judiPDEfull ##################################################

@testset "judiPDEfull Unit Test with $(nsrc) sources" for nsrc=[1, 2]

    info = example_info(nsrc=nsrc)
    model = example_model()
    rec_geometry = example_rec_geometry(nsrc=nsrc)
    src_geometry = example_src_geometry(nsrc=nsrc)
    PDE = judiModeling(info, model, src_geometry, rec_geometry; options=Options())

    @test isequal(typeof(PDE), judiPDEfull{Float32, Float32})
    @test isequal(PDE.recGeometry, rec_geometry)
    @test isequal(PDE.srcGeometry, src_geometry)

    @test test_transpose(PDE)
    @test test_getindex(PDE, nsrc)

end

######################################################## judiJacobian ##################################################

@testset "judiJacobian Unit Test with $(nsrc) sources" for nsrc=[1, 2]

    info = example_info(nsrc=nsrc)
    model = example_model()
    rec_geometry = example_rec_geometry(nsrc=nsrc)
    src_geometry = example_src_geometry(nsrc=nsrc)
    wavelet = randn(Float32, src_geometry.nt[1])
    PDE = judiModeling(info, model, src_geometry, rec_geometry; options=Options())
    q = judiVector(src_geometry, wavelet)

    J = judiJacobian(PDE, q)

    @test isequal(typeof(J), judiJacobian{Float32, Float32})
    @test isequal(J.recGeometry, rec_geometry)
    @test isequal(J.source.geometry, src_geometry)
    @test all(isequal(J.source.data[i], wavelet) for i=1:nsrc)
    @test isequal(size(J)[2], prod(model.n))
    @test test_transpose(J)

    # get index
    J_sub = J[1]
    @test isequal(J_sub.info.nsrc, 1)
    @test isequal(J_sub.model, J.model)
    @test isequal(size(J_sub)[1], Int(size(J)[1]/nsrc))

    inds = nsrc > 1 ? (1:nsrc) : 1
    J_sub = J[inds]
    @test isequal(J_sub.info.nsrc, nsrc)
    @test isequal(J_sub.model, J.model)
    @test isequal(size(J_sub), size(J))

end

######################################################## judiJacobian ##################################################

@testset "judiJacobianExtended Unit Test with $(nsrc) sources" for nsrc=[1, 2]

    info = example_info(nsrc=nsrc)
    model = example_model()
    rec_geometry = example_rec_geometry(nsrc=nsrc)
    wavelet = randn(Float32, rec_geometry.nt[1])
    # Setup operators
    Pr = judiProjection(info, rec_geometry)
    F = judiModeling(info, model)
    Pw = judiLRWF(info, wavelet)
    w = judiWeights(randn(Float32, model.n); nsrc=nsrc)
    J = judiJacobian(Pr*F*Pw', w)

    @test isequal(typeof(J), judiJacobianExQ{Float32, Float32})
    @test isequal(J.recGeometry, rec_geometry)
    for i=1:nsrc
        @test isapprox(J.wavelet[i], wavelet)
        @test isapprox(J.weights[i], w.weights[i])
    end
    @test isequal(size(J)[2], prod(model.n))
    @test test_transpose(J)

    # get index
    J_sub = J[1]
    @test isequal(J_sub.info.nsrc, 1)
    @test isequal(J_sub.model, J.model)
    @test isapprox(J_sub.weights, J.weights[1:1])
    @test isequal(size(J_sub)[1], Int(size(J)[1]/nsrc))

    inds = nsrc > 1 ? (1:nsrc) : 1
    J_sub = J[inds]
    @test isequal(J_sub.info.nsrc, nsrc)
    @test isequal(J_sub.model, J.model)
    @test isapprox(J_sub.weights, J.weights[1:nsrc])
    @test isequal(size(J_sub), size(J))

    # Test Pw alone
    P1 = subsample(Pw, 1)
    @test isapprox(P1.wavelet, Pw[1].wavelet)
    @test isapprox(conj(Pw).wavelet, Pw.wavelet)
    @test size(conj(Pw)) == size(Pw)
    @test size(transpose(Pw)) == size(Pw)[end:-1:1]
    @test isapprox(transpose(Pw).wavelet, Pw.wavelet)

    P1 = adjoint(Pw) * w
    @test isequal(typeof(P1), judiExtendedSource{Float32})

end

####################################################### judiProjection #################################################

@testset "judiProjection Unit Test with $(nsrc) sources" for nsrc=[1, 2]

    info = example_info(nsrc=nsrc)
    rec_geometry = example_rec_geometry(nsrc=nsrc)
    src_geometry = example_src_geometry(nsrc=nsrc)
    wavelet = randn(Float32, src_geometry.nt[1])
    q = judiVector(src_geometry, wavelet)

    Pr = judiProjection(info, rec_geometry)
    Ps = judiProjection(info, src_geometry)

    @test isequal(typeof(Pr), judiProjection{Float32, Float32})
    @test isequal(typeof(Ps), judiProjection{Float32, Float32})
    @test isequal(Pr.geometry, rec_geometry)
    @test isequal(Ps.geometry, src_geometry)

    @test test_transpose(Pr)
    @test test_transpose(Ps)

    Pr_sub = Pr[1]
    @test isequal(Pr_sub.info.nsrc, 1)
    @test isequal(size(Pr_sub), convert(Tuple{Int64, Int64}, size(Pr) ./ nsrc))

    inds = nsrc > 1 ? (1:nsrc) : 1
    Pr_sub = Pr[inds]
    @test isequal(Pr_sub.info.nsrc, nsrc)
    @test isequal(size(Pr_sub), size(Pr))

    RHS = transpose(Ps)*q
    @test isequal(typeof(RHS), judiRHS{Float32})
    @test isequal(RHS.geometry, q.geometry)

end

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
function gemvtest(n, iter)
    A = rand(n,n)
    x = rand(n)
    z = similar(x)
    for i = 1:iter
        z = A * x
    end
    z
end

problemsizes =[(2, 10^6, "tiny"), (2^4, 10^5, "small"), (2^6, 10^4, "medium"), (2^8, 10^3, "large"), (2^10, 10^2, "huge")]
testdata = [(gemvtest, "gemv", "matrix-vector multiplication", problemsizes)]
include("../perfgeneric.jl")

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
module TestHMList

using Reagents
using Test
using ..Utils: concurrently, random_sleep

include("../../../examples/hmlist.jl")

function test_12()
    list = HMList{Int}()
    @test !(1 in list)
    push!(list, 1)
    @test 1 in list
    @test !(2 in list)
    push!(list, 2)
    @test 1 in list
    @test 2 in list
    delete!(list, 2)
    @test 1 in list
    @test !(2 in list)
    delete!(list, 1)
    @test !(1 in list)
    @test !(2 in list)
end

function test_1_to_9()
    list = HMList{Int}()
    for i in 1:9
        push!(list, i)
    end
    for i in 1:9
        @test i in list
    end
    for i in 1:9
        push!(list, i)
    end
    for i in 1:9
        @test i in list
    end
    for i in 1:9
        delete!(list, i)
    end
    for i in 1:9
        @test !(i in list)
    end
end

function test_spawn()
    list = HMList{Int}()
    iterated = Threads.Atomic{Int}(0)
    maxiter = 100_000

    function notdone()
        n = Threads.atomic_add!(iterated, 1)
        if mod(n, 1000) == 0
            random_sleep(false)
            # @show 11 in list
            # mod(n, 10000) == 0 && @show n
        end
        return n < maxiter
    end

    persistent = 10:10:90
    pushdelete = persistent .+ 1

    for i in persistent
        push!(list, i)
    end

    function pusher()
        while notdone()
            for i in pushdelete
                push!(list, i)
            end
        end
    end
    function deleter()
        while notdone()
            for i in pushdelete
                delete!(list, i)
            end
        end
    end
    function always_contains!(ref, xs)
        function repeat_always_contains!()
            while notdone()
                for i in xs
                    ref[] &= i in list
                    if !ref[]
                        iterated[] = maxiter
                        break
                    end
                end
            end
        end
    end
    function never_contains!(ref, xs)
        function repeat_never_contains!()
            while notdone()
                for i in xs
                    ref[] &= !(i in list)
                    if !ref[]
                        iterated[] = maxiter
                        break
                    end
                end
            end
        end
    end

    a1 = Ref(true)
    a2 = Ref(true)
    n1 = Ref(true)
    n2 = Ref(true)

    concurrently(
        pusher,
        pusher,
        deleter,
        deleter,
        always_contains!(a1, persistent),
        always_contains!(a2, persistent),
        never_contains!(n1, pushdelete .+ 1),
        never_contains!(n2, pushdelete .+ 1),
        spawn = true,
    )
    @test a1[] === a2[] === n1[] === n2[] === true
end

end  # module


text:
function acrobot_dynamics!(Ẋ, X, U)
    m1 = 1.0; #mass of link 1
    m2 = 1.0; #mass of link 2
    L1 = 1.0; #length of link 1
    L2 = 1.0; #length of link 2
    Lc1 = 0.5; #distance from pivot to center of link 1
    Lc2 = 0.5; #distance from pivot to center of link 2
    J1 = (1.0/3.0)*m1*L1*L1; #moment of inertia of link 1 about pivot
    J2 = (1.0/3.0)*m2*L2*L2; #moment of inertia of link 2 about pivot
    g = 9.81;

    q = X[1:2];
    qd = X[3:4];

    if isfinite(q[2])
        s = sin(q[2])
        c = cos(q[2])
    else
        s = Inf
        c = Inf
    end

    s1 = sin(q[1]);
    s12 = sin(q[1]+q[2]);
    s2 = sin(q[2]);
    c2 = cos(q[2]);

    H = [J1+J2+m2*L1*L1+2*m2*L1*Lc2*c2  J2+m2*L1*Lc2*c2;  J2+m2*L1*Lc2*c2  J2];
    C = [-2*m2*L1*Lc2*s2*qd[2]  -m2*L1*Lc2*s2*qd[2];  m2*L1*Lc2*s2*qd[1]  0];
    G = [(m1*Lc1 + m2*L1)*g*s1+m2*g*L2*s12;  m2*g*L2*s12];
    B = [0; 1];

    qdd = -H\(C*qd + G - B*U');

    Ẋ[1:2] = qd
    Ẋ[3:4] = qdd
    return nothing
end

n = 4
m = 1
model = Model(acrobot_dynamics!,n,m)

x0 = [0.0;0.0;0.0;0.0]
xf = [pi;0.0;0.0;0.0]

Q = 1.0*Diagonal(I,n)
Qf = 1000.0*Diagonal(I,n)
R = 0.1*Diagonal(I,m)

dt = 0.05
tf = 10.0
N = tf/dt+1

obj_uncon = LQRObjective(Q, R, Qf, tf, x0, xf)
opts = SolverOptions()
opts.verbose = true;
opts.square_root = true;

solver = Solver(model,obj_uncon,dt=dt,integration=:rk3,opts=opts);
U0 = rand(solver.model.m,solver.N-1);
@time results, stats = solve(solver,U0);

plot(results.X,label="")
plot(results.U)

u_min = -5
u_max = 5
obj_con = ConstrainedObjective(obj_uncon,u_min=u_min,u_max=u_max,use_xf_equality_constraint=false)

solver = Solver(model,obj_con,dt=dt,integration=:rk3,opts=opts);
U0 = rand(solver.model.m,solver.N-1);
@time results, stats = solve(solver,U0);

plot(results.X,label="")
plot(results.U)

## dircol
dircol_options = Dict("tol"=>solver.opts.cost_tolerance,"constr_viol_tol"=>solver.opts.constraint_tolerance)
X0 = rollout(solver,U0)
@time results_dircol, stats_dircol = TrajectoryOptimization.solve_dircol(solver, X0, U0, options=dircol_options)

plot(results_dircol.X',label="")
plot(results_dircol.U[:,1:solver.N-1]')

## Square root
x0 = [0.0;0.0;0.0;0.0]
xf = [pi;0.0;0.0;0.0]

Q = 1.0*Diagonal(I,n)
Q[1,1] = 1e4
Q[2,2] = 1e-12
Q[3,3] = 1e4
Q[4,4] = 1e-12
Qf = 1000.0*Diagonal(I,n)
Qf[1,1] = 1e6
Qf[2,2] = 1e-6
Qf[3,3] = 1.0
Qf[4,4] = 1e-6
R = 1e-6*Diagonal(I,m)

dt = 0.05
tf = 10.0
N = tf/dt+1

obj_uncon = LQRObjective(Q, R, Qf, tf, x0, xf)
opts = SolverOptions()
opts.verbose = true;
opts.square_root = false;

solver = Solver(model,obj_uncon,dt=dt,integration=:rk3,opts=opts);
U0 = rand(solver.model.m,solver.N-1);
@time results, stats = solve(solver,U0);

plot(results.X,label="")
plot(results.U)

norm(results.X[end] - solver.obj.xf)

u_min = -5
u_max = 5
obj_con = ConstrainedObjective(obj_uncon,u_min=u_min,u_max=u_max,use_xf_equality_constraint=false)

opts.constraint_tolerance = 1e-12
opts.penalty_max = 1e12
opts.penalty_scaling = 100.0
solver = Solver(model,obj_con,dt=dt,integration=:rk3,opts=opts);
solver.opts.square_root = false
solver_sqrt = Solver(model,obj_con,dt=dt,integration=:rk3,opts=opts);
solver_sqrt.opts.square_root = true
U0 = rand(solver.model.m,solver.N-1);
@time results, stats = solve(solver,U0);
@time results_sqrt, stats_sqrt = solve(solver_sqrt,U0);

plot(stats["c_max"])
plot!(stats_sqrt["c_max"])

plot(stats["S condition number"])
plot!(stats_sqrt["S condition number"])


text:
using DataStructures: OrderedDict

@testset "Datasets" begin
    data = [rand(4, 5, 12), rand(4, 5, 12), rand(4, 5)]
    axlist1 = [
        RangeAxis("XVals", 1.0:4.0),
        CategoricalAxis("YVals", [1, 2, 3, 4, 5]),
        RangeAxis("Time", Date(2001, 1, 15):Month(1):Date(2001, 12, 15)),
    ]
    axlist2 = [RangeAxis("XVals", 1.0:4.0), CategoricalAxis("YVals", [1, 2, 3, 4, 5])]
    props = [Dict("att$i" => i) for i = 1:3]
    c1, c2, c3 = (
        YAXArray(axlist1, data[1], props[1]),
        YAXArray(axlist1, data[2], props[2]),
        YAXArray(axlist2, data[3], props[3]),
    )
    ds = Dataset(avar = c1, something = c2, smaller = c3)
    @testset "Basic functions" begin
        b = IOBuffer()
        show(b, ds)
        s = split(String(take!(b)), "\n")
        s2 = """
        YAXArray Dataset
        Dimensions:
        XVals               Axis with 4 Elements from 1.0 to 4.0
        YVals               Axis with 5 elements: 1 2 3 4 5
        Time                Axis with 12 Elements from 2001-01-15 to 2001-12-15
        Variables: avar something smaller """
        s2 = split(s2, "\n")
        #     @test s[[1,2,6]] == s2[[1,2,6]]
        #     @test all(i->in(i,s2), s[3:5])
        for n in [:avar, :something, :smaller, :XVals, :Time, :YVals]
            @test n in propertynames(ds)
            @test n in propertynames(ds, true)
        end
        @test :axes ∉ propertynames(ds)
        @test :cubes ∉ propertynames(ds)
        @test :axes ∈ propertynames(ds, true)
        # Test getproperty
        @test all(i -> in(i, values(ds.axes)), axlist1)
        @test collect(keys(ds.cubes)) == [:avar, :something, :smaller]
        @test collect(values(ds.cubes)) == [c1, c2, c3]
        @test ds.avar === c1
        @test ds.something === c2
        @test ds.smaller === c3
        @test ds[:avar] === c1
        ds2 = ds[[:avar, :smaller]]
        @test collect(keys(ds2.cubes)) == [:avar, :smaller]
        @test collect(values(ds2.cubes)) == [c1, c3]
        @test YAXArrays.Datasets.fuzzyfind("hal", ["hallo", "bla", "something"]) == 1
        ds3 = ds[["av", "some"]]
        @test collect(keys(ds3.cubes)) == [:av, :some]
        @test collect(values(ds3.cubes)) == [c1, c2]
        @test ds["avar"] === c1
    end
    @testset "Dataset interface" begin
        struct MockDataset
            vars::Any
            dims::Any
            attrs::Any
            path::Any
        end
        Base.getindex(d::MockDataset, i) = d.vars[i]
        Base.haskey(d::MockDataset, i) = haskey(d.vars, i)
        YAXArrayBase.get_varnames(d::MockDataset) = (keys(d.vars)...,)
        YAXArrayBase.get_var_dims(d::MockDataset, name) = d.dims[name]
        YAXArrayBase.get_var_attrs(d::MockDataset, name) = d.attrs[name]
        YAXArrayBase.allow_missings(d::MockDataset) = !occursin("nomissings", d.path)
        function YAXArrayBase.create_empty(::Type{MockDataset}, path)
            mkpath(dirname(path))
            open(_ -> nothing, path, "w")
            MockDataset(Dict(), Dict(), Dict(), path)
        end
        function YAXArrayBase.add_var(ds::MockDataset, T, name, s, dimlist, atts; kwargs...)
            data = Array{T}(undef, s...)
            ds.vars[name] = data
            ds.dims[name] = dimlist
            ds.attrs[name] = atts
            data
        end
        YAXArrayBase.backendlist[:mock] = MockDataset
        push!(YAXArrayBase.backendregex, r".mock$" => MockDataset)
        data1, data2, data3, d1, d2, d3 =
            (rand(12, 5, 2), rand(12, 5), rand(12, 5, 2), 1:12, 0.1:0.1:0.5, ["One", "Two"])
        att1 = Dict("att1" => 5, "_ARRAY_OFFSET" => (2, 0, 0))
        att2 = Dict("att2" => 6, "_ARRAY_OFFSET" => (2, 0))
        attd1 = Dict(
            "_ARRAY_OFFSET" => 2,
            "units" => "days since 2001-01-01",
            "calendar" => "gregorian",
        )
        attd2 = Dict("attd" => "d")
        attd3 = Dict("attd" => "d")
        function MockDataset(p)
            MockDataset(
                OrderedDict(
                    "Var1" => data1,
                    "Var2" => data2,
                    "Var3" => data3,
                    "time" => d1,
                    "d2" => d2,
                    "d3" => d3,
                ),
                Dict(
                    "Var1" => ("time", "d2", "d3"),
                    "Var2" => ("time", "d2"),
                    "Var3" => ("time", "d2", "d3"),
                    "time" => ("time",),
                    "d2" => ["d2"],
                    "d3" => ["d3"],
                ),
                Dict(
                    "Var1" => att1,
                    "Var2" => att2,
                    "Var3" => att1,
                    "time" => attd1,
                    "d2" => attd2,
                    "d3" => attd3,
                ),
                p,
            )
        end
        m = MockDataset("testpath.mock")
        @testset "collectdims" begin
            dcollect = YAXArrays.Datasets.collectdims(m)
            @test dcollect["time"].ax isa RangeAxis
            @test YAXArrays.Cubes.Axes.axname(dcollect["time"].ax) == "time"
            @test dcollect["time"].ax.values ==
                  DateTime(2001, 1, 4):Day(1):DateTime(2001, 1, 13)
            @test dcollect["time"].offs == 2
            @test dcollect["d2"].ax isa RangeAxis
            @test YAXArrays.Cubes.Axes.axname(dcollect["d2"].ax) == "d2"
            @test dcollect["d2"].ax.values == 0.1:0.1:0.5
            @test dcollect["d2"].offs == 0
            @test dcollect["d3"].ax isa CategoricalAxis
            @test YAXArrays.Cubes.Axes.axname(dcollect["d3"].ax) == "d3"
            @test dcollect["d3"].ax.values == ["One", "Two"]
            @test dcollect["d3"].offs == 0
            a1 = [0.1, 0.2, 0.3, 0.4]
            a2 = [0.1, 0.21, 0.3, 0.4]
            @test YAXArrays.Datasets.testrange(a1) == 0.1:0.1:0.4
            @test YAXArrays.Datasets.testrange(a2) isa Array
            @test YAXArrays.Datasets.testrange(a2) == [0.1, 0.21, 0.3, 0.4]
        end
        @testset "open_dataset" begin
            ds = open_dataset("test.mock")
            @test size(ds.Var1) == (10, 5, 2)
            @test size(ds.Var2) == (10, 5)
            @test all(in(keys(ds.axes)), (:time, :d2, :d3))
            ar = Cube(ds)
            @test ar isa YAXArray
            @test size(ar) == (10, 5, 2, 2)
            @test YAXArrays.Cubes.Axes.axname.(ar.axes) == ["time", "d2", "d3", "Variable"]
            @test ar.axes[4].values == ["Var1", "Var3"]
        end
        @testset "Dataset creation" begin
            al = [
                RangeAxis("Time", Date(2001):Month(1):Date(2001, 12, 31)),
                CategoricalAxis("Variable", ["A", "B"]),
                RangeAxis("Xvals", 1:10),
            ]
            # Basic
            newds, newds2 = YAXArrays.Datasets.createdataset(MockDataset, al)
            @test YAXArrays.Cubes.axsym.(newds2.axes) == [:Time, :Xvals, :Variable]
            @test newds2.axes[1].values == Date(2001):Month(1):Date(2001, 12, 31)
            @test newds2.axes[3].values == ["A", "B"]
            @test newds2.axes[2].values == 1:10
            @test newds2.data isa YAXArrays.Cubes.DiskArrayTools.DiskArrayStack
            # A bit more advanced
            fn = string(tempname(), ".mock")
            newds, newds2 = YAXArrays.Datasets.createdataset(
                MockDataset,
                al,
                path = fn,
                persist = false,
                chunksize = (4, 2, 4),
                chunkoffset = (2, 0, 3),
                properties = Dict("att1" => 5),
                datasetaxis = "A",
            )
            @test size(newds.data) == (12, 2, 10)
            @test size(newds.data.parent) == (14, 2, 13)
            @test eltype(newds.data) <: Union{Float32,Missing}
            @test newds.properties["att1"] == 5
            @test isfile(fn)
            newds = nothing
            newds2 = nothing
            # Without missings
            fn = string(tempname(), "nomissings.mock")
            newds = YAXArrays.Datasets.createdataset(
                MockDataset,
                al,
                path = fn,
                datasetaxis = "A",
            )
        end
    end
end

@testset "Saving and loading between different backends" begin
    using NetCDF, Zarr
    x = rand(10, 5)
    ax1 = CategoricalAxis("Ax1", string.(1:10))
    ax2 = RangeAxis("Ax2", 1:5)
    p = tempname()
    savecube(YAXArray([ax1, ax2], x), p, backend = :zarr)
    cube1 = Cube(p)
    @test cube1.Ax1 == ax1
    @test cube1.Ax2 == ax2
    @test eltype(cube1.Ax2.values) <: Int64
    @test cube1.data == x
    p2 = string(tempname(), ".nc")
    savecube(cube1, p2, backend = :netcdf)
    cube2 = Cube(p2)
    @test cube2.Ax1 == ax1
    @test cube2.Ax2 == ax2
    @test cube2.data == x
    @test eltype(cube2.Ax2.values) <: Int64
end


text:
#
# Proquint.jl by P. Bayer, MIT license
#

const conson = "bdfghjklmnprstvz"
const vowels = "aiou"
const d = Dict{Char,Int}()
for i in eachindex(conson)
    d[conson[i]] = i-1
end
for i in eachindex(vowels)
    d[vowels[i]] = i-1
end
d['x'] = 0

c2i(c) = d[c]
i2v(i) = vowels[i+1]
i2c(i) = conson[i+1]

function quint(i)
    q = fill('x',5)
    q[5] = i2c(i & 0xF)
    q[4] = i2v((i>>4) & 0x3)
    q[3] = i2c((i>>6) & 0xF)
    q[2] = i2v((i>>10) & 0x3)
    q[1] = i2c((i>>12) & 0xF)
    return String(q)
end

function squint(i)
    s = ((i>>(12,10,6,4,0)[j]) & (0xF,0x3,0xF,0x3,0xF)[j] for j in 1:5)
    q = fill('x',5); push = false; k = 1
    for j in s
        j != 0 && (push = true)
        push && (q[k] = isodd(k) ? i2c(j) : i2v(j))
        k += 1
    end
    str = String(filter(!=('x'), q))
    return isempty(str) ? "x" : str
end

function uint(q::T) where T<:AbstractString
    res = 0
    j = 5 - length(q)
    for i in eachindex(q)
        res += c2i(q[i]) << (12,10,6,4,0)[i+j]
    end
    return UInt(res)
end

"""
    uint2quint(x::T, sep="-"; short=false) where T<:Integer

Convert an integer value into a corresponding proquint string identifier.

# Arguments
- `x`: an integer argument,
- `sep="-"`: a separator string,
- `short=false`: determines whether a short string should be created.

# Example
```julia
julia> uint2quint(0x7f000001)
"lusab-babad"

julia> uint2quint(0x7f000001, short=true)
"lusab-d"
```
"""
function uint2quint(x::T, sep="-"; short=false) where T<:Integer
    f = short ? squint : quint
    qs = sizeof(x) ÷ 2
    @assert qs ≥ 1 "size of x=$(repr(x)) must be at least 2 bytes!"
    s = ""
    for q in 1:qs
        s = f(x)*s
        x >>= 16
        q == qs || (s = sep*s)
    end
    return s
end

"""
    quint2uint(qs::String, U::Type{T}=UInt) where T<:Unsigned

Convert a proquint string identifier `qs` into a corresponding integer value.

# Arguments
- `U`: unsigned integer type e.g. `UInt16`, `UInt32`, `UInt64` or `UInt128`

# Example
```julia
julia> quint2uint("lusab-babad", UInt32)
0x7f000001

julia> quint2uint("lusab-d", UInt32)
0x7f000001
```
"""
function quint2uint(qs::String, U::Type{T}=UInt) where T<:Unsigned
    qa = split(qs, "-") |> reverse
    l = length(qa)
    res = 0
    for i in 1:l
        res += uint(qa[i]) << ((i-1)*16)
    end
    return U(res)
end


text:
# ForwardDiff integration

byhand = [:exp2, :log2, :exp10, :log10, :abs]

for f in libdevice
  if haskey(ForwardDiff.DiffRules.DEFINED_DIFFRULES, (:Base,f,1))
    f ∈ byhand && continue
    diffrule = ForwardDiff.DiffRules.DEFINED_DIFFRULES[(:Base,f,1)]
    ForwardDiff.DiffRules.DEFINED_DIFFRULES[(:CUDAnative,f,1)] =
      (args...) -> replace_device(diffrule(args...))
    eval(ForwardDiff.unary_dual_definition(:CUDAnative, f))
  end
end

# byhand: exp2
ForwardDiff.DiffRules.DEFINED_DIFFRULES[(:CUDAnative, :exp2, 1)] = x ->
  :((CuArrays.cufunc(exp2))(x) * (CuArrays.cufunc(log))(oftype(x, 2)))
eval(ForwardDiff.unary_dual_definition(:CUDAnative, :exp2))

# byhand: log2
ForwardDiff.DiffRules.DEFINED_DIFFRULES[(:CUDAnative, :log2, 1)] = x ->
   :(inv(x) / (CuArrays.cufunc(log))(oftype(x, 2)))
eval(ForwardDiff.unary_dual_definition(:CUDAnative, :log2))

# byhand: exp10
ForwardDiff.DiffRules.DEFINED_DIFFRULES[(:CUDAnative, :exp10, 1)] = x ->
  :((CuArrays.cufunc(exp10))(x) * (CuArrays.cufunc(log))(oftype(x, 10)))
eval(ForwardDiff.unary_dual_definition(:CUDAnative, :exp10))

# byhand: log10
ForwardDiff.DiffRules.DEFINED_DIFFRULES[(:CUDAnative, :log10, 1)] = x ->
   :(inv(x) / (CuArrays.cufunc(log))(oftype(x, 10)))
eval(ForwardDiff.unary_dual_definition(:CUDAnative, :log10))

# byhand: abs
ForwardDiff.DiffRules.DEFINED_DIFFRULES[(:CUDAnative, :abs, 1)] = x ->
   :(signbit(x) ? -one(x) : one(x))
eval(ForwardDiff.unary_dual_definition(:CUDAnative, :abs))

ForwardDiff.DiffRules.DEFINED_DIFFRULES[(:CUDAnative, :pow, 2)] = (x, y) ->
  replace_device.(ForwardDiff.DiffRules.DEFINED_DIFFRULES[(:Base, :^, 2)](x, y))

@eval begin
  ForwardDiff.@define_binary_dual_op(
      CUDAnative.pow,
      begin
        vx = ForwardDiff.value(x)
        vy = ForwardDiff.value(y)
        expv = (CUDAnative.pow)(vx, vy)

        powval = vy * CUDAnative.pow(vx , vy - Int32(1))

        py = ForwardDiff.partials(y)
        px = ForwardDiff.partials(x)

        cond = all(py.values) do x
          x == zero(x)
        end

        if cond
          logval = one(expv)
        else
          logval = expv * CUDAnative.log(vx)
        end

        new_partials = powval * px + logval * py
        return ForwardDiff.Dual{Txy}(expv, new_partials)
      end,
      begin
        v = ForwardDiff.value(x)
        expv = (CUDAnative.pow)(v, y)
        if y == zero(y)
          new_partials = zero(ForwardDiff.partials(x))
        else
          new_partials = ForwardDiff.partials(x) * y * (CUDAnative.pow)(v, y - Int32(1))
        end
        return ForwardDiff.Dual{Tx}(expv, new_partials)
      end,
      begin
        v = ForwardDiff.value(y)
        expv = (CUDAnative.pow)(x, v)
        deriv = expv*CUDAnative.log(x)
        return ForwardDiff.Dual{Ty}(expv, deriv * ForwardDiff.partials(y))
      end
    )
end


text:
import ContentHashes
const C = ContentHashes

using Test

@testset "ContentHashes.jl" begin
    struct T{S}
        x::S
    end
    @testset "primatives" begin
        @test C.hash(1) === C.hash(1)
        @test C.hash("string"^10) === C.hash("string"^10)
        @test C.hash([]) === C.hash([])
        @test C.hash(()) === C.hash(())
        @test C.hash(Dict(1 => 2, 2 => 1)) === C.hash(Dict(2 => 1, 1 => 2))
        @test C.hash(Matrix{Any}(undef, 5, 5)) === C.hash(Matrix{Any}(undef, 5, 5))
        @test C.hash(ntuple(n -> T(undef), 50)) === C.hash(ntuple(n -> T(undef), 50))

        @test C.hash(1, 1) === C.hash(1, 1)
        @test C.hash(1, 1) !== C.hash(1, 2)
    end
    @testset "anonymous functions" begin
        f = x -> T(@fastmath sin(x[1:5] + 1 / 2))
        g = x -> T(@fastmath sin(x[1:5] + 1 / 2))

        @test hash(f) !== hash(g)
        @test C.hash(f) === C.hash(g)
    end
    @testset "circular references" begin
        a = []
        push!(a, a)
        b = []
        push!(b, b)
        @test C.hash(a) === C.hash(b)

        push!(a, a)
        push!(b, b)
        @test C.hash(a) === C.hash(b)
        push!(a, T(Ref([])))
        push!(b, T(Ref([])))
        @test C.hash(a) === C.hash(b)

        push!(a[1], a)
        push!(b[1], b)
        @test C.hash(a) === C.hash(b)
    end
end


text:
#=

# Interfacing with the C API

## Structures and variable preservation

As the API is written in C, there are a lot of pointers to deal with and handling them is not always an easy task. With a little practice, one can figure out how to wrap function calls with `cconvert` and `unsafe_convert` provided by Julia. Those functions provide automatic conversions and `ccall` GC-roots `cconvert`ed variables to ensure that pointers will point to valid memory (by explicitly telling the compiler not to garbage-collect nor optimize away the original variable).

However, the situation gets a lot more complicated when you deal with pointers as type fields. We will look at a naive example that show how difficult it can get for a Julia developer not used to calling C code. If we wanted to create a `VkInstance`, we might be tempted to do:

=#

using Vulkan.core

function create_instance(app_name, engine_name)
    app_info = VkApplicationInfo(
        VK_STRUCTURE_TYPE_APPLICATION_INFO, # sType
        C_NULL, # pNext
        pointer(app_name), # application name
        1, # application version
        pointer(engine_name), # engine name
        0, # engine version
        VK_VERSION_1_2, # requested API version
    )
    create_info = InstanceCreateInfo(
        VK_STRUCTURE_TYPE_INSTANCE_CREATE_INFO, # sType
        C_NULL, # pNext
        0, # flags
        Base.unsafe_convert(Ptr{VkApplicationInfo}, (Ref(app_info))), # application info
        0, # layer count
        C_NULL, # layers (none requested)
        0, # extension count
        C_NULL, # extensions (none requested)
    )
    p_instance = Ref{VkInstance}()

    GC.@preserve app_info begin
        vkCreateInstance(
            Ref(create_info),
            C_NULL, # custom allocator (we choose the default one provided by Vulkan)
            p_instance,
        )
    end

    p_instance[]
end

## instance = create_instance("AppName", "NoEngine") # very likely to segfault

#=

which will probably result in a segmentation fault. Why?

Two causes may lead to such a result:
1. `app_name` and `engine_name` may never be allocated if the compiler decides not to, so there is no guarantee that `pointer(app_name)` and `pointer(engine_name)` will point to anything valid. Additionally, even if those variables were allocated with valid pointer addresses at some point, they can be garbage collected **at any time**, including before the call to `vkCreateInstance`.
2. `app_info` is not what should be preserved. It cannot be converted to a pointer, but a `Ref` to it can. Therefore it is the reference that needs to be `GC.@preserve`d, not `app_info`. So, `Ref(app_info)` must be assigned to a variable, and replace `app_info` in the call to `GC.@preserve`.

Basically, it all comes down to having to preserve everything you take a pointer of. And, if you need to create an intermediary object when converting a variable to a pointer, you need to preserve it too. For example, take of an array of `String`s, that need to be converted as a `Ptr{Cstring}`. You first need to create an array of `Cstring`s, then convert that array to a pointer. The `String`s and the `Cstring` array need to be preseved.

This is exactly what `cconvert` and `unsafe_convert` are for. `cconvert` converts a variable to a type that can be converted to the desired (possibly pointer) type using `unsafe_convert`. In addition of chaining both conversions, `ccall` also preserves the `cconvert`ed variable, so that the unsafe conversion becomes safe.

Because we cannot use `ccall` in this case, we need to `cconvert` any argument that will be transformed to a pointer, and store the result as long as the desired struct may be used. Then, `unsafe_convert` can be called on this result, to get the desired (pointer) type necessary to construct the API struct.

There are several possibilities for preserving what we may call "pointer dependencies". One of them is to reference them inside a global variable, such as a `Dict`, and deleting them once we no longer need it. This has the severe disadvantage of requiring to explicitly manage every dependency, along with large performance issues. Another possibility, which we have taken in this wrapper, is to create a new structure that will store both the API structure and the required dependencies. That way, we can safely rely on the GC for preserving what we need just when we need it.

Therefore, every API structure is wrapped inside another one (without the "Vk" prefix), as follows:

=#

abstract type VulkanStruct{has_deps} end

struct InstanceCreateInfo <: VulkanStruct{true}
    vks::VkInstanceCreateInfo # API struct
    deps::Vector{Any}         # contains all required dependencies
end

#=

and every structure exposes a convenient constructor that works perfectly with `String`s and mutable `AbstractArray`s. No manual `Ref`s/`cconvert`/`unsafe_convert` needed.

We hope that the additional `Vector{Any}` will not introduce too much overhead. In the future, this might be changed to a `NTuple{N, Any}` or a `StaticArrays.SVector{N, Any}`. We could also have stored dependencies as additional fields, but this does not scale well with nested structs. It would either require putting an additional field for each dependency (be it direct, or indirect dependencies coming from a pointer to another struct), possibly defining other structures that hold dependencies to avoid having a large number of fields, inducing additional compilation time.

!!! tip
    `cconvert`/`unsafe_convert` were extended on wrapper types so that, when using an API function directly, [`ccall`](https://docs.julialang.org/en/v1/base/c/#ccall) will convert a struct to its API-compatible version.

## Version numbers

Version numbers are encoded in unsigned integers (often `UInt32`) in C. Most APIs have a specific encoding (and decoding) scheme. Ideally we want to use version numbers from Julia without having to look for this encoding or decoding function, and this has been automated for the Vulkan API.

=#


text:
# ---
# title: 967. Numbers With Same Consecutive Differences
# id: problem967
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search, Breadth-first Search
# link: <https://leetcode.com/problems/numbers-with-same-consecutive-differences/description/>
# hidden: true
# ---
#
# Return all **non-negative** integers of length `n` such that the absolute
# difference between every two consecutive digits is `k`.
#
# Note that **every** number in the answer **must not** have leading zeros
# **except** for the number `0` itself. For example, `01` has one leading zero
# and is invalid, but `0` is valid.
#
# You may return the answer in **any order**.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 3, k = 7
#     Output: [181,292,707,818,929]
#     Explanation: Note that 070 is not a valid number, because it has leading zeroes.
#
#
# **Example 2:**
#
#
#
#     Input: n = 2, k = 1
#     Output: [10,12,21,23,32,34,43,45,54,56,65,67,76,78,87,89,98]
#
#
# **Example 3:**
#
#
#
#     Input: n = 2, k = 0
#     Output: [11,22,33,44,55,66,77,88,99]
#
#
# **Example 4:**
#
#
#
#     Input: n = 2, k = 1
#     Output: [10,12,21,23,32,34,43,45,54,56,65,67,76,78,87,89,98]
#
#
# **Example 5:**
#
#
#
#     Input: n = 2, k = 2
#     Output: [13,20,24,31,35,42,46,53,57,64,68,75,79,86,97]
#
#
#
#
# **Constraints:**
#
#   * `2 <= n <= 9`
#   * `0 <= k <= 9`
#
#
## @lc code=start
using LeetCode

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


text:
export emst

using mlpack._Internal.io

import mlpack_jll
const emstLibrary = mlpack_jll.libmlpack_julia_emst

# Call the C binding of the mlpack emst binding.
function emst_mlpackMain()
  success = ccall((:emst, emstLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module emst_internal
  import ..emstLibrary

end # module

"""
    emst(input; [leaf_size, naive, verbose])

This program can compute the Euclidean minimum spanning tree of a set of input
points using the dual-tree Boruvka algorithm.

The set to calculate the minimum spanning tree of is specified with the `input`
parameter, and the output may be saved with the `output` output parameter.

The `leaf_size` parameter controls the leaf size of the kd-tree that is used to
calculate the minimum spanning tree, and if the `naive` option is given, then
brute-force search is used (this is typically much slower in low dimensions).
The leaf size does not affect the results, but it may have some effect on the
runtime of the algorithm.

For example, the minimum spanning tree of the input dataset `data` can be
calculated with a leaf size of 20 and stored as `spanning_tree` using the
following command:

```julia
julia> using CSV
julia> data = CSV.read("data.csv")
julia> spanning_tree = emst(data; leaf_size=20)
```

The output matrix is a three-dimensional matrix, where each row indicates an
edge.  The first dimension corresponds to the lesser index of the edge; the
second dimension corresponds to the greater index of the edge; and the third
column corresponds to the distance between the two points.

# Arguments

 - `input::Array{Float64, 2}`: Input data matrix.
 - `leaf_size::Int`: Leaf size in the kd-tree.  One-element leaves give
      the empirically best performance, but at the cost of greater memory
      requirements.  Default value `1`.

 - `naive::Bool`: Compute the MST using O(n^2) naive algorithm.  Default
      value `false`.

 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `output::Array{Float64, 2}`: Output data.  Stored as an edge list.

"""
function emst(input;
              leaf_size::Union{Int, Missing} = missing,
              naive::Union{Bool, Missing} = missing,
              verbose::Union{Bool, Missing} = missing,
              points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, emstLibrary), Nothing, ());

  IORestoreSettings("Fast Euclidean Minimum Spanning Tree")

  # Process each input argument before calling mlpackMain().
  IOSetParamMat("input", input, points_are_rows)
  if !ismissing(leaf_size)
    IOSetParam("leaf_size", convert(Int, leaf_size))
  end
  if !ismissing(naive)
    IOSetParam("naive", convert(Bool, naive))
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output")
  # Call the program.
  emst_mlpackMain()

  return IOGetParamMat("output", points_are_rows)
end