text:
const MAX_M = 1000

struct PreComp{T}
    BINS::SparseArrays.SparseMatrixCSC{T,Int64}
    LAGUERRE::SparseArrays.SparseMatrixCSC{T,Int64}
    LAGUERRE2::SparseArrays.SparseMatrixCSC{T,Int64}
    FACTS::Vector{T}
end

function bigfloat_precomputations(m)
    setprecision(2048)
    m = big(m)
    BINS = zeros(BigInt,(m,m))
    FACTS = zeros(BigInt,m)
    LAGUERRE = zeros(BigFloat,(m,m))
    LAGUERRE2 = deepcopy(LAGUERRE)
    for i in 1:m
        FACTS[i] = factorial(i-1)
    end
    for i in 1:m, j in 1:m
        BINS[j,i] = binomial(i-1,j-1)
        LAGUERRE[j,i] = BINS[j,i]/FACTS[j]*(-big(2))^(j-1)
        LAGUERRE2[j,i] = LAGUERRE[j,i] * sqrt(big(2)) * FACTS[j]
    end
    LAGUERRE = SparseArrays.sparse(LAGUERRE)
    LAGUERRE2 = SparseArrays.sparse(LAGUERRE2)
    BINS = SparseArrays.sparse(BigFloat.(BINS))
    FACTS = BigFloat.(FACTS)
    PreComp(BINS,LAGUERRE,LAGUERRE2,FACTS)
end

const big_P = bigfloat_precomputations(MAX_M)

const Precomp_Dict = Dict{Tuple{Int64,DataType},Any}((MAX_M,BigFloat)=>big_P)

"""
    get_precomp(type,m)

Query and eventualy compute and store the precomputations for laguerre things.
"""
function get_precomp(type::DataType,m::T) where {T <: Int}
    # The getter should :
    if (m,type) in keys(Precomp_Dict)
        return Precomp_Dict[(m,type)]::PreComp{type}
    end
    # if this is not the case, we will just construct it :
    if m > MAX_M
        error("You ask for too much")
    end
    P = big_P
    Precomp_Dict[(m,type)] = PreComp(type.(P.BINS[1:m,1:m]),type.(P.LAGUERRE[1:m,1:m]),type.(P.LAGUERRE2[1:m,1:m]),type.(P.FACTS[1:m]))
    return Precomp_Dict[(m,type)]::PreComp{type}
end


text:
using GAMSFiles, QuadDIRECT, ProgressMeter

struct ProblemInfo
    smoothness::String
    convexity::String
    nvars::Int
    minval::Float64
end

gamsdir = joinpath(@__DIR__, "gams")
if !isdir(gamsdir)
    mkdir(gamsdir)
end
# Download the problems studied by L. M. Rios and N. V. Sahinidis, "Derivative-free optimization:
# a review of algorithms and comparison of software implementations."
# Journal of Global Optimization,  July 2013, Volume 56, Issue 3, pp 1247–1293.
site = "http://thales.cheme.cmu.edu/dfo/comparison"
gamsdirs = ("convex_nonsmooth_gams",
            "convex_smooth_gams",
            "nonconvex_nonsmooth_gams",
            "nonconvex_smooth_gams")
cd(gamsdir) do
    for dir in gamsdirs
        subdir = joinpath(gamsdir, dir)
        if !isdir(subdir)
            pathfile = subdir*".zip"
            download(joinpath(site, dir*".zip"), pathfile)
            mkdir(subdir)
            run(`unzip $pathfile -d $subdir`)
        end
    end
end
problemdatafile = joinpath(gamsdir, "problemdata.zip")
if !isfile(problemdatafile)
    download(joinpath(site, "problemdata.zip"), problemdatafile)
    cd(gamsdir) do
        mkdir("problemdata")
        run(`unzip $problemdatafile -d problemdata`)
    end
end
answerfile = joinpath(gamsdir, "answers.tsv")
if !isfile(answerfile)
    error("must parse the answer file using tips from https://eagereyes.org/data/scrape-tables-using-google-docs using URL ", joinpath(site, "comp.html"))
end
answers_array = readdlm(answerfile, '\t')
@assert(answers_array[1,2] == "problem")
@assert(answers_array[1,3] == "smoothness")
@assert(answers_array[1,4] == "convexity")
@assert(answers_array[1,5] == "variables")
@assert(answers_array[1,7] == "solution")
answers = Dict{String,ProblemInfo}()
for i = 2:size(answers_array, 1)
    prob = answers_array[i,2]
    info = ProblemInfo(answers_array[i,3], answers_array[i,4], answers_array[i,5], answers_array[i,7])
    answers[prob] = info
end

function read_problemdata(file)
    lower, upper, x0 = open(file) do io
        values = split(read(io, String))
        k = 1
        nvars = parse(Int, values[1])
        lower = Vector{Float64}(undef, nvars)
        upper = similar(lower)
        x0 = similar(lower)
        for i = 1:nvars
            lower[i] = parse(Float64, values[k+=1])
        end
        for i = 1:nvars
            upper[i] = parse(Float64, values[k+=1])
        end
        for i = 1:nvars
            x0[i] = parse(Float64, values[k+=1])
        end
        @assert(k == length(values))
        return lower, upper, x0
    end
    return lower, upper, x0
end

stopping_criterion(v) = max(1.01*v, v+0.01)
stopping_criterion(info::ProblemInfo) = stopping_criterion(info.minval)

# Debugging utility: run a single GAMS file
function gams1(filename::AbstractString; kwargs...)
    dirname, fname = splitdir(filename)
    basename, ext = splitext(fname)
    lower, upper, x0 = read_problemdata(joinpath(gamsdir, "problemdata", basename)*".problem.data")
    gams = parsegams(filename)
    modex, axs = parsegams(Module, basename, gams)
    mod = eval(modex)
    f = getfield(mod, :objective)
    fwrap = LoggedFunction(f)
    root, x0 = @eval analyze($fwrap, $x0, $lower, $upper; rtol=0, fvalue=stopping_criterion($(answers[basename])), $(kwargs...))
    return root, x0, fwrap
end

failures = String[]
successes = Dict{String,Any}()

# Because of https://github.com/JuliaLang/julia/issues/25927
excludes = ["nonconvex_smooth_gams/nnls.gms",
            "nonconvex_smooth_gams/chebyqad.gms"]

mode = "optimize"   # any choice other than "optimize" tests that the objective can be evaluated

terminating = false
for dir in gamsdirs
    terminating && break
    subdir = joinpath(gamsdir, dir)
    cd(subdir) do
        files = readdir()
        @showprogress 1 "Working on $subdir" for file in files
            endswith(file, ".gms") || (println("skipping ", file); continue)
            fullfile = joinpath(subdir, file)
            exclude = false
            for excl in excludes
                if endswith(fullfile, excl)
                    exclude = true
                    println("excluding ", fullfile)
                end
            end
            exclude && continue
            basename, ext = splitext(file)
            lower, upper, x0 = read_problemdata(joinpath(gamsdir, "problemdata", basename)*".problem.data")
            local x0g
            try
                modex, axs = parsegams(Module, file)
                x0g = axs[1]  # initialization provided by the file
                if isdefined(Main, Symbol(basename))
                    mod = getfield(Main, Symbol(basename))
                else
                    mod = eval(modex)
                end
                f = getfield(mod, :objective)
                # # FIXME: we need a better way of setting the splits
                if mode == "optimize"
                    fwrap = LoggedFunction(f)
                    root, x0 = @eval analyze($fwrap, $x0, $lower, $upper; rtol=0, fvalue=stopping_criterion($(answers[basename])))
                    successes[file] = fwrap.values
                else
                    val = Base.invokelatest(f, x0)
                    if !isfinite(val)
                        # Also try any initialization present in the GAMS file
                        val = Base.invokelatest(f, x0g)
                        if isfinite(val)
                            println(fullfile, ": finite result from file initialization but not problemdata")
                        else
                            println("\nnot finite on ", fullfile)
                            println()
                            push!(failures, fullfile)
                        end
                    end
                    if isfinite(val)
                        successes[fullfile] = val
                    end
                end
            catch ex
                if ex isa InterruptException
                    terminating = true
                    break
                end
                push!(failures, fullfile)
                println("\nfailed on ", fullfile)
                println("problemdata axes ", axes(x0), ", and from the file it's ", axes(x0g))
                showerror(STDOUT, ex)
                println()
            end
        end
    end
end

ext: jl


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

module XoshiroSimd
# Getting the xoroshiro RNG to reliably vectorize is somewhat of a hassle without Simd.jl.
import ..Random: TaskLocalRNG, rand, rand!, Xoshiro, CloseOpen01, UnsafeView,
                 SamplerType, SamplerTrivial
using Base: BitInteger_types
using Core.Intrinsics: llvmcall

# Vector-width. Influences random stream.
xoshiroWidth() = Val(8)
# Simd threshold. Influences random stream.
simdThreshold(::Type{T}) where T = 64
simdThreshold(::Type{Bool}) = 640

@inline _rotl45(x::UInt64) = (x<<45)|(x>>19)
@inline _shl17(x::UInt64) = x<<17
@inline _rotl23(x::UInt64) = (x<<23)|(x>>41)
@inline _plus(x::UInt64,y::UInt64) = x+y
@inline _xor(x::UInt64,y::UInt64) = xor(x,y)
@inline _and(x::UInt64, y::UInt64) = x & y
@inline _or(x::UInt64, y::UInt64) = x | y
@inline _lshr(x, y::Int32) = _lshr(x, y % Int64)
@inline _lshr(x::UInt64, y::Int64) = llvmcall("""
    %res = lshr i64 %0, %1
    ret i64 %res
    """,
    UInt64,
    Tuple{UInt64, Int64},
    x, y)

@inline _bits2float(x::UInt64, ::Type{Float64}) = reinterpret(UInt64, Float64(x >>> 11) * 0x1.0p-53)
@inline function _bits2float(x::UInt64, ::Type{Float32})
    #=
    # this implementation uses more high bits, but is harder to vectorize
    x = x >>> 16  # discard low 16 bits
    u = Float32(x >>> 24) * Float32(0x1.0p-24)
    l = Float32(x & 0x00ffffff) * Float32(0x1.0p-24)
    =#
    ui = (x>>>32) % UInt32
    li = x % UInt32
    u = Float32(ui >>> 8) * Float32(0x1.0p-24)
    l = Float32(li >>> 8) * Float32(0x1.0p-24)
    (UInt64(reinterpret(UInt32, u)) << 32) | UInt64(reinterpret(UInt32, l))
end

# required operations. These could be written more concisely with `ntuple`, but the compiler
# sometimes refuses to properly vectorize.
for N in [4,8,16]
    let code, s, fshl = "llvm.fshl.v$(N)i64",
        VT = :(NTuple{$N, VecElement{UInt64}})

        s = ntuple(_->VecElement(UInt64(45)), N)
        @eval @inline _rotl45(x::$VT) = ccall($fshl, llvmcall, $VT, ($VT, $VT, $VT), x, x, $s)

        s = ntuple(_->VecElement(UInt64(23)), N)
        @eval @inline _rotl23(x::$VT) = ccall($fshl, llvmcall, $VT, ($VT, $VT, $VT), x, x, $s)

        code = """
        %lshiftOp = shufflevector <1 x i64> <i64 17>, <1 x i64> undef, <$N x i32> zeroinitializer
        %res = shl <$N x i64> %0, %lshiftOp
        ret <$N x i64> %res
        """
        @eval @inline _shl17(x::$VT) = llvmcall($code, $VT, Tuple{$VT}, x)

        code = """
        %res = add <$N x i64> %1, %0
        ret <$N x i64> %res
        """
        @eval @inline _plus(x::$VT, y::$VT) = llvmcall($code, $VT, Tuple{$VT, $VT}, x, y)

        code = """
        %res = xor <$N x i64> %1, %0
        ret <$N x i64> %res
        """
        @eval @inline _xor(x::$VT, y::$VT) = llvmcall($code, $VT, Tuple{$VT, $VT}, x, y)

        code = """
        %res = and <$N x i64> %1, %0
        ret <$N x i64> %res
        """
        @eval @inline _and(x::$VT, y::$VT) = llvmcall($code, $VT, Tuple{$VT, $VT}, x, y)

        code = """
        %res = or <$N x i64> %1, %0
        ret <$N x i64> %res
        """
        @eval @inline _or(x::$VT, y::$VT) = llvmcall($code, $VT, Tuple{$VT, $VT}, x, y)

        code = """
        %tmp = insertelement <1 x i64> undef, i64 %1, i32 0
        %shift = shufflevector <1 x i64> %tmp, <1 x i64> %tmp, <$N x i32> zeroinitializer
        %res = lshr <$N x i64> %0, %shift
        ret <$N x i64> %res
        """
        @eval @inline _lshr(x::$VT, y::Int64) = llvmcall($code, $VT, Tuple{$VT, Int64}, x, y)

        code = """
        %shiftamt = shufflevector <1 x i64> <i64 11>, <1 x i64> undef, <$N x i32> zeroinitializer
        %sh = lshr <$N x i64> %0, %shiftamt
        %f = uitofp <$N x i64> %sh to <$N x double>
        %scale = shufflevector <1 x double> <double 0x3ca0000000000000>, <1 x double> undef, <$N x i32> zeroinitializer
        %m = fmul <$N x double> %f, %scale
        %i = bitcast <$N x double> %m to <$N x i64>
        ret <$N x i64> %i
        """
        @eval @inline _bits2float(x::$VT, ::Type{Float64}) = llvmcall($code, $VT, Tuple{$VT}, x)

        code = """
        %as32 = bitcast <$N x i64> %0 to <$(2N) x i32>
        %shiftamt = shufflevector <1 x i32> <i32 8>, <1 x i32> undef, <$(2N) x i32> zeroinitializer
        %sh = lshr <$(2N) x i32> %as32, %shiftamt
        %f = uitofp <$(2N) x i32> %sh to <$(2N) x float>
        %scale = shufflevector <1 x float> <float 0x3e70000000000000>, <1 x float> undef, <$(2N) x i32> zeroinitializer
        %m = fmul <$(2N) x float> %f, %scale
        %i = bitcast <$(2N) x float> %m to <$N x i64>
        ret <$N x i64> %i
        """
        @eval @inline _bits2float(x::$VT, ::Type{Float32}) = llvmcall($code, $VT, Tuple{$VT}, x)
    end
end

function forkRand(rng::Union{TaskLocalRNG, Xoshiro}, ::Val{N}) where N
    # constants have nothing up their sleeve. For more discussion, cf rng_split in task.c
    # 0x02011ce34bce797f == hash(UInt(1))|0x01
    # 0x5a94851fb48a6e05 == hash(UInt(2))|0x01
    # 0x3688cf5d48899fa7 == hash(UInt(3))|0x01
    # 0x867b4bb4c42e5661 == hash(UInt(4))|0x01
    s0 = ntuple(i->VecElement(0x02011ce34bce797f * rand(rng, UInt64)), Val(N))
    s1 = ntuple(i->VecElement(0x5a94851fb48a6e05 * rand(rng, UInt64)), Val(N))
    s2 = ntuple(i->VecElement(0x3688cf5d48899fa7 * rand(rng, UInt64)), Val(N))
    s3 = ntuple(i->VecElement(0x867b4bb4c42e5661 * rand(rng, UInt64)), Val(N))
    (s0, s1, s2, s3)
end

_id(x, T) = x

@inline function xoshiro_bulk(rng::Union{TaskLocalRNG, Xoshiro}, dst::Ptr{UInt8}, len::Int, T::Union{Type{UInt8}, Type{Bool}, Type{Float32}, Type{Float64}}, ::Val{N}, f::F = _id) where {N, F}
    if len >= simdThreshold(T)
        written = xoshiro_bulk_simd(rng, dst, len, T, Val(N), f)
        len -= written
        dst += written
    end
    if len != 0
        xoshiro_bulk_nosimd(rng, dst, len, T, f)
    end
    nothing
end

@noinline function xoshiro_bulk_nosimd(rng::Union{TaskLocalRNG, Xoshiro}, dst::Ptr{UInt8}, len::Int, ::Type{T}, f::F) where {T, F}
    if rng isa TaskLocalRNG
        task = current_task()
        s0, s1, s2, s3 = task.rngState0, task.rngState1, task.rngState2, task.rngState3
    else
        (; s0, s1, s2, s3) = rng::Xoshiro
    end

    i = 0
    while i+8 <= len
        res = _plus(_rotl23(_plus(s0,s3)),s0)
        unsafe_store!(reinterpret(Ptr{UInt64}, dst + i), f(res, T))
        t = _shl17(s1)
        s2 = _xor(s2, s0)
        s3 = _xor(s3, s1)
        s1 = _xor(s1, s2)
        s0 = _xor(s0, s3)
        s2 = _xor(s2, t)
        s3 = _rotl45(s3)
        i += 8
    end
    if i < len
        res = _plus(_rotl23(_plus(s0,s3)),s0)
        t = _shl17(s1)
        s2 = _xor(s2, s0)
        s3 = _xor(s3, s1)
        s1 = _xor(s1, s2)
        s0 = _xor(s0, s3)
        s2 = _xor(s2, t)
        s3 = _rotl45(s3)
        ref = Ref(f(res, T))
        # TODO: This may make the random-stream dependent on system endianness
        ccall(:memcpy, Ptr{Cvoid}, (Ptr{UInt8}, Ptr{UInt64}, Csize_t), dst+i, ref, len-i)
    end
    if rng isa TaskLocalRNG
        task.rngState0, task.rngState1, task.rngState2, task.rngState3 = s0, s1, s2, s3
    else
       rng.s0, rng.s1, rng.s2, rng.s3 =  s0, s1, s2, s3
    end
    nothing
end

@noinline function xoshiro_bulk_nosimd(rng::Union{TaskLocalRNG, Xoshiro}, dst::Ptr{UInt8}, len::Int, ::Type{Bool}, f)
    if rng isa TaskLocalRNG
        task = current_task()
        s0, s1, s2, s3 = task.rngState0, task.rngState1, task.rngState2, task.rngState3
    else
        (; s0, s1, s2, s3) = rng::Xoshiro
    end

    i = 0
    while i+8 <= len
        res = _plus(_rotl23(_plus(s0,s3)),s0)
        shift = 0
        while i+8 <= len && shift < 8
            resLoc = _and(_lshr(res, shift), 0x0101010101010101)
            unsafe_store!(reinterpret(Ptr{UInt64}, dst + i), resLoc)
            i += 8
            shift += 1
        end

        t = _shl17(s1)
        s2 = _xor(s2, s0)
        s3 = _xor(s3, s1)
        s1 = _xor(s1, s2)
        s0 = _xor(s0, s3)
        s2 = _xor(s2, t)
        s3 = _rotl45(s3)
    end
    if i < len
        # we may overgenerate some bytes here, if len mod 64 <= 56 and len mod 8 != 0
        res = _plus(_rotl23(_plus(s0,s3)),s0)
        resLoc = _and(res, 0x0101010101010101)
        ref = Ref(resLoc)
        ccall(:memcpy, Ptr{Cvoid}, (Ptr{UInt8}, Ptr{UInt64}, Csize_t), dst+i, ref, len-i)
        t = _shl17(s1)
        s2 = _xor(s2, s0)
        s3 = _xor(s3, s1)
        s1 = _xor(s1, s2)
        s0 = _xor(s0, s3)
        s2 = _xor(s2, t)
        s3 = _rotl45(s3)
    end
    if rng isa TaskLocalRNG
        task.rngState0, task.rngState1, task.rngState2, task.rngState3 = s0, s1, s2, s3
    else
        rng.s0, rng.s1, rng.s2, rng.s3 = s0, s1, s2, s3
    end
    nothing
end

@noinline function xoshiro_bulk_simd(rng::Union{TaskLocalRNG, Xoshiro}, dst::Ptr{UInt8}, len::Int, ::Type{T}, ::Val{N}, f::F) where {T,N,F}
    s0, s1, s2, s3 = forkRand(rng, Val(N))

    i = 0
    while i + 8*N <= len
        res = _plus(_rotl23(_plus(s0,s3)),s0)
        t = _shl17(s1)
        s2 = _xor(s2, s0)
        s3 = _xor(s3, s1)
        s1 = _xor(s1, s2)
        s0 = _xor(s0, s3)
        s2 = _xor(s2, t)
        s3 = _rotl45(s3)
        unsafe_store!(reinterpret(Ptr{NTuple{N,VecElement{UInt64}}}, dst + i), f(res, T))
        i += 8*N
    end
    return i
end

@noinline function xoshiro_bulk_simd(rng::Union{TaskLocalRNG, Xoshiro}, dst::Ptr{UInt8}, len::Int, ::Type{Bool}, ::Val{N}, f) where {N}
    s0, s1, s2, s3 = forkRand(rng, Val(N))
    msk = ntuple(i->VecElement(0x0101010101010101), Val(N))
    i = 0
    while i + 64*N <= len
        res = _plus(_rotl23(_plus(s0,s3)),s0)
        t = _shl17(s1)
        s2 = _xor(s2, s0)
        s3 = _xor(s3, s1)
        s1 = _xor(s1, s2)
        s0 = _xor(s0, s3)
        s2 = _xor(s2, t)
        s3 = _rotl45(s3)
        for k=0:7
            tmp = _lshr(res, k)
            toWrite = _and(tmp, msk)
            unsafe_store!(reinterpret(Ptr{NTuple{N,VecElement{UInt64}}}, dst + i + k*N*8), toWrite)
        end
        i += 64*N
    end
    return i
end

function rand!(rng::Union{TaskLocalRNG, Xoshiro}, dst::Array{Float32}, ::SamplerTrivial{CloseOpen01{Float32}})
    GC.@preserve dst xoshiro_bulk(rng, convert(Ptr{UInt8}, pointer(dst)), length(dst)*4, Float32, xoshiroWidth(), _bits2float)
    dst
end

function rand!(rng::Union{TaskLocalRNG, Xoshiro}, dst::Array{Float64}, ::SamplerTrivial{CloseOpen01{Float64}})
    GC.@preserve dst xoshiro_bulk(rng, convert(Ptr{UInt8}, pointer(dst)), length(dst)*8, Float64, xoshiroWidth(), _bits2float)
    dst
end

for T in BitInteger_types
    @eval function rand!(rng::Union{TaskLocalRNG, Xoshiro}, dst::Union{Array{$T}, UnsafeView{$T}}, ::SamplerType{$T})
        GC.@preserve dst xoshiro_bulk(rng, convert(Ptr{UInt8}, pointer(dst)), length(dst)*sizeof($T), UInt8, xoshiroWidth())
        dst
    end
end

function rand!(rng::Union{TaskLocalRNG, Xoshiro}, dst::Array{Bool}, ::SamplerType{Bool})
    GC.@preserve dst xoshiro_bulk(rng, convert(Ptr{UInt8}, pointer(dst)), length(dst), Bool, xoshiroWidth())
    dst
end

end # module

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
# See https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method

function sqrt(x)
    t = 1
    for i in 1:10
        t = (t + x / t) / 2
    end
    t
end

expressions = [
    :(sqrt(2))
    :(sqrt(81))
    :(sqrt(62232491515607091882574410635924603070626544377175485625797))
]

for expr in expressions
    println("[] :julia-answer \"\"\"", escape_string(string(expr)), " = ", eval(expr), "\"\"\".")
end

max_stars_repo_name: pmoura/eye

max_issues_repo_name: pmoura/eye

max_forks_repo_name: pmoura/eye


text:
###########################################################################
#
#  Sequential Monte-Carlo
#
#  Takes MCMCTasks and performs Sequential Monte-Carlo sampling
#   ( a kind of population MC )
#
#  The last of the MCMCTasks should be the target distribution of interest,
#   others can either be the same target or tempered versions (see ref).
#
#  Particles are an arbitrary number n of initial values on which the tasks
#   will be run in succession, the number of tasks (target distributions) is not
#   related to the number of particles.
#
#  Ref : "On population-based simulation for static inference" A.Jasra -
#          D.Stephens - C.Holmes, p 271
#
###########################################################################

export SeqMC

println("Loading SeqMC(steps, burnin, trigger) runner")

immutable SeqMC <: MCMCRunner
  steps::Int
  burnin::Int
  trigger::Float64

  function SeqMC(steps::Int, burnin::Int, trigger::Float64)
	  @assert burnin >= 0 "Burnin rounds ($burnin) should be >= 0"
	  @assert steps > burnin "Steps ($steps) should be > to burnin ($burnin)"

    new(steps, burnin, trigger)
  end
end

SeqMC(; steps::Int=1, burnin::Int=0, trigger::Float64=1e-10) = SeqMC(steps, burnin, trigger)

function run_seqmc(targets::Array{MCMCTask}; particles::Vector{Vector{Float64}} = [[randn()] for i in 1:100])
	local ntargets = length(targets)
	local npart = length(particles)
	local tsize = targets[end].model.size
  local steps = targets[end].runner.steps
  local burnin = targets[end].runner.burnin
  local trigger = targets[end].runner.trigger

  @assert all(map(t->t.model.size, targets) .== tsize) "Models do not have the same parameter vector size"

	tic() # start timer

	# initialize all tasks
	map(t -> consume(t.task), targets)

	# initialize samples and weights arrays
	samples = fill(NaN, tsize, (steps-burnin)*npart)
	weights = fill(NaN, (steps-burnin)*npart)

	local pars = deepcopy(particles)
	local logW = zeros(npart)  # log of particle weights
	local logtarget = zeros(npart) # loglik of previous target distrib

	for i in 1:steps  # i = 1

		for t in targets  # t = targets[1]
			# mutate each particle with task t
			for n in 1:npart  #  n = 1
				MCMC.reset(t, pars[n])  # force pars of task #t to particle #n
				sample = consume(t.task)
				pars[n], plogtarget, ll0 = sample.ppars, sample.plogtarget, sample.logtarget
				logW[n] += ll0 - logtarget[n]
				logtarget[n] = plogtarget
			end

			# resample if likelihood variance of particles is too low
			#  TODO : improve, make user-settable, ..
			local W = exp(logW)
			if var(W) < trigger
				cp = cumsum(W) / sum(W)
				rs = fill(0, npart)
				for n in 1:npart  #  n = 1
					l = rand()
					rs[n] = findfirst(p-> (p>=l), cp)
				end
				pars = pars[rs]
				logW = zeros(npart)
				logtarget = logtarget[rs]
				# println("resampled !")
			end
		end

		println("iter $i, var $(var(exp(logW)))")
		logtarget = zeros(npart)

		if i > burnin # store betas of all particles in the result chain
			pos = (i-burnin-1) * npart
			for n in 1:npart  #  n = 1
				samples[:, pos+n] = pars[n]
				weights[pos+n] = exp(logW[n])
			end
		end
	end

	# generate column names for the samples DataFrame
	cn = ASCIIString[]
	for (k,v) in targets[end].model.pmap
	  if length(v.dims) == 0 # scalar
	    push!(cn, string(k))
	  elseif length(v.dims) == 1 # vector
	    cn = vcat(cn, ASCIIString[ "$k.$i" for i in 1:v.dims[1] ])
	  elseif length(v.dims) == 2 # matrix
	    cn = vcat(cn, ASCIIString[ "$k.$i.$j" for i in 1:v.dims[1], j in 1:v.dims[2] ])
	  end
	end

	# create Chain
	MCMCChain((burnin+1):1:((steps-burnin)*npart),
	  DataFrame(samples', cn),
		DataFrame(),  # TODO, store gradient here, needs to be passed by newprop
		{"weigths" => weights, "particle" => rep([1:npart],(steps-burnin))},
		targets,
		toq())
end

# TODO: check that all elements of array contain MCMCTasks of the same type
function resume_seqmc(targets::Array{MCMCTask}; steps::Int=100)
	run(MCMCTask[targets[i].model * targets[i].sampler * SeqMC(steps=steps, trigger=targets[i].runner.trigger)
		for i in 1:length(targets)])
end


text:
mutable struct WealthInSpace <: AbstractAgent
    id::Int
    pos::NTuple{2,Int}
    wealth::Int
end

"""
``` julia
wealth_distribution(;
    dims = (25, 25),
    wealth = 1,
    M = 1000
)
```
Same as in [Wealth distribution model](@ref).
"""
function wealth_distribution(; dims = (25, 25), wealth = 1, M = 1000)
    space = GridSpace(dims, periodic = true)
    model = ABM(WealthInSpace, space; scheduler = schedule_randomly)
    for i in 1:M # add agents in random positions
        add_agent!(model, wealth)
    end
    return model, wealth_distribution_agent_step!, dummystep
end

function wealth_distribution_agent_step!(agent, model)
    agent.wealth == 0 && return # do nothing
    neighboring_positions = collect(nearby_positions(agent.pos, model))
    push!(neighboring_positions, agent.pos) # also consider current position
    rpos = rand(model.rng, neighboring_positions) # the position that we will exchange with
    available_ids = ids_in_position(rpos, model)
    if length(available_ids) > 0
        random_neighbor_agent = model[rand(model.rng, available_ids)]
        agent.wealth -= 1
        random_neighbor_agent.wealth += 1
    end
end


text:
@test cleanzipcode("01234") == "01234"
@test cleanzipcode("1234") == "01234"
@test cleanzipcode(" 1234 ") == "01234"
@test cleanzipcode(" 1234-9999") == "01234"
@test isna(cleanzipcode("a"))
@test cleanzipcode("a", enforcestring = true) == "NA"

@test cleanzipcode(1234) == "01234"
@test cleanzipcode(0123456789, returnNA=false) == "123456789"

zipvec1 = ["1234", 1234, "1234-1234"]
@test cleanzipcode(zipvec1, enforcestring=true) == ["01234", "01234", "01234"]

cleaner1 = removewhitespace()
cleaner2 = removezipsuffix()
cleaner3 = padleftzeros()

@test ismatch(cleaner1.pattern, " 01234 ")
@test cleaner1.cleanfun(" 01234 ") == "01234"
@test ismatch(cleaner2.pattern, "01234-9999")
@test cleaner2.cleanfun("01234-9999") == "01234"
@test ismatch(cleaner3.pattern, "1234")
@test cleaner3.cleanfun("1234") == "01234"


text:
module TestFeatures
using SimString
using Test

@testset "Test feature extraction" begin
    char_ngram_res = SimString.extract_features(CharacterNGrams(3, " "), "prepress")
    @test char_ngram_res[5] == ("pre", 2)

    # Unicode text test
    @test SimString.extract_features(CharacterNGrams(2, " "), "∀∃😄🍕")[3] == ("😄🍕", 1)

    word_ngram_res = SimString.extract_features(WordNGrams(2, " ", " "), "You are a really really really cool dude 😄🍕")
    @test word_ngram_res[5] == (["really", "really"], 2)
    @test word_ngram_res[8] == (["dude", "😄🍕"], 1)
end

@testset "Test padding" begin
    @test SimString.pad_string(["one", "word"], " ") == [" ", "one", "word", " "]
    @test SimString.pad_string("one word", " ") == " one word "
end

@testset "Test cummulative_ngram_count" begin

end

end  # module


text:
module DoloTests

using Test
using Dolo, DataStructures

tests = length(ARGS) > 0 ? ARGS : [
                                   "model_types",
                                   "model_import",
                                   "test_calibration",
                                   "test_algos",
                                   "test_perfect_foresight",
                                   "test_features",
                                   "test_linter",
                                   "test_minilang",
                                   "test_decision_rules",
                                   "test_discretize_process"
                                   ]

for t in tests
    include("$(t).jl")
end

end  # module


text:
using ResumableFunctions
using Test

@resumable function test_for(a::Int=0; b::Int=a+1) :: Int
  for i in 1:10
    @yield a
    a, b = b, a+b
  end
end

@testset "test_for" begin
@test collect(test_for(4)) == [4, 5, 9, 14, 23, 37, 60, 97, 157, 254]
end

@resumable function test_try(io)
  try
    a = 1
    @yield a
    a = 2
    c = @yield a
    println(io,c)
  catch except
    println(io,except)
  finally
    println(io,"Always")
  end
end

struct SpecialException <: Exception end
@testset "test_try" begin
io = IOBuffer()
try_me = test_try(io)
try_me()
try_me(SpecialException())
@test String(take!(copy(io))) == "SpecialException()\nAlways\n"

io = IOBuffer()
try_me = test_try(io)
try_me()
try_me()
try_me("Hello")
@test String(take!(copy(io))) == "Hello\nAlways\n"

io = IOBuffer()
try_me = test_try(io)
try_me()
try_me()
try_me(SpecialException())
@test_throws ErrorException try_me()
@test String(take!(copy(io))) == "SpecialException()\nAlways\n"
end #test_try

@resumable function (test_where1(a::N) :: N) where {N<:Number}
  b = a + one(N)
  for i in 1:10
    @yield a
    a, b = b, a+b
  end
end

@resumable function (test_where2(a::N=4; b::N=a + one(N)) :: N) where N
  for i in 1:10
    @yield a
    a, b = b, a+b
  end
end

@testset "test_where" begin
@test collect(test_where1(4.0)) == [4.0, 5.0, 9.0, 14.0, 23.0, 37.0, 60.0, 97.0, 157.0, 254.0]
@test collect(test_where2(4)) == [4, 5, 9, 14, 23, 37, 60, 97, 157, 254]
end

@resumable function test_varargs(a...)
  for (i, e) in enumerate(a)
    @yield e
  end
end

@testset "test_varargs" begin
@test collect(test_varargs(1, 2, 3)) == [1, 2, 3]
end

@resumable function test_let()
  for u in [[(1,2),(3,4)], [(5,6),(7,8)]]
    for i in 1:2
      let i=i
        val = [a[i] for a in u]
      end
      @yield val
    end
  end
end

@testset "test_let" begin
@test collect(test_let()) == [[1,3],[2,4],[5,7],[6,8]]
end

@resumable function test_return_value()
  return 1
end

@testset "test_return_value" begin
  @test collect(test_return_value()) == []
end


text:
using LinearManifoldCluster
import Random
using Statistics
using MultivariateStats
using LinearAlgebra

import LaTeXStrings: latexstring
using Plots
upscale = 2
Plots.scalefontsizes(upscale)
default(size=(600*upscale,400*upscale))
pgfplots()
# gr()
# clibrary(:colorcet)
# theme(:wong2)

const SEED = 921369837
Random.seed!(SEED)

R(θ) = [cos(θ) -sin(θ); sin(θ) cos(θ)]

function elliptical_bound(S::AbstractMatrix{<:Real}, χ::Real)
    F = eigen(Symmetric(S))
    λ = F.values
    FI = sortperm(λ, rev=true)
    EV = F.vectors[:,FI[1]]
    ϕ = atan(EV[2], EV[1])
    if ϕ < 0
        ϕ += 2π
    end

    θ = range(0.0, 2π, step=0.1)
    a, b = χ*sqrt(λ[FI[1]]), χ*sqrt(λ[FI[2]])
    x, y  = a.*cos.(θ), b.*sin.(θ)
    ellipse = R(ϕ) * hcat(x, y)'
    return ellipse
end

# generate
n = 200          # number of points in generated manifold
N = 2            # space dimensionality
M = [1]         # number of manifolds
τ = .2 => .8   # translation vector bounds
E = [[.12]]      # bound on an extend from manifold for all manifolds)
Φ = Vector{Float64}[[.9]]   # bounds on manifold points (per manifold wrt dims)
X = generate(n,N,M,τ,Φ,E)[1][1]

# rotate
RX = R(π/5)*X
RX .+=  [0.5; -0.13]

# perform PCA
M = fit(PCA, RX)
P = projection(M)
O = mean(M)

e = 0.1
L1 = [O+P[:,1] O-P[:,1]]
L2 = [O+P[:,2]*e O-P[:,2]*e]

plot(RX[1,:], RX[2,:], seriestype=:scatter, aspect_ratio=:equal,
     xlims=(0.0,1.0), ylims=(0.2,0.85), markersize=2,
     legend=:bottomright, lab=latexstring("LM cluster, \$C_{\\Lambda,\\theta}\$"), palette=:viridis)
plot!(L1[1,:], L1[2,:], seriestype=:line, linewidth=2, lab=latexstring("linear manifold, \$\\Lambda\$"))
plot!(L2[1,:], L2[2,:], seriestype=:line, linewidth=2, linestyle=:dash, lab="orthoganal subspace")

# boundary
C = cov(RX.-O, dims=2)

for χ in [2.0, 2.5, 3.0]
    B = elliptical_bound(C, χ).+O
    plot!(B[1,:], B[2,:], seriestype=:path, linewidth=1,
          lab=latexstring("geom. LM cluster, \$C_{\\Lambda,\\chi=$χ}\$"))
end

Plots.pdf("linear-manifold-cluster-1d")


text:
"""
    basis(elem::Line,N,r)

Computes the generalized Vandermonde matrix V of degree N (along with the
derivative matrix Vr) at points r.
"""
function basis(elem::Line, N, r)
    V1D, Vr1D = ntuple(x->zeros(length(r), N+1),2)
    for j = 1:N+1
        V1D[:,j] .= jacobiP(r[:], 0, 0, j-1)
        Vr1D[:,j] .= grad_jacobiP(r[:], 0, 0, j-1)
    end
    return V1D, Vr1D
end

"""
    nodes(elem::Line,N)

Computes interpolation nodes of degree N.
"""
nodes(elem::Line, N) = first(gauss_lobatto_quad(0, 0, N))

"""
    equi_nodes(elem::Line,N)

Computes equally spaced nodes of degree N.
"""
equi_nodes(elem::Line, N) = collect(LinRange(-1, 1, N+1))

"""
    quad_nodes(elem::Line,N)

Computes (N+1)-point Gauss quadrature rule (exact for degree 2N+1 polynomials)

"""
quad_nodes(elem::Line, N) = gauss_quad(0, 0, N)


text:
<<<<<<< Updated upstream
function run_Hermite4(ps, params_nbody, filepath)
    @unpack Δt, t_end, ϵ, α, save_interval = params_nbody
=======
function run_Hermite4(ps, params_sim, filepath)
    @unpack Δt, t_end, ϵ, α, save_interval = params_sim
>>>>>>> Stashed changes

    times = (0:Δt:t_end)
    initialize!(ps)

    @show Δt
    @show t_end
    println(length(times))

    open(filepath, "w") do io

        E₀ = getTotalEnergy(ps)
        for (i, t) in enumerate(times)
            (i-1)%save_interval == 0 && save_snapshot_txt(io, t, ps)

            predict!(ps, Δt)
            evaluate!(ps, ϵ)
            collect!(ps, Δt)
            evaluate!(ps, ϵ)

            prepare!(ps, Δt)
        end
        E₁ = getTotalEnergy(ps)
        ΔE = (E₁ - E₀)/E₀
        @show E₀
        @show E₁
        @show ΔE
    end
end

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

<<<<<<< Updated upstream
function run_Hermite4_test(ps, params_nbody, filepath)
    open(filepath, "w") do io

        @unpack η, t_end, save_interval = params_nbody
=======
function run_Hermite4_test(ps, params_sim, filepath)
    open(filepath, "w") do io

        @unpack η, t_end, save_interval = params_sim
>>>>>>> Stashed changes

        initialize!(ps)
        Δt = η * minimum(norm(p.a)/norm(p.a¹) for p in ps)
        t = 0.
        i = 0

        E₀ = getTotalEnergy(ps)
        while t < t_end
            i%save_interval == 0 && save_snapshot_txt(io, t, ps)
            collide(ps) != (0, 0) && break

            Hermite4_shared!(ps, params_sim, Δt)
            Δt = get_Δt_Aarseth(ps, params_sim)
            t += Δt
            i += 1
        end
        E₁ = getTotalEnergy(ps)
        ΔE = (E₁ - E₀)/E₀
        @show E₀
        @show E₁
        @show ΔE
        println("Collision: ", collide(ps))
    end
end

<<<<<<< Updated upstream
function Hermite4_shared!(ps, params_nbody, Δt)
    # @unpack ϵ, α = params_nbody
=======
function Hermite4_shared!(ps, params_sim, Δt)
    # @unpack ϵ, α = params_sim
>>>>>>> Stashed changes

    predict!(ps, Δt)
    evaluate!(ps)
    collect!(ps, Δt)
    evaluate!(ps)

    prepare!(ps, Δt)
end

#############################################################
#                   Predict & collect
#############################################################

function predict!(ps, Δt)
    Δt² = Δt * Δt
    Δt³ = Δt * Δt²

    for p in ps
        r, v, a, a¹ = SVector{3}(p.r), SVector{3}(p.v), SVector{3}(p.a), SVector{3}(p.a¹)

        p.⁺r .= r + Δt*v + Δt²/2*a + Δt³/6*a¹
        p.⁺v .= v + Δt*a + Δt²/2*a¹
    end
end

function evaluate!(ps, ϵ)
    for p in ps
        p.⁺a  .= 0
        p.⁺a¹ .= 0
    end

    for i in eachindex(ps)
        for j in eachindex(ps)
            i == j && continue
            ps[j].m == 0 && continue

            r = SVector{3}(ps[j].⁺r) - SVector{3}(ps[i].⁺r)
            v = SVector{3}(ps[j].⁺v) - SVector{3}(ps[i].⁺v)

            r² = norm(r)^2 + ϵ^2
            r⁻³ = r²^(-3/2)
            r⁻⁵ = r²^(-5/2)

            ps[i].⁺a  .+= G * ps[j].m * r * r⁻³
            ps[i].⁺a¹ .+= G * ps[j].m * (v*r⁻³ - 3(v ⋅ r)*r*r⁻⁵)
        end
    end
end

evaluate!(ps) = evaluate!(ps, 0)

function collect!(ps, Δt, α)
    Δt² = Δt * Δt
    Δt³ = Δt * Δt²
    Δt⁴ = Δt * Δt³
    Δt⁵ = Δt * Δt⁴

    for p in ps
         a,  a¹ = SVector{3}(p.a),  SVector{3}(p.a¹)
        ⁺a, ⁺a¹ = SVector{3}(p.⁺a), SVector{3}(p.⁺a¹)

        p.a² .= ( -6(a - ⁺a) - Δt*(4a¹ + 2⁺a¹) ) / Δt²
        p.a³ .= ( 12(a - ⁺a) + 6Δt*(a¹ +  ⁺a¹) ) / Δt³

        a², a³ = SVector{3}(p.a²), SVector{3}(p.a³)

        p.⁺r .+= Δt⁴/24*a² + α*Δt⁵/120*a³
        p.⁺v .+= Δt³/ 6*a² +   Δt⁴/ 24*a³
    end
end

collect!(ps, Δt) = collect!(ps, Δt, 7/6)

function prepare!(ps, Δt)
    for p in ps
        @. p.r  = p.⁺r
        @. p.v  = p.⁺v
        @. p.a  = p.⁺a
        @. p.a¹ = p.⁺a¹
<<<<<<< Updated upstream

        # @. p.⁺a³ = p.a³
        @. p.a² += Δt*p.a³  # for calculating time step
=======

        # @. p.⁺a³ = p.a³
        @. p.a² += Δt*p.a³
>>>>>>> Stashed changes
    end
end

#############################################################
#                         Time step
#############################################################

<<<<<<< Updated upstream
get_Δt_Aarseth(ps, params_nbody) = minimum(get_Δt_Aarseth(p, params_nbody) for p in ps)

function get_Δt_Aarseth(p::AbstractParticle, params_nbody)
    @unpack η = params_nbody
=======
get_Δt_Aarseth(ps, params_sim) = minimum(get_Δt_Aarseth(p, params_sim) for p in ps)

function get_Δt_Aarseth(p::AbstractParticle, params_sim)
    @unpack η = params_sim
>>>>>>> Stashed changes

    a  = norm(SVector{3}(p.a ))
    a¹ = norm(SVector{3}(p.a¹))
    a² = norm(SVector{3}(p.a²))
    a³ = norm(SVector{3}(p.a³))

    Δt = η * sqrt((a*a² + a¹*a¹) / (a¹*a³ + a²*a²))
end

max_stars_repo_name: MasanoriKanamaru/Astroshaper

max_issues_repo_name: MasanoriKanamaru/Astroshaper

max_forks_repo_name: MasanoriKanamaru/Astroshaper


text:
#!/usr/bin/env julia

# LD_LIBRARY_PATH=/usr/lib/llvm-3.7/lib

using Clang.wrap_c

ctx = wrap_c.init(
        headers=["/usr/include/expat.h"],
        common_file="common.jl",
        output_file="output.jl",
        header_library="libexpat",
        header_wrapped=(header, cursorname) -> startswith(cursorname, "/usr/include/expat"))
wrap_c.run(ctx)

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
using Gencode, Preprocessing

function gmshcall(pde, filename, nd, elemtype)

opts = "-format msh3";

# find gmsh executable
gmsh = Preprocessing.findexec(pde.gmsh, pde.version);

print("Gmsh mesh generator...\n");
mystr = gmsh * " " * filename * ".geo -" * string(nd) * " " * opts;
run(Gencode.string2cmd(mystr));

tm = readlines(filename * ".msh");

nlines = length(tm);
j = 0;
for i = 1:nlines
    if tm[i] == "\$Nodes"
        j = i;
        break;
    end
end
np = parse(Int64,tm[j+1]);
p = zeros(nd,np);
i = j;
for ii = 1:np
    j = i+1+ii;
    pii = parse.(Float64,split(tm[j]));
    p[:,ii] = pii[2:(1+nd)];
end
i = j+1;
for ii = i:nlines
    if tm[ii] == "\$Elements"
        j = ii;
        break;
    end
end
ne = parse(Int64,tm[j+1]);
if nd==1 # line
    nve = 2;
    wcase = 1;
elseif nd==2 && elemtype==0 # tri
    nve = 3;
    wcase = 2;
elseif nd==2 && elemtype==1 # quad
    nve = 4;
    wcase = 3;
elseif nd==3 && elemtype==0 # tet
    nve = 4;
    wcase = 4;
elseif nd==3 && elemtype==1 # hex
    nve = 8;
    wcase = 5;
end
t = zeros(Int,nve,ne);
m = 0;
i = j;
for ii = 1:ne
    j = i+1+ii;
    tii = parse.(Int,split(tm[j]));
    if (tii[2] == wcase)
        m = m + 1;
        t[:,m] = tii[5:end];
    end
end
t = t[:,1:m];

return p, t

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    _reverse_mode(d::NLPEvaluator, x)

Run reverse-mode automatic differentiation on `d` given the primal solution `x`.

This function updates many of the data-structures inside `d` in-place.

At a high level, reverse-mode AD has two phases:

In Phase I, we evalute the problem in `d` at the primal solution `x`, and
stores the primal solution of each expression in the tree and the first-order
partial derivative information for each node with respect to its arguments.

Because the nodes in our data structure are topologically sorted, we can make a
single pass through the tree by iterating backwards through the vector of stored
nodes.

In Phase II, we propagate the partial derivative information back down the tree
to find the derivative of each function with respect to the input.

Because the nodes in our data structure are topologically sorted, we can make a
single pass through the tree by iterating forwards through the vector of stored
nodes.
"""
function _reverse_mode(d::NLPEvaluator, x)
    if d.last_x == x
        # Fail fast if the primal solution has not changed since last call.
        return
    end
    # Phase I
    for k in d.subexpression_order
        d.subexpression_forward_values[k] =
            _forward_eval(d.subexpressions[k], d, x)
    end
    if d.objective !== nothing
        _forward_eval(d.objective::_FunctionStorage, d, x)
    end
    for con in d.constraints
        _forward_eval(con, d, x)
    end
    # Phase II
    for k in d.subexpression_order
        _reverse_eval(d.subexpressions[k])
    end
    if d.objective !== nothing
        _reverse_eval(d.objective::_FunctionStorage)
    end
    for con in d.constraints
        _reverse_eval(con)
    end
    copyto!(d.last_x, x)
    return
end

"""
    _forward_eval(
        f::Union{_FunctionStorage,_SubexpressionStorage},
        d::NLPEvaluator,
        x::AbstractVector{T},
    ) where {T}

Forward-mode evaluation of an expression tree given in `f`.

 * This function assumes that the values of all sub-expressions have already
   been computed and are stored in `d.subexpression_forward_values`.
 * `f.partials_storage[k]` is the partial derivative of `nodes[k].parent` with
   respect to the value of node `k`. It's efficient to compute this at the same
   time as the value of the parent because we use it in reverse mode and in dual
   forward mode. Note that `partials_storage`` makes a subtle assumption that we
   have a tree instead of a general DAG. If we have a DAG, then need to
   associate storage with each edge of the DAG.
"""
function _forward_eval(
    # !!! warning
    #     This Union depends upon _FunctionStorage and _SubexpressionStorage
    #     having similarly named fields.
    f::Union{_FunctionStorage,_SubexpressionStorage},
    d::NLPEvaluator,
    x::AbstractVector{T},
)::T where {T}
    @assert length(f.forward_storage) >= length(f.nodes)
    @assert length(f.partials_storage) >= length(f.nodes)
    operators = d.data.operators
    # f.nodes is already in order such that parents always appear before
    # children, so a backwards pass through f.nodes is a forward pass through
    # the tree.
    children_arr = SparseArrays.rowvals(f.adj)
    for k in length(f.nodes):-1:1
        node = f.nodes[k]
        f.partials_storage[k] = zero(T)
        if node.type == Nonlinear.NODE_VARIABLE
            f.forward_storage[k] = x[node.index]
            # This should never happen, because we will have replaced these by now.
            # elseif node.type == Nonlinear.NODE_MOI_VARIABLE
            #     f.forward_storage[k] = x[node.index]
        elseif node.type == Nonlinear.NODE_VALUE
            f.forward_storage[k] = f.const_values[node.index]
        elseif node.type == Nonlinear.NODE_SUBEXPRESSION
            f.forward_storage[k] = d.subexpression_forward_values[node.index]
        elseif node.type == Nonlinear.NODE_PARAMETER
            f.forward_storage[k] = d.data.parameters[node.index]
        elseif node.type == Nonlinear.NODE_CALL_MULTIVARIATE
            children_indices = SparseArrays.nzrange(f.adj, k)
            N = length(children_indices)
            # TODO(odow);
            # With appropriate benchmarking, the special-cased if-statements can
            # be removed in favor of the generic user-defined function case.
            if node.index == 1 # :+
                tmp_sum = zero(T)
                for c_idx in children_indices
                    @inbounds ix = children_arr[c_idx]
                    @inbounds f.partials_storage[ix] = one(T)
                    @inbounds tmp_sum += f.forward_storage[ix]
                end
                f.forward_storage[k] = tmp_sum
            elseif node.index == 2 # :-
                @assert N == 2
                child1 = first(children_indices)
                @inbounds ix1 = children_arr[child1]
                @inbounds ix2 = children_arr[child1+1]
                @inbounds tmp_sub = f.forward_storage[ix1]
                @inbounds tmp_sub -= f.forward_storage[ix2]
                @inbounds f.partials_storage[ix1] = one(T)
                @inbounds f.partials_storage[ix2] = -one(T)
                f.forward_storage[k] = tmp_sub
            elseif node.index == 3 # :*
                tmp_prod = one(T)
                for c_idx in children_indices
                    @inbounds tmp_prod *= f.forward_storage[children_arr[c_idx]]
                end
                if tmp_prod == zero(T) || N <= 2
                    # This is inefficient if there are a lot of children.
                    # 2 is chosen as a limit because (x*y)/y does not always
                    # equal x for floating-point numbers. This can produce
                    # unexpected error in partials. There's still an error when
                    # multiplying three or more terms, but users are less likely
                    # to complain about it.
                    for c_idx in children_indices
                        prod_others = one(T)
                        for c_idx2 in children_indices
                            (c_idx == c_idx2) && continue
                            ix = children_arr[c_idx2]
                            prod_others *= f.forward_storage[ix]
                        end
                        f.partials_storage[children_arr[c_idx]] = prod_others
                    end
                else
                    # Compute all-minus-one partial derivatives by dividing from
                    # the total product.
                    for c_idx in children_indices
                        ix = children_arr[c_idx]
                        f.partials_storage[ix] =
                            tmp_prod / f.forward_storage[ix]
                    end
                end
                @inbounds f.forward_storage[k] = tmp_prod
            elseif node.index == 4 # :^
                @assert N == 2
                idx1 = first(children_indices)
                idx2 = last(children_indices)
                @inbounds ix1 = children_arr[idx1]
                @inbounds ix2 = children_arr[idx2]
                @inbounds base = f.forward_storage[ix1]
                @inbounds exponent = f.forward_storage[ix2]
                if exponent == 2
                    @inbounds f.forward_storage[k] = base * base
                    @inbounds f.partials_storage[ix1] = 2 * base
                elseif exponent == 1
                    @inbounds f.forward_storage[k] = base
                    @inbounds f.partials_storage[ix1] = 1.0
                else
                    f.forward_storage[k] = pow(base, exponent)
                    f.partials_storage[ix1] = exponent * pow(base, exponent - 1)
                end
                f.partials_storage[ix2] = f.forward_storage[k] * log(base)
            elseif node.index == 5 # :/
                @assert N == 2
                idx1 = first(children_indices)
                idx2 = last(children_indices)
                @inbounds ix1 = children_arr[idx1]
                @inbounds ix2 = children_arr[idx2]
                @inbounds numerator = f.forward_storage[ix1]
                @inbounds denominator = f.forward_storage[ix2]
                recip_denominator = 1 / denominator
                @inbounds f.partials_storage[ix1] = recip_denominator
                f.partials_storage[ix2] =
                    -numerator * recip_denominator * recip_denominator
                f.forward_storage[k] = numerator * recip_denominator
            elseif node.index == 6 # ifelse
                @assert N == 3
                idx1 = first(children_indices)
                @inbounds condition = f.forward_storage[children_arr[idx1]]
                @inbounds lhs = f.forward_storage[children_arr[idx1+1]]
                @inbounds rhs = f.forward_storage[children_arr[idx1+2]]
                @inbounds f.partials_storage[children_arr[idx1+1]] =
                    condition == 1
                @inbounds f.partials_storage[children_arr[idx1+2]] =
                    !(condition == 1)
                f.forward_storage[k] = ifelse(condition == 1, lhs, rhs)
            else
                f_input = _UnsafeVectorView(d.jac_storage, N)
                ∇f = _UnsafeVectorView(d.user_output_buffer, N)
                for (r, i) in enumerate(children_indices)
                    f_input[r] = f.forward_storage[children_arr[i]]
                    ∇f[r] = 0.0
                end
                f.forward_storage[k] = Nonlinear.eval_multivariate_function(
                    operators,
                    operators.multivariate_operators[node.index],
                    f_input,
                )
                Nonlinear.eval_multivariate_gradient(
                    operators,
                    operators.multivariate_operators[node.index],
                    ∇f,
                    f_input,
                )
                for (r, i) in enumerate(children_indices)
                    f.partials_storage[children_arr[i]] = ∇f[r]
                end
            end
        elseif node.type == Nonlinear.NODE_CALL_UNIVARIATE
            child_idx = children_arr[f.adj.colptr[k]]
            f.forward_storage[k] = Nonlinear.eval_univariate_function(
                operators,
                operators.univariate_operators[node.index],
                f.forward_storage[child_idx],
            )
            f.partials_storage[child_idx] = Nonlinear.eval_univariate_gradient(
                operators,
                operators.univariate_operators[node.index],
                f.forward_storage[child_idx],
            )
        elseif node.type == Nonlinear.NODE_COMPARISON
            children_idx = SparseArrays.nzrange(f.adj, k)
            result = true
            f.partials_storage[children_arr[children_idx[1]]] = zero(T)
            for r in 2:length(children_idx)
                lhs = children_arr[children_idx[r-1]]
                rhs = children_arr[children_idx[r]]
                result &= Nonlinear.eval_comparison_function(
                    operators,
                    operators.comparison_operators[node.index],
                    f.forward_storage[lhs],
                    f.forward_storage[rhs],
                )
                f.partials_storage[rhs] = zero(T)
            end
            f.forward_storage[k] = result
        else
            @assert node.type == Nonlinear.NODE_LOGIC
            children_idx = SparseArrays.nzrange(f.adj, k)
            lhs = children_arr[children_idx[1]]
            rhs = children_arr[children_idx[2]]
            f.forward_storage[k] = Nonlinear.eval_logic_function(
                operators,
                operators.logic_operators[node.index],
                f.forward_storage[lhs] == 1,
                f.forward_storage[rhs] == 1,
            )
            f.partials_storage[lhs] = zero(T)
            f.partials_storage[rhs] = zero(T)
        end
    end
    return f.forward_storage[1]
end

"""
    _reverse_eval(f::Union{_FunctionStorage,_SubexpressionStorage})

Reverse-mode evaluation of an expression tree given in `f`.

 * This function assumes `f.partials_storage` is already updated.
 * This function assumes that `f.reverse_storage` has been initialized with 0.0.
"""
function _reverse_eval(
    # !!! warning
    #     This Union depends upon _FunctionStorage and _SubexpressionStorage
    #     having similarly named fields.
    f::Union{_FunctionStorage,_SubexpressionStorage},
)
    @assert length(f.reverse_storage) >= length(f.nodes)
    @assert length(f.partials_storage) >= length(f.nodes)
    # f.nodes is already in order such that parents always appear before
    # children so a forward pass through nodes is a backwards pass through the
    # tree.
    f.reverse_storage[1] = one(Float64)
    for k in 2:length(f.nodes)
        node = f.nodes[k]
        if node.type == Nonlinear.NODE_VALUE ||
           node.type == Nonlinear.NODE_LOGIC ||
           node.type == Nonlinear.NODE_COMPARISON ||
           node.type == Nonlinear.NODE_PARAMETER
            continue
        end
        rev_parent = f.reverse_storage[node.parent]
        partial = f.partials_storage[k]
        f.reverse_storage[k] = ifelse(
            rev_parent == 0.0 && !isfinite(partial),
            rev_parent,
            rev_parent * partial,
        )
    end
    return
end

"""
    _extract_reverse_pass(
        g::Vector{T},
        d::NLPEvaluator,
        f::Union{_FunctionStorage,_SubexpressionStorage},
    ) where {T}

Fill the gradient vector `g` with the values from the reverse pass. Assumes you
have already called `_reverse_eval_all(d, x)`.
"""
function _extract_reverse_pass(
    g::Vector{T},
    d::NLPEvaluator,
    f::Union{_FunctionStorage,_SubexpressionStorage},
) where {T}
    for i in f.dependent_subexpressions
        d.subexpression_reverse_values[i] = 0.0
    end
    _extract_reverse_pass_inner(g, f, d.subexpression_reverse_values, 1.0)
    for i in length(f.dependent_subexpressions):-1:1
        k = f.dependent_subexpressions[i]
        _extract_reverse_pass_inner(
            g,
            d.subexpressions[k],
            d.subexpression_reverse_values,
            d.subexpression_reverse_values[k],
        )
    end
    return
end

function _extract_reverse_pass_inner(
    output::AbstractVector{T},
    # !!! warning
    #     This Union depends upon _FunctionStorage and _SubexpressionStorage
    #     having similarly named fields.
    f::Union{_FunctionStorage,_SubexpressionStorage},
    subexpressions::AbstractVector{T},
    scale::T,
) where {T}
    @assert length(f.reverse_storage) >= length(f.nodes)
    for (k, node) in enumerate(f.nodes)
        if node.type == Nonlinear.NODE_VARIABLE
            output[node.index] += scale * f.reverse_storage[k]
        elseif node.type == Nonlinear.NODE_SUBEXPRESSION
            subexpressions[node.index] += scale * f.reverse_storage[k]
        end
    end
    return
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
module SorbetHelperDashFallingBlock

using ..Ahorn, Maple

@mapdef Entity "SorbetHelper/DashFallingBlock" DashFallingBlock(
    x::Integer,
    y::Integer,
    width::Integer=Maple.defaultBlockWidth,
    height::Integer=Maple.defaultBlockHeight,
    tiletype::String="3",
    shakeSfx::String="event:/game/general/fallblock_shake",
    impactSfx::String="event:/game/general/fallblock_impact",
    fallOnTouch::Bool=false,
    climbFall::Bool=true,
    fallOnStaticMover::Bool=false,
    depth::Integer=-9000,
    allowWavedash::Bool=false,
    dashCornerCorrection::Bool=false,
)

const placements = Ahorn.PlacementDict(
    "Dash Falling Block (Sorbet Helper)" => Ahorn.EntityPlacement(
        DashFallingBlock,
        "rectangle",
        Dict{String, Any}(),
        Ahorn.tileEntityFinalizer
    ),
)

Ahorn.editingOptions(entity::DashFallingBlock) = Dict{String, Any}(
    "tiletype" => Ahorn.tiletypeEditingOptions(),
    "shakeSfx" => Dict{String, String}(
        "Default" => "event:/game/general/fallblock_shake",
        "Snow" => "event:/game/01_forsaken_city/fallblock_ice_shake",
        "Wood" => "event:/game/03_resort/fallblock_wood_shake",
        "Reflection" => "event:/game/06_reflection/fallblock_boss_shake"
    ),
    "impactSfx" => Dict{String, String}(
        "Default" => "event:/game/general/fallblock_impact",
        "Snow" => "event:/game/01_forsaken_city/fallblock_ice_impact",
        "Wood" => "event:/game/03_resort/fallblock_wood_impact",
        "Reflection" => "event:/game/06_reflection/fallblock_boss_impact"
    )
)

Ahorn.minimumSize(entity::DashFallingBlock) = 8, 8
Ahorn.resizable(entity::DashFallingBlock) = true, true

Ahorn.selection(entity::DashFallingBlock) = Ahorn.getEntityRectangle(entity)

Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::DashFallingBlock, room::Maple.Room) = Ahorn.drawTileEntity(ctx, room, entity)

Ahorn.editingOrder(entity::DashFallingBlock) = String["x", "y", "width", "height", "shakeSfx", "impactSfx", "tiletype", "depth", "fallOnTouch", "climbFall", "fallOnStaticMover", "allowWavedash", "dashCornerCorrection"]

end

max_stars_repo_name: Kitty-Cats/SorbetHelper

max_issues_repo_name: Kitty-Cats/SorbetHelper

max_forks_repo_name: Kitty-Cats/SorbetHelper


text:
struct GeneDataView{G <: AbstractGene} <: AbstractArray{G, 1}
    parent::Vector{Record{G}}
    indices::Vector{UInt}
    property::Symbol
end

function Base.getproperty(gene::G, name::Symbol) where {G <: AbstractGene}
    if name in propertynames(gene)
        return parent(gene).genedata[index(gene), name]
    else
        return missing
    end
end

function Base.getproperty(gv::GeneDataView{G}, name::Symbol) where {G <: AbstractGene}
    return getfield(gv, name)
end

function Base.getproperty(genes::AbstractArray{G, 1}, name::Symbol) where {G <: AbstractGene}
    GeneDataView(parent.(genes), getfield.(genes, Ref(:index)), name)
end

function Base.setproperty!(gene::G, name::Symbol, x::T) where {G <: AbstractGene, T}
    if hasproperty(parent(gene).genedata, name)
        parent(gene).genedata[index(gene), name] = x
    else
        s = size(parent(gene).genedata, 1)
        parent(gene).genedata[!, name] = Vector{Union{Missing, T}}(missing, s)
        parent(gene).genedata[index(gene), name] = x
    end
    return x
end

Base.size(gv::GeneDataView) = size(gv.indices)
Base.getindex(gv::GeneDataView, i::Int) = getproperty(gv.parent[i].genes[gv.indices[i]], gv.property)
Base.getindex(gv::GeneDataView, I::AbstractArray) = GeneDataView(gv.parent[I], gv.indices[I], gv.property)
Base.setindex!(gv::GeneDataView, v, i::Int) = (Base.setproperty!(gv.parent[i].genes[gv.indices[i]], gv.property, v))
Base.similar(gv::GeneDataView{G}) where {G <: AbstractGene} = GeneDataView(gv.parent, gv.indices, gv.property)
Base.copy(gv::GeneDataView{G}) where {G <: AbstractGene} = GeneDataView(gv.parent, gv.indices, gv.property)
Base.view(gv::GeneDataView{G}, I) where {G <: AbstractGene} = GeneDataView(gv.parent[I], gv.indices[I], gv.property)

function Base.fill!(gv::GeneDataView{G}, x) where {G <: AbstractGene}
    chrs = unique(gv.parent)
    for chr in chrs
        if hasproperty(chr.genedata, gv.property)
            xT = convert(eltype(chr.genedata[!, gv.property]), x)
        else
            xT = x
        end
    end
    for i in eachindex(gv.indices)
        @inbounds gv[i] = x
    end
    gv
end


text:
# Draw some text to a PNG image.

using Compose

lines = "hello and goodbye\nFooSubBarSup\nA Third Line\nFooSubBarSupPooh\nA Fifth Line\nA Sixth Line"
c = compose(context(),
            (context(), text(0px, 250px, lines), fontsize(8pt), fill("tomato")),
            (context(), text(100px, 250px, lines), fontsize(24pt), fill("bisque")))
draw(SVG("text-fontfallback.svg", 400px, 400px), c)

import Cairo, Fontconfig

draw(SVG("text-pango.svg", 400px, 400px), c)
draw(PNG("text.png", 400px, 400px, dpi=192), c)
draw(PDF("text.pdf", 400px, 400px), c)


text:
function nutsepsilon(x::GeneralNode, logfgrad::Function, delta::Float64, target::Float64)
    n = size(x)[1] - 1
    log_target = log(target)
    # molifier is necessary!
    # blv = get_branchlength_vector(x)
    # set_branchlength_vector!(x, molifier.(blv, delta))

    logf0, gr = logfgrad(x)

    r0 = randn(n)
    x0 = Tree_HMC_State(deepcopy(x), r0, gr, logf0)
    x1 = transfer(x0)
    H0 = hamiltonian(x0)
    epsilon = 1.0
    _ = refraction!(x0, epsilon, logfgrad, delta)
    Hp = hamiltonian(x0)

    prob = Hp - H0
    direction = prob > log_target ? 1 : -1

    while direction == 1 ? prob > log_target : prob < log_target
        epsilon = direction == 1 ? 2 * epsilon : 0.5 * epsilon

        x2 = transfer(x1)
        _ = refraction!(x2, epsilon, logfgrad, delta)

        Hp = hamiltonian(x2)

        prob = Hp - H0
    end

    epsilon
end

function nutsepsilon(x::Vector{<:Real}, logfgrad::Function, target::Float64)

    logf0, gr = logfgrad(x)
    target = log(target)

    r0 = randn(size(x))
    x0 = Array_HMC_State(x, r0, gr, logf0)
    x1 = transfer(x0)
    H0 = hamiltonian(x0)
    epsilon = 1.0
    leapfrog!(x0, epsilon, logfgrad)
    Hp = hamiltonian(x0)

    prob = Hp - H0

    direction = prob > target ? 1 : -1
    restart = false
    while direction == 1 ? prob > target : prob < target
        epsilon = direction == 1 ? 2 * epsilon : 0.5 * epsilon
        if isapprox(epsilon, 0)
            throw("The estimated stepsize is too small. ϵ: $epsilon")
        elseif epsilon > 1e7
            throw("The estimated stepsize is too large. ϵ: $epsilon")
        end
        x2 = transfer(x1)
        leapfrog!(x2, epsilon, logfgrad)

        Hp = hamiltonian(x2)

        prob = Hp - H0

    end

    if restart
        @warn "restarting initialization of ϵ"
        epsilon = nutsepsilon(x, logfgrad, exp(target))
    end
    epsilon
end

@inline function scale_fac(x::T, delta::T)::T where {T<:Real}
    x < delta ? x/delta : 1.0
end

@inline function molifier(x::T, delta::Real)::T where {T<:Real}
    x >= delta ? x : (x^2 + delta^2) / (2.0 * delta)
end # function

@inline function abs_adapter(x::R)::Float64 where R <:Real
    x / (1+abs(x))
end

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
#
# Code adapted from: [Nanosoldier.jl](https://github.com/JuliaCI/Nanosoldier.jl)
# See: src/jobs/BenchmarkJob.jl
#
# Reporting should probably be moved into BenchmarkTooks.jl or into a package separate from
# Nanosoldier.jl
#
# Note: Only minimal changes have been done to the reporting code.

const REGRESS_MARK = ":x:"
const IMPROVE_MARK = ":white_check_mark:"

function printreport(io::IO, judged)
    # print result table #
    #--------------------#

    iscomparisonjob = true
    tablegroup = judged

    println(io, """
                ## Results
                *Note: If Chrome is your browser, I strongly recommend installing the [Wide GitHub](https://chrome.google.com/webstore/detail/wide-github/kaalofacklcidaampbokdplbklpeldpj?hl=en)
                extension, which makes the result table easier to read.*
                Below is a table of this job's results, obtained by running the benchmarks found in
                [JuliaCI/BaseBenchmarks.jl](https://github.com/JuliaCI/BaseBenchmarks.jl). The values
                listed in the `ID` column have the structure `[parent_group, child_group, ..., key]`,
                and can be used to index into the BaseBenchmarks suite to retrieve the corresponding
                benchmarks.
                The percentages accompanying time and memory values in the below table are noise tolerances. The "true"
                time/memory value for a given benchmark is expected to fall within this percentage of the reported value.
                """)

    if iscomparisonjob
        print(io, """
                  A ratio greater than `1.0` denotes a possible regression (marked with $(REGRESS_MARK)), while a ratio less
                  than `1.0` denotes a possible improvement (marked with $(IMPROVE_MARK)). Only significant results - results
                  that indicate possible regressions or improvements - are shown below (thus, an empty table means that all
                  benchmark results remained invariant between builds).
                  | ID | time ratio | memory ratio |
                  |----|------------|--------------|
                  """)
    else
        print(io, """
                  | ID | time | GC time | memory | allocations |
                  |----|------|---------|--------|-------------|
                  """)
    end

    entries = BenchmarkTools.leaves(tablegroup)

    try
        entries = entries[sortperm(map(x -> string(first(x)), entries))]
    end

    for (ids, t) in entries
        if !(iscomparisonjob) || BenchmarkTools.isregression(t) || BenchmarkTools.isimprovement(t)
            println(io, resultrow(ids, t))
        end
    end

    println(io)

    # print list of executed benchmarks #
    #-----------------------------------#
    println(io, """
                ## Benchmark Group List
                Here's a list of all the benchmark groups executed by this job:
                """)

    for id in unique(map(pair -> pair[1][1:end-1], entries))
        println(io, "- `", idrepr(id), "`")
    end

    println(io)

    # # print build version info #
    # #--------------------------#
    # print(io, """
    #           ## Version Info
    #           #### Primary Build
    #           ```
    #           $(build.vinfo)
    #           ```
    #           """)

    # if hasagainstbuild
    #     println(io)
    #     print(io, """
    #               #### Comparison Build
    #               ```
    #               $(get(job.against).vinfo)
    #               ```
    #               """)
    # end
    return nothing
end

idrepr(id) = (str = repr(id); str[searchindex(str, '['):end])

intpercent(p) = string(ceil(Int, p * 100), "%")

resultrow(ids, t::BenchmarkTools.Trial) = resultrow(ids, minimum(t))

function resultrow(ids, t::BenchmarkTools.TrialEstimate)
    t_tol = intpercent(BenchmarkTools.params(t).time_tolerance)
    m_tol = intpercent(BenchmarkTools.params(t).memory_tolerance)
    timestr = string(BenchmarkTools.prettytime(BenchmarkTools.time(t)), " (", t_tol, ")")
    memstr = string(BenchmarkTools.prettymemory(BenchmarkTools.memory(t)), " (", m_tol, ")")
    gcstr = BenchmarkTools.prettytime(BenchmarkTools.gctime(t))
    allocstr = string(BenchmarkTools.allocs(t))
    return "| `$(idrepr(ids))` | $(timestr) | $(gcstr) | $(memstr) | $(allocstr) |"
end

function resultrow(ids, t::BenchmarkTools.TrialJudgement)
    t_tol = intpercent(BenchmarkTools.params(t).time_tolerance)
    m_tol = intpercent(BenchmarkTools.params(t).memory_tolerance)
    t_ratio = @sprintf("%.2f", BenchmarkTools.time(BenchmarkTools.ratio(t)))
    m_ratio =  @sprintf("%.2f", BenchmarkTools.memory(BenchmarkTools.ratio(t)))
    t_mark = resultmark(BenchmarkTools.time(t))
    m_mark = resultmark(BenchmarkTools.memory(t))
    timestr = "$(t_ratio) ($(t_tol)) $(t_mark)"
    memstr = "$(m_ratio) ($(m_tol)) $(m_mark)"
    return "| `$(idrepr(ids))` | $(timestr) | $(memstr) |"
end

resultmark(sym::Symbol) = sym == :regression ? REGRESS_MARK : (sym == :improvement ? IMPROVE_MARK : "")


text:
using HomotopyContinuation

@testset "whitney umbrella" begin
    @var x y z
    V = HomotopyOpt.ConstraintVariety([x,y,z], [x^2-y^2*z], 3, 2, 100);
    p0 = V.samples[1]
    objective(x) = sin(x[1])+cos(x[2])+cos(sin(x[3]))

    println("gaussnewtonstep")
    @time resultminimum = HomotopyOpt.findminima(p0, 1e-4, V, objective; maxseconds=120, whichstep="gaussnewtonstep", initialstepsize=0.5);
    @test(resultminimum.converged==true)
    @test(resultminimum.lastpointisminimum==true)

    println("EDStep")
    @time resultminimum = HomotopyOpt.findminima(p0, 1e-4, V, objective; maxseconds=120, whichstep="EDStep", initialstepsize=0.5);
    @test(resultminimum.converged==true)
    @test(resultminimum.lastpointisminimum==true)

    println("twostep")
    @time resultminimum = HomotopyOpt.findminima(p0, 1e-4, V, objective; maxseconds=120, whichstep="twostep", initialstepsize=0.5);
    @test(resultminimum.converged==true)
    @test(resultminimum.lastpointisminimum==true)
end

@testset "space circle Test" begin
    @var x y z
    V = HomotopyOpt.ConstraintVariety([x,y,z], [x^2+y^2+z^2-1, z], 3, 1, 100);
    p0 = V.samples[1]
    objective(x) = exp(x[1]+x[2]+x[3])

    println("gaussnewtonstep")
    @time resultminimum = HomotopyOpt.findminima(p0, 1e-4, V, objective; maxseconds=120, whichstep="gaussnewtonstep", initialstepsize=0.5);
    @test(resultminimum.converged==true)
    @test(resultminimum.lastpointisminimum==true)

    println("EDStep")
    @time resultminimum = HomotopyOpt.findminima(p0, 1e-4, V, objective; maxseconds=120, whichstep="EDStep", initialstepsize=0.5);
    @test(resultminimum.converged==true)
    @test(resultminimum.lastpointisminimum==true)

    println("twostep")
    @time resultminimum = HomotopyOpt.findminima(p0, 1e-4, V, objective; maxseconds=120, whichstep="twostep", initialstepsize=0.5);
    @test(resultminimum.converged==true)
    @test(resultminimum.lastpointisminimum==true)
end


text:
import Base: copy, copyto!, ==
import Random: rand, seed!
import RandomNumbers: gen_seed, seed_type

# Generic functions

@inline seed!(r::AbstractPCG{StateType}, seed::Integer=gen_seed(StateType)) where
    StateType <: PCGUInt = pcg_seed!(r, seed % StateType)
@inline seed!(r::PCGStateSetseq{StateType}, seed::NTuple{2, Integer}=gen_seed(StateType, 2)) where
    StateType <: PCGUInt = pcg_seed!(r, seed[1] % StateType, seed[2] % StateType)

@inline seed_type(::Type{PCGStateOneseq{T, T1, T2}}) where {T, T1, T2} = T
@inline seed_type(::Type{PCGStateMCG{   T, T1, T2}}) where {T, T1, T2} = T
@inline seed_type(::Type{PCGStateUnique{T, T1, T2}}) where {T, T1, T2} = T
@inline seed_type(::Type{PCGStateSetseq{T, T1, T2}}) where {T, T1, T2} = NTuple{2, T}

function copyto!(dest::T, src::T) where T <: AbstractPCG
    dest.state = src.state
    dest
end
function copyto!(dest::T, src::T) where T <: PCGStateSetseq
    dest.state = src.state
    dest.inc = src.inc
    dest
end

copy(src::T) where T <: AbstractPCG = copyto!(T(), src)

==(r1::T, r2::T) where T <: AbstractPCG = r1.state == r2.state
==(r1::T, r2::T) where T <: PCGStateSetseq = r1.state == r2.state && r1.inc == r2.inc

@inline function rand(r::AbstractPCG{StateType, MethodType, OutputType}, ::Type{OutputType}) where
        {StateType <: Union{pcg_uints[1:end-1]...}, MethodType <: PCGMethod, OutputType <: PCGUInt}
    old_state = r.state
    pcg_step!(r)
    pcg_output(old_state, MethodType)
end

@inline function rand(r::AbstractPCG{UInt128, MethodType, OutputType}, ::Type{OutputType}) where
        {MethodType <: PCGMethod, OutputType <: PCGUInt}
    pcg_step!(r)
    pcg_output(r.state, MethodType)
end

"""
```julia
bounded_rand(r, bound)
```

Producing a random number less than a given `bound` in the output type.
"""
@inline function bounded_rand(s::AbstractPCG{StateType, MethodType, OutputType}, bound::OutputType) where
        {StateType <: PCGUInt, MethodType <: PCGMethod, OutputType <: PCGUInt}
    threshold = (-bound) % bound
    r = rand(s, OutputType)
    while r < threshold
        r = rand(s, OutputType)
    end
    r % bound
end

"""
```julia
advance!(r, Δ)
```

Advance a PCG object `r` for `Δ` steps.

# Examples
```jldoctest
julia> r = PCGStateSetseq(UInt64, PCG_RXS_M_XS, (123, 321))
PCGStateSetseq{UInt64,Val{:RXS_M_XS},UInt64}(0x45389f8b27528b29, 0x0000000000000283)

julia> A = rand(r, UInt64, 2);

julia> p = rand(r);

julia> r
PCGStateSetseq{UInt64,Val{:RXS_M_XS},UInt64}(0x9b1fc763ae0ad702, 0x0000000000000283)

julia> advance!(r, -3)
PCGStateSetseq{UInt64,Val{:RXS_M_XS},UInt64}(0x45389f8b27528b29, 0x0000000000000283)

julia> @Test.test A == rand(r, UInt64, 2)
Test Passed

julia> @Test.test p == rand(r)
Test Passed
```
"""
@inline function advance!(r::AbstractPCG{StateType}, Δ::Integer) where StateType <: PCGUInt
    pcg_advance!(r, Δ % StateType)
    r
end

# Constructors.
for (pcg_type_t, uint_type, method_symbol, output_type) in PCG_LIST
    pcg_type = Symbol("PCGState$pcg_type_t")
    method = Val{method_symbol}

    if pcg_type_t != :Setseq
        @eval function $pcg_type(output_type::Type{$output_type}, method::Type{$method},
                seed::Integer=gen_seed($uint_type))
            s = $pcg_type{$uint_type, method, output_type}()
            seed!(s, seed)
            s
        end
    else
        @eval function $pcg_type(output_type::Type{$output_type}, method::Type{$method},
                seed::NTuple{2, Integer}=gen_seed($uint_type, 2))
            s = $pcg_type{$uint_type, method, output_type}()
            seed!(s, seed)
            s
        end
    end
end

"""
```julia
PCGStateOneseq{StateType<:PCGUInt, MethodType<:PCGMethod, OutputType<:PCGUInt} <:
    AbstractPCG{StateType, MethodType, OutputType}
PCGStateOneseq([seed])
PCGStateOneseq(output_type[, seed])
PCGStateOneseq(method[, seed])
PCGStateOneseq(output_type[, method, seed])
```

PCG generator with *single streams*, where all instances use the same fixed constant, thus the RNG always
somewhere in same sequence.

`seed` is an `Integer` which will be automatically converted to the state type.

`output_type` is the type of the PCG's output. If missing it is set to `UInt64`.

`method` is one of the [`PCGMethod`](@ref). If missing it is set to `PCG_XSH_RS`.

See [`PCG_LIST`](@ref) for the available parameter combinations.
"""
PCGStateOneseq(seed::Integer=gen_seed(UInt128)) = PCGStateOneseq(UInt64, PCG_XSH_RS, seed)
PCGStateOneseq(::Type{T}, seed::Integer=gen_seed(UInt128)) where
    T <: PCGUInt = PCGStateOneseq(T, PCG_XSH_RS, seed)
PCGStateOneseq(::Type{T}, seed::Integer=gen_seed(UInt128)) where
    T <: PCGMethod = PCGStateOneseq(UInt64, T, seed)

"""
```julia
PCGStateMCG{StateType<:PCGUInt, MethodType<:PCGMethod, OutputType<:PCGUInt} <:
    AbstractPCG{StateType, MethodType, OutputType}
PCGStateMCG([seed])
PCGStateMCG(output_type[, seed])
PCGStateMCG(method[, seed])
PCGStateMCG(output_type[, method, seed])
```

PCG generator with *MCG*, where the increment is zero, resulting in a single stream and reduced period.

`seed` is an `Integer` which will be automatically converted to the state type.

`output_type` is the type of the PCG's output. If missing it is set to `UInt64`.

`method` is one of the [`PCGMethod`](@ref). If missing it is set to `PCG_XSH_RS`.

See [`PCG_LIST`](@ref) for the available parameter combinations.
"""
PCGStateMCG(seed::Integer=gen_seed(UInt128)) = PCGStateMCG(UInt64, PCG_XSH_RS, seed)
PCGStateMCG(::Type{T}, seed::Integer=gen_seed(UInt128)) where
    T <: PCGUInt = PCGStateMCG(T, PCG_XSH_RS, seed)
PCGStateMCG(::Type{T}, seed::Integer=gen_seed(UInt128)) where
    T <: PCGMethod = PCGStateMCG(UInt64, T, seed)

"""
```julia
PCGStateSetseq{StateType<:PCGUInt, MethodType<:PCGMethod, OutputType<:PCGUInt} <:
    AbstractPCG{StateType, MethodType, OutputType}
PCGStateSetseq([seed])
PCGStateSetseq(output_type[, seed])
PCGStateSetseq(method[, seed])
PCGStateSetseq(output_type[, method, seed])
```

PCG generator with *specific streams*, where the constant can be changed at any time, selecting a different
random sequence.

`seed` is a `Tuple` of two `Integer`s which will both be automatically converted to the state type.

`output_type` is the type of the PCG's output. If missing it is set to `UInt64`.

`method` is one of the [`PCGMethod`](@ref). If missing it is set to `PCG_XSH_RS`.

See [`PCG_LIST`](@ref) for the available parameter combinations.
"""
PCGStateSetseq(seed::NTuple{2, Integer}=gen_seed(UInt128, 2)) = PCGStateSetseq(UInt64, PCG_XSH_RS, seed)
PCGStateSetseq(::Type{T}, seed::NTuple{2, Integer}=gen_seed(UInt128, 2)) where
    T <: PCGUInt = PCGStateSetseq(T, PCG_XSH_RS, seed)
PCGStateSetseq(::Type{T}, seed::NTuple{2, Integer}=gen_seed(UInt128, 2)) where
    T <: PCGMethod = PCGStateSetseq(UInt64, T, seed)

"""
```julia
PCGStateUnique{StateType<:PCGUInt, MethodType<:PCGMethod, OutputType<:PCGUInt} <:
    AbstractPCG{StateType, MethodType, OutputType}
PCGStateUnique([seed])
PCGStateUnique(output_type[, seed])
PCGStateUnique(method[, seed])
PCGStateUnique(output_type[, method, seed])
```

PCG generator with *unique streams*, where the constant is based on the memory address of the object, thus
every RNG has its own unique sequence.

`seed` is an `Integer` which will be automatically converted to the state type.

`output_type` is the type of the PCG's output. If missing it is set to `UInt64`.

`method` is one of the [`PCGMethod`](@ref). If missing it is set to `PCG_XSH_RS`.

See [`PCG_LIST`](@ref) for the available parameter combinations.
"""
PCGStateUnique(seed::Integer=gen_seed(UInt128)) = PCGStateUnique(UInt64, PCG_XSH_RS, seed)
PCGStateUnique(::Type{T}, seed::Integer=gen_seed(UInt128)) where
    T <: PCGUInt = PCGStateUnique(T, PCG_XSH_RS, seed)
PCGStateUnique(::Type{T}, seed::Integer=gen_seed(UInt128)) where
    T <: PCGMethod = PCGStateUnique(UInt64, T, seed)


text:
using Microbiome
using Test
using Microbiome.SparseArrays
using Microbiome.Tables
using Microbiome.Dictionaries
import Microbiome.MultivariateStats: MDS

@testset "Samples and Features" begin
    @testset "MicriobiomeSamples and metadata" begin
        ms = MicrobiomeSample("sample")
        @test name(ms) == "sample"
        @test isempty(metadata(ms))
        @test metadata(ms) isa Dictionary
        @test_throws Dictionaries.IndexError ms.thing = "metadata"
        @test_throws Dictionaries.IndexError ms[:thing] = "metadata"

        @test insert!(ms, :thing, "metadata") isa MicrobiomeSample
        @test ms.thing == ms[:thing] == "metadata"
        @test let
            ms.thing = "other metadata"
            @test ms.thing == ms[:thing] == "other metadata"
            ms[:thing] = "metadata"
            ms.thing == ms[:thing] == "metadata"
        end
        @test_throws MethodError ms["thing"] = "metadata"
        @test_throws Dictionaries.IndexError insert!(ms, :thing, "still other metadata")
        @test_throws Dictionaries.IndexError delete!(ms, :thing2)
        @test let
            delete!(ms, :thing)
            unset!(ms, :thing2)
            true
        end
        @test let
            set!(ms, :thing2, "metadata2")
            ms.thing2 == "metadata2"
        end
    end

    @testset "Taxa" begin
        _clades = (:domain, :kingdom, :phylum, :class, :order, :family, :genus, :species, :subspecies, :strain)
        txm = Taxon("taxon", missing)
        @test txm === Taxon("taxon")
        @test ismissing(clade(txm))
        @test !hasclade(txm)

        for (i, c) in enumerate(_clades)
            tx = Taxon("taxon", c)
            @test clade(tx) == c
            @test tx === Taxon("taxon", i-1)
            @test tx !== txm
            @test tx == txm
        end

        @test_throws ErrorException Taxon("taxon", :invalid)
        @test_throws ErrorException Taxon("taxon", 10)
        @test let tx = Taxon("taxon", :kingdom)
            hasclade(tx)
        end
    end

    @testset "Gene Functions" begin
        gfm = GeneFunction("gene", missing)
        @test name(gfm) == "gene"
        @test gfm === GeneFunction("gene")
        @test ismissing(taxon(gfm))
        @test !hastaxon(gfm)

        gf1 = GeneFunction("gene", Taxon("sp1", :species))
        gf2 = GeneFunction("gene", Taxon("sp1"))
        @test name(gf1) == "gene"

        @test gf1 == gf2 == gfm
        @test gf1 !== gf2
        @test gf1 !== gfm
        @test hastaxon(gf1)
        @test !ismissing(taxon(gf1))
        @test taxon(gf1) == taxon(gf2)
    end
end

@testset "Profiles" begin
    _clades = (:domain, :kingdom, :phylum, :class, :order, :family, :genus, :species, :subspecies)
    mss = [MicrobiomeSample("sample$i") for i in 1:5]
    txs = [Taxon("taxon$i", _clades[i]) for i in 1:9]
    push!(txs, Taxon("taxon10", missing))

    mat = spzeros(10,5)
    for i in 1:5; mat[i,i] = 1.; end
    for i in 1:5; mat[i+5,i] = 0.6; end

    comm = CommunityProfile(mat, txs, mss)

    @testset "Profile operations" begin
        @test CommunityProfile{Float64, Taxon, MicrobiomeSample}(mat, txs, mss) isa CommunityProfile
        @test nsamples(comm) == 5
        @test nfeatures(comm) == 10
        @test size(comm) == (10, 5)
        @test profiletype(comm) == Taxon
        @test clades(comm)[1:9] == [_clades...]
        @test sampletotals(comm) == [1.6 1.6 1.6 1.6 1.6]
        @test featuretotals(comm) == reshape([1.0, 1.0, 1.0, 1.0, 1.0, 0.6, 0.6, 0.6, 0.6, 0.6], 10, 1)
        @test features(comm) == txs
        @test samples(comm) == mss

        @test present(0.1)
        @test !present(0.1, 0.2)
        @test_throws DomainError present(-0.1)
        @test_throws DomainError present(0., -0.1)

        @test prevalence([0.0, 0.1, 0.2, 0.3]) ≈ 0.75
        @test prevalence([0.0, 0.1, 0.2, 0.3], 0.15) ≈ 0.5

        @test all(≈(0.2), prevalence(comm))
        @test all(prevalence(comm, 0.7) .≈ [0.2, 0.2, 0.2, 0.2, 0.2, 0, 0, 0, 0, 0])

        @test let c2 = deepcopy(comm)
            relativeabundance!(c2)
            @test all(≈(0.625), featuretotals(c2)[1:5])
            all(≈(0.375), featuretotals(c2)[6:end])
        end
        @test let c2 = deepcopy(comm)
            relativeabundance!(c2, kind=:percent)
            @test all(≈(62.5), featuretotals(c2)[1:5])
            all(≈(37.5), featuretotals(c2)[6:end])
        end

        @test_throws ArgumentError relativeabundance!(comm, kind=:invalid)
    end

    @testset "Indexing and Tables integration" begin
        for i in 1:5
            @test abundances(comm[:, "sample$i"]) == mat[:, [i]]
            @test abundances(comm["taxon$i", :]) == mat[[i], :]
        end

        for (i, col) in enumerate(Tables.columns(comm))
            if i == 1
                @test col == name.(txs)
            else
                @test col ==  mat[:, [i-1]]
            end
        end

        for (i, row) in enumerate(Tables.rows(comm))
            @test row == (; :features => name(txs[i]), (Symbol("sample$(j)") => mat[i, j] for j in 1:5)...)
        end

        tbl = Tables.columntable(comm)
        @test tbl.features == featurenames(comm)
        @test keys(tbl) == (:features, Symbol.(samplenames(comm))...)
        @test tbl.sample1 == abundances(comm[:, "sample1"])
    end

    @testset "Diversity" begin
        R = 10
        s1 = collect(0:R-1) # high diversity
        s2 = [i % 2 == 0 ? s1[i] : 0 for i in eachindex(s1)] # low diversity
        s3 = ones(R) # uniform
        s4 = [10, zeros(R-1)...] # no diversity

        @test shannon(s1) > shannon(s2)
        @test shannon(s3) ≈ log(R)
        @test shannon(s4) ≈ 0.0
        for s in (s1, s2, s3, s4)
            shannon(s ./ sum(s)) ≈ shannon(s)
        end

        @test ginisimpson(s1) > ginisimpson(s2)
        @test ginisimpson(s3) ≈ 1.0 - 1/R
        @test ginisimpson(s4) ≈ 0.0
        for s in (s1, s2, s3, s4)
            ginisimpson(s ./ sum(s)) ≈ ginisimpson(s)
        end

        @test size(shannon(comm)) == (1, 5)
        @test size(ginisimpson(comm)) == (1, 5)

        @test let c2 = deepcopy(comm)
            shannon!(c2)
            @test all(s-> haskey(s, :shannon), samples(c2))
            @test_throws Dictionaries.IndexError shannon!(c2)
            shannon!(c2, overwrite=true)
            true
        end
        @test let c2 = deepcopy(comm)
            ginisimpson!(c2)
            @test all(s-> haskey(s, :ginisimpson), samples(c2))
            @test_throws Dictionaries.IndexError ginisimpson!(c2)
            ginisimpson!(c2, overwrite=true)
            true
        end

        @test all(braycurtis(comm) .== [0 1 1 1 1; 1 0 1 1 1; 1 1 0 1 1; 1 1 1 0 1; 1 1 1 1 0])
        @test all(jaccard(comm) .== [0 1 1 1 1; 1 0 1 1 1; 1 1 0 1 1; 1 1 1 0 1; 1 1 1 1 0])
        @test all(hellinger(comm) .== [0 1 1 1 1; 1 0 1 1 1; 1 1 0 1 1; 1 1 1 0 1; 1 1 1 1 0])
        @test pcoa(comm) isa MDS
    end

end


text:
# ---
# title: 1298. Maximum Candies You Can Get from Boxes
# id: problem1298
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Breadth-first Search
# link: <https://leetcode.com/problems/maximum-candies-you-can-get-from-boxes/description/>
# hidden: true
# ---
#
# Given `n` boxes, each box is given in the format `[status, candies, keys,
# containedBoxes]` where:
#
#   * `status[i]`: an integer which is **1** if `box[i]` is open and **0** if `box[i]` is closed.
#   * `candies[i]`: an integer representing the number of candies in `box[i]`.
#   * `keys[i]`: an array contains the indices of the boxes you can open with the key in `box[i]`.
#   * `containedBoxes[i]`: an array contains the indices of the boxes found in `box[i]`.
#
# You will start with some boxes given in `initialBoxes` array. You can take all
# the candies in any open box and you can use the keys in it to open new boxes
# and you also can use the boxes you find in it.
#
# Return _the maximum number of candies_ you can get following the rules above.
#
#
#
# **Example 1:**
#
#
#
#     Input: status = [1,0,1,0], candies = [7,5,4,100], keys = [[],[],[1],[]], containedBoxes = [[1,2],[3],[],[]], initialBoxes = [0]
#     Output: 16
#     Explanation: You will be initially given box 0. You will find 7 candies in it and boxes 1 and 2. Box 1 is closed and you don't have a key for it so you will open box 2. You will find 4 candies and a key to box 1 in box 2.
#     In box 1, you will find 5 candies and box 3 but you will not find a key to box 3 so box 3 will remain closed.
#     Total number of candies collected = 7 + 4 + 5 = 16 candy.
#
#
# **Example 2:**
#
#
#
#     Input: status = [1,0,0,0,0,0], candies = [1,1,1,1,1,1], keys = [[1,2,3,4,5],[],[],[],[],[]], containedBoxes = [[1,2,3,4,5],[],[],[],[],[]], initialBoxes = [0]
#     Output: 6
#     Explanation: You have initially box 0. Opening it you can find boxes 1,2,3,4 and 5 and their keys. The total number of candies will be 6.
#
#
# **Example 3:**
#
#
#
#     Input: status = [1,1,1], candies = [100,1,100], keys = [[],[0,2],[]], containedBoxes = [[],[],[]], initialBoxes = [1]
#     Output: 1
#
#
# **Example 4:**
#
#
#
#     Input: status = [1], candies = [100], keys = [[]], containedBoxes = [[]], initialBoxes = []
#     Output: 0
#
#
# **Example 5:**
#
#
#
#     Input: status = [1,1,1], candies = [2,3,2], keys = [[],[],[]], containedBoxes = [[],[],[]], initialBoxes = [2,1,0]
#     Output: 7
#
#
#
#
# **Constraints:**
#
#   * `1 <= status.length <= 1000`
#   * `status.length == candies.length == keys.length == containedBoxes.length == n`
#   * `status[i]` is `0` or `1`.
#   * `1 <= candies[i] <= 1000`
#   * `0 <= keys[i].length <= status.length`
#   * `0 <= keys[i][j] < status.length`
#   * All values in `keys[i]` are unique.
#   * `0 <= containedBoxes[i].length <= status.length`
#   * `0 <= containedBoxes[i][j] < status.length`
#   * All values in `containedBoxes[i]` are unique.
#   * Each box is contained in one box at most.
#   * `0 <= initialBoxes.length <= status.length`
#   * `0 <= initialBoxes[i] < status.length`
#
#
## @lc code=start
using LeetCode

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


text:
mutable struct RbfCovarianceCache
    last_theta::AbstractArray{Float64, 1}
    D2::AbstractArray{Float64, 2}
    K::AbstractArray{Float64, 2}
end

"""
    SquaredExponentialIsoKernel <: AbstractGPKernel

Squared exponential kernel with uniform length scale across all dimensions, ``l``.
```math
\\begin{aligned}
K(r) & = \\sigma_f^2 e^{-r/2} \\\\
r_{ij} & = \\sum_{k=1}^d\\frac{(x_{ik}-z_{jk})^2}{l^2}
\\end{aligned}
```

``r_{ij}`` are computed by [`scaled_squared_distance`](@ref)

# Hyperparameters
Hyperparameters vector for this kernel must contain two elements, in the following order:
1. ``\\sigma_f``: the signal standard deviation
2. ``l``: the length scale
"""
struct SquaredExponentialIsoKernel <: AbstractGPKernel
    cache::RbfCovarianceCache
end

function SquaredExponentialIsoKernel()
    SquaredExponentialIsoKernel(RbfCovarianceCache(Array{Float64}(undef, 0), Array{Float64}(undef, 0, 0), Array{Float64}(undef, 0, 0)))
end

"""
    SquaredExponentialArdKernel <: AbstractGPKernel

Squared exponential kernel with distinct length scale for each dimention, ``l_k``.
```math
\\begin{aligned}
K(r) & = \\sigma_f^2 e^{-r/2} \\\\
r_{ij} & = \\sum_{k=1}^d\\frac{(x_{ik}-z_{jk})^2}{l_k^2}
\\end{aligned}
```

``r_{ij}`` are computed by [`scaled_squared_distance`](@ref)

# Hyperparameters
The length of hyperparameters array for this kernel depends on the dimensionality
of the data. Assuming each data point is a vector in a ``d``-dimensional space,
this kernel needs ``d+1`` hyperparameters, in the following order:
1. ``\\sigma_f``: the signal standard deviation
2. ``l_1, \\ldots, l_d``: the length scales for each dimension
"""
struct SquaredExponentialArdKernel <: AbstractGPKernel
    cache::RbfCovarianceCache
end

function SquaredExponentialArdKernel()
    SquaredExponentialArdKernel(RbfCovarianceCache(Array{Float64}(undef, 0), Array{Float64}(undef, 0, 0), Array{Float64}(undef, 0, 0)))
end

function get_hyperparameters_size(kernel::SquaredExponentialIsoKernel, training_data::AbstractArray{Float64, 2})
    return 2
end

function get_hyperparameters_size(kernel::SquaredExponentialArdKernel, training_data::AbstractArray{Float64, 2})
    return size(training_data, 2) + 1
end

function rbf_covariance_grad_common(cache::RbfCovarianceCache, log_theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2}, R::AbstractArray{Float64, 2})
    if log_theta == cache.last_theta && size(cache.K, 1) == size(x, 1) && size(cache.D2, 1) == size(x, 1)
        K = cache.K
        D2 = cache.D2
    else
        sigma_f = exp(log_theta[1] * 2)
        D2 = scaled_squared_distance(log_theta[2:end], x, x)
        K = sigma_f .* exp.(-D2 ./ 2)
        cache.last_theta = copy(log_theta)
        cache.D2 = D2
        cache.K = K
    end
    KR = K .* R
    d_sigma_f = 2 * sum(KR)
    if length(log_theta) == 2
        d_ell = KR[:]' * D2[:]
    else
        d_ell = -scaled_squared_distance_grad(log_theta[2:end], x, x, KR)/2
    end
    return [d_sigma_f; d_ell]
end

function rbf_covariance_common(log_theta::AbstractArray{Float64, 1},
        x1::AbstractArray{Float64, 2}, x2::AbstractArray{Float64, 2})
    sigma_f = exp(log_theta[1] * 2)
    D2 = scaled_squared_distance(log_theta[2:end], x1, x2)
    K = sigma_f .* exp.(-D2 ./ 2)
end

function covariance(ker::SquaredExponentialIsoKernel, log_theta::AbstractArray{Float64, 1},
        x1::AbstractArray{Float64, 2}, x2::AbstractArray{Float64, 2})
    rbf_covariance_common(log_theta, x1, x2)
end

function covariance(ker::SquaredExponentialArdKernel, log_theta::AbstractArray{Float64, 1},
        x1::AbstractArray{Float64, 2}, x2::AbstractArray{Float64, 2})
    rbf_covariance_common(log_theta, x1, x2)
end

function rbf_covariance_common_training(cache::RbfCovarianceCache,
        log_theta::AbstractArray{Float64, 1}, x::AbstractArray{Float64, 2})
    if log_theta != cache.last_theta || size(cache.D2, 1) != size(x, 1) || size(cache.K, 1) != size(x, 1)
        cache.last_theta = copy(log_theta)
        sigma_f = exp(log_theta[1] * 2)
        D2 = scaled_squared_distance(log_theta[2:end], x, x)
        cache.D2 = D2
        cache.K = sigma_f .* exp.(-D2 ./ 2)
    end
    cache.K
end

function covariance_training(ker::SquaredExponentialIsoKernel, log_theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2})
    rbf_covariance_common_training(ker.cache, log_theta, x)
end

function covariance_training(ker::SquaredExponentialArdKernel, log_theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2})
    rbf_covariance_common_training(ker.cache, log_theta, x)
end

function covariance_diagonal(ker::SquaredExponentialIsoKernel, log_theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2})
    fill(exp(log_theta[1] * 2), (size(x, 1), 1))
end

function covariance_diagonal(ker::SquaredExponentialArdKernel, log_theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2})
    fill(exp(log_theta[1] * 2), (size(x, 1), 1))
end

function covariance_grad(ker::SquaredExponentialIsoKernel, log_theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2}, R::AbstractArray{Float64, 2})
    rbf_covariance_grad_common(ker.cache, log_theta, x, R)
end

function covariance_grad(ker::SquaredExponentialArdKernel, log_theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2}, R::AbstractArray{Float64, 2})
    rbf_covariance_grad_common(ker.cache, log_theta, x, R)
end