text:
@testset "Data types and missings" begin
    # TODO: Convert these test_group things to nested testsets
    #test_group("missings")
    @test ismissing(3 == missing)
    @test ismissing(missing == 3)
    @test ismissing(missing == missing)

    #test_group("DataVector creation")
    dvint = @data [1, 2, missing, 4]
    dvint2 = DataArray(collect(5:8))
    dvint3 = convert(DataArray, 5:8)
    dvflt = @data [1.0, 2, missing, 4]
    dvstr = @data ["one", "two", missing, "four"]
    # FIXME: triggers a segfault on Julia 0.6.0
    # dvdict = DataArray(Dict, 4) # for issue DataFrames#199
    dvany = convert(DataArray{Any, 1}, dvint)

    @test isa(dvint, DataVector{Int})
    @test isa(dvint2, DataVector{Int})
    @test isa(dvint3, DataVector{Int})
    @test isa(dvflt, DataVector{Float64})
    @test isa(dvstr, DataVector{String})
    @test_throws ArgumentError DataArray([5:8], falses(2))

    #test_group("PooledDataVector creation")
    pdvstr = @pdata ["one", "one", "two", "two", missing, "one", "one"]
    @test isa(pdvstr, PooledDataVector{String})
    @test isequal(PooledDataArray(pdvstr), pdvstr)

    #test_group("PooledDataVector creation with predetermined pool")
    pdvpp = PooledDataArray([1, 2, 2, 3], [1, 2, 3, 4])
    @test isequal(pdvpp.pool, [1, 2, 3, 4])
    @test string(pdvpp) == "[1, 2, 2, 3]"
    pdvpp = PooledDataArray([1, 2, 2, 3, 2, 1], [1, 2, 3, 4])
    @test isequal(pdvpp.pool, [1, 2, 3, 4])
    @test string(pdvpp) == "[1, 2, 2, 3, 2, 1]"
    pdvpp = PooledDataArray(["one", "two", "two"], ["one", "two", "three"])
    @test isequal(convert(DataArray, pdvpp), @data(["one", "two", "two"]))
    @test isequal(levels(pdvpp), @data(["one", "two", "three"]))
    @test isequal(pdvpp.pool, ["one", "two", "three"])
    @test string(pdvpp) == "[one, two, two]"
    @test string(PooledDataArray(["one", "two", "four"],
                                   ["one", "two", "three"])) ==
            "[one, two, missing]"

    #test_group("PooledDataVector utf8 support")
    pdvpp = PooledDataArray([String("hello")], [false])
    @test isa(pdvpp[1], String)
    pdvpp = PooledDataArray([String("hello")])
    @test isa(pdvpp[1], String)

    #test_group("DataVector access")
    @test dvint[1] == 1
    @test ismissing(dvint[3])
    @test isequal(dvflt[3:4], @data([missing, 4.0]))
    @test isequal(dvint[[true, false, true, false]], @data([1, missing]))
    @test isequal(dvstr[[1, 2, 1, 4]], @data(["one", "two", "one", "four"]))
    # Indexing produces #undef?
    # @test isequal(dvstr[[1, 2, 1, 3]], DataVector["one", "two", "one", missing])

    #test_group("PooledDataVector access")
    @test pdvstr[1] == "one"
    @test ismissing(pdvstr[5])
    @test isequal(pdvstr[1:3], @data(["one", "one", "two"]))
    @test isequal(pdvstr[[true, false, true, false, true, false, true]],
                    @pdata(["one", "two", missing, "one"]))
    @test isequal(pdvstr[[1, 3, 1, 2]], @data(["one", "two", "one", "one"]))

    #test_group("DataVector methods")
    @test size(dvint) == (4,)
    @test length(dvint) == 4
    @test sum(ismissing.(dvint)) == 1
    @test eltype(dvint) == Union{Int,Missing}

    #test_group("PooledDataVector methods")
    @test size(pdvstr) == (7,)
    @test length(pdvstr) == 7
    @test sum(ismissing.(pdvstr)) == 1
    @test eltype(pdvstr) == Union{String,Missing}

    #test_group("DataVector operations")
    @test isequal(dvint .+ 1, DataArray([2, 3, 4, 5], [false, false, true, false]))
    @test isequal(dvint .* 2, @data([2, 4, missing, 8]))
    @test isequal(dvint .== 2, @data([false, true, missing, false]))
    @test isequal(dvint .> 1, @data([false, true, missing, true]))

    #test_group("PooledDataVector operations")
    # @test isequal(pdvstr .== "two", PooledDataVector[false, false, true, true, missing, false, false])

    #test_group("DataVector to something else")
    @test collect(skipmissing(dvint)) == [1, 2, 4]
    @test all(convert(Vector, dvint, 0) .== [1, 2, 0, 4])
    @test all(convert(Vector, dvany, 0) .== [1, 2, 0, 4])
    utf8three = convert(String, "three")
    asciithree = convert(String, "three")
    @test all(convert(Vector, dvstr, utf8three) .== ["one", "two", "three", "four"])
    @test all(convert(Vector, dvstr, asciithree) .== ["one", "two", "three", "four"])
    @test all(convert(Vector{Int}, dvint2) .== [5:8;])
    @test all([i + 1 for i in dvint2] .== [6:9;])
    #@test all([length(x)::Int for x in dvstr] == [3, 3, 1, 4])
    # Julia 0.6 and 0.7 differ in ordering of Unions
    @test repr(dvint) in ("Union{$Int, Missings.Missing}[1, 2, missing, 4]",
                          "Union{Missings.Missing, $Int}[1, 2, missing, 4]",
                          "Union{$Int, Missing}[1, 2, missing, 4]",
                          "Union{Missing, $Int}[1, 2, missing, 4]")

    #test_group("PooledDataVector to something else")
    @test collect(skipmissing(pdvstr)) == ["one", "one", "two", "two", "one", "one"]
    @test all(convert(Vector, pdvstr, "nine") .== ["one", "one", "two", "two", "nine", "one", "one"])
    #@test all([length(i)::Int for i in pdvstr] .== [3, 3, 3, 3, 1, 3, 3])
    @test string(pdvstr[1:3]) == "[one, one, two]"

    #test_group("DataVector Filter and Replace")
    @test collect(skipmissing(dvint)) == [1, 2, 4]
    @test isequal(convert(Vector, dvint, 7), [1, 2, 7, 4])
    @test sum(skipmissing(dvint)) == 7
    @test sum(convert(Vector, dvint, 7)) == 14

    #test_group("PooledDataVector Filter and Replace")
    @test reduce(string, "", skipmissing(pdvstr)) == "oneonetwotwooneone"
    @test reduce(string, "", convert(Vector, pdvstr, "!")) == "oneonetwotwo!oneone"

    #test_group("DataVector assignment")
    assigntest = @data [1, 2, missing, 4]
    assigntest[1] = 8
    @test isequal(assigntest, (@data [8, 2, missing, 4]))
    assigntest[1:2] = 9
    @test isequal(assigntest, (@data [9, 9, missing, 4]))
    assigntest[[1,3]] = 10
    @test isequal(assigntest, (@data [10, 9, 10, 4]))
    assigntest[[true, false, true, true]] = 11
    @test isequal(assigntest, (@data [11, 9, 11, 11]))
    assigntest[1:2] = [12, 13]
    @test isequal(assigntest, (@data [12, 13, 11, 11]))
    assigntest[[1, 4]] = [14, 15]
    @test isequal(assigntest, (@data [14, 13, 11, 15]))
    assigntest[[true, false, true, false]] = [16, 17]
    @test isequal(assigntest, (@data [16, 13, 17, 15]))
    assigntest[1] = missing
    @test isequal(assigntest, (@data [missing, 13, 17, 15]))
    assigntest[[1, 2]] = missing
    @test isequal(assigntest, (@data [missing, missing, 17, 15]))
    assigntest[[true, false, true, false]] = missing
    @test isequal(assigntest, (@data [missing, missing, missing, 15]))
    assigntest[1] = 1
    assigntest[2:4] = missing
    @test isequal(assigntest, (@data [1, missing, missing, missing]))

    #test_group("PooledDataVector assignment")
    ret = (pdvstr[2] = "three")
    @test ret == "three"
    @test pdvstr[2] == "three"
    ret = pdvstr[[1,2]] = "two"
    @test ret == "two"
    @test pdvstr[2] == "two"
    pdvstr2 = @pdata ["one", "one", "two", "two"]
    ret = (pdvstr2[[true, false, true, false]] = "three")
    @test ret == "three"
    @test pdvstr2[1] == "three"
    ret = (pdvstr2[[false, true, false, true]] = ["four", "five"])
    @test isequal(ret, ["four", "five"])
    @test isequal(pdvstr2[3:4], (@data ["three", "five"]))
    pdvstr2 = @pdata ["one", "one", "two", "two"]
    ret = (pdvstr2[2:3] = "three")
    @test ret == "three"
    @test isequal(pdvstr2[3:4], (@data ["three", "two"]))
    ret = (pdvstr2[2:3] = ["four", "five"])
    @test ret == ["four", "five"]
    @test isequal(pdvstr2[1:2], (@data ["one", "four"]))
    pdvstr2 = @pdata ["one", "one", "two", "two", "three"]
    @test ismissing(begin pdvstr2[1] = missing end)
    @test all(ismissing(begin pdvstr2[[1, 2]] = missing end))
    @test all(ismissing(begin pdvstr2[[false, false, true, false, false]] = missing end))
    @test all(ismissing(begin pdvstr2[4:5] = missing end))
    @test all(ismissing.(pdvstr2))

    #test_group("PooledDataVector replace!")
    pdvstr2 = @pdata ["one", "one", "two", "two", "three"]
    @test replace!(pdvstr2, "two", "four") == "four"
    @test replace!(pdvstr2, "three", "four") == "four"
    @test ismissing.(replace!(pdvstr2, "one", missing))
    @test replace!(pdvstr2, missing, "five") == "five"
    @test isequal(pdvstr2, (@data ["five", "five", "four", "four", "four"]))
end


text:
# Bond Helpers
type FixedRateBondHelper <: BondHelper
  price::Quote
  bond::FixedRateBond
end
value(b::FixedRateBondHelper) = b.price.value
maturity_date(b::FixedRateBondHelper) = maturity_date(b.bond)

# bond helper functions
function implied_quote{B <: BondHelper}(bond_h::B, clean::Bool = true)
  bond = bond_h.bond
  recalculate!(bond)
  settlement_value = bond.settlementValue
  return clean ? clean_price(bond, settlement_value, settlement_date(bond)) : dirty_price(bond, settlement_value, settlement_date(bond))
end

function clone(fixedRateBondHelper::FixedRateBondHelper, ts::TermStructure)
  # first we have to clone the PE
  newPE = clone(fixedRateBondHelper.bond.pricingEngine, ts)
  # then we have to clone the bond
  newBond = clone(fixedRateBondHelper.bond, newPE)

  return FixedRateBondHelper(fixedRateBondHelper.price, newBond)
end

update_termstructure(bondHelper::FixedRateBondHelper, ts::TermStructure) = clone(bondHelper, ts)

ext: jl


text:
function update_particles_nofluid_nocl!(pList,wList,param,data)
    pUarr,pVarr = compute_particle_velocity_nofluids_nocl(pList,wList,param,data)
    for ti=1:data.n_particles
        pList[ti].pos.x += pUarr[ti]*param.Δt
        pList[ti].pos.y += pVarr[ti]*param.Δt
    end

    nothing
end

# does not use cell lists
function compute_particle_velocity_nofluids_nocl(pList,wList,param,data)
    compute_gravity_force!(data.gfX,data.gfY,data.n_particles,param.G)
    compute_cohesion_force_backup!(data.cfX,data.cfY,pList,param.sC,param.ϵ)
    compute_adhesion_force!(data.afX,data.afY,pList,wList,data.pointOnWall,param.sA,param.ϵ)

    pUarr = data.gfX + data.cfX + data.afX
    pVarr = data.gfY + data.cfY + data.afY

    return pUarr,pVarr
end

function update_particles_nofluid!(pList,wList,param,data,cl)
    pUarr,pVarr = compute_particle_velocity_nofluids(pList,wList,param,data,cl)
    for ti=1:data.n_particles
        pList[ti].pos.x += pUarr[ti]*param.Δt
        pList[ti].pos.y += pVarr[ti]*param.Δt
    end

    nothing
end

# does not use cell lists
function compute_particle_velocity_nofluids(pList,wList,param,data,cl)
    compute_gravity_force!(data.gfX,data.gfY,data.n_particles,param.G)
    compute_cohesion_force_cl!(data.cfX,data.cfY,cl,pList,param.sC,param.ϵ)
    compute_adhesion_force!(data.afX,data.afY,pList,wList,data.pointOnWall,param.sA,param.ϵ)

    pUarr = data.gfX + data.cfX + data.afX
    pVarr = data.gfY + data.cfY + data.afY

    return pUarr,pVarr
end


text:
using MyCalculus
using Test

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


text:
# This file includes the testset for links

@testset "LinkTestSet" begin
    @info "Running LinkTestSet ..."

    # Link construction.
    l = Link(5)
    @test eltype(l) == Float64
    @test eltype(l.channel) == Float64
    @test l.channel.sz_max == 0
    @test length(l.buffer) == 5
    @test mode(l.buffer) == Cyclic
    @test !iswritable(l)
    @test !isreadable(l)

    # More on Buffer construction
    l = Link{Int}(5)
    @test size(l.buffer) == (5,)
    @test eltype(l) == Int
    l = Link{Bool}()
    @test size(l.buffer) == (64,)
    l = Link()
    @test eltype(l) == Float64
    @test size(l.buffer) == (64,)

    # Putting values to link
    l = Link()
    t = @async while true
        take!(l) === NaN && break
    end
    vals = collect(1:5)
    for i = 1 : length(vals)
        put!(l, vals[i])
        @test l.buffer[1] == vals[i]
    end
    close(l)
    @test istaskdone(t)
    @test !isopen(l.channel)

    # Taking values from the link
    l = Link()
    vals = collect(1 : 10)
    t = launch(l, vals)
    val = take!(l)
    @test val == 1.
    for i = 2 : 10
        @test take!(l) == vals[i]
    end
    close(l)
    wait(t)
    @test istaskdone(t)

    @info "Done LinkTestSet."
end  # testset


text:
"""
    LazyNode{Name, D, C} <: AbstractMillNode

Data node storing data of type `D` in a lazy manner and optional metadata of type `C`.

Source of data or its type is specified in `Name`.

See also: [`AbstractMillNode`](@ref), [`LazyModel`](@ref), [`Mill.unpack2mill`](@ref).
"""
struct LazyNode{Name, D, C} <: AbstractMillNode
    data::D
    metadata::C
end

"""
    LazyNode([Name::Symbol], d, m=nothing)
    LazyNode{Name}(d, m=nothing)

Construct a new [`LazyNode`](@ref) with name `Name`, data `d`, and metadata `m`.

# Examples
```jldoctest
julia> LazyNode(:Codons, ["GGGCGGCGA", "CCTCGCGGG"])
LazyNode{Codons} 	# 2 obs, 98 bytes
```

See also: [`AbstractMillNode`](@ref), [`LazyModel`](@ref), [`Mill.unpack2mill`](@ref).
"""
LazyNode(Name::Symbol, d::T, m::C=nothing) where {T, C} = LazyNode{Name, T, C}(d, m)
LazyNode{Name}(d::T, m::C=nothing) where {Name, T, C} = LazyNode{Name, T, C}(d, m)

"""
    Mill.unpack2mill(x::LazyNode)

Return a representation of [`LazyNode`](@ref) `x` using `Mill.jl` structures. Every custom
[`LazyNode`](@ref) should have a special method as it is used in [`LazyModel`](@ref).

# Examples
```jldoctest unpack2mill; output=false
function Mill.unpack2mill(ds::LazyNode{:Sentence})
    s = split.(ds.data, " ")
    x = NGramMatrix(reduce(vcat, s))
    BagNode(ArrayNode(x), Mill.length2bags(length.(s)))
end
# output

```
```jldoctest unpack2mill
julia> LazyNode{:Sentence}(["foo bar", "baz"]) |> Mill.unpack2mill
BagNode 	# 2 obs, 120 bytes
  └── ArrayNode(2053×3 NGramMatrix with Int64 elements) 	# 3 obs, 274 bytes
```

See also: [`LazyNode`](@ref), [`LazyModel`](@ref).
"""
function unpack2mill end

Base.ndims(x::LazyNode) = Colon()
StatsBase.nobs(a::LazyNode) = length(a.data)
StatsBase.nobs(a::LazyNode, ::Type{ObsDim.Last}) = nobs(a.data)

mapdata(f, x::LazyNode{N}) where N = LazyNode{N}(mapdata(f, x.data), x.metadata)

dropmeta(x::LazyNode{N}) where N = LazyNode{N}(x.data)

function Base.reduce(::typeof(catobs), as::Vector{<: LazyNode{N}}) where N
    LazyNode{N}(reduce(catobs, data.(as)), reduce(catobs, metadata.(as)))
end

Base.getindex(x::LazyNode{N, T}, i::VecOrRange{<:Int}) where {N, T} = LazyNode{N}(subset(x.data, i))

Base.hash(e::LazyNode, h::UInt) = hash((e.data), h)
(e1::LazyNode == e2::LazyNode) = e1.data == e2.data
Base.isequal(e1::LazyNode, e2::LazyNode) = isequal(e1.data, e2.data)


text:
# --------- Defaults & Helpers
_default_matsize_stresstest() = 2048 # linear size of matrices; 2048 x 2048 should be efficiently implemented in CUBLAS
_default_use_mem_stresstest() = 0.9 # use 90% of available (i.e. free) memory

function _estimate_t_matmul(C, A, B)
    t_matmul = device!(device(C)) do
        nsamples = 5
        CUDA.@elapsed mul!(C, A, B) # warmup
        t_matmul = Inf
        for _ in 1:nsamples
            t_sample = CUDA.@elapsed for _ in 1:5
                mul!(C, A, B)
            end
            t_matmul = min(t_sample / 5, t_matmul)
        end
        # CUDA.@sync mul!(C, A, B) # warmup
        # t_matmul = Inf
        # for _ in 1:nsamples
        #     t_sample = @elapsed CUDA.@sync for _ in 1:5 mul!(C, A, B) end
        #     t_matmul = min(t_sample/5, t_matmul)
        # end
        return t_matmul
    end
    return t_matmul
end

# --------- AbstractStressTest
abstract type AbstractStressTest{T} end

_eltype(::Type{AbstractStressTest{T}}) where {T} = T
_eltype(::AbstractStressTest{T}) where {T} = T

# --------- StressTestBatched
"""
GPU stress test (matrix multiplications) in which we try to run for a given time period.
We try to keep the CUDA stream continously busy with matmuls at any point in time.
Concretely, we submit batches of matmuls and, after half of them, we record a CUDA event.
On the host, after submitting a batch, we (non-blockingly) synchronize on, i.e. wait for, the CUDA event
and, if we haven't exceeded the desired duration already, submit another batch.
"""
struct StressTestBatched{T,S} <: AbstractStressTest{T}
    device::CuDevice
    C::S
    A::S
    B::S
    duration::Int64
    size::Int64
    batchsize::Int64
    batch_duration::Int64

    function StressTestBatched(
        device::CuDevice=device();
        duration=60,
        batch_duration=nothing,
        dtype=Float32,
        verbose=true,
        size=_default_matsize_stresstest(),
    )
        C, A, B = device!(device) do
            C = CUDA.zeros(dtype, size, size)
            A = CUDA.rand(dtype, size, size)
            B = CUDA.rand(dtype, size, size)
            return C, A, B
        end
        if isnothing(batch_duration)
            batch_duration = ceil(Int, duration / 10) # reasonable default?
        end
        t_matmul = _estimate_t_matmul(C, A, B)
        batchsize = ceil(Int, batch_duration / t_matmul)
        verbose && @info(
            "Chose batch size such that a batch takes approx. $(batch_duration) seconds.",
            t_matmul,
            batchsize,
        )
        return new{dtype,typeof(A)}(
            device, C, A, B, duration, size, batchsize, batch_duration
        )
    end
end

function (st::StressTestBatched)(; verbose=false)
    verbose && @info(
        "Running StressTestBatched{$(_eltype(st))} on Julia thread $(Threads.threadid()) and $(st.device)."
    )
    counter = device!(st.device) do
        C, A, B, duration, batchsize = st.C, st.A, st.B, st.duration, st.batchsize
        event_halfbatch = CuEvent()
        i_record = Int(batchsize ÷ 2)
        counter = 0
        t = time()
        CUDA.@sync @inbounds while ((time() - t) < duration)
            for i in 1:batchsize
                mul!(C, A, B)
                counter += 1
                i == i_record && CUDA.record(event_halfbatch)
            end
            synchronize(event_halfbatch)
            # alternative:
            # while !CUDA.isdone(event_halfbatch)
            #     nothing
            #     # println("not done.")
            # end
            # println("Batch half done.")
        end
        return counter
    end
    verbose && @info("Ran $counter iterations on $(st.device).")
    return nothing
end

# --------- StressTestEnforced
"""
GPU stress test (matrix multiplications) in which we run almost precisely for a given time period (duration is enforced).
"""
struct StressTestEnforced{T,S} <: AbstractStressTest{T}
    device::CuDevice
    C::S
    A::S
    B::S
    enforced_duration::Int64
    size::Int64

    function StressTestEnforced(
        device::CuDevice=device();
        enforced_duration=5,
        dtype=Float32,
        verbose=true,
        size=_default_matsize_stresstest(),
    )
        C, A, B = device!(device) do
            C = CUDA.zeros(dtype, size, size)
            A = CUDA.rand(dtype, size, size)
            B = CUDA.rand(dtype, size, size)
            return C, A, B
        end
        return new{dtype,typeof(A)}(device, C, A, B, enforced_duration, size)
    end
end

function (st::StressTestEnforced)(; verbose=false)
    verbose && @info(
        "Running StressTestEnforced{$(_eltype(st))} on Julia thread $(Threads.threadid()) and $(st.device)."
    )
    counter = 0
    device!(st.device) do
        C, A, B, enforced_duration = st.C, st.A, st.B, st.enforced_duration
        t = time()
        @inbounds while (time() - t) < enforced_duration
            CUDA.@sync mul!(C, A, B)
            counter += 1
        end
    end
    verbose && @info("Ran $counter iterations on $(st.device).")
    return nothing
end

# --------- StressTestStoreResults (gpu-burn-like)
"""
GPU stress test (matrix multiplications) in which we store all matmul results and try to
run as many iterations as possible for a certain memory limit (default: 90% of free memory).

This stress test is somewhat inspired by [gpu-burn](https://github.com/wilicc/gpu-burn) by Ville Timonen.
"""
struct StressTestStoreResults{T,S} <: AbstractStressTest{T}
    device::CuDevice
    C::Vector{S}
    A::S
    B::S
    niters::Int64
    size::Int64

    function StressTestStoreResults(
        device::CuDevice=device();
        mem=nothing,
        dtype=Float32,
        verbose=true,
        size=_default_matsize_stresstest(),
    )
        max_use_bytes = device!(() -> _memarg2maxusebytes(mem), device)
        mat_sizeof = sizeof(dtype) * size * size
        niters = Int((max_use_bytes - 2 * mat_sizeof) ÷ mat_sizeof) # we substract sizeof A and B
        use_bytes = niters * mat_sizeof
        if verbose
            mem_free, mem_total = device!(CUDA.Mem.info, device)
            @info(
                "Allocating $(bytes(use_bytes)) of memory on $(device) for results. ($(bytes(mem_free)) free, $(bytes(mem_total)) total)"
            )
            @info("Will be able to run $niters iterations.")
        end
        if use_bytes < 3 * mat_sizeof
            error("Low memory for result. Aborting.")
        end
        C, A, B = device!(device) do
            C = [CUDA.zeros(dtype, size, size) for _ in 1:niters]
            A = CUDA.rand(dtype, size, size)
            B = CUDA.rand(dtype, size, size)
            return C, A, B
        end
        return new{dtype,typeof(A)}(device, C, A, B, niters, size)
    end
end

function (st::StressTestStoreResults)(; verbose=false)
    verbose && @info(
        "Running StressTestStoreResults{$(_eltype(st))} on Julia thread $(Threads.threadid()) and $(st.device)."
    )
    device!(st.device) do
        C, A, B, niters = st.C, st.A, st.B, st.niters
        @inbounds for i in 1:niters
            mul!(C[i], A, B)
        end
    end
    return nothing
end

_containertype(::Type{StressTestStoreResults{T,S}}) where {T,S} = S
_containertype(::StressTestStoreResults{T,S}) where {T,S} = S

_memarg2maxusebytes(mem::Nothing) = CUDA.available_memory() * _default_use_mem_stresstest()
_memarg2maxusebytes(mem::UnitPrefixedBytes) = bytes(mem)
function _memarg2maxusebytes(mem::Real)
    if 0 < mem < 1
        return CUDA.available_memory() * mem
    else
        throw(
            ArgumentError(
                "Provided `mem` argument is < 0 or > 1 and thus not a valid fraction."
            ),
        )
    end
end

# --------- StressTestFixedIter
"""
GPU stress test (matrix multiplications) in which we run for a given number of iteration,
or try to run for a given time period (with potentially high uncertainty!).
In the latter case, we estimate how long a synced matmul takes and set `niter` accordingly.
"""
struct StressTestFixedIter{T,S} <: AbstractStressTest{T}
    device::CuDevice
    C::S
    A::S
    B::S
    size::Int64
    niter::Int64

    function StressTestFixedIter(
        device::CuDevice=device();
        approx_duration=60,
        niter=nothing,
        dtype=Float32,
        verbose=true,
        size=_default_matsize_stresstest(),
    )
        C, A, B = device!(device) do
            C = CUDA.zeros(dtype, size, size)
            A = CUDA.rand(dtype, size, size)
            B = CUDA.rand(dtype, size, size)
            return C, A, B
        end
        if isnothing(niter)
            t_matmul = _estimate_t_matmul(C, A, B)
            niter = ceil(Int, approx_duration / t_matmul)
            verbose && @info(
                "Estimated required number of iterations for desired duration.",
                t_matmul,
                niter
            )
        end
        return new{dtype,typeof(A)}(device, C, A, B, size, niter)
    end
end

function (st::StressTestFixedIter)(; verbose=false)
    verbose && @info(
        "Running StressTestFixedIter{$(_eltype(st))} on Julia thread $(Threads.threadid()) and $(st.device)."
    )
    device!(st.device) do
        C, A, B, niter = st.C, st.A, st.B, st.niter
        CUDA.@sync @inbounds for _ in 1:niter
            mul!(C, A, B)
        end
    end
    verbose && @info("Ran $(st.niter) iterations on $(st.device).")
    return nothing
end


text:
# Computes A^(⊗p) = A ⊗ A ⊗ ... ⊗ A
# where ⊗ is the Kronecker product
# Since p > 0, B will be kron'd at least once so it will be a Matrix{T}
# This is not easy to guess for Julia inference so I annotate the return type
function (kronpow(A::MT, p)::MT) where MT <: AbstractMatrix
    @assert p > 0
    B = 1
    C = A
    bitmap = 1
    while bitmap <= p
        if (bitmap & p) != 0
            B = kron(B, C)
        end
        C = kron(C, C)
        bitmap <<= 1
    end
    B
end


text:
abstract type Examples end
abstract type JuliaExamples <: Examples end
abstract type PSYExamples <: Examples end
abstract type PSIExamples <: Examples end
abstract type PSDExamples <: Examples end

const PACKAGE_DIR = dirname(dirname(pathof(SIIPExamples)))
const SCRIPT_DIR = joinpath(PACKAGE_DIR, "script")
const TEST_DIR = joinpath(PACKAGE_DIR, "test")
const NB_DIR = joinpath(PACKAGE_DIR, "notebook")

const JULIA_EX_FOLDER = "1_introduction"
const PSY_EX_FOLDER = "2_PowerSystems_examples"
const PSI_EX_FOLDER = "3_PowerSimulations_examples"
const PSD_EX_FOLDER = "4_PowerSimulationsDynamics_examples"

get_dir(::Type{JuliaExamples}) = JULIA_EX_FOLDER
get_dir(::Type{PSYExamples}) = PSY_EX_FOLDER
get_dir(::Type{PSIExamples}) = PSI_EX_FOLDER
get_dir(::Type{PSDExamples}) = PSD_EX_FOLDER


text:
module Utils

        using DataValues
        using Missings

        export make_type_block, counter, basiccensor, averagegain, zeroormissing
        export assignrand, distribute_proportionately

        function Vec( n :: Integer, initial :: Real ) :: Vector
            zeros(n) .+ initial
        end

        function make_type_block( r::UnitRange, thistype::DataType=Int) :: Dict{Int64,DataType}
                d = Dict{Int64, DataType}()
                for i in r
                        d[i]=thistype
                end;
                d
        end

        function assignrand( props :: Vector, payments :: Vector ) :: Real
                l = length( payments )
                @assert length( props ) == l
                r = rand()
                for j in 1:l
                    println( "j=$j rand=$r props[$j]")
                    if r <= props[j]
                        return payments[j]
                    end
                end
                @assert "black" == "white"
        end

        function Base.strip( m :: Missing )
        end

        function counter( v, what :: Real = 0.0 )::Integer
                c = 0
                for x in v
                        if x > what
                                c += 1
                        end
                end
                c
        end

        """
        returns the string converted to a form suitable to be used as (e.g.) a Symbol,
        with leading/trailing blanks removed, forced to lowercase, and with various
        characters replaced with '_' (at most '_' in a run).
        """
        function basiccensor( s :: AbstractString ) :: AbstractString
                s = strip( lowercase( s ))
                s = replace( s, r"[ \-,\t–]" => "_" )
                s = replace( s, r"[=\:\)\('’‘]" => "" )
                s = replace( s,  r"[\";:\.\?\*”“]" => "" )
                s = replace( s,  r"_$"=> "" )
                s = replace( s,  r"^_"=> "" )
                s = replace( s,  r"^_"=> "" )
                s = replace( s,  r"\/"=> "_or_" )
                s = replace( s,  r"\&"=> "_and_" )
                s = replace( s,  r"\+"=> "_plus_" )
                s = replace( s,  r"_\$+$"=> "" )
                if occursin( r"^[\d].*", s )
                        s = string("v_", s ) # leading digit
                end
                s = replace( s,  r"__+"=> "_" )
                s = replace( s, r"^_" => "" )
                s = replace( s, r"_$" => "" )
                return s
        end

        function averagegain( v :: Vector )
            n = 0.0
            c = 0.0
            for x in v
                if x > 0
                    n += 1
                    c += x
                end
            end
            println( "c=$c n=$n")
            c/n
        end

        function zeroormissing( d :: DataValue )
            return (d == NA ) || ( d == 0 )
        end

        function zeroormissing( d :: Real )
            return ( d == 0 )
        end

        function zeroormissing( t :: Tuple )
            d = zeroormissing.( t )
            return any( d->d, d )
        end

        """
        distribute `total_to_add` proportionately between the elements of v
        """
        function distribute_proportionately( v::AbstractArray, total_to_add ::Number )::AbstractArray
            total = sum( skipmissing( v ))
            m  = total_to_add/total
            return v .* (1+m)
        end
end


text:
"""
This file contains all the functions in ClimateIsca.jl that resorts the Isca GCM output and resaves
into NetCDF files grouped according to parameter/variables and resaved by year.

"""

function iscaresort(parvec,iroot::Dict)

    @info "$(Dates.now()) - Resorting the following Isca GCM output variables ..."
    pinfo = iscaparload(parvec);
    for ii = 1 : size(pinfo,1); @info "$(Dates.now()) - $(ii)) $(pinfo[ii,3]) | $(pinfo[ii,4])" end

    fol,fname = iscadinfo(iroot["raw"]); ndir = length(fol);

    pf  = ncread(fname,"pfull"); ph = ncread(fname,"phalf"); np = length(pf);
    lon = ncread(fname,"lon"); nlon = length(lon);
    lat = ncread(fname,"lat"); nlat = length(lat);
    tinfo,nt = iscatime(fname); dim = [nlon,nlat,nt];

    ncattr = Dict("fullpre"=>pf,"halfpre"=>ph,"nlevels"=>np,"dimensions"=>dim,
                  "longitude"=>lon,"latitude"=>lat,
                  "timestep"=>tinfo["dt"],"timespan"=>tinfo["span"]);

    for yr = 1 : ndir; cd(joinpath(iroot["raw"],fol[yr]))

        for par in parvec; pinfo = iscaparameter(par)

            pdata = iscaparncread(fname,pinfo,ncattr);
                    iscaparncsave(pdata,pinfo,ncattr,yr,iroot["data"]);

        end

    end

    cd(iroot["isca"]); @save "isca_info.jld2" ncattr;

end


text:
# Engheim's "functional" approach is mostly about using closures
# But just because he doesnt want to use struct based types, doesn't mean you have to not
# use multiple dispatch, which he does in his supposedly oo solution
# here is a cleaner functional implementation which preserves the use of closures

# Since the alphabet we are using is hardcoded in the function, the length can be hardcoded as well
# no need for a const
function encrypt(ch::Char, shift::Int)
    ch in 'A':'Z' ? 'A' + mod((ch - 'A') + shift, 26) :
    ch in 'a':'z' ? 'a' + mod((ch - 'a') + shift, 26) :
    #= otherwise =# ch
end

# if we dont plan to allow single char encyption we dont need this funcion
function decrypt(ch::Char, shift::Int)
    # we reuse Engheim's nice trick
    encrypt(ch,-shift)
end

# once you define the encrypt/decrypt pair for the cipher
# these functions handle strings for you
function encrypt(shift::Int)
    message -> map(message) do ch
        encrypt(ch, shift)
    end
end

function decrypt(shift::Int)
    message -> map(message) do ch
        encrypt(ch, -shift)
    end
end

# Example
# As Engheim notes the mdisadvantage of using a closure in this case:
# “You don’t have any guarantees that the decryption function is using the same shift
# value as the encryption function.” But in a real world application we can define
# a const that preserves this value across associated encrypt/decrypt
# julia> const s1 = 1
# 1
# There is also an extra step in having to define two closures, whereas in the so-called "OO"
# solution you only have to instanciate the cipher once. In that sense, preserving state in a
# struct does seem a better solution for this case
# julia> cipher = encrypt(s1)
# #5 (generic function with 1 method)

# julia> decipher = decrypt(s1)
# #9 (generic function with 1 method)

# julia> message = "ABC"
# "ABC"

# julia> emessage = cipher(message)
# "BCD"

# julia> decipher(emessage)
# "ABC"

# Substitution cipher
using Random
# Engheim's solution only works because he requires uppercase passwords
# this solution is a bit more complex

# a nice helper function to create the shuffles
function makeshuffle()::Tuple{Dict{Char, Char}, Dict{Char, Char}}
    lowdict = Dict(zip('a':'z', shuffle('a':'z')))
    updict = Dict(zip('A':'Z', shuffle('A':'Z')))
    return updict, lowdict
end

function encrypt(ch::Char, updict::Dict{Char, Char}, lowdict::Dict{Char, Char})
    ch in 'A':'Z' ? get(updict, ch, ch) :
    ch in 'a':'z' ? get(lowdict, ch, ch) :
    #= otherwise =# ch
end

# Once again this function is not necessary if we don't need single char decryption
function decrypt(ch::Char, rupdict::Dict{Char, Char}, rlowdict::Dict{Char, Char})
    ch in 'A':'Z' ? get(rupdict, ch, ch) :
    ch in 'a':'z' ? get(rlowdict, ch, ch) :
    #= otherwise =# ch
end

# Closures
function encrypt(updict::Dict{Char, Char}, lowdict::Dict{Char, Char})
    message -> map(message) do ch
        encrypt(ch, updict, lowdict)
    end
end

function decrypt(updict::Dict{Char, Char}, lowdict::Dict{Char, Char})
    rupdict = Dict(values(updict) .=> keys(updict))
    rlowdict = Dict(values(lowdict) .=> keys(lowdict))
    message -> map(message) do ch
        encrypt(ch, rupdict, rlowdict)
    end
end

# Example
# julia> dicts = makeshuffle()
# (Dict('E' => 'M', 'Z' => 'A', 'X' => 'O', 'B' => 'F', 'C' => 'L', 'D' => 'K', 'A' => 'Y', 'R' => 'P', 'G' => 'X', 'F' => 'J'…), Dict('n' => 't', 'f' => 'c', 'w' => 's', 'd' => 'k', 'e' => 'd', 'o' => 'w', 'h' => 'i', 'j' => 'b', 'i' => 'v', 'k' => 'm'…))

# julia> cipher = encrypt(dicts[1], dicts[2])
# #41 (generic function with 1 method)

# julia> decipher = decrypt(dicts[1], dicts[2])
# #45 (generic function with 1 method)

# julia> emessage = cipher("Google sucks!")
# "Xwwqhd ueamu!"

# julia> decipher(emessage)
# "Google sucks!"

# Password Keeper
mutable struct Vault
    passwords::Dict{String, String}
    encrypt::Function
    decrypt::Function
end

function Vault(encrypter, decrypter)
    Vault(Dict{String, String}(), encrypter, decrypter)
end

function addlogin!(vault::Vault, login::AbstractString, password::AbstractString)
    vault.passwords[login] = vault.encrypt(password)
end

function getpassword(vault::Vault, login::AbstractString)
    vault.decrypt(vault.passwords[login])
end

# Example
# Caeser
# julia> const s1 = 23
# 23

# julia> vault = Vault(encrypt(s1), decrypt(s1))
# Vault(Dict{String, String}(), var"#1#3"{Int64}(23), var"#5#7"{Int64}(23))

# julia> addlogin!(vault, "google", "BING")
# "YFKD"

# julia> getpassword(vault, "google")
# "BING"

# julia> addlogin!(vault, "amazon", "Secret!")
# "Pbzobq!"

# julia> getpassword(vault, "amazon")
# "Secret!"

# Substitution
# julia> using Random

# julia> dicts = makeshuffle()
# (Dict('E' => 'X', 'Z' => 'A', 'X' => 'C', 'B' => 'W', 'C' => 'P', 'D' => 'U', 'A' => 'K', 'R' => 'G', 'G' => 'S', 'F' => 'Y'…), Dict('n' => 'q', 'f' => 'y', 'w' => 'j', 'd' => 'v', 'e' => 'c', 'o' => 'h', 'h' => 'w', 'j' => 'd', 'i' => 'p', 'k' => 'f'…))

# julia> vaults = Vault(encrypt(dicts[1], dicts[2]), decrypt(dicts[1], dicts[2]))
# Vault(Dict{String, String}(), var"#9#11"{Dict{Char, Char}, Dict{Char, Char}}(Dict('E' => 'X', 'Z' => 'A', 'X' => 'C', 'B' => 'W', 'C' => 'P', 'D' => 'U', 'A' => 'K', 'R' => 'G', 'G' => 'S', 'F' => 'Y'…), Dict('n' => 'q', 'f' => 'y', 'w' => 'j', 'd' => 'v', 'e' => 'c', 'o' => 'h', 'h' => 'w', 'j' => 'd', 'i' => 'p', 'k' => 'f'…)), var"#13#15"{Dict{Char, Char}, Dict{Char, Char}}(Dict('n' => 'l', 'f' => 'k', 'w' => 'h', 'd' => 'j', 'e' => 'p', 'o' => 'q', 'j' => 'w', 'y' => 'f', 'h' => 'o', 'r' => 'r'…), Dict('E' => 'L', 'Z' => 'H', 'X' => 'E', 'C' => 'X', 'B' => 'S', 'D' => 'O', 'A' => 'Z', 'R' => 'P', 'G' => 'R', 'Q' => 'K'…)))

# julia> addlogin!(vaults, "apple", "My Super Secret!")
# "Hk Bxecr Bcsrci!"

# julia> getpassword(vaults, "apple")
# "My Super Secret!"

max_stars_repo_name: FourMInfo/JuliaforBeginners

max_issues_repo_name: FourMInfo/JuliaforBeginners

max_forks_repo_name: FourMInfo/JuliaforBeginners


text:
function absorb end
function PID end
function TID end
function between end

function decompose(f::FormulaTerm, data::AbstractDataFrame, contrasts::Dict{Symbol})
    data = dropmissing(data[termvars(f)], disallowmissing = true)
    rhs = isa(f.rhs, Tuple) ? collect(f.rhs) : [f.rhs]
    absorbed = findall(t -> isa(t, FunctionTerm{typeof(absorb)}), rhs)
    if isempty(absorbed)
        absorbed = (Vector{Symbol}(), Vector{Vector{Vector{Int}}}())
    elseif length(absorbed) == 1
        absorbed = rhs[absorbed[1]]
        if isa(absorbed, AbstractVector)
            absorbed = unique(flatten(termvars(t) for t ∈ absorbed))
        else
            absorbed = termvars(absorbed)
        end
        absorbed = (absorbed,
                    [ [ findall(isequal(l), data[dim]) for l ∈ levels(data[dim]) ] for dim ∈ absorbed ])
    else
        throw(ArgumentError("There can only be at most one absorb term"))
    end
    effect = findall(t -> isa(t, FunctionTerm{typeof(between)}), rhs)
    if isempty(effect)
        effect = (Vector{Symbol}(), Vector{Vector{Vector{Int}}}())
    elseif length(effect) == 1
        effect = termvars(f.rhs[effect[1]])
        effect = (effect,
                  [ findall(isequal(l), data[effect[1]]) for l ∈ levels(data[effect[1]]) ])
    else
        throw(ArgumentError("There can only be at most one effect for the between estimator"))
    end
    pid = findall(t -> isa(t, FunctionTerm{typeof(PID)}), rhs)
    if isempty(pid)
        pid = (Vector{Symbol}(), Vector{Vector{Vector{Int}}}())
    elseif length(pid) == 1
        pid = termvars(f.rhs[pid[1]])[1]
        pid = ([pid], [ findall(isequal(l), data[pid]) for l ∈ levels(data[pid]) ])
    else
        throw(ArgumentError("There can only be one panel identifier"))
    end
    tid = findall(t -> isa(t, FunctionTerm{typeof(TID)}), rhs)
    if isempty(tid)
        tid = (Vector{Symbol}(), Vector{Vector{Vector{Int}}}())
    elseif length(tid) == 1
        tid = termvars(f.rhs[tid[1]])[1]
        tid = ([tid], [ findall(isequal(l), data[tid]) for l ∈ levels(data[tid]) ])
    else
        throw(ArgumentError("There can only be one time identifier"))
    end
    wts = findall(t -> isa(t, FunctionTerm{typeof(weights)}), rhs)
    if isempty(wts)
        wts = FrequencyWeights(zeros(0))
    elseif length(wts) == 1
        wts = data[termvars(f.rhs[wts[1]])[1]]
    else
        throw(ArgumentError("There can only be one weight variable"))
    end
    exo_rhs = filter(t -> !(isa(t, FormulaTerm) ||
                     isa(t, FunctionTerm{<:Union{<:typeof(absorb),
                                                 <:typeof(PID),
                                                 <:typeof(TID),
                                                 <:typeof(weights),
                                                 <:typeof(between)}})),
                    rhs)
    exogenous = FormulaTerm(f.lhs,
                            isempty(exo_rhs) ? Tuple([ConstantTerm(1)]) : Tuple(exo_rhs))
    @assert length(terms(exogenous.lhs)) == 1 "Response should be a single variable"
    iv = findall(t -> isa(t, FormulaTerm), rhs)
    if length(iv) == 0
        iv = @eval(@formula(0 ~ 0))
    elseif length(iv) == 1
        iv = rhs[iv[1]]
    else
        throw(ArgumentError("Formula syntax is response ~ exogenous + (endogenous ~ instruments)"))
    end
    sc = schema(data, contrasts)
    f = apply_schema(f, sc)
    exogenous = apply_schema(exogenous, sc, EconometricsModel)
    iv = apply_schema(iv, sc)
    if isa(exogenous.lhs, CategoricalTerm)
        y = data[exogenous.lhs.sym]
        X = modelcols(exogenous.rhs, data)
        if isordered(y)
            X = X[:,2:end]
        end
    else
        y, X = modelcols(exogenous, data)
    end
    wts = isempty(wts) ? FrequencyWeights(Fill(1, size(y, 1))) : wts
    if !isa(iv.lhs, InterceptTerm)
        z, Z = modelcols(iv, data)
        z = isa(z, Tuple) ? hcat(z...) : z
    else
        z = zeros(0, 0)
        Z = zeros(0, 0)
    end
    data, f, exogenous, iv, absorbed, pid, tid, wts, effect, y, X, z, Z
end


text:
struct Biomass{V} <: Measurement
    val::V
end

Base.rand(::Type{Biomass}, names::Vector{Symbol}, dist::Distribution, x::Int, y::Int) =
    begin
        dict = Dict()

        for name in names
            dict[name] = rand(Biomass, dist, x, y)
        end
        return dict
    end

Base.rand(::Type{Biomass}, x::Int, y::Int) = Base.rand(Biomass, 0.5, x, y)
Base.rand(::Type{Biomass}, dist::Distribution, x::Int, y::Int) = begin
    z = zeros(Biomass, x, y)
    for i in eachindex(z)
        z[i] = rand(dist)
    end
    return z
end

Base.zero(::Biomass{T}) where {T} = Biomass{T}(zero(T))
Base.zero(::Type{Biomass}) = Biomass{Float64}(0)
Base.zero(::Type{<:Biomass{T}}) where {T} = Biomass(zero(T))

Base.zero(::NamedTuple{T,V}) where {T,V<:Tuple{Biomass}} = Biomass{Float32(0.0)}

Base.show(io::IO, t::Biomass{T}) where {T} = Base.show(io, T(t.val))

Base.convert(::Type{Biomass}, v::V) where {V<:Number} = Biomass{V}(v)

Base.isless(a::Biomass, b::Biomass) = isless(a.val, b.val)
Base.isless(a::Number, b::Biomass) = isless(a, b.val)
Base.isless(a::Biomass, b::Number) = isless(a.val, b)

Base.oneunit(::Type{<:Biomass{T}}) where {T} = Biomass(one(T))
Base.one(::Type{<:Biomass{T}}) where {T} = Biomass(one(T))
Base.one(::Biomass{T}) where {T} = Biomass{T}(one(T))

Base.:*(x::Biomass, v::Number) = Biomass(x.val * v)
Base.:*(v::Number, x::Biomass) = Biomass(x.val * v)

Base.:/(x::Biomass, v::Biomass) = Biomass(x.val / v.val)
Base.:/(x::Biomass, v::Number) = Biomass(x.val / v)
Base.:/(x::Number, v::Biomass) = Biomass(x / v.val)

Base.:+(x1::Number, x2::Biomass) = Biomass(x1 + x2.val)
Base.:+(x1::Biomass, x2::Number) = Biomass(x1.val + x2)

Base.:-(x1::Number, x2::Biomass) = Biomass(x1 - x2.val)
Base.:-(x1::Biomass, x2::Number) = Biomass(x1.val - x2)

Base.:^(x1::Biomass, x2::Number) = x1.val^x2
Base.:^(x1::Number, x2::Biomass) = x1^x2.val

Base.convert(::Type{Float64}, t::Biomass) = Float64(t.val)

Base.:+(x1::Biomass, x2::Biomass) = Biomass(x1.val + x2.val)
Base.:-(x1::Biomass, x2::Biomass) = Biomass(x1.val - x2.val)
Base.:*(v::Biomass, x::Biomass) = Biomass(x.val * v.val)

DynamicGrids.to_rgb(a::Any, v::Biomass) = DynamicGrids.to_rgb(a, v.val)


text:
function test_eftRecip{T<:AbstractFloat}(x::T)
    hi,lo = eftRecip(x)
    b = big(1)/big(x)
    tst = Float64(b - (big(hi)+big(lo)))
    eps(abs(lo))/2 >= abs(tst)
end


text:
# library.jl
## https://developer.spotify.com/documentation/web-api/reference/library/

## https://developer.spotify.com/documentation/web-api/reference/library/get-users-saved-tracks/
@doc """
# Get a User's Saved Tracks
**Summary**: Get a list of the songs saved in the current Spotify user's 'Your Music' library.\n

`limit` _Optional_: The maximum number of objects to return. Default 20. Minimum 1. Maximum 50. \n
`offset` _Optional_: The index of the first object to return. Default 0.\n
`market` _Optional_: An ISO 3166-1 alpha-2 country code.

[Reference](https://developer.spotify.com/documentation/web-api/reference/library/get-users-saved-tracks/)
""" ->
function library_get_saved_tracks(limit=20, offset=0, market="US")
    return spotify_request("me/tracks?market=$market&offset=$offset&limit=$limit")
end

## https://developer.spotify.com/documentation/web-api/reference/library/get-users-saved-shows/

@doc """
# Get a User's Saved Shows
**Summary**: Get a list of the songs saved in the current Spotify user's 'Your Music' library.\n

`limit` _Optional_: The maximum number of objects to return. Default 20. Minimum 1. Maximum 50. \n
`offset` _Optional_: The index of the first object to return. Default 0.\n

[Reference](https://developer.spotify.com/documentation/web-api/reference/library/get-users-saved-shows/)
""" ->
function library_get_saved_shows(limit=20, offset=0)
    return spotify_request("me/shows?limit=$limit&offset=$offset")
end

## https://developer.spotify.com/documentation/web-api/reference/library/get-users-saved-albums/
@doc """
# Get a User's Saved Albums
**Summary**: Get a list of the albums saved in the current Spotify user's 'Your Music' library.\n

`limit` _Optional_: The maximum number of objects to return. Default 20. Minimum 1. Maximum 50. \n
`offset` _Optional_: The index of the first object to return. Default 0.\n
`market` _Optional_: An ISO 3166-1 alpha-2 country code.

[Reference](https://developer.spotify.com/documentation/web-api/reference/library/get-users-saved-albums/)
""" ->
function library_get_saved_albums(limit=20, offset=0, market="US")
    return spotify_request("me/albums?limit=$limit&offset=$offset&market=$market")
end

## https://developer.spotify.com/documentation/web-api/reference/library/check-users-saved-tracks/
@doc """
# Check User's Saved Tracks
**Summary**: Check if one or more tracks is already saved in the current Spotify user's 'Your Music' library.\n

`track_ids` _Required_: A comma-separated list of the Spotify IDs for the tracks. Maximum 50.\n

[Reference](https://developer.spotify.com/documentation/web-api/reference/library/check-users-saved-tracks/)
""" ->
function library_check_saved_tracks(track_ids)
    return  spotify_request("me/tracks/contains?ids=$track_ids")
end

## https://developer.spotify.com/documentation/web-api/reference/library/check-users-saved-shows/
@doc """
# Check User's Saved Shows
**Summary**: Check if one or more shows is already saved in the current Spotify user's library.\n

`show_ids` _Required_: A comma separated list of the Spotify IDs for the shows. Maximum 50\n

[Reference](https://developer.spotify.com/documentation/web-api/reference/library/check-users-saved-shows/)
""" ->
function library_check_saved_shows(show_ids)
    return spotify_request("me/shows/contains?ids=$show_ids")
end

## https://developer.spotify.com/documentation/web-api/reference/library/check-users-saved-albums/
@doc """
# Check User's Saved Albums
**Summary**: Check if one or more albums is already saved in the current Spotify user's 'Your Music' library.\n

`album_ids` _Required_: A comma separated list of the Spotify IDs for the albums. Max 50.\n

[Required](https://developer.spotify.com/documentation/web-api/reference/library/check-users-saved-albums/)
""" ->
function library_check_saved_albums(album_ids)
    return spotify_request("me/albums/contains?ids=$album_ids")
end

#### DELETE ####

## https://developer.spotify.com/documentation/web-api/reference/library/remove-albums-user/
@doc """
# Remove Albums for Current User
**Summary**: Remove one or more albums for the current user's 'Your Music' library.\n

`album_ids` _Required_: A comma-separated list of the Spotify IDs. Maximum 50.\n

[Reference](https://developer.spotify.com/documentation/web-api/reference/library/remove-albums-user/)
""" ->
function library_remove_albums(album_ids)
    return spotify_request("me/albums?ids=$album_ids", method="DELETE")
end

## https://developer.spotify.com/documentation/web-api/reference/library/remove-shows-user/
@doc """
# Remove Shows for Current User
**Summary**: Remove one or more shows for the current user's library.\n

`show_ids` _Required_: A comma-separated list of the Spotify IDs. Maximum 50.\n

[Reference](https://developer.spotify.com/documentation/web-api/reference/library/remove-shows-user/)
""" ->
function library_remove_shows(show_ids)
    return spotify_request("me/shows?ids=$show_ids", method="DELETE")
end

## https://developer.spotify.com/documentation/web-api/reference/library/remove-tracks-user/
@doc """
# Remove Tracks for Current User
**Summary**: Remove one or more tracks for the current user's 'Your Music' library.\n

`track_ids` _Required_: A comma-separated list of the Spotify IDs. Maximum 50.\n

[Reference](https://developer.spotify.com/documentation/web-api/reference/library/remove-tracks-user/)
""" ->
function library_remove_tracks(track_ids)
    return spotify_request("me/tracks?ids=$track_ids", method="DELETE")
end

#### PUT #####

## https://developer.spotify.com/documentation/web-api/reference/library/save-albums-user/
@doc """
# Save Albums for Current User
** Summary**: Save one or more albums to the current user's 'Your Music' library.\n

`album_ids` _Required_: A comma-separated list of Spotify IDs. Maximum 50. \n

[Reference](https://developer.spotify.com/documentation/web-api/reference/library/save-albums-user/)
""" ->
function library_save_album(album_ids)
    return spotify_request("me/albums?ids=$album_ids", method="PUT")
end

## https://developer.spotify.com/documentation/web-api/reference/library/save-shows-user/
@doc """
# Save Shows for Current User
** Summary**: Save one or more shows to the current user's library.\n

`shows_ids` _Required_: A comma-separated list of Spotify IDs. Maximum 50. \n

[Reference](https://developer.spotify.com/documentation/web-api/reference/library/save-shows-user/)
""" ->
function library_save_show(shows_ids)
    return spotify_request("me/shows?ids=$show_ids", method="PUT")
end

## https://developer.spotify.com/documentation/web-api/reference/library/save-tracks-user/
@doc """
# Save Tracks for Current User
** Summary**: Save one or more tracks to the current user's 'Your Music' library.\n

`track_ids` _Required_: A comma-separated list of Spotify IDs. Maximum 50. \n

[Reference](https://developer.spotify.com/documentation/web-api/reference/library/save-tracks-user/)
""" ->
function library_save_track(track_ids)
    return spotify_request("me/tracks?ids=$track_ids", method="PUT")
end


text:
using QuantEcon

"""
The data for a consumption problem, including some default values.
"""
struct ConsumptionProblem{TF<:AbstractFloat}
    beta::TF
    y::Vector{TF}
    b0::TF
    P::Matrix{TF}
end

"""
Parameters
----------

beta : discount factor
P    : 2x2 transition matrix
y    : Array containing the two income levels
b0   : debt in period 0 (= state_1 debt level)
"""
function ConsumptionProblem(;
                 beta = 0.96,
                 y = [2.0, 1.5],
                 b0 = 3.0,
                 P = [0.8 0.2;
                      0.4 0.6])

    ConsumptionProblem(beta, y, b0, P)
end

"""
Computes endogenous values for the complete market case.

Parameters
----------

cp : instance of ConsumptionProblem

Returns
-------

    c_bar : constant consumption
    b1    : rolled over b0
    b2    : debt in state_2

associated with the price system

    Q = beta * P

"""

function consumption_complete(cp::ConsumptionProblem)

    beta, P, y, b0 = cp.beta, cp.P, cp.y, cp.b0  # Unpack

    y1, y2 = y          # extract income levels
    b1 = b0             # b1 is known to be equal to b0
    Q = beta * P        # assumed price system

    # Using equation (7) calculate b2
    b2 = (y2 - y1 - (Q[1, 1] - Q[2, 1] - 1) * b1)/(Q[1, 2] + 1 - Q[2, 2])

    # Using equation (5) calculae c_bar
    c_bar = y1 - b0 + ([b1 b2] * Q[1, :] )[1]

    return c_bar, b1, b2
end

"""
Computes endogenous values for the incomplete market case.

Parameters
----------

cp : instance of ConsumptionProblem
N_simul : Integer

"""
function consumption_incomplete(cp::ConsumptionProblem;
                                N_simul::Integer=150)

    beta, P, y, b0 = cp.beta, cp.P, cp.y, cp.b0  # Unpack
    # For the simulation define a quantecon MC class
    mc = MarkovChain(P)

    # Useful variables
    y = y''
    v = inv(eye(2) - beta * P) * y

    # Simulat state path
    s_path = simulate(mc, N_simul, init = 1)

    # Store consumption and debt path
    b_path, c_path = ones(N_simul + 1), ones(N_simul)
    b_path[1] = b0

    # Optimal decisions from (12) and (13)
    db = ((1 - beta) * v - y) / beta

    for (i, s) in enumerate(s_path)
        c_path[i] = (1 - beta) * (v - b_path[i] * ones(2, 1))[s, 1]
        b_path[i + 1] = b_path[i] + db[s, 1]
    end

    return c_path, b_path[1:end-1], y[s_path], s_path
end


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

mutable struct IoK8sApiAutoscalingV2beta2MetricValueStatus <: SwaggerModel
    averageUtilization::Any # spec type: Union{ Nothing, Int32 } # spec name: averageUtilization
    averageValue::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApiResourceQuantity } # spec name: averageValue
    value::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApiResourceQuantity } # spec name: value

    function IoK8sApiAutoscalingV2beta2MetricValueStatus(;averageUtilization=nothing, averageValue=nothing, value=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV2beta2MetricValueStatus, Symbol("averageUtilization"), averageUtilization)
        setfield!(o, Symbol("averageUtilization"), averageUtilization)
        validate_property(IoK8sApiAutoscalingV2beta2MetricValueStatus, Symbol("averageValue"), averageValue)
        setfield!(o, Symbol("averageValue"), averageValue)
        validate_property(IoK8sApiAutoscalingV2beta2MetricValueStatus, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        o
    end
end # type IoK8sApiAutoscalingV2beta2MetricValueStatus

const _property_map_IoK8sApiAutoscalingV2beta2MetricValueStatus = Dict{Symbol,Symbol}(Symbol("averageUtilization")=>Symbol("averageUtilization"), Symbol("averageValue")=>Symbol("averageValue"), Symbol("value")=>Symbol("value"))
const _property_types_IoK8sApiAutoscalingV2beta2MetricValueStatus = Dict{Symbol,String}(Symbol("averageUtilization")=>"Int32", Symbol("averageValue")=>"IoK8sApimachineryPkgApiResourceQuantity", Symbol("value")=>"IoK8sApimachineryPkgApiResourceQuantity")
Base.propertynames(::Type{ IoK8sApiAutoscalingV2beta2MetricValueStatus }) = collect(keys(_property_map_IoK8sApiAutoscalingV2beta2MetricValueStatus))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV2beta2MetricValueStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAutoscalingV2beta2MetricValueStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV2beta2MetricValueStatus }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV2beta2MetricValueStatus[property_name]

function check_required(o::IoK8sApiAutoscalingV2beta2MetricValueStatus)
    true
end

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


text:
function pos_datadep_register()
    register(DataDep("POS Model Weights",
        """
        The weights for POS Sequence Labelling Model.
        """,
        "https://github.com/JuliaText/TextAnalysis.jl/releases/download/v0.6.0/pos_weights.tar.xz",
        "b02e891ea913be6834ff67d6ecf2ddae6754d55509bb3d9c078dbfc7eed27988";
        post_fetch_method = function(fn)
            unpack(fn)
            dir = "weights"
            innerfiles = readdir(dir)
            mv.(joinpath.(dir, innerfiles), innerfiles)
            rm(dir)
        end
    ))

    register(DataDep("POS Model Dicts",
        """
        The character and words dict for POS Sequence Labelling Model.
        """,
        "https://github.com/JuliaText/TextAnalysis.jl/releases/download/v0.6.0/pos_model_dicts.tar.xz",
        "4d7fe8238ff0cfb92d195dfa745b4ed08f916d4707e3dbe27a1b3144c9282f41";
        post_fetch_method = function(fn)
            unpack(fn)
            dir = "model_dicts"
            innerfiles = readdir(dir)
            mv.(joinpath.(dir, innerfiles), innerfiles)
            rm(dir)
        end
    ))
end


text:
reload("src/OBC.jl")
using HDF5, JLD
using Stats
using Distributions

loc = expanduser("~/largeresearch/seq-data/rnaseq.jld")
@load loc luad_norm luad_norm_sampids luad_norm_geneids

blas_set_num_threads(1)

d = 10.0
iters = 10000
burn = 2000
thin = 50
D = 2
rounds = 1
numsamps = 100

srand(123)

samps = sample(1:size(luad_norm,2), numsamps)
genes = sample(1:size(luad_norm,1), D*rounds)

i=1
data = luad_norm[genes[(i-1)*D+1:i*D],samps]'

######################################################################
# Fit Model
######################################################################

#samplers = MPM.AMWGRecord[]
#for i=1:5
    #push!(samplers, MPM.mpm_model(data,burn=burn,thin=thin,d=d,usepriors=true))
    #sample(samplers[end],iters)
#end

dnostics = OBC.gelman_rubin(samplers)

#for i=1:5
    #sample(samplers[i],iters)
#end

dnostics2 = OBC.gelman_rubin(samplers)


text:
function second_order_cross(X)
    n = size(X, 1)
    k = size(X, 2)
    extra_k = div((1 + k) * k ,2)
    X_2 = zeros(n, extra_k)
    count = 0
    for i in 1:k
        for j in 1:k
            if i <= j
                count = count + 1
                X_2[:, count] = X[:, i] .* X[:, j]
            end
        end
    end
    return X_2
end

function third_order_cross(X)
    n = size(X, 1)
    k = size(X, 2)
    extra_k = div(k * (1 + k) * (2 + k),  6)
    X_3 = zeros(n, extra_k)
    count = 0
    for i in 1:k
        for j in 1:k
            for z in 1:k
                if (i <= j) && (j <= z)
                    count = count + 1
                    X_3[:, count] = X[:, i] .* X[:, j] .* X[:, z]
                end
            end
        end
    end
    return X_3
end

max_stars_repo_name: jungangl/forecast-machine-learning

max_issues_repo_name: jungangl/forecast-machine-learning

max_forks_repo_name: jungangl/forecast-machine-learning


text:
using OrdinaryDiffEq, DelayDiffEq, DiffEqDevTools, DiffEqBase, Test
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
using DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_2Dlinear, prob_ode_linear

using DiffEqProblemLibrary.DDEProblemLibrary: importddeproblems; importddeproblems()
using DiffEqProblemLibrary.DDEProblemLibrary: prob_dde_constant_1delay_ip

## Setup Tests

prob = prob_ode_linear

probs = Vector{ODEProblem}(undef, 2)
probs[1] = prob
probs[2] = prob_ode_2Dlinear

tspan = [0,1]
tspans = Vector{Vector{Int64}}(undef, 2)
tspans[1] = tspan; tspans[2] = tspan
abstols = 1 ./ 10 .^(3:10)
reltols = 1 ./ 10 .^(3:10)

setups = [Dict(:alg=>RK4());Dict(:alg=>Euler());Dict(:alg=>BS3());
          Dict(:alg=>Midpoint());Dict(:alg=>BS5());Dict(:alg=>DP5())]

t1 = @elapsed sol = solve(prob,RK4(),dt=1/2^(4))
t2 = @elapsed sol2 = solve(prob,setups[1][:alg],dt=1/2^(4))

@test (sol2[end] == sol[end])

## Shootout Tests

println("Shootout Tests")

shoot = Shootout(prob,setups,dt=1/2^(4))

#show(shoot)
#println(shoot)
shoot[end]

set = ShootoutSet(probs,setups;dt=1/2^(4))

#println(set[1])
#println(set[:])
set[end]
set[1][:]

## WorkPrecision Tests
println("WorkPrecision Tests")
println("Test DP5")
wp = WorkPrecision(prob,DP5(),abstols,reltols;name="Dormand-Prince 4/5")

wp[1]
wp[:]
wp[end]
#show(wp)

wp_set = WorkPrecisionSet(prob,abstols,reltols,setups;dt=1/2^4,numruns=2)

wp_set[1]
wp_set[:]
wp_set[end]
#println(wp_set)
#show(wp_set)
@test (minimum(diff(wp_set[2].errors).==0)) # The errors for a fixed timestep method should be constant

prob = prob_ode_2Dlinear

abstols = 1 ./10 .^(3:7)
reltols = 1 ./10 .^(0:4)

setups = [Dict(:alg=>DP5())
          Dict(:alg=>Tsit5(),
		  	   :abstols => 1 ./10 .^(4:7),
			   :reltols => 1 ./10 .^(1:4))]

sol = solve(prob,Vern7(),abstol=1/10^14,reltol=1/10^14)
test_sol1 = TestSolution(sol)
println("Test DP5 and Tsit5")
wp = WorkPrecisionSet(prob,abstols,reltols,setups;save_everystep=false)

@test length(wp[1]) == 5
@test length(wp[2]) == 4

function lotka(du,u,p,t)
  du[1] = 1.5 * u[1] - u[1]*u[2]
  du[2] = -3 * u[2] + u[1]*u[2]
end

prob = ODEProblem(lotka,[1.0;1.0],(0.0,10.0))

abstols = 1 ./10 .^(6:9)
reltols = 1 ./10 .^(3:6)
sol = solve(prob,Vern7(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)

setups = [Dict(:alg=>DP5())
          Dict(:alg=>Tsit5())
          Dict(:alg=>Vern6())
          ]
println("Test DP5, Tsit5, and Vern6")
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=20,maxiters=10000)

# Dual Problem

probs = [prob,prob_ode_2Dlinear]
setups = [Dict(:alg=>DP5(),:prob_choice => 1)
          Dict(:alg=>Tsit5(), :prob_choice => 2)
          Dict(:alg=>Vern6(), :prob_choice => 2)
          ]
println("Test DP5, Tsit5, and Vern6")
wp = WorkPrecisionSet(probs,abstols,reltols,setups;appxsol=[test_sol,nothing],save_everystep=false,numruns=20,maxiters=10000)
wp = WorkPrecisionSet(probs,abstols,reltols,setups;appxsol=[test_sol,test_sol1],save_everystep=false,numruns=20,maxiters=10000)

# DDE problem
prob = prob_dde_constant_1delay_ip

abstols = 1 ./10 .^(7:10)
reltols = 1 ./10 .^(4:7)
sol = solve(prob, MethodOfSteps(Vern9(), fpsolve = NLFunctional(; max_iter = 1000)); reltol=1e-8, abstol=1e-8)
test_sol = TestSolution(sol)

setups = [Dict(:alg => MethodOfSteps(BS3()))
          Dict(:alg => MethodOfSteps(Tsit5()))]
println("Test MethodOfSteps BS3 and Tsit5")
#Travis compile time issue
#wp = WorkPrecisionSet(prob, abstols, reltols, setups; appxsol = test_sol)
println("DDE Done")

__f(u,p,t) = 1.01*u
u0=1/2
tspan = (0.0,1.0)
prob = ODEProblem(__f,u0,tspan)
sol = solve(prob,Rodas4(),reltol=1e-8,abstol=1e-8)
setups =  [Dict(:alg=>RadauIIA5()),
				  Dict(:alg=>RosShamp4())]
shoot = Shootout(prob,setups;appxsol = TestSolution(sol))


text:
@testset "vertical processes" begin
    @test all(
        isapprox.(
            Wflow.rainfall_interception_gash(3.0, 0.11, 0.24, 18.0, 1.5, 4.0),
            (13.568000000000001, 4.0, 0.432, 1.5),
        ),
    )
    @test all(
        isapprox.(
            Wflow.rainfall_interception_modrut(8.6, 3.8, 1.5, 0.45, 2.8),
            (4.343, 3.87, 0.387, 0.0, 3.8, 2.043),
        ),
    )
    @test all(
        isapprox.(
            Wflow.acttransp_unsat_sbm(
                300.0,
                55.0,
                400.0,
                3.6,
                1.85,
                10.5,
                300.0,
                0.6,
                0.1,
                10.0,
                false,
            ),
            (55.0, 1.85, 3.6),
        ),
    )
    @test all(
        isapprox.(
            Wflow.infiltration(27.5, 0.2, 0.038, 8.9, 50.0, 5.0, 23.5, false, false),
            (23.5, 19.14814814814815, 4.351851851851852, 22.0, 5.5, 0.5),
        ),
    )
    @test all(
        isapprox.(
            Wflow.unsatzone_flow_layer(43.5, 256.0, 135.0, 12.6),
            (43.49983744545384, 0.00016255454615829025),
        ),
    )
    @test all(
        isapprox.(
            Wflow.unsatzone_flow_sbm(67.0, 950.0, 600.0, 250.0, 200.0, 0.6, 0.1),
            (19.142857142857146, 47.857142857142854),
        ),
    )
    @test all(
        isapprox.(
            Wflow.snowpack_hbv(201.5, 15.0, 30.1, 0.54, 2.0, 0.0, 0.0, 2.5, 0.10),
            (207.073, 20.707300000000004, 1.35, 18.819699999999997, 6.923),
        ),
    )
    @test Wflow.scurve(2.0, 0.0, 3.0, 2.5) ≈ 0.3325863502664285
    @test all(
        isapprox.(
            Wflow.glacier_hbv(0.35, 500.0, 9.5, 5.0, 0.0, 3.4, 0.2, Second(Day(1))),
            (8.835, 1.9, 484.9, 17.0),
        ),
    )
end


text:
##############################################################################
##
## Some things about Strings

# length: number of characters
# ncodeunits: Return the number of code units in a string (aking to index of vector).
# Not all such indices  are valid – they may not be the start of a character,.
# sizeof:  Size, in bytes, of the string str. Equal to the number of code units in str
# multiplied by the size, in bytes, of one code unit in str.

# lastindex: Return the last index of a collection
# nextinds(s, i):  return the index of the start of the character whose encoding starts after index i
# nextind(s, 0, N): return the index of the Nth character of s (or, if there are
# less than N characters, return ncodeunits(str) + (N - length(s))

##############################################################################

# This type allows to compute length once and for all
struct StringWithLength{T <: AbstractString} <: AbstractString
    s::T
    l::Int
end
string_with_length(s::AbstractString) = StringWithLength(s, length(s))
# Not really needed but avoid multi-encapsulation
string_with_length(s::StringWithLength) = s
Base.length(s::StringWithLength) = s.l
Base.iterate(s::StringWithLength) = iterate(s.s)
Base.iterate(s::StringWithLength, i::Integer) = iterate(s.s, i)
Base.nextind(s::StringWithLength, i::Int, n::Int = 1) = nextind(s.s, i, n)
Base.ncodeunits(s::StringWithLength) = ncodeunits(s.s)
Base.isvalid(s::StringWithLength, i::Int) = isvalid(s.s, i)
function reorder(s1::AbstractString, s2::AbstractString)
    s1 = string_with_length(s1)
    s2 = string_with_length(s2)
    (length(s1) <= length(s2)) ? (s1, s2) : (s2, s1)
end

function reorder(s1, s2)
    (length(s1) <= length(s2)) ? (s1, s2) : (s2, s1)
end

function common_prefix(s1, s2)
    l = 0
    for (ch1, ch2) in zip(s1, s2)
        ch1 != ch2 && break
        l += 1
    end
    return l
end


text:
mutable struct BezierPathSegment <: AbstractArray{Luxor.Point, 1}
    p1::Point
    cp1::Point
    cp2::Point
    p2::Point
end

Base.size(bps::BezierPathSegment) = (4, )
Base.length(bps::BezierPathSegment) = 4
Base.eltype(::Type{BezierPathSegment}) = Luxor.Point
Base.getindex(bps::BezierPathSegment, i::T) where T <: Integer = [bps.p1, bps.cp1, bps.cp2, bps.p2][i]
Base.setindex!(bp::BezierPathSegment, v, i::T) where T <: Integer = [bps.p1, bps.cp1, bps.cp2, bps.p2][i] = v
Base.IndexStyle(::Type{<:BezierPathSegment}) = IndexLinear()

# state is integer referring to one of the four points
function Base.iterate(bps::BezierPathSegment)
    return (bps.p1, 2)
end

function Base.iterate(bps::BezierPathSegment, s)
    if (s > length(bps))
        return
    end
    return bps[s], s+1
end

# function Base.show(io::IO, bps::BezierPathSegment)
#     println(io, "p1  $(bps.p1)  cp1 $(bps.cp1)")
#     println(io, "cp2  $(bps.cp2)  p2 $(bps.p2)")
# end

mutable struct BezierPath <: AbstractArray{BezierPathSegment, 1}
    segments::Vector{BezierPathSegment}
end

BezierPath(bps::BezierPathSegment) = BezierPath([bps])
BezierPath() = BezierPath([])

Base.size(bp::BezierPath) = size(bp.segments)
Base.length(bp::BezierPath) = length(bp.segments)
Base.eltype(::Type{BezierPath}) = BezierPathSegment
Base.getindex(bp::BezierPath, i::T) where T <: Integer = bp.segments[i]
Base.setindex!(bp::BezierPath, v, i::T) where T <: Integer = bp.segments[i] = v
Base.IndexStyle(::Type{<:BezierPath}) = IndexLinear()

Base.push!(bp::BezierPath, v::BezierPathSegment) = Base.push!(bp.segments, v)

# Base.show(io::IO, bp::BezierPath) =
#     print(io, "BezierPath with $(length(bp.segments)) segments")

"""
    bezier(t, A::Point, A1::Point, B1::Point, B::Point)

Return the result of evaluating the Bezier cubic curve function, `t` going from
0 to 1, starting at A, finishing at B, control points A1 (controlling A), and B1
(controlling B).
"""
bezier(t, A::Point, A1::Point, B1::Point, B::Point) =
      A * (1.0 - t)^3 +
      (A1 * 3t * (1.0 - t)^2) +
      (B1 * 3t^2 * (1.0 - t)) +
      (B * t^3)

"""
  bezier′(t, A::Point, A1::Point, B1::Point, B::Point)

Return the first derivative of the Bezier function.
"""
bezier′(t, A, A1, B1, B) = 3(1.0-t)^2 * (A1-A) + 6(1.0-t) * t * (B1-A1) + 3t^2 * (B-B1)

"""
    bezier′′(t, A::Point, A1::Point, B1::Point, B::Point)

Return the second derivative of Bezier function.
"""
bezier′′(t, A, A1, B1, B) = 6(1-t) * (B1-2A1+A) + 6t * (B-2B1+A1)

"""
    beziercurvature(t, A::Point, A1::Point, B1::Point, B::Point)

Return the curvature of the Bezier curve at `t` ([0-1]), given start and end
points A and B, and control points A1 and B1. The value (kappa) will typically
be a value between -0.001 and 0.001 for points with coordinates in the 100-500
range.

κ(t) is the curvature of the curve at point t, which for a parametric planar
curve is:

```math
\\begin{equation}
\\kappa = \\frac{\\mid \\dot{x}\\ddot{y}-\\dot{y}\\ddot{x}\\mid}
    {(\\dot{x}^2 + \\dot{y}^2)^{\\frac{3}{2}}}
\end{equation}
```

The radius of curvature, or the radius of an osculating circle at a point, is
1/κ(t). Values of 1/κ will typically be in the range -1000 to 1000 for points
with coordinates in the 100-500 range.

TODO Fix overshoot...

...The value of kappa can sometimes collapse near 0, returning NaN (and Inf for
radius of curvature).
"""
function beziercurvature(t, A::Point, A1::Point, B1::Point, B::Point)
    x′ = bezier′(t, A, A1, B1, B).x
    y′ = bezier′(t, A, A1, B1, B).y
    x′′ = bezier′′(t, A, A1, B1, B).x
    y′′ = bezier′′(t, A, A1, B1, B).y
    return ((x′ * y′′) - (y′ * x′′)) / (x′^2 + y′^2)^(3/2)
end

"""
    findbeziercontrolpoints(previouspt::Point,
        pt1::Point,
        pt2::Point,
        nextpt::Point;
        smooth_value=0.5)

Find the Bézier control points for the line between `pt1` and `pt2`, where the
point before `pt1` is `previouspt` and the next point after `pt2` is `nextpt`.
"""
function findbeziercontrolpoints(previouspt::Point, point1::Point, point2::Point, nextpt::Point;
    smoothing=0.5)

    xc1 = (previouspt.x + point1.x)/2.0 ; yc1 = (previouspt.y + point1.y)/2.0
    xc2 = (point1.x + point2.x)/2.0     ; yc2 = (point1.y + point2.y)/2.0
    xc3 = (point2.x + nextpt.x)/2.0     ; yc3 = (point2.y + nextpt.y)/2.0

    l1, l2, l3 = distance(previouspt, point1), distance(point2, point1), distance(nextpt, point2)
    k1 = l1 / (l1 + l2)   ; k2 = l2 / (l2 + l3)
    xm1 = xc1 + (xc2-xc1) * k1  ; ym1 = yc1 + (yc2-yc1) * k1
    xm2 = xc2 + (xc3-xc2) * k2  ; ym2 = yc2 + (yc3-yc2) * k2
    c1x = xm1 + (xc2-xm1) * smoothing + point1.x-xm1
    c1y = ym1 + (yc2-ym1) * smoothing + point1.y-ym1
    c2x = xm2 + (xc2-xm2) * smoothing + point2.x-xm2
    c2y = ym2 + (yc2-ym2) * smoothing + point2.y-ym2
    return (Point(c1x, c1y),  Point(c2x, c2y))
end

"""
    makebezierpath(pgon::Array{Point, 1}; smoothing=1.0)

Return a Bézier path (a BezierPath) that represents a polygon (an array of points). The Bézier
path is an array of segments (tuples of 4 points); each segment contains the
four points that make up a section of the entire Bézier path. `smoothing`
determines how closely the curve follows the polygon. A value of 0 returns a
straight-sided path; as values move above 1 the paths deviate further from
the original polygon's edges.
"""
function makebezierpath(pgon::Array{Point, 1}; smoothing=1.0)
    lpg = length(pgon)
    newpath = BezierPath()
    for i in 1:lpg
        cp1, cp2 = findbeziercontrolpoints(
            pgon[mod1(i - 1, lpg)],
            pgon[mod1(i,     lpg)],
            pgon[mod1(i + 1, lpg)],
            pgon[mod1(i + 2, lpg)], smoothing=smoothing)
        push!(newpath, BezierPathSegment(pgon[mod1(i, lpg)], cp1, cp2, pgon[mod1(i + 1, lpg)]))
    end
    return newpath
end

"""
    drawbezierpath(bezierpath::BezierPath, action=:none;
        close=true)

Draw the Bézier path, and apply the action, such as `:none`, `:stroke`, `:fill`,
etc. By default the path is closed.
"""
function drawbezierpath(bezierpath::BezierPath, action=:none; close=true)
    move(bezierpath[1].p1)
    for i in 1:length(bezierpath) - 1
        c = bezierpath[i]
        curve(c.cp1, c.cp2, c.p2)
    end
    if close == true
        c = bezierpath[end]
        curve(c.cp1, c.cp2, c.p2)
    end
    do_action(action)
end

"""
    drawbezierpath(bps::BezierPathSegment, action=:none;
        close=false)

Draw the Bézier path segment, and apply the action, such as `:none`, `:stroke`,
`:fill`, etc. By default the path is open.
"""
function drawbezierpath(bps::BezierPathSegment,
    action=:none; close=false)
    move(bps.p1)
    curve(bps.cp1, bps.cp2, bps.p2)
    if close == true
        closepath()
    end
    do_action(action)
end

drawbezierpath(bpsa::Array{BezierPathSegment, 1}, action=:none) =
   foreach(b -> drawbezierpath(b, action), bpsa)

"""
    beziertopoly(bpseg::BezierPathSegment; steps=10)

Convert a Bezier segment to a polygon (an array of points).
"""
function beziertopoly(bpseg::BezierPathSegment; steps=10)
    thepolygon = Point[]
    anchor1, control1, control2, anchor2 = bpseg.p1, bpseg.cp1, bpseg.cp2, bpseg.p2
    dx1 = control1.x - anchor1.x
    dy1 = control1.y - anchor1.y
    dx2 = control2.x - control1.x
    dy2 = control2.y - control1.y
    dx3 = anchor2.x - control2.x
    dy3 = anchor2.y - control2.y

    ss1 = 1.0/(steps + 1)
    ss2 = ss1 * ss1
    ss3 = ss2 * ss1

    p1 = 3.0 * ss1; p2 = 3.0 * ss2; p4 = 6.0 * ss2; p5 = 6.0 * ss3

    t1x = anchor1.x - control1.x * 2.0 + control2.x
    t1y = anchor1.y - control1.y * 2.0 + control2.y

    t2x = (control1.x - control2.x) * 3.0 - anchor1.x + anchor2.x
    t2y = (control1.y - control2.y) * 3.0 - anchor1.y + anchor2.y

    fx = anchor1.x ; fy = anchor1.y

    dfx = (control1.x - anchor1.x) * p1 + t1x * p2 + t2x * ss3
    dfy = (control1.y - anchor1.y) * p1 + t1y * p2 + t2y * ss3

    d2fx = t1x * p4 + t2x * p5 ; d2fy = t1y * p4 + t2y * p5
    d3fx = t2x * p5            ; d3fy = t2y * p5

    n = steps
    push!(thepolygon, Point(anchor1.x, anchor1.y))
    while n > 0
        fx   += dfx
        fy   += dfy
        dfx  += d2fx
        dfy  += d2fy
        d2fx += d3fx
        d2fy += d3fy
        push!(thepolygon, Point(fx, fy))
        n -= 1
    end
    # should we add the last point? no, do that downstream if required:
    # push!(thepolygon, Point(x4, y4))
    return thepolygon
end

"""
    bezierpathtopoly(bezierpath::BezierPath; steps=10)

Convert a Bezier path (an array of Bezier segments, where each segment is a
tuple of four points: anchor1, control1, control2, anchor2) to a polygon.

To make a Bezier path, use `makebezierpath()` on a polygon.

The `steps` optional keyword determines how many line sections are used for each
path.
"""
function bezierpathtopoly(bezierpath::BezierPath; steps=10)
    resultpoly = Point[]
    for bp in bezierpath
        p = beziertopoly(bp, steps=steps)
        append!(resultpoly, p)
    end
    if resultpoly[1] == resultpoly[end]
        pop!(resultpoly)
    end
    return resultpoly
end

"""
    pathtobezierpaths(
        ; flat=true)

Convert the current path (which may consist of one or more paths) to an array of
Bezier paths. Each Bezier path is, in turn, an array of path segments. Each path
segment is a tuple of four points. A straight line is converted to a Bezier segment
in which the control points are set to be the same as the end points.

If `flat` is true, use `getpathflat()` rather than `getpath()`.

# Example

This code draws the Bezier segments and shows the control points as "handles", like
a vector-editing program might.

```
@svg begin
    fontface("MyanmarMN-Bold")
    st = "goo"
    thefontsize = 100
    fontsize(thefontsize)
    sethue("red")
    fontsize(thefontsize)
    textpath(st)
    nbps = pathtobezierpaths()
    for nbp in nbps
        setline(.15)
        sethue("grey50")
        drawbezierpath(nbp, :stroke)
        for p in nbp
            sethue("red")
            circle(p[2], 0.16, :fill)
            circle(p[3], 0.16, :fill)
            line(p[2], p[1], :stroke)
            line(p[3], p[4], :stroke)
            if p[1] != p[4]
                sethue("black")
                circle(p[1], 0.26, :fill)
                circle(p[4], 0.26, :fill)
            end
        end
    end
end
```

"""
function pathtobezierpaths(; flat=true)
    flat ? (originalpath = getpath()) : (originalpath = getpathflat())
    # to store all the Bezier paths
    result = BezierPath[]
    # to store one Bezier path
    newbezpath = BezierPath()
    # to store a Bezier segment
    currentpos = Point(0, 0)
    if length(originalpath) > 0
        for e in originalpath
            if e.element_type == Cairo.CAIRO_PATH_MOVE_TO                # 0
                currentpos = Point(e.points[1], e.points[2])
            elseif e.element_type == Cairo.CAIRO_PATH_LINE_TO            # 1
                # add a straight line segment in the form of a bezsegment
                bezsegment = BezierPathSegment(currentpos,
                              currentpos,
                              Point(e.points[1], e.points[2]),
                              Point(e.points[1], e.points[2]))
                push!(newbezpath, bezsegment)
                currentpos = Point(e.points[1], e.points[2])
            elseif e.element_type == Cairo.CAIRO_PATH_CURVE_TO           # 2
                bezsegment = BezierPathSegment(currentpos,
                              Point(e.points[1], e.points[2]),
                              Point(e.points[3], e.points[4]),
                              Point(e.points[5], e.points[6]))
                currentpos = Point(e.points[5], e.points[6])
                push!(newbezpath, bezsegment)
            elseif e.element_type == Cairo.CAIRO_PATH_CLOSE_PATH         # 3
                # finish this path and store it away
                if length(newbezpath) > 0
                    push!(result, newbezpath)
                end
                # make new path
                newbezpath = BezierPath()
            else
                error("unknown CairoPathEntry " * repr(e.element_type))
                error("unknown CairoPathEntry " * repr(e.points))
            end
        end
        # tidy up, not everything gets closed
        if length(newbezpath) > 0
            push!(result, newbezpath)
        end
    end
    return result
end

# More Bezier madness

function _solvexy(a, b, c, d, e, f)
# solve linear equation ai + bj = c and di + ej = f
    j = (c - a / d * f) / (b - a * e / d)
    i = (c - (b * j)) / a
    return i, j
end

# basis functions
_b0(t) = (1 - t) ^ 3
_b1(t) = t * (1 - t) * (1 - t) * 3
_b2(t) = (1 - t) * t * t * 3
_b3(t) = (t ^ 3)

"""
    bezierfrompoints(startpoint::Point, pointonline1::Point,
        pointonline2::Point, endpoint::Point)

Given four points, return the Bezier curve that passes through all four points,
starting at `startpoint` and ending at `endpoint`. The two middle points of the
returned BezierPathSegment are the two control points that make the curve pass
through the two middle points supplied.
"""
function bezierfrompoints(startpoint::Point, pointonline1::Point, pointonline2::Point, endpoint::Point)
    # chord lengths
    c1 = sqrt((pointonline1.x - startpoint.x) * (pointonline1.x - startpoint.x) + (pointonline1.y - startpoint.y) * (pointonline1.y - startpoint.y))
    c2 = sqrt((pointonline2.x - pointonline1.x) * (pointonline2.x - pointonline1.x) + (pointonline2.y - pointonline1.y) * (pointonline2.y - pointonline1.y))
    c3 = sqrt((endpoint.x - pointonline2.x) * (endpoint.x - pointonline2.x) + (endpoint.y - pointonline2.y) * (endpoint.y - pointonline2.y))
    # best guess for t
    t1 = c1 / (c1 + c2 + c3)
    t2 = (c1 + c2) / (c1 + c2 + c3)
    (x1, x2) = _solvexy(_b1(t1), _b2(t1), pointonline1.x - (startpoint.x * _b0(t1)) - (endpoint.x * _b3(t1)), _b1(t2), _b2(t2), pointonline2.x - (startpoint.x * _b0(t2)) - (endpoint.x * _b3(t2)))
    (y1, y2) = _solvexy(_b1(t1), _b2(t1), pointonline1.y - (startpoint.y * _b0(t1)) - (endpoint.y * _b3(t1)), _b1(t2), _b2(t2), pointonline2.y - (startpoint.y * _b0(t2)) - (endpoint.y * _b3(t2)))
    return BezierPathSegment(startpoint, Point(x1, y1), Point(x2, y2), endpoint)
end

"""
    bezierfrompoints(ptslist::Array{Point, 1})

Given four points, return the Bezier curve that passes through all four points.
"""
bezierfrompoints(ptslist::Array{Point, 1}) = bezierfrompoints(ptslist...)

"""
    bezierstroke(point1, point2, width=0.0)

Return a BezierPath, a stroked version of a straight line between two points.

It wil have 2 or 6 Bezier path segments that define a brush or pen shape. If
width is 0, the brush shape starts and ends at a point. Otherwise the brush shape
starts and ends with the thick end.

To draw it, use eg `drawbezierpath(..., :fill)`.
"""
function bezierstroke(p1::Point, p2::Point, width=0.0)
    bezpath = BezierPath()
    if isapprox(width, 0.0)
        # simple wavy stroke starting ending at point
        push!(bezpath, BezierPathSegment(p1, p1, p2, p2))
        push!(bezpath, BezierPathSegment(p2, p2, p1, p1))
        result = setbezierhandles.(bezpath, angles=[0.1, -0.1],handles=[0.3, 0.3])
    else
        # stroke with broad opening and closing
        # make the paths, then give them some control power
        cp1 = perpendicular(p1, p2, -width)
        seg = BezierPathSegment(p1, p1, cp1, cp1)
        adjustedseg = setbezierhandles(seg, angles=[0.1, -0.1], handles=[0.3, 0.3])
        push!(bezpath, adjustedseg)

        cp2 = perpendicular(p2, p1, width)
        seg = BezierPathSegment(cp1, cp1, cp2, cp2)
        adjustedseg = setbezierhandles(seg, angles=[0.1, -0.1], handles=[0.3, 0.3])
        push!(bezpath, adjustedseg)

        seg = BezierPathSegment(cp2, cp2, p2, p2)
        adjustedseg = setbezierhandles(seg, angles=[0.1, -0.1], handles=[0.3, 0.3])
        push!(bezpath, adjustedseg)

        # TODO the segments should be collinear, perhaps?
        cp3 = perpendicular(p2, p1, -width)
        seg = BezierPathSegment(p2, p2, cp3, cp3)
        adjustedseg = setbezierhandles(seg, angles=[0.1, -0.1], handles=[0.3, 0.3])
        push!(bezpath, adjustedseg)

        cp4 = perpendicular(p1, p2, width)
        seg = BezierPathSegment(cp3, cp3, cp4, cp4)
        adjustedseg = setbezierhandles(seg, angles=[0.1, -0.1], handles=[0.3, 0.3])
        push!(bezpath, adjustedseg)

        seg = BezierPathSegment(cp4, cp4, p1, p1)
        adjustedseg = setbezierhandles(seg, angles=[0.1, -0.1], handles=[0.3, 0.3])
        push!(bezpath, adjustedseg)
        result = bezpath
    end
    return result
end

"""
    setbezierhandles(bps::BezierPathSegment;
            angles  = [0.05, -0.1],
            handles = [0.3, 0.3])

Return a new Bezier path segment with new locations for the Bezier control
points in the Bezier path segment `bps`.

`angles` are the two angles that the "handles" make with the line direciton.

` handles` are the lengths of the "handles". 0.3 is a typical value.
"""
function setbezierhandles(bps::BezierPathSegment;
            angles=[0.05, -0.1],
            handles=[0.3, 0.3])
    d = distance(bps[1], bps[4])
    s = slope(bps[1], bps[4])
    bezhandle1 = bps[1] + polar(d * handles[1], s - angles[1])
    s = slope(bps[4], bps[1])
    bezhandle2 = bps[4] + polar(d * handles[2], s + angles[2])
    return BezierPathSegment(bps[1], bezhandle1, bezhandle2, bps[4])
end

"""
    setbezierhandles(bezpath::BezierPath;
            angles=[0 .05, -0.1],
            handles=[0.3, 0.3])

Return a new Bezierpath with new locations for the Bezier control points in
every Bezier path segment of the BezierPath in `bezpath`.

`angles` are the two angles that the "handles" make with the line direciton.

` handles` are the lengths of the "handles". 0.3 is a typical value.
"""
function setbezierhandles(bezpath::BezierPath;
        angles=[0.05, -0.1],
        handles=[0.4, 0.4])
    newbezpath = BezierPath()
    for bps in bezpath
        push!(newbezpath, setbezierhandles(bps, angles=angles, handles=handles))
    end
    return newbezpath
end

"""
    shiftbezierhandles(bps::BezierPathSegment;
        angles=[0.1, -0.1], handles=[1.1, 1.1])

Return a new BezierPathSegment that modifies the Bezier path in `bps` by moving
the control handles. The values in `angles` increase the angle of the handles;
the values in `handles` modifies the lengths: 1 preserves the length, 0.5 halves
the length of the  handles, 2 doubles them.
"""
function shiftbezierhandles(bps::BezierPathSegment;
        angles  = [0.1, -0.1],
        handles = [0.1, 0.1])
    p1, cp1, cp2, p2 = bps
    # find slope of curve at the end points
    spt1 = bezier′(0.0, p1, cp1, cp2, p2)
    spt2 = bezier′(1.0, p1, cp1, cp2, p2)
    # angle of handle to curve
    s1 = slope(p1, cp1)
    s2 = slope(p2, cp2)
    # length of handles
    l1 = distance(p1, cp1)
    l2 = distance(p2, cp2)
    # new angle
    s1 += angles[1]
    s2 += angles[2]
    # new handle lengths
    l1 = l1 * handles[1]
    l2 = l2 * handles[2]
    newcontrolpoint1 = p1 + polar(l1, s1)
    newcontrolpoint2 = p2 + polar(l2, s2)
    return BezierPathSegment(p1, newcontrolpoint1, newcontrolpoint2, p2)
end

"""
    brush(pt1, pt2, width=10;
        strokes=10,
        minwidth=0.01,
        maxwidth=0.03,
        twist = -1,
        lowhandle  = 0.3,
        highhandle = 0.7,
        randomopacity = true,
        tidystart = false,
        action = :fill,
        strokefunction = (nbpb) -> nbpb))

Draw a composite brush stroke made up of some randomized individual filled
Bezier paths.

`strokefunction` allows a function to process a BezierPathSegment or do other things before it's drawn.

!!! note

    There is a lot of randomness in this function. Results are unpredictable.
"""
function brush(pt1, pt2, width=10;
        strokes=10,
        minwidth=0.01,
        maxwidth=0.03,
        twist = -1,
        lowhandle  = 0.3,
        highhandle = 0.7,
        randomopacity = true,
        tidystart = false,
        action = :fill,
        strokefunction = (nbpb) -> nbpb)
    @layer begin
        sl = slope(pt1, pt2)
        n = distance(pt1, pt2)
        translate(pt1)
        rotate(sl - pi/2)
        widthsteps = (maxwidth - minwidth)/10
        for j in 1:strokes
            shiftedline = [O + (0, 0), O + (0, n)]
            if tidystart
                shiftedline .+= Point(rand(-width/2:width/2), 0)
            else
                shiftedline .+= Point(rand(-width/2:width/2), rand(-width/2:width/2))
            end
            pbp = bezierstroke(shiftedline[1],
                shiftedline[2],
                rand(Bool) ? 0.0 : rand(minwidth:widthsteps:maxwidth)
                )
            for bps in pbp
                randomopacity ? setopacity(rand(0.3:0.1:0.9)) : setopacity(1.0)
                nbpb = shiftbezierhandles(bps,
                    angles  = [twist * rand(minwidth:widthsteps:maxwidth), -twist * rand(minwidth:widthsteps:maxwidth)],
                    handles = [rand(lowhandle:0.01:highhandle, 2)...]
                    )
                strokefunction(nbpb)
                drawbezierpath(nbpb, action, close=false)
            end
        end
    end
end


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

## initdefs.jl - initialization and runtime management definitions

"""
    PROGRAM_FILE

A string containing the script name passed to Julia from the command line. Note that the
script name remains unchanged from within included files. Alternatively see
[`@__FILE__`](@ref).
"""
global PROGRAM_FILE = ""

"""
    ARGS

An array of the command line arguments passed to Julia, as strings.
"""
const ARGS = String[]

"""
    exit(code=0)

Stop the program with an exit code. The default exit code is zero, indicating that the
program completed successfully. In an interactive session, `exit()` can be called with
the keyboard shortcut `^D`.
"""
exit(n) = ccall(:jl_exit, Cvoid, (Int32,), n)
exit() = exit(0)

const roottask = current_task()

is_interactive = false

"""
    isinteractive() -> Bool

Determine whether Julia is running an interactive session.
"""
isinteractive() = (is_interactive::Bool)

## package depots (registries, packages, environments) ##

const DEPOT_PATH = String[]

function init_depot_path(BINDIR::String = Sys.BINDIR)
    if haskey(ENV, "JULIA_DEPOT_PATH")
        depots = split(ENV["JULIA_DEPOT_PATH"], Sys.iswindows() ? ';' : ':')
        append!(empty!(DEPOT_PATH), map(expanduser, depots))
    else
        push!(empty!(DEPOT_PATH), joinpath(homedir(), ".julia"))
        push!(DEPOT_PATH, abspath(BINDIR, "..", "local", "share", "julia"))
        push!(DEPOT_PATH, abspath(BINDIR, "..", "share", "julia"))
    end
end

## LOAD_PATH ##

# split on `:` (or `;` on Windows)
# first empty entry is replaced with DEFAULT_LOAD_PATH, the rest are skipped
# entries starting with `@` are named environments:
#  - the first three `#`s in a named environment are replaced with version numbers
#  - `@stdlib` is a special name for the standard library and expands to its path

# if you want a current env setup, use direnv and
# have your .envrc do something like this:
#
#   export JULIA_LOAD_PATH="$(pwd):$JULIA_LOAD_PATH"
#
# this will inherit an existing JULIA_LOAD_PATH value or if there is none, leave
# a trailing empty entry in JULIA_LOAD_PATH which will be replaced with defaults.

const DEFAULT_LOAD_PATH = ["@v#.#", "@stdlib"]

"""
    LOAD_PATH

An array of paths for `using` and `import` statements to consdier as project
environments or package directories when loading code. See Code Loading.
"""
const LOAD_PATH = copy(DEFAULT_LOAD_PATH)

function current_env(dir::AbstractString = pwd())
    # look for project file in current dir and parents
    home = homedir()
    while true
        for proj in project_names
            file = joinpath(dir, proj)
            isfile_casesensitive(file) && return file
        end
        # bail at home directory or top of git repo
        (dir == home || ispath(joinpath(dir, ".git"))) && break
        old, dir = dir, dirname(dir)
        dir == old && break
    end
end

function parse_load_path(str::String)
    envs = String[]
    isempty(str) && return envs
    first_empty = true
    for env in split(str, Sys.iswindows() ? ';' : ':')
        if isempty(env)
            first_empty && append!(envs, DEFAULT_LOAD_PATH)
            first_empty = false
        elseif env == "@"
            dir = current_env()
            dir !== nothing && push!(envs, dir)
        else
            push!(envs, env)
        end
    end
    return envs
end

function init_load_path(BINDIR::String = Sys.BINDIR)
    if Base.creating_sysimg
        load_path = ["@stdlib"]
    elseif haskey(ENV, "JULIA_LOAD_PATH")
        load_path = parse_load_path(ENV["JULIA_LOAD_PATH"])
    else
        load_path = DEFAULT_LOAD_PATH
    end
    append!(empty!(LOAD_PATH), load_path)
end

const atexit_hooks = []

"""
    atexit(f)

Register a zero-argument function `f()` to be called at process exit. `atexit()` hooks are
called in last in first out (LIFO) order and run before object finalizers.
"""
atexit(f::Function) = (pushfirst!(atexit_hooks, f); nothing)

function _atexit()
    for f in atexit_hooks
        try
            f()
        catch err
            show(stderr, err)
            println(stderr)
        end
    end
end

max_stars_repo_name: perrutquist/julia

max_issues_repo_name: perrutquist/julia

max_forks_repo_name: perrutquist/julia


text:
using Statistics, StatsBase, Random

function MNmodel(θ, rndseed=1234)
    Random.seed!(rndseed)
    n = 1000
    μ1, μ2, σ1, σ2, prob = θ
    d1=randn(n).*σ1 .+ μ1
    d2=randn(n).*(σ1+σ2) .+ (μ1 - μ2) # second component lower mean and higher variance
    ps=rand(n).<prob
    data=zeros(n)
    data[ps].=d1[ps]
    data[.!ps].=d2[.!ps]
    return data
end

function auxstat(θ, reps)
    data = [MNmodel(θ, rand(1:Int64(1e12))) for i = 1:reps]  # reps draws of data
    auxstat.(data)
end

function auxstat(data)
    r = 0.0 : 0.1 : 1.0
    sqrt(Float64(size(data,1)))*vcat(mean(data), std(data), skewness(data), kurtosis(data), quantile.(Ref(data),r))
end

function TrueParameters()
    [1.0, 1.0, 0.2, 1.8, 0.4] # first component N(1,0.2) second component N(0,2)
end

function PriorSupport()
    lb = [0.0, 0.0, 0.0, 0.0, 0.05] # there is always at least 5% prob for each component
    ub = [3.0, 3.0, 1.0, 3.0, 0.95]
    lb,ub
end

function InSupport(θ)
    lb,ub = PriorSupport()
    all(θ .>= lb) & all(θ .<= ub)
end

function PriorMean()
    lb,ub = PriorSupport()
    (ub - lb) ./ 2.0
end

function PriorDraw()
    lb, ub = PriorSupport()
    θ = (ub-lb).*rand(size(lb,1)) + lb
end

function Prior(θ)
    lb,ub = PriorSupport()
    a = 0.0
    if (all(θ .>= lb) & all(θ .<= ub))
        a = 1.0
    end
    return a
end

max_stars_repo_name: mcreel/ILSNM

max_issues_repo_name: mcreel/ILSNM

max_forks_repo_name: mcreel/ILSNM


text:
using OrdinaryDiffEq, ParameterizedFunctions, Plots, DiffEqDevTools, Sundials, LSODA
gr() #gr(fmt=:png)
f = @ode_def Orego begin
  dy1 = p1*(y2+y1*(1-p2*y1-y2))
  dy2 = (y3-(1+y1)*y2)/p1
  dy3 = p3*(y1-y3)
end p1 p2 p3

p = [77.27,8.375e-6,0.161]
prob = ODEProblem(f,[1.0,2.0,3.0],(0.0,30.0),p)
sol = solve(prob,Rodas5(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)

# group 1
setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>lsoda())]

# High Tolerance
abstols = 1 ./ 10 .^ (5:8)
reltols = 1 ./ 10 .^ (1:4);

wp = WorkPrecisionSet(prob,abstols,reltols,setups;dense = false,verbose=false,
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
DIAGRAMS["Orego-HighTol-Group1"] = plot(wp)

# group 2
setups = [Dict(:alg=>Rodas4()),
          Dict(:alg=>Rodas5()),
          Dict(:alg=>Rodas5P()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>RadauIIA5()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>lsoda())]

# Low Tolerance
abstols = 1 ./ 10 .^ (7:13)
reltols = 1 ./ 10 .^ (4:10)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;verbose=false,
                      dense=false,appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
DIAGRAMS["Orego-LowTol-Group2"] = plot(wp)


text:
function cnt_total(img::AbstractArray, periodic::Bool)
    cnt_total(size(img), periodic)
end

function cnt_total(img::CuArray, periodic::Bool)
    cu(cnt_total(size(img), periodic))
end

function cnt_total(img_size::Tuple{Vararg{Int}}, periodic::Bool)
    if periodic
        result = *(img_size...)
    else
        result = zeros(Int, img_size)
        for ix in CartesianIndices(img_size)
            dir = Tuple(ix) .- 1
            result[ix] = *((img_size .- dir)...)
        end
    end
    result
end

function gradient_norm(img, kernelfactor=KernelFactors.sobel)
    dimgs = imgradients(img, kernelfactor)

    map((x...) -> norm(x), dimgs...)
end

function gradient_norm(img::CuArray, kernelfactor=KernelFactors.sobel)
    imgCPU = Array(img)
    cu(gradient_norm(imgCPU, kernelfactor))
end

"""
    cross_correlation!(f, g)

Compute cross correlation using FFT for xᵢ ≥ 0
and periodic boundary conditions.
Store result in `f`.
"""
function cross_correlation!(f, g)
    fft!(f)
    fft!(g)
    @. f = conj(f) * g
    ifft!(f)
end

"""
    self_correlation!(f)

Compute self correlation using FFT for xᵢ ≥ 0
and periodic boundary conditions.
Store result in `f`.
"""
function self_correlation!(f)
    fft!(f)
    @. f = abs2(f)
    ifft!(f)
end

function size_slice(img::AbstractArray{T, 3}, dim::Int) where T
    if dim == 1
        size(img)[2:3]
    elseif dim == 2
        (size(img, 1), size(img, 3))
    elseif dim == 3
        size(img)[1:2]
    else
        error("dim is not in {1, 2, 3}")
    end
end

function size_slice(img::AbstractArray{T, 2}, dim::Int) where T
    if dim == 1
        (size(img, 2),)
    elseif dim == 2
        (size(img, 1),)
    else
        error("dim is not in {1, 2}")
    end
end

function size_slice(img::AbstractVector, dim::Int)
    if dim == 1
        ()
    else
        error("dim is not 1")
    end
end


text:
abstract type ValueHistory end
abstract type UnivalueHistory{I} <: ValueHistory end
abstract type MultivalueHistory <: ValueHistory end

Base.push!(history::UnivalueHistory, iteration, value) =
    throw(ArgumentError("The specified arguments are of incompatible type"))

# length(history::ValueHistory) =  error()
# enumerate(history::ValueHistory) = error()
# get(history::ValueHistory) = error()
# first(history::ValueHistory) = error()
# last(history::ValueHistory) = error()

"""
    drop!(hist, it)

Deletes all values stored for iterations greater than `it`. If
`it>first(last(hist))` then it does nothing.
"""
function drop!(history::UnivalueHistory, it)
    iter, vals = get(history)
    n = findlast(x->x<=it, iter)
    isnothing(n) && return last(iter)

    resize!(history, n)
    return history.lastiter
end


text:
mutable struct LinearAllocator
    buffer::BufferBlock{MemoryBlock}
    last_offset::Int
    base_ptr::Ptr{Cvoid}
end

device(allocator::LinearAllocator) = device(allocator.buffer)

function LinearAllocator(device, size)
    buffer = BufferBlock(device, size; usage = Vk.BUFFER_USAGE_SHADER_DEVICE_ADDRESS_BIT)
    allocate!(buffer, MEMORY_DOMAIN_DEVICE_HOST)
    LinearAllocator(buffer, 0, map(memory(buffer)))
end

device_address(la::LinearAllocator) = device_address(la.buffer)

_collect_data(data::AbstractArray) = collect(data)
_collect_data(data::AbstractString) = string(data)
_collect_data(data) = data
Base.copyto!(la::LinearAllocator, data, alignment = 8) = _copyto!(la, _collect_data(data), alignment)

function _copyto!(la::LinearAllocator, data, alignment) where {T}
    offset = get_offset(la, data, alignment)
    ptrcopy!(la.base_ptr + offset, data)
    size = sizeof(data)
    la.last_offset = offset + size
    @view la.buffer[offset:offset + size]
end

function get_offset(la::LinearAllocator, data, alignment)
    offset = alignment * cld(la.last_offset, alignment)
    offset + sizeof(data) ≤ size(la.buffer) || error("Data does not fit in memory (available: $(size(la.buffer) - offset), requested: $(sizeof(data))).")
    offset
end

function reset!(la::LinearAllocator)
    la.last_offset = 0
end


text:
module TestRecording

using EventTracker
using EventTracker: EVENT_INTERVAL, EVENT_POINT
using EventTrackerBase: LocationProxy, EventTrackerBase
using Test

function test_location_hack()
    swb = EventTrackerBase
    @test swb._check_location1() ===  swb._check_location1()::LocationProxy
    @test swb._check_location2() ===  swb._check_location2()::LocationProxy
    @test swb._check_location3() ===  swb._check_location3()::LocationProxy
    @test Base.PkgId(swb._check_location1()._module) === Base.PkgId(swb)
    @test swb._check_location1().tag === swb.NOTAG
    @test swb._check_location1().line.line + 1 == swb._check_location2().line.line
end

function record_once()
    ans = @recordinterval :record_once begin
        y = 1
    end
    return (ans, y)
end

function test_record_once()
    @test record_once() == (1, 1)
end

function no_tag()
    ans = @recordinterval y = 1
    return (ans, y)
end

function test_no_tag()
    @test no_tag() == (1, 1)
end

record_point_once() = @recordpoint

function test_record_point_once()
    @test record_point_once() === nothing
end

function test_record_in_spawn()
    tasks = map(1:10) do i
        Threads.@spawn begin
            ans = @recordinterval :record_in_spawn begin
                y = i * 10
            end
            return (ans, y)
        end
    end
    results = fetch.(tasks)
    @test first.(results) == (1:10) .* 10
    @test last.(results) == (1:10) .* 10
end

end  # module