text:
export JacobiWeight

"""
    JacobiWeight(β,α,s::Space)

weights a space `s` by a Jacobi weight, which on `-1..1`
is `(1+x)^β*(1-x)^α`.
For other domains, the weight is inferred by mapping to `-1..1`.
"""
struct JacobiWeight{S,DD,RR,T<:Real} <: WeightSpace{S,DD,RR}
    β::T
    α::T
    space::S
    function JacobiWeight{S,DD,RR,T}(β::T,α::T,space::S) where {S<:Space,DD,RR,T}
        if isa(space,JacobiWeight)
            new(β+space.β,α+space.α,space.space) else
            new(β,α,space)
        end
    end
end

JacobiWeight{S,DD,RR,T}(β,α,space::Space) where {S,DD,RR,T} =
    JacobiWeight{S,DD,RR,T}(convert(T,β)::T, convert(T,α)::T, convert(S,space)::S)

JacobiWeight(a::Number, b::Number, d::Space) =
    JacobiWeight{typeof(d),domaintype(d),rangetype(d),promote_type(eltype(a),eltype(b))}(a,b,d)
JacobiWeight(β::Number, α::Number, d::JacobiWeight) =  JacobiWeight(β+d.β,α+d.α,d.space)
JacobiWeight(a::Number, b::Number, d::IntervalOrSegment) = JacobiWeight(a,b,Space(d))
JacobiWeight(a::Number, b::Number, d) = JacobiWeight(a,b,Space(d))
JacobiWeight(a::Number, b::Number) = JacobiWeight(a,b,Chebyshev())

JacobiWeight(a::Number, b::Number,s::PiecewiseSpace) = PiecewiseSpace(JacobiWeight(a,b,components(s)))

Fun(::typeof(identity), S::JacobiWeight) =
    isapproxinteger(S.β) && isapproxinteger(S.α) ? Fun(x->x,S) : Fun(identity,domain(S))

order(S::JacobiWeight{Ultraspherical{Int,D,R},D,R}) where {D,R} = order(S.space)

spacescompatible(A::JacobiWeight,B::JacobiWeight) =
    A.β ≈ B.β && A.α ≈ B.α && spacescompatible(A.space,B.space)
spacescompatible(A::JacobiWeight,B::Space{DD,RR}) where {DD<:IntervalOrSegment,RR<:Real} =
    spacescompatible(A,JacobiWeight(0,0,B))
spacescompatible(B::Space{DD,RR},A::JacobiWeight) where {DD<:IntervalOrSegment,RR<:Real} =
    spacescompatible(A,JacobiWeight(0,0,B))

transformtimes(f::Fun{JW1},g::Fun{JW2}) where {JW1<:JacobiWeight,JW2<:JacobiWeight}=
            Fun(JacobiWeight(f.space.β+g.space.β,f.space.α+g.space.α,f.space.space),
                coefficients(transformtimes(Fun(f.space.space,f.coefficients),
                                            Fun(g.space.space,g.coefficients))))
transformtimes(f::Fun{JW},g::Fun) where {JW<:JacobiWeight} =
    Fun(f.space,coefficients(transformtimes(Fun(f.space.space,f.coefficients),g)))
transformtimes(f::Fun,g::Fun{JW}) where {JW<:JacobiWeight} =
    Fun(g.space,coefficients(transformtimes(Fun(g.space.space,g.coefficients),f)))

jacobiweight(β,α,x) = -1 ≤ x ≤ 1 ? (1+x)^β*(1-x)^α : zero(x)
jacobiweight(β,α,d::Domain) = Fun(JacobiWeight(β,α,ConstantSpace(d)),[1.])
jacobiweight(β,α) = jacobiweight(β,α,ChebyshevInterval())

weight(sp::JacobiWeight,x) = jacobiweight(sp.β,sp.α,tocanonical(sp,x))
dimension(sp::JacobiWeight) = dimension(sp.space)

Base.first(f::Fun{JW}) where {JW<:JacobiWeight} = space(f).β>0 ? zero(cfstype(f)) : f(leftendpoint(domain(f)))
Base.last(f::Fun{JW}) where {JW<:JacobiWeight} = space(f).α>0 ? zero(cfstype(f)) : f(rightendpoint(domain(f)))

setdomain(sp::JacobiWeight,d::Domain)=JacobiWeight(sp.β,sp.α,setdomain(sp.space,d))

# we assume that points avoids singularities

##TODO: paradigm for same space
function coefficients(f::AbstractVector,sp1::JacobiWeight{SJ1,DD},sp2::JacobiWeight{SJ2,DD}) where {SJ1,SJ2,DD<:IntervalOrSegment}
    β,α=sp1.β,sp1.α
    c,d=sp2.β,sp2.α

    if isapprox(c,β) && isapprox(d,α)
        # remove wrapper spaces and then convert
        coefficients(f,sp1.space,sp2.space)
    else
        # go back to default
        defaultcoefficients(f,sp1,sp2)
    end
end
coefficients(f::AbstractVector,sp::JacobiWeight{SJ,DD},
             S2::SubSpace{S,IT,DD,RR}) where {SJ,S,IT,DD<:IntervalOrSegment,RR<:Real} = subspace_coefficients(f,sp,S2)
coefficients(f::AbstractVector,S2::SubSpace{S,IT,DD,RR},
             sp::JacobiWeight{SJ,DD}) where {SJ,S,IT,DD<:IntervalOrSegment,RR<:Real} = subspace_coefficients(f,S2,sp)
#TODO: it could be possible that we want to JacobiWeight a SumSpace....
coefficients(f::AbstractVector,sp::JacobiWeight{SJ,DD},S2::SumSpace{SV,DD,RR}) where {SJ,SV,DD<:IntervalOrSegment,RR<:Real} =
    sumspacecoefficients(f,sp,S2)
coefficients(f::AbstractVector,sp::JacobiWeight{SJ,Segment{Vec{2,TT}}},S2::TensorSpace{SV,TTT,DD}) where {SJ,TT,SV,TTT,DD<:Domain2d} =
    coefficients(f,sp,JacobiWeight(0,0,S2))

coefficients(f::AbstractVector,sp::JacobiWeight{SJ,DD},S2::Space{DD,RR}) where {SJ,DD<:IntervalOrSegment,RR<:Real} =
    coefficients(f,sp,JacobiWeight(0,0,S2))
coefficients(f::AbstractVector,sp::ConstantSpace{DD},ts::JacobiWeight{SJ,DD}) where {SJ,DD<:IntervalOrSegment} =
    f.coefficients[1]*ones(ts).coefficients
coefficients(f::AbstractVector,S2::Space{DD,RR},sp::JacobiWeight{SJ,DD}) where {SJ,DD<:IntervalOrSegment,RR<:Real} =
    coefficients(f,JacobiWeight(0,0,S2),sp)

"""
`increase_jacobi_parameter(f)` multiplies by `1-x^2` on the unit interval.
`increase_jacobi_parameter(-1,f)` multiplies by `1+x` on the unit interval.
`increase_jacobi_parameter(+1,f)` multiplies by `1-x` on the unit interval.
On other domains this is accomplished by mapping to the unit interval.
"""
increase_jacobi_parameter(f) = Fun(f,JacobiWeight(f.space.β+1,f.space.α+1,space(f).space))
increase_jacobi_parameter(s,f) = s==-1 ? Fun(f,JacobiWeight(f.space.β+1,f.space.α,space(f).space)) :
                                       Fun(f,JacobiWeight(f.space.β,f.space.α+1,space(f).space))

function canonicalspace(S::JacobiWeight)
    if isapprox(S.β,0) && isapprox(S.α,0)
        canonicalspace(S.space)
    else
        #TODO: promote singularities?
        JacobiWeight(S.β,S.α,canonicalspace(S.space))
    end
end

function union_rule(A::ConstantSpace,B::JacobiWeight{P}) where P<:PolynomialSpace
    # we can convert to a space that contains contants provided
    # that the parameters are integers
    # when the parameters are -1 we keep them
    if isapproxinteger(B.β) && isapproxinteger(B.α)
        JacobiWeight(min(B.β,0.),min(B.α,0.),B.space)
    else
        NoSpace()
    end
end

## Algebra

function /(c::Number,f::Fun{JW}) where JW<:JacobiWeight
    g=c/Fun(space(f).space,f.coefficients)
    Fun(JacobiWeight(-f.space.β,-f.space.α,space(g)),g.coefficients)
end

function ^(f::Fun{JW}, k::AbstractFloat) where JW<:JacobiWeight
    S=space(f)
    g=Fun(S.space,coefficients(f))^k
    Fun(JacobiWeight(k*S.β,k*S.α,space(g)),coefficients(g))
end

function *(f::Fun{JW1},g::Fun{JW2}) where {JW1<:JacobiWeight,JW2<:JacobiWeight}
    @assert domainscompatible(f,g)
    fβ,fα=f.space.β,f.space.α
    gβ,gα=g.space.β,g.space.α
    m=(Fun(space(f).space,f.coefficients).*Fun(space(g).space,g.coefficients))
    if isapprox(fβ+gβ,0)&&isapprox(fα+gα,0)
        m
    else
        Fun(JacobiWeight(fβ+gβ,fα+gα,space(m)),m.coefficients)
    end
end

/(f::Fun{JW1},g::Fun{JW2}) where {JW1<:JacobiWeight,JW2<:JacobiWeight}=f*(1/g)

# O(min(m,n)) Ultraspherical conjugated inner product

function conjugatedinnerproduct(sp::Ultraspherical,u::AbstractVector{S},v::AbstractVector{V}) where {S,V}
    λ=order(sp)
    if λ==1
        mn = min(length(u),length(v))
        if mn > 0
            return dotu(u[1:mn],v[1:mn])*π/2
        else
            return zero(promote_type(eltype(u),eltype(v)))
        end
    else
        T,mn = promote_type(S,V),min(length(u),length(v))
        if mn > 1
            wi = sqrt(convert(T,π))*gamma(λ+one(T)/2)/gamma(λ+one(T))
            ret = u[1]*wi*v[1]
            for i=2:mn
              wi *= (i-2one(T)+2λ)/(i-one(T)+λ)*(i-2one(T)+λ)/(i-one(T))
              ret += u[i]*wi*v[i]
            end
            return ret
        elseif mn > 0
            wi = sqrt(convert(T,π))*gamma(λ+one(T)/2)/gamma(λ+one(T))
            return u[1]*wi*v[1]
        else
            return zero(promote_type(eltype(u),eltype(v)))
        end
    end
end

function conjugatedinnerproduct(::Chebyshev,u::AbstractVector,v::AbstractVector)
    mn = min(length(u),length(v))
    if mn > 1
        return (2u[1]*v[1]+dotu(u[2:mn],v[2:mn]))*π/2
    elseif mn > 0
        return u[1]*v[1]*π
    else
        return zero(promote_type(eltype(u),eltype(v)))
    end
end

function bilinearform(f::Fun{JacobiWeight{Ultraspherical{LT,D,R},D,R,TT}},g::Fun{Ultraspherical{LT,D,R}}) where {LT,D,R,TT}
    d = domain(f)
    @assert d == domain(g)
    λ = order(space(f).space)
    if order(space(g)) == λ && f.space.β == f.space.α == λ-0.5
        return complexlength(d)/2*conjugatedinnerproduct(Ultraspherical(λ,d),f.coefficients,g.coefficients)
    else
        return defaultbilinearform(f,g)
    end
end

function bilinearform(f::Fun{Ultraspherical{LT,D,R}},
                    g::Fun{JacobiWeight{Ultraspherical{LT,D,R},D,R,TT}}) where {LT,D,R,TT}
    d = domain(f)
    @assert d == domain(g)
    λ = order(space(f))
    if order(space(g).space) == λ && g.space.β == g.space.α == λ-0.5
        return complexlength(d)/2*conjugatedinnerproduct(Ultraspherical(λ,d),f.coefficients,g.coefficients)
    else
        return defaultbilinearform(f,g)
    end
end

function bilinearform(f::Fun{JacobiWeight{Ultraspherical{LT,D,R},D,R,TT1}},
                    g::Fun{JacobiWeight{Ultraspherical{LT,D,R},D,R,TT2}}) where {LT,D,R,TT1,TT2}
    d = domain(f)
    @assert d == domain(g)
    λ = order(space(f).space)
    if order(space(g).space) == λ && f.space.β+g.space.β == f.space.α+g.space.α == λ-0.5
        return complexlength(domain(f))/2*conjugatedinnerproduct(Ultraspherical(λ,d),f.coefficients,g.coefficients)
    else
        return defaultbilinearform(f,g)
    end
end

function linebilinearform(f::Fun{JacobiWeight{Ultraspherical{LT,D,R},D,R,TT}},g::Fun{Ultraspherical{LT,D,R}}) where {LT,D,R,TT}
    d = domain(f)
    @assert d == domain(g)
    λ = order(space(f).space)
    if order(space(g)) == λ && f.space.β == f.space.α == λ-0.5
        return arclength(d)/2*conjugatedinnerproduct(Ultraspherical(λ,d),f.coefficients,g.coefficients)
    else
        return defaultlinebilinearform(f,g)
    end
end

function linebilinearform(f::Fun{Ultraspherical{LT,D,R}},g::Fun{JacobiWeight{Ultraspherical{LT,D,R},D,R,TT}}) where {LT,D,R,TT}
    d = domain(f)
    @assert d == domain(g)
    λ = order(space(f))
    if order(space(g).space) == λ &&  g.space.β == g.space.α == λ-0.5
        return arclength(d)/2*conjugatedinnerproduct(Ultraspherical(λ,d),f.coefficients,g.coefficients)
    else
        return defaultlinebilinearform(f,g)
    end
end

function linebilinearform(f::Fun{JacobiWeight{Ultraspherical{LT,D,R},D,R,TT1}},g::Fun{JacobiWeight{Ultraspherical{LT,D,R},D,R,TT2}}) where {LT,D,R,TT1,TT2}
    d = domain(f)
    @assert d == domain(g)
    λ = order(space(f).space)
    if order(space(g).space) == λ &&  f.space.β+g.space.β == f.space.α+g.space.α == λ-0.5
        return arclength(d)/2*conjugatedinnerproduct(Ultraspherical(λ,d),f.coefficients,g.coefficients)
    else
        return defaultlinebilinearform(f,g)
    end
end


text:
@testset "reference counting tests" begin
  let x = K(1)
    @test auto_r0(K_new, x) do p; xr(p) == 1 end
    @test xr(kpointer(x)) == 0
  end
end

ext: jl


text:
using Underscores
using Utils

struct Password
    counts::UnitRange{UInt8}
    character::Char
    password::String
end

function load()
    reg = r"(\d+)-(\d+) (\w): (\w+)"
    re = Password[]
    open(joinpath(@__DIR__, "02_passwords.txt"), "r") do f
        for line in eachline(f)
            f, t, c, p = match(reg, line).captures
            push!(re, Password(parse(UInt8, f):parse(UInt8, t), c[1], p))
        end
    end
    re
end

function isvalid(pw::Password)
    (@_ count(_ == pw.character, pw.password)) in pw.counts
end

function isvalid_positions(pw::Password)
    (pw.password[pw.counts.start] == pw.character) ⊻
    (pw.password[pw.counts.stop] == pw.character)
end

function count_valid(which)
    @_ count(which(_), pws)
end

const pws = load()
println("Valid passwords variant 1: $(count_valid(isvalid))/$(length(pws))")
println("Valid passwords variant 2: $(count_valid(isvalid_positions))/$(length(pws))")

using Test

@testset "Password problem" begin
    only_fst2 = Password(2:4, 'a', "bbbbbbbbbaa")
    only_fst3 = Password(2:4, 'a', "bbbbbbbbbaaa")
    only_fst4 = Password(2:4, 'a', "bbbbbbbbbaaaa")
    only_snd1 = Password(2:4, 'a', "babcbbbb")
    only_snd2 = Password(2:4, 'a', "bcbabbbb")

    @testset "variant 1" begin
        @test isvalid(only_fst2)
        @test isvalid(only_fst3)
        @test isvalid(only_fst4)
        @test !isvalid(only_snd1)
        @test !isvalid(only_snd2)
    end

    @testset "variant 2" begin
        @test !isvalid_positions(only_fst2)
        @test !isvalid_positions(only_fst3)
        @test !isvalid_positions(only_fst4)
        @test isvalid_positions(only_snd1)
        @test isvalid_positions(only_snd2)
    end

    @testset "solutions" begin
        @test count_valid(isvalid) == 474
        @test count_valid(isvalid_positions) == 745
    end
end

max_stars_repo_name: CmdQ/AoC


text:
using Test, MPI

using MPIPreferences

# load test packages to trigger precompilation
using DoubleFloats
if get(ENV, "JULIA_MPI_TEST_ARRAYTYPE", "") == "CuArray"
    import CUDA
    CUDA.version()
    CUDA.precompile_runtime()
    ArrayType = CUDA.CuArray
else
    ArrayType = Array
end

nprocs_str = get(ENV, "JULIA_MPI_TEST_NPROCS", "")
nprocs = nprocs_str == "" ? clamp(Sys.CPU_THREADS, 2, 4) : parse(Int, nprocs_str)

@info "Running MPI tests" ArrayType nprocs MPIPreferences.abi MPIPreferences.binary

if haskey(ENV,"JULIA_MPI_TEST_BINARY")
    @test ENV["JULIA_MPI_TEST_BINARY"] == MPIPreferences.binary
end
if haskey(ENV,"JULIA_MPI_TEST_ABI")
    @test ENV["JULIA_MPI_TEST_ABI"] == MPIPreferences.abi
end

if Sys.isunix()
    # This test doesn't need to be run with mpiexec.  `mpiexecjl` is currently
    # available only on Unix systems
    include("mpiexecjl.jl")
end

testdir = @__DIR__
istest(f) = endswith(f, ".jl") && startswith(f, "test_")
testfiles = sort(filter(istest, readdir(testdir)))

@testset "$f" for f in testfiles
    mpiexec() do mpirun
        cmd(n=nprocs) = `$mpirun -n $n $(Base.julia_cmd()) $(joinpath(testdir, f))`
        if f == "test_spawn.jl"
            # Some command as the others, but always use a single process
            run(cmd(1))
        elseif f == "test_threads.jl"
            withenv("JULIA_NUM_THREADS" => "4") do
                run(cmd())
            end
        elseif f == "test_error.jl"
            r = run(ignorestatus(cmd()))
            @test !success(r)
        elseif f == "test_errorhandler.jl" && (MPI.identify_implementation()[1] == MPI.UnknownMPI ||
            # Fujitsu MPI is known to not work with custom error handlers
            startswith(MPI.MPI_LIBRARY_VERSION_STRING, "FUJITSU MPI"))
            try
                run(cmd())
            catch e
                @error """
                       $(f) tests failed.  This may due to the fact this implementation of MPI doesn't support custom error handlers.
                       See the full error message for more details.  Some messages may have been written above.
                       """ exception=(e, catch_backtrace())
                @test_broken false
            end
        else
            # MPI_Reduce with MPICH 3.4.2 on macOS when root != 0 and
            # when recvbuf == C_NULL segfaults
            # <https://github.com/pmodels/mpich/issues/5700>
            if get(ENV, "JULIA_MPI_TEST_DISABLE_REDUCE_ON_APPLE", "") != "" && Sys.isapple() && f == "test_reduce.jl"
                return
            end
            run(cmd())
        end
        @test true
    end
end


text:
"""
    loadArrayFromFile(path, type)

Load binary data from a file `path` and interpret the data as an array of type `type`.

# Examples
```julia-repl
julia> data = loadArrayFromFile("data.bin", Float32)
```
"""
function loadArrayFromFile(path::String, T::DataType)::Vector{T}
    if !isfile(path)
        throw(ErrorException("file not found."))
    end

    fileSize = filesize(path);
    numElems = fileSize÷sizeof(T);
    if fileSize != numElems*sizeof(T)
        println("Warn: file size is not a multiple of element size.\nLoading truncated data.")
    end

    array = Vector{T}(undef, numElems);
    read!(path, array);
    array
end


text:
using Julog

include("population.jl")
include("zebra.jl")


text:
module Helmholtz #hide
# # Helmholtz equation
# In this example, we use the SLEPc to find the eigenvalues of the following Helmholtz equation:
# ``u'' + \omega^2 u = 0`` associated to Dirichlet boundary conditions on the domain ``[0,1]``. Hence
# the theoritical eigenvalues are ``\omega = k \pi`` with ``k \in \mathbb{Z}^*``; and the associated
# eigenvectors are ``u(x) = \sin(k\pix)``.
# A centered finite difference scheme is used for the spatial discretization.
#
# The equation is written in matrix form ``Au = \alpha Bu`` where ``\alpha = \omega^2``.
#
# To run this example, simplfy excute `mpirun -n your_favourite_integer julia helmholtz_FD.jl`
#
# In this example, PETSc/SLEPc legacy method names are used. For more fancy names, check the next example.
#
# Note that the way we achieve things in the document can be highly improved and the purpose of this example
# is only demonstrate some method calls to give an overview.
#
# Start by importing both `PetscWrap`, for the distributed matrices, and `SlepcWrap` for the eigenvalues.
using PetscWrap
using SlepcWrap

# Number of mesh points and mesh step
n = 21
Δx = 1. / (n - 1)

# Initialize SLEPc. Either without arguments, calling `SlepcInitialize()` or using "command-line" arguments.
# To do so, either provide the arguments as one string, for instance
# `SlepcInitialize("-eps_max_it 100 -eps_tol 1e-5")` or provide each argument in
# separate strings : `PetscInitialize(["-eps_max_it", "100", "-eps_tol", "1e-5")`.
# Here we ask for the five closest eigenvalues to ``0``, using a non-zero pivot for the LU factorization and a
# "shift-inverse" process.
SlepcInitialize("-eps_target 0 -eps_nev 5 -st_pc_factor_shift_type NONZERO -st_type sinvert")

# Create the problem matrices, set sizes and apply "command-line" options. Note that we should
# set the number of preallocated non-zeros to increase performance.
A = MatCreate()
B = MatCreate()
MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n)
MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, n, n)
MatSetFromOptions(A)
MatSetFromOptions(B)
MatSetUp(A)
MatSetUp(B)

# Get rows handled by the local processor
A_rstart, A_rend = MatGetOwnershipRange(A)
B_rstart, B_rend = MatGetOwnershipRange(B)

# Fill matrix A  with second order derivative central scheme
for i in A_rstart:A_rend-1
    if (i == 0)
        MatSetValues(A, [0], [0, 1], [-2., 1] / Δx^2, INSERT_VALUES)
    elseif (i == n-1)
        MatSetValues(A, [n-1], [n-2, n-1], [1., -2.] / Δx^2, INSERT_VALUES)
    else
        MatSetValues(A, [i], i-1:i+1, [1., -2., 1.] / Δx^2, INSERT_VALUES)
    end
end

# Fill matrix B with identity matrix
for i in B_rstart:B_rend-1
    MatSetValue(B, i, i, -1., INSERT_VALUES)
end

# Set boundary conditions : u(0) = 0 and u(1) = 0. Only the processor
# handling the corresponding rows are playing a role here.
(A_rstart == 0) && MatSetValues(A, [0], [0,1], [1., 0.], INSERT_VALUES)
(B_rstart == 0) && MatSetValue(B, 0, 0, 0., INSERT_VALUES)

(A_rend == n) && MatSetValues(A, [n-1], [n-2,n-1], [0., 1.], INSERT_VALUES)
(B_rend == n) && MatSetValue(B, n-1, n-1, 0., INSERT_VALUES)

# Assemble the matrices
MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)
MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)
MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)
MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)

# Now we set up the eigenvalue solver
eps = EPSCreate()
EPSSetOperators(eps, A, B)
EPSSetFromOptions(eps)
EPSSetUp(eps)

# Then we solve
EPSSolve(eps)

# And finally we can inspect the solution. Let's first get the number of converged eigenvalues:
nconv = EPSGetConverged(eps)

# Then we can get/display these eigenvalues (more precisely their square root, i.e ``\simeq \omega``)
for ieig in 0:nconv - 1
    vpr, vpi = EPSGetEigenvalue(eps, ieig)
    @show √(vpr), √(vpi)
end

# We can also play with eigen vectors. First, create two Petsc vectors to allocate memory
vecr, veci = MatCreateVecs(A)

# Then loop over the eigen pairs and retrieve eigenvectors
for ieig in 0:nconv-1
    vpr, vpi, vecpr, vecpi = EPSGetEigenpair(eps, ieig, vecr, veci)

    ## At this point, you can call VecGetArray to obtain a Julia array (see PetscWrap examples).
    ## If you are on one processor, you can even plot the solution to check that you have a sinus
    ## solution. On multiple processors, this would require to "gather" the solution on one processor only.
end

# Finally, let's free the memory
MatDestroy(A)
MatDestroy(B)
EPSDestroy(eps)

# And call finalize when you're done
SlepcFinalize()

end #hide


text:
module Quadratures

    using ..Utils

    export Quadrature
    export RiemannQuadratureLeft, RiemannQuadratureRight
    export GaussLegendreQuadrature, LobattoLegendreQuadrature,
           GaussChebyshevQuadrature, LobattoChebyshevQuadrature,
           ClenshawCurtisQuadrature
    export quadrature, weights

    include("quadratures/quadrature.jl")

end


text:
"""
$(TYPEDEF)

Introduce direct observations on all dimensions of a Pose2 variable:

Example:
--------
```julia
PriorPose2( MvNormal([10; 10; pi/6.0], Matrix(Diagonal([0.1;0.1;0.05].^2))) )
```
"""
mutable struct PriorPose2{T} <: IncrementalInference.AbstractPrior  where {T <: IncrementalInference.SamplableBelief}
  Z::T
  # empty constructor
  PriorPose2{T}() where T = new{T}()
  # regular constructor
  PriorPose2{T}(x::T) where {T <: IncrementalInference.SamplableBelief}  = new{T}(x)
end
# convencience and default object helper
PriorPose2(x::T) where {T <: IncrementalInference.SamplableBelief} = PriorPose2{T}(x)

function getSample(cfo::CalcFactor{<:PriorPose2}, N::Int=1)
  return (rand(cfo.factor.Z,N), )
end

function (s::CalcFactor{<:PriorPose2})(meas,
                                       wXi)
  #
  iXihat = SE2(meas) \ SE2(wXi)
  return se2vee(iXihat)
end

## Serialization support

"""
$(TYPEDEF)
"""
mutable struct PackedPriorPose2  <: IncrementalInference.PackedInferenceType
    str::String
    PackedPriorPose2() = new()
    PackedPriorPose2(x::String) = new(x)
end
function convert(::Type{PackedPriorPose2}, d::PriorPose2)
  return PackedPriorPose2(convert(PackedSamplableBelief, d.Z))
end
function convert(::Type{PriorPose2}, d::PackedPriorPose2)
  distr = convert(SamplableBelief, d.str)
  return PriorPose2{typeof(distr)}(distr)
end

## NOTE likely deprecated comparitors, see DFG compareFields, compareAll instead
function compare(a::PriorPose2,b::PriorPose2; tol::Float64=1e-10)
  compareDensity(a.Z, b.Z)
end


text:
# This file was generated, do not modify it. # hide
r  = range(model, :(linear_regressor.lambda), lower=1e-2, upper=100_000, scale=:log10)
tm = TunedModel(model=model, ranges=r, tuning=Grid(resolution=50),
                resampling=CV(nfolds=3, rng=4141), measure=rms)
mtm = machine(tm, Xc, y)
fit!(mtm, rows=train)

best_mdl = fitted_params(mtm).best_model
round(best_mdl.linear_regressor.lambda, sigdigits=4)


text:
"""
    DMMethod

Abstract type for nesting the various methods that can be used to perform Diebold-Mariano tests. Subtypes are: \n
    DMHAC
    DMBoot \n
The subtypes have entries in the help (?) menu.
"""
abstract type DMMethod end

"""
    DMHAC(data ; alpha::Float64=0.05, kernelfunction, bandwidth::Int=-1)

Method type for doing a Diebold-Mariano test using a HAC variance estimator. A
constructor that takes the data and three keyword arguments is provided. Relevant
keyword arguments are: \n
    alpha <- The confidence level for the test \n
    kernelfunction <- Kernel function to use in HAC variance estimator. Valid values are
        KernelEpanechnikov(), KernelGaussian(), KernelUniform(), KernelBartlett(). You can
        also :epanechnikov, :gaussian, :uniform, :bartlett, or String variants of these symbols. \n
    bandwidth <- Bandwidth to use in HAC variance estimator (set less than or equal to -1 to estimate bandwidth using Politis (2003) "Adaptive Bandwidth Choice") \n
"""
struct DMHAC{Tk<:KernelFunction} <: DMMethod
    alpha::Float64
    kernelfunction::Tk
    bandwidth::Int
    function DMHAC(alpha::Float64, kernelfunction::Tk, bandwidth::Int) where {Tk<:KernelFunction}
        !(0.0 < alpha < 0.5) && error("Confidence level set to $(alpha) which is not on the (0, 0.5) interval")
        new{Tk}(alpha, kernelfunction, bandwidth)
    end
end
DMHAC( ; alpha::Float64=0.05, kernelfunction=KernelEpanechnikov(), bandwidth::Int=-1) = DMHAC(alpha, get_kernel_type(kernelfunction), bandwidth)
DMHAC(data ; alpha::Float64=0.05, kernelfunction=KernelEpanechnikov(), bandwidth::Int=-1) = DMHAC(alpha, get_kernel_type(kernelfunction), bandwidth) #Superfluous but included for consistency

"""
    DMBoot(data ; alpha, bootinput_kwargs...)

Method type for doing a Diebold-Mariano test using a dependent bootstrap. A constructor
that takes the data and several keyword arguments is provided. Relevant keyword arguments are: \n
    alpha=0.05 <- The confidence level to use in the test \n
    blocklength=0.0 <- Block length to use with the bootstrap. Default of 0.0 implies
        the block length will be optimally estimated from the data use the method deemed
        most appropriate by the DependentBootstrap package (typically the selection procedure
        of Patton, Politis, and White (2009)). \n
    numresample=1000 <- Number of resamples to use when bootstrapping. \n
    bootmethod=:stationary <- Bootstrap methodology to use, where the default is the
        stationary bootstrap of Politis and Romano (1993) \n
See the DependentBootstrap package docs for more info on bootstrap input keyword arguments.
"""
struct DMBoot <: DMMethod
    alpha::Float64
    bootinput::BootInput
    function DMBoot(alpha::Float64, bootinput::BootInput)
        !(0.0 < alpha < 0.5) && error("Confidence level set to $(alpha) which is not on the (0, 0.5) interval")
        bootinput.flevel1 != mean && error("flevel1 field in BootInput must be set equal to mean. Please use the keyword constructor for DMBoot as this automatically enforces this behaviour: $(bootinput.flevel1)")
        new(alpha, bootinput)
    end
end
DMBoot(data ; alpha::Float64=0.05, kwargs...) = DMBoot(alpha, BootInput(data ; flevel1=mean, kwargs...))

"""
    DMTest(rejH0::Int, pvalue::Float64, bestinput::Int, teststat::Float64, dmmethod::DMMethod)

Output type for a Diebold-Mariano test. A description of the fields follows: \n
    rejH0 <- true if the null is rejected, false otherwise
    pvalue <- p-value from the test
    bestinput <- 1 if forecast 1 is more accurate, and 2 if forecast 2 is more accurate. See ?dm for definition of forecast 1 versus 2.
    teststat <- If dmmethod == :hac then is the mean of loss difference scaled by HAC variance
                If dmmethod == :boot then is the mean of loss difference
    dmmethod <- Diebold-Mariano method used in the test. See ?DMMethod for more detail.
"""
struct DMTest
    rejH0::Bool
    pvalue::Float64
    bestinput::Int
    teststat::Float64
    dmmethod::DMMethod
    DMTest(rejH0::Bool, pvalue::Float64, bestinput::Int, teststat::Float64, dmmethod::DMMethod) = new(rejH0, pvalue, bestinput, teststat, dmmethod)
end
function Base.show(io::IO, a::DMTest)
    println(io, "Diebold-Mariano test results:")
    println(io, "    Reject null of equal predictive ability: $(a.rejH0)")
    println(io, "    p-value: $(a.pvalue)")
    println(io, "    Preferred forecast: $(a.bestinput)")
end

"""
    dm(lossdiff::Vector{<:Number}, dmmethod ; kwargs...)

This function implements the test proposed in Diebold, Mariano (1995) "Comparing Predictive Accuracy". \n
Let x_1 denote forecast 1, x_2 denote forecast 2, and let y denote the forecast target. Let L(., .) denote a loss function.
Then the first argument lossdiff is assumed to be a vector created by the following operation: \n
L(x_1, y) - L(x_2, y) \n
The second argument, dmmethod, can be an explicit method type, currently DMHAC or DMBoot,
(see ?DMHAC and ?DMBoot for more detail), in which case the keyword arguments are not needed. \n
Alternatively, dmmethod can be set to the Symbol :hac or :boot, depending on whether the user
wants to use the HAC method or the bootstrap method. In this instance, the keyword arguments
provided will be passed on to DMHAC or DMBoot constructors (see ?DMHAC and ?DMBoot for more detail). \n
Finally, the second argument can be omitted entirely, in which case the method will default to
the default bootstrap method, which is the stationary bootstrap of Politis and Romano (1993) with
block length estimated followed Patton, Politis, and White (2009). \n
The output of a Diebold-Mariano test is of type DMTest. Use ?DMTest for more information.
"""
function dm(lossDiff::Vector{<:Number}, method::DMHAC)::DMTest
	length(lossDiff) < 2 && error("Input data vector has length of $(length(lossDiff))")
	m = mean(lossDiff)
	(v, _) = hacvariance(lossDiff, kf=method.kernelfunction, bw=method.bandwidth)
	testStat = m / sqrt(v / length(lossDiff))
	pVal = pvaluelocal(Normal(), testStat, tail=:both)
	pVal > method.alpha ? (rejH0 = false) : (rejH0 = true)
    testStat > 0 ? (bestInput = 2) : (bestInput = 1)
	return DMTest(rejH0, pVal, bestInput, testStat, method)
end
function dm(lossDiff::Vector{<:Number}, method::DMBoot)::DMTest
    length(lossDiff) < 2 && error("Input data vector has length of $(length(lossDiff))")
    method.bootinput.flevel1 != mean && error("flevel1 field in BootInput must be set equal to mean. Please use the keyword constructor for DMBoot as this automatically enforces this behaviour: $(bootinput.flevel1)")
    statVec = dbootlevel1(lossDiff, method.bootinput)
    pVal = pvaluelocal(statVec, 0.0, tail=:both, as=false)
    pVal > method.alpha ? (rejH0 = false) : (rejH0 = true)
    testStat = mean(lossDiff)
    testStat > 0 ? (bestInput = 2) : (bestInput = 1)
	return DMTest(rejH0, pVal, bestInput, testStat, method)
end
#Keyword constructor
function dm(lossDiff::Vector{<:Number}, dmmethod::Symbol ; kwargs...)::DMTest
    dmmethod == :hac && return dm(lossDiff, DMHAC(lossDiff ; kwargs...))
    (dmmethod == :boot || dmmethod == :bootstrap) && return dm(lossDiff, DMBoot(lossDiff ; kwargs...))
    error("Invalid dmmethod: $(dmmethod)")
end
dm(lossDiff::Vector{<:Number})::DMTest = dm(lossDiff, DMBoot(lossDiff))


text:
## Update all hyperparameters for the full batch GP models ##
function  update_hyperparameters!(model::Union{GP,VGP,VStP})
    Jnn = kernelderivativematrix.([model.X],model.kernel)
    f_l,f_v,f_μ₀ = hyperparameter_gradient_function(model)
    grads_l = map(compute_hyperparameter_gradient,model.kernel,fill(f_l,model.nPrior),Jnn,1:model.nPrior)
    grads_v = map(f_v,model.kernel,1:model.nPrior)
    grads_μ₀ = map(f_μ₀,1:model.nPrior)

    apply_gradients_lengthscale!.(model.kernel,grads_l) #Send the derivative of the matrix to the specific gradient of the model
    apply_gradients_variance!.(model.kernel,grads_v) #Send the derivative of the matrix to the specific gradient of the model
    update!.(model.μ₀,grads_μ₀)

    model.inference.HyperParametersUpdated = true
end

## Update all hyperparameters for the sparse variational GP models ##
function update_hyperparameters!(model::SVGP{T}) where {T<:Real}
    Jmm = kernelderivativematrix.(model.Z,model.kernel)
    Jnm = kernelderivativematrix.([model.inference.x],model.Z,model.kernel)
    Jnn = kernelderivativediagmatrix.([model.inference.x],model.kernel)
    f_l,f_v,f_μ₀ = hyperparameter_gradient_function(model)
    grads_l = map(compute_hyperparameter_gradient,model.kernel,fill(f_l,model.nPrior),Jmm,Jnm,Jnn,collect(1:model.nPrior))
    grads_v = map(f_v,model.kernel,1:model.nPrior)
    if !isnothing(model.Zoptimizer)
        Z_gradients = inducingpoints_gradient(model) #Compute the gradient given the inducing points location
        model.Z += GradDescent.update(model.Zoptimizer,Z_gradients) #Apply the gradients on the location
    end
    grads_μ₀ = map(f_μ₀,1:model.nPrior)
    apply_gradients_lengthscale!.(model.kernel,grads_l)
    apply_gradients_variance!.(model.kernel,grads_v)
    update!.(model.μ₀,grads_μ₀)
    model.inference.HyperParametersUpdated = true
end

## Return the derivative of the KL divergence between the posterior and the GP prior ##
function hyperparameter_KL_gradient(J::AbstractMatrix{T},A::AbstractMatrix{T}) where {T<:Real}
    return 0.5*opt_trace(J,A)
end

## Return functions computing gradients of the ELBO given the kernel hyperparameters for a non-sparse model ##
function hyperparameter_gradient_function(model::VGP{T}) where {T<:Real}
    A = ([Diagonal{T}(I,model.nFeatures)].-model.invKnn.*(model.Σ.+(model.µ.-model.μ₀).*transpose.(model.μ.-model.μ₀))).*model.invKnn
    if model.IndependentPriors
        return (function(Jnn,index)
                    return -hyperparameter_KL_gradient(Jnn,A[index])
                end,
                function(kernel,index)
                    return -1.0/getvariance(kernel)*hyperparameter_KL_gradient(model.Knn[index],A[index])
                end,
                function(index)
                    return -model.invKnn[index]*(model.μ₀[index]-model.μ[index])
                end)
    else
        return (function(Jnn,index)
            return -sum(hyperparameter_KL_gradient.([Jnn],A))
                end,
                function(kernel,index)
                    return -1.0/getvariance(kernel)*sum(hyperparameter_KL_gradient.(model.Knn,A))
                end,
                function(index)
                    return -sum(model.invKnn.*(model.μ₀.-model.μ))
                end)
    end
end

function hyperparameter_local_gradient_function(model::VGP{T}) where {T<:Real}
    A = ([Diagonal{T}(I,model.nFeatures)].-model.likelihood.invK.*(model.likelihood.Σ.+(model.likelihood.µ.-model.likelihood.μ₀).*transpose.(model.likelihood.μ.-model.likelihood.μ₀))).*model.likelihood.invK
    return (function(Jnn,index)
                return -hyperparameter_KL_gradient(Jnn,A[index])
            end,
            function(kernel,index)
                return -1.0/getvariance(kernel)*hyperparameter_KL_gradient(model.Knn[index],A[index])
            end,
            function(index)
                return -model.likelihood.invK[index]*(model.likelihood.μ₀[index]-model.likelihood.μ[index])
            end)
end

## Return functions computing gradients of the ELBO given the kernel hyperparameters for a non-sparse model ##
function hyperparameter_gradient_function(model::SVGP{T}) where {T<:Real}
    A = ([Diagonal{T}(I,model.nFeatures)].-model.invKmm.*(model.Σ.+model.µ.*transpose.(model.μ))).*model.invKmm
    ι = Matrix{T}(undef,model.inference.nSamplesUsed,model.nFeatures) #Empty container to save data allocation
    κΣ = model.κ.*model.Σ
    if model.IndependentPriors
        return (function(Jmm::Symmetric{T,Matrix{T}},Jnm::Matrix{T},Jnn::Vector{T},index::Int)
                    return (hyperparameter_expec_gradient(model,ι,κΣ[index],Jmm,Jnm,Jnn,index)-hyperparameter_KL_gradient(Jmm,A[index]))
                end,
                function(kernel::Kernel{T},index::Int)
                    return 1.0/getvariance(kernel)*(
                            - model.inference.ρ*dot(∇E_Σ(model,index),model.K̃[index])
                            - hyperparameter_KL_gradient(model.Kmm[index],A[index]))
                end,
                function(index)
                    return -model.invKmm[index]*(model.μ₀[index]-model.μ[index])
                end)
    else
        return (function(Jmm::Symmetric{T,Matrix{T}},Jnm::Matrix{T},Jnn::Vector{T},index::Int)
                    return  (hyperparameter_expec_gradient(model,ι,κΣ,Jmm,Jnm,Jnn)
                           - sum(hyperparameter_KL_gradient.([Jmm],A)))
                end,
                function(kernel::Kernel{T},index::Int)
                    return 1.0/getvariance(kernel)*(model.inference.ρ*sum(
                            -dot(∇E_Σ(model,i),model.K̃[1]) for i in 1:model.nLatent)
                            - sum(hyperparameter_KL_gradient.(model.Kmm,A)))
                end,
                function(index)
                    return -sum(model.invKmm.*(model.μ₀.-model.μ))
                end)
    end
end

function hyperparameter_gradient_function(model::VStP{T}) where {T<:Real}
    A = ([Diagonal{T}(I,model.nFeatures)].-invK(model).*(model.Σ.+(model.µ.-model.μ₀).*transpose.(model.μ.-model.μ₀))).*invK(model)
    if model.IndependentPriors
        return (function(Jnn,index)
                    return -hyperparameter_KL_gradient(Jnn,A[index])
                end,
                function(kernel,index)
                    return -1.0/getvariance(kernel)*hyperparameter_KL_gradient(model.Knn[index].*model.χ[index],A[index])
                end,
                function(index)
                    return -invK(model,index)*(model.μ₀[index]-model.μ[index])
                end)
    else
        return (function(Jnn,index)
            return -sum(hyperparameter_KL_gradient.([Jnn],A))
                end,
                function(kernel,index)
                    return -1.0/getvariance(kernel)*sum(hyperparameter_KL_gradient.(model.Knn[index].*model.χ[index],A) for index in 1:model.nPrior)
                end,
                function(index)
                    return -sum(invK(model).*(model.μ₀.-model.μ))
                end)
    end
end

## Gradient with respect to hyperparameter with independent priors for analytical VI ##
function hyperparameter_expec_gradient(model::SVGP{T,<:Likelihood,<:AnalyticVI},ι::Matrix{T},κΣ::Matrix{T},Jmm::Symmetric{T,Matrix{T}},Jnm::Matrix{T},Jnn::Vector{T},index::Integer) where {T<:Real}
    mul!(ι,(Jnm-model.κ[index]*Jmm),model.invKmm[index])
    Jnn .-= opt_diag(ι,model.Knm[index]) + opt_diag(model.κ[index],Jnm)
    dμ = dot(∇E_μ(model,index),ι*model.μ[index])
    dΣ =  -dot(∇E_Σ(model,index),Jnn)
    dΣ += -dot(∇E_Σ(model,index),2.0*(opt_diag(ι*model.Σ[index],model.κ[index])))
    dΣ += -dot(∇E_Σ(model,index),2.0*(ι*model.μ[index]).*(model.κ[index]*model.μ[index]))
    return model.inference.ρ*(dμ+dΣ)
end

## Gradient with respect to hyperparameter with independent priors for numerical VI ##
function hyperparameter_expec_gradient(model::SVGP{T,<:Likelihood,<:NumericalVI},ι::Matrix{T},κΣ::Matrix{T},Jmm::Symmetric{T,Matrix{T}},Jnm::Matrix{T},Jnn::Vector{T},index::Integer) where {T<:Real}
    mul!(ι,(Jnm-model.κ[index]*Jmm),model.invKmm[index])
    Jnn .-= opt_diag(ι,model.Knm[index]) + opt_diag(model.κ[index],Jnm)
    dμ = dot(∇E_μ(model,index),ι*model.μ[index])
    dΣ = dot(∇E_Σ(model,index),Jnn+2.0*(opt_diag(ι*model.Σ[index],model.κ[index])))
    return model.inference.ρ*(dμ+dΣ)
end

## Gradient with respect to hyperparameter with shared priors for analytical VI ##
function hyperparameter_expec_gradient(model::SVGP{T,<:Likelihood,<:AnalyticVI},ι::Matrix{T},κΣ::Vector{Matrix{T}},Jmm::Symmetric{T,Matrix{T}},Jnm::Matrix{T},Jnn::Vector{T}) where {T<:Real}
    mul!(ι,(Jnm-model.κ[1]*Jmm),model.invKmm[1])
    Jnn .-= opt_diag(ι,model.Knm[1]) + opt_diag(model.κ[1],Jnm)
    dμ = sum(dot(∇E_μ(model,i),ι*model.μ[i]) for i in 1:model.nLatent)
    dΣ = -0.5*sum(dot(∇E_Σ(model,i),Jnn+2.0*opt_diag(ι,κΣ[i])) for i in 1:model.nLatent)
    dΣ += -0.5*sum(dot(∇E_Σ(model,i),2.0*(ι*model.μ[i]).*(model.κ[1]*model.μ[i])) for i in 1:model.nLatent)
    return model.inference.ρ*(dμ+dΣ)
end

## Gradient with respect to hyperparameter with shared priors for numerical VI ##
function hyperparameter_expec_gradient(model::SVGP{T,<:Likelihood,<:NumericalVI},ι::Matrix{T},κΣ::Vector{Matrix{T}},Jmm::Symmetric{T,Matrix{T}},Jnm::Matrix{T},Jnn::Vector{T}) where {T<:Real}
    mul!(ι,(Jnm-model.κ[1]*Jmm),model.invKmm[1])
    Jnn .-= opt_diag(ι,model.Knm[1]) + opt_diag(model.κ[1],Jnm)
    dμ = sum(dot(∇E_μ(model,i),ι*model.μ[i]) for i in 1:model.nLatent)
    dΣ = sum(dot(∇E_Σ(model,i),Jnn+2.0*opt_diag(ι,κΣ[i])) for i in 1:model.nLatent)
    return model.inference.ρ*(dμ+dΣ)
end

## Return a function computing the gradient of the ELBO given the inducing point locations ##
function inducingpoints_gradient(model::SVGP{T}) where {T<:Real}
    if model.IndependentPriors
        gradients_inducing_points = zero(model.Z[1])
        A = ([I].-model.invKmm.*(model.Σ.+model.µ.*transpose.(model.μ))).*model.invKmm
        #preallocation
        ι = Matrix{T}(undef,model.nSamplesUsed,model.nInducingPoints)
        for k in 1:model.nPrior
            for i in 1:model.nInducingPoints #Iterate over the points
                Jnm,Jmm = computeIndPointsJ(model,i) #TODO
                for j in 1:model.nDim #iterate over the dimensions
                    @views mul!(ι,(Jnm[j,:,:]-model.κ[k]*Jmm[j,:,:]),model.invKmm[k])
                    @views gradients_inducing_points[c][i,j] =  hyperparameter_expec_gradient(model,ι,κΣ,Jmm[j,:,:],Jnm[j,:,:],zeros(T,model.nSamplesUsed))-hyperparameter_KL_gradient(Jmm[j,:,:],A[k])
                end
            end
        end
        return gradients_inducing_points
    else
        @warn "Inducing points for shared prior not implemented yet"
        gradients_inducing_points = zero(model.inducingPoints[1]) #TODO
    end
end

## Apply gradients on mean prior ##
function apply_gradients_mean_prior!(μ::AbstractVector{<:Real},grad_μ::AbstractVector{<:Real},opt::Optimizer)
    μ .+= update(opt,grad_μ)
end


text:
# https://github.com/JuliaDebug/Debugger.jl/blob/4cf99c662ab89da0fe7380c1e81461e2428e8b00/src/limitio.jl

mutable struct LimitIO{IO_t <: IO} <: IO
    io::IO_t
    maxbytes::Int
    n::Int
end
LimitIO(io::IO, maxbytes) = LimitIO(io, maxbytes, 0)

struct LimitIOException <: Exception end

function Base.write(io::LimitIO, v::UInt8)
    io.n > io.maxbytes && throw(LimitIOException())
    io.n += write(io.io, v)
end

function sprintlimited(args...; func = show, limit::Int = 30, ellipsis::AbstractString = "…", color = false)
    io = IOBuffer()
    ioctx = IOContext(LimitIO(io, limit - length(ellipsis)), :limit => true, :color => color, :displaysize => (30, 64))

    try
        Base.invokelatest(func, ioctx, args...)
    catch err
        if err isa LimitIOException
            print(io, ellipsis)
        else
            rethrow(err)
        end
    end

    str = filter(isvalid, String(take!(io)))

    return color ? str : remove_ansi_control_chars(str)
end

function strlimit(str; limit::Int = 30, ellipsis::AbstractString = "…")
    will_append = length(str) > limit

    io = IOBuffer()
    i = 1
    for c in str
        will_append && i > limit - length(ellipsis) && break
        isvalid(c) || continue

        print(io, c)
        i += 1
    end
    will_append && print(io, ellipsis)

    return String(take!(io))
end

# https://stackoverflow.com/a/33925425/12113178

function remove_ansi_control_chars(str::String)
    replace(str, r"(\x9B|\x1B\[)[0-?]*[ -\/]*[@-~]" => "")
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule GAP_pkg_datastructures_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("GAP_pkg_datastructures")
JLLWrappers.@generate_main_file("GAP_pkg_datastructures", UUID("66bcec33-4fa3-551a-b2da-ca256b6829d9"))
end  # module GAP_pkg_datastructures_jll


text:
module TestIterators

using Phylo
using Compat.Test
using Compat

@testset "$TreeType" for TreeType in [Nonultrametric, Ultrametric]
    ntip = 10
    nur = rand(TreeType(ntip))

    nni = nodenameiter(nur)
    @test Compat.IteratorEltype(nni) == Base.HasEltype()
    @test eltype(nni) == nodenametype(nur)
    @test Compat.IteratorSize(nni) == Base.HasLength()
    @test length(nni) == 2 * ntip - 1

    ni = nodeiter(nur)
    @test Compat.IteratorEltype(ni) == Base.HasEltype()
    @test eltype(ni) == nodetype(nur)
    @test Compat.IteratorSize(ni) == Base.HasLength()
    @test length(ni) == 2 * ntip - 1
    @test Set(ni) == Set(map(n -> getnode(nur, n), nni))

    bni = branchnameiter(nur)
    @test Compat.IteratorEltype(bni) == Base.HasEltype()
    @test eltype(bni) == branchnametype(nur)
    @test Compat.IteratorSize(bni) == Base.HasLength()
    @test length(bni) == ntip * 2 - 2

    bi = branchiter(nur)
    @test Compat.IteratorEltype(bi) == Base.HasEltype()
    @test eltype(bi) == branchtype(nur)
    @test Compat.IteratorSize(bi) == Base.HasLength()
    @test length(bi) == 2 * ntip - 2
    @test Set(bi) == Set(map(b -> getbranch(nur, b), bni))
end

@testset "$TreeType with filter" for TreeType in [Nonultrametric, Ultrametric]
    ntip = 10
    nur = rand(TreeType(ntip))

    nni = nodenamefilter(isleaf, nur)
    @test Compat.IteratorEltype(nni) == Base.HasEltype()
    @test eltype(nni) == nodenametype(nur)
    @test Compat.IteratorSize(nni) == Base.HasLength()
    @test length(nni) == ntip
    leaves = collect(nni)
    @test leaves == getleafnames(nur)

    ni = nodefilter(isleaf, nur)
    @test Compat.IteratorEltype(ni) == Base.HasEltype()
    @test eltype(ni) == nodetype(nur)
    @test Compat.IteratorSize(ni) == Base.HasLength()
    @test length(ni) == ntip
    @test all(isleaf, collect(ni))

    bni = branchnamefilter(nur) do branch
        return isleaf(nur, dst(branch))
    end
    @test Compat.IteratorEltype(bni) == Base.HasEltype()
    @test eltype(bni) == branchnametype(nur)
    @test Compat.IteratorSize(bni) == Base.HasLength()
    @test length(bni) == ntip

    bi = branchfilter(nur) do branch
        return isleaf(nur, dst(branch))
    end
    @test Compat.IteratorEltype(bi) == Base.HasEltype()
    @test eltype(bi) == branchtype(nur)
    @test Compat.IteratorSize(bi) == Base.HasLength()
    @test length(bi) == ntip
    bileaves = map(dst, bi)
    @test Set(bileaves) == Set(leaves)
end

end


text:
abstract type ObjectiveFunc end
abstract type Objective1D  <: ObjectiveFunc end
abstract type Objective2D  <: ObjectiveFunc end

Base.string(o::ObjectiveFunc) = o.name

#Test functions
mutable struct DistributionFunc <: Objective1D
    name::String
    d::Distribution
    σ::Float64
    xlim::Tuple{Float64,Float64}
    g_max::Float64
    g_min::Float64
    x_max::Float64

    function DistributionFunc(name::String, d::Distribution, σ::Float64; xmin=-1.0, xmax=1.0, n=1000)
        o = new(name, d, σ, (xmin,xmax))
        xs = linspace(xmin, xmax, n)
        o.g_max, i_max = findmax(truth.(o, xs))
        o.x_max = xs[i_max]
        o.g_min = minimum(truth.(o, xs))
        o
    end
end
function NarrowBump(d=Normal(-0.25,0.05); σ=0.3)
    DistributionFunc("NarrowBump", d, σ)
end
function WideBump(d=Normal(0.15,0.25); σ=0.3)
    DistributionFunc("WideBump", d, σ)
end
function TwoBumps(d=MixtureModel([Normal(-0.5,0.12), Normal(0.1,0.3)],
                                  [0.35,0.65]); σ=0.3)
    DistributionFunc("TwoBumps", d, σ)
end
function ThreeBumps(d=MixtureModel([Normal(-0.5,0.15),Normal(0.0,0.3),Normal(0.5,0.07)],
                                   [0.27,0.57,0.16]); σ=0.3)
    DistributionFunc("ThreeBumps", d, σ)
end
truth(o::DistributionFunc, x) = pdf(o.d, x)
function (o::DistributionFunc)(x, rng::AbstractRNG=Base.GLOBAL_RNG)
    truth(o,x) + o.σ*randn(rng)
end

mutable struct SincFunc <: Objective1D
    name::String
    a::Float64
    b::Float64
    c::Float64
    σ::Float64
    xlim::Tuple{Float64,Float64}
    g_max::Float64
    x_max::Float64
    g_min::Float64

    function SincFunc(name::String, a::Float64, b::Float64, c::Float64, σ::Float64; xmin=-1.0, xmax=1.0, n=1000)
        o = new(name, a, b, c, σ, (xmin,xmax))
        xs = linspace(xmin, xmax, n)
        o.g_max, i_max = findmax(truth.(o, xs))
        o.x_max = xs[i_max]
        o.g_min = minimum(truth.(o, xs))
        o
    end
end
function Sinc(; a=20.0, b=0.3, c=0.5, σ=0.3)
    SincFunc("Sinc", a,b,c,σ)
end
sinc_truth(a::Float64, b::Float64, c::Float64, x::Float64) = sinc(a * (x-b) + c)
truth(o::SincFunc, x) = sinc_truth(o.a, o.b, o.c, x)
function (o::SincFunc)(x, rng::AbstractRNG=Base.GLOBAL_RNG)
    truth(o,x) + o.σ*randn(rng)
end
function circle(x,y,r)
    c = Shape(Plots.partialcircle(0,2π,20,r))
    translate!(c, x, y)
    c
end
@recipe function plot(o::Objective1D; xmin=-1.0, xmax=1.0, n=100)
    xlim := o.xlim
    ylim := (o.g_min - 1.0, o.g_max + 1.0)
    xlabel := "x"
    ylabel := "f(x)"
    @series begin
        seriestype := :path
        label := "true mean"
        xs = linspace(xmin, xmax, n)
        xs, truth.(o, xs)
    end
end
@recipe function plot(o::Objective1D, rng::AbstractRNG; xmin=-1.0, xmax=1.0, n=100)
    xlim := o.xlim
    ylim := (o.g_min - 1.0, o.g_max + 1.0)
    xlabel := "x"
    ylabel := "f(x)"
    @series begin
        seriestype := :path
        label := "noisy sample"
        xs = linspace(xmin, xmax, n)
        xs, o.(xs, rng)
    end
end

Base.maximum(o::Objective1D) = o.g_max
Base.minimum(o::Objective1D) = o.g_min

#####
# 2D
mutable struct DistributionFunc2D <: Objective2D
    name::String
    d::Distribution
    σ::Float64
    xlim::Tuple{Float64,Float64}
    ylim::Tuple{Float64,Float64}
    g_max::Float64
    g_min::Float64
    x_max::Vector{Float64}

    function DistributionFunc2D(name::String, d::Distribution, σ::Float64;
                                xlim=(-2.0,2.0), ylim=(-2.0,2.0), n=1000)
        o = new(name, d, σ, xlim, ylim)
        Xs = [[x,y] for x in linspace(xlim..., n), y in linspace(ylim..., n)]
        o.g_max, i_max = findmax(truth.(o,Xs))
        o.x_max = Xs[i_max]
        o.g_min = minimum(truth.(o,Xs))
        o
    end
end
function NarrowBump2D(d=MvNormal([0.4,0.4],[0.2,0.2]); σ=2.0)
    DistributionFunc2D("NarrowBump2D", d, σ)
end
function WideBump2D(d=MvNormal([0.25,0.25],[0.2,0.7]); σ=0.5)
    DistributionFunc2D("WideBump2D", d, σ)
end
truth(o::DistributionFunc2D, x) = pdf(o.d, x)
function (o::DistributionFunc2D)(x, rng::AbstractRNG=Base.GLOBAL_RNG)
    truth(o,x) + o.σ*randn(rng)
end

mutable struct Rosenbrock <: Objective2D
    name::String
    σ::Float64
    xlim::Tuple{Float64,Float64}
    ylim::Tuple{Float64,Float64}
    g_max::Float64
    g_min::Float64
    x_max::Vector{Float64}

    function Rosenbrock(σ::Float64=10.0; xlim=(-2.0,2.0), ylim=(-2.0,2.0), n=1000)
        o = new("Rosenbrock", σ, xlim, ylim)
        Xs = [[x,y] for x in linspace(xlim..., n), y in linspace(ylim..., n)]
        o.g_max, i_max = findmax(truth.(o,Xs))
        o.x_max = Xs[i_max]
        o.g_min = minimum(truth.(o,Xs))
        o
    end
end
truth(o::Rosenbrock, x) = -rosenbrock(x)
function (o::Rosenbrock)(x, rng::AbstractRNG=Base.GLOBAL_RNG)
    truth(o,x) + o.σ*randn(rng)
end

@recipe function plot(o::Objective2D; n=100)
    xlim := o.xlim
    ylim := o.ylim
    zlim := (o.g_min - 1.0, o.g_max + 1.0)
    xlabel := "x"
    ylabel := "y"
    zlabel := "f(x,y)"
    @series begin
        seriestype := :heatmap
        label := "true mean"
        xs = linspace(o.xlim..., n)
        ys = linspace(o.ylim..., n)
        xs, ys, (x,y)->truth(o, [x,y])
    end
end
@recipe function plot(o::Objective2D, rng::AbstractRNG; n=100)
    xlim := o.xlim
    ylim := o.ylim
    zlim := (o.g_min - 1.0, o.g_max + 1.0)
    xlabel := "x"
    ylabel := "y"
    zlabel := "f(x,y)"
    @series begin
        seriestype := :heatmap
        label := "noisy sample"
        xs = linspace(o.xlim..., n)
        ys = linspace(o.ylim..., n)
        xs, ys, (x,y)->o([x,y],rng)
    end
end

Base.maximum(o::Objective2D) = o.g_max
Base.minimum(o::Objective2D) = o.g_min


text:
testOnAllContexts("Synthesizer") do ctx
    println("\tconcept: select indices by chance")
    const depth = 2
    const count = 4
    s = range(ctx, UInt32, [count])
    values = moddims(range(ctx, Float32, [count * depth]), count, depth)
    dist = randu(ctx, Float32, count, depth)
    rnd = randu(ctx, Float32, count, depth)
    test = rnd ./ maxOf(dist, 0.000001f0)
    max, idx = imax(test, 1)
    idx2 = idx .* count .+ s
    result = values[idx2]
    @test dims(result, 0) == count
    @test dims(result, 1) == 1
    @test host(result) == [0.0f0, 1.0f0, 6.0f0, 3.0f0]

    println("\tSynthesizer + scalar and set rule")

    syn = Synthesizer(ctx)
    genes = Vector{Float32}()
    firstRule = rule = SetRule(["A", "B", "C"], 3)
    define!(syn, rule)
    @test resultSize(rule) == 3
    @test dnaSize(rule) == 9

    push!(genes, 0.1f0)
    push!(genes, 0.1f0)
    push!(genes, 0.1f0)
    push!(genes, 0.1f0)
    push!(genes, 0.1f0)
    push!(genes, 0.1f0)
    push!(genes, 0.9f0)
    push!(genes, 0.9f0)
    push!(genes, 0.9f0)

    rule = ScalarRule(count = 3)
    define!(syn, rule)
    @test resultSize(rule) == 3
    @test dnaSize(rule) == 3

    push!(genes, 0.5f0)
    push!(genes, 0.6f0)
    push!(genes, 0.7f0)

    rule = ScalarRule(count = 2, round = true, min = 1f0, max = 100f0)
    define!(syn, rule)
    @test resultSize(rule) == 2
    @test dnaSize(rule) == 2

    push!(genes, 0.5f0)
    push!(genes, 0.6f0)

    rule = ScalarRule(count = 2, round = true, min = 1, max = 100, variationCount = 2)
    define!(syn, rule)
    @test resultSize(rule) == 2
    @test dnaSize(rule) == 2 * 2 * 2

    push!(genes, 0.001f0)
    push!(genes, 0.01f0)
    push!(genes, 0.99f0)
    push!(genes, 0.59f0)

    push!(genes, 0.1f0)
    push!(genes, 0.19f0)
    push!(genes, 0.12f0)
    push!(genes, 0.2f0)

    dna = DNA(ctx, genes)
    result = decode(syn, dna)

    @test typeof(result) == Vector{Any}

    @test dnaSize!(syn) == length(dna)
    @test resultSize!(syn) == length(result)

    i = 0

    @test result[i += 1] in ["A", "B", "C"]
    @test result[i += 1] in ["A", "B", "C"]
    @test result[i += 1] in ["A", "B", "C"]

    @test result[i += 1] == 0.5f0
    @test result[i += 1] == 0.6f0
    @test result[i += 1] == 0.7f0

    @test result[i += 1] == 51f0
    @test result[i += 1] == 60f0

    @test result[i += 1] in [1, 2, 59, 99]
    @test result[i += 1] in [1, 2, 59, 99]

    println("\tdecode at")
    caResult = decodeAsContextArray(syn, dna)

    @test typeof(caResult) != Vector{Any}

    @test resultSize!(syn) == dims(caResult, 0)
    @test length(size(caResult)) == 1
    @test resultSize!(syn) ==  size(caResult)[1]

    decoded = decodeAt(syn, caResult, [DecodeRule(firstRule), DecodeRule(2, true), DecodeRule(syn.rules[3], false)])
    @test typeof(decoded) == Vector{Any}

    result1 = decoded[1]
    result2 = decoded[2]
    result3 = decoded[3]

    @test isa(result1, Vector{Any})
    @test length(result1) == 3
    for item in result1
        @test item in ["A", "B", "C"]
    end

    result2H = host(result2)
    @test length(result2H) == 3
    @test isa(result2H, Vector{Float32})
    @test result2H[1] == 0.5f0
    @test result2H[2] == 0.6f0
    @test result2H[3] == 0.7f0

    @test isa(result3, Vector{Any})
    @test length(result3) == 2
    @test result3[1] == 51f0
    @test result3[2] == 60f0

    @test_throws BoundsError decodeAt(syn, caResult, [DecodeRule(42)])

    @test_throws ErrorException decodeAt(syn, caResult, [DecodeRule(ScalarRule(count = 3))])
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
module Utils
    using Dates
#---------------------Parse utils-----------------------------------------------------#
"""
Searches a directory for all files containing the key.
"""
searchdir(path::String, key) = joinpath.((path,), filter(x -> occursin(key, x), readdir(path)))
export searchdir

"""
Parse an array of strings into an array of a type.
"""
function parse_string_array(T::Type, array)
    return map(x -> (v = tryparse(T, x); v == nothing ? 0.0 : v), array)
end
export parse_string_array

"""
Parse a line for occurrences of type T.
"""
parse_line(T::Type, line::String) = parse_string_array(T, split(line))
export parse_line

"""
Splits a line using arguments, then strips spaces from the splits.
"""
strip_split(line, args...) = strip.(split(line, args...))
export strip_split

function replace_multiple(str, replacements::Pair{String,String}...)
    tstr = deepcopy(str)
    for r in replacements
        tstr = replace(tstr, r)
    end
    return tstr
end
export replace_multiple

function cut_after(line, c)
    t = findfirst(isequal(c), line)
    if t == nothing
        return line
    elseif t == 1
        return ""
    else
        return line[1:t-1]
    end
end
export cut_after

#--------------------------------------------------------------------------------------#
function fort2julia(f_type)
    f_type = lowercase(f_type)
    if f_type == "real"
        return Float32
    elseif f_type == "real(kind=dp)"
        return Float64
    elseif f_type == "complex(kind=dp)"
        return Complex{Float64}
    elseif occursin("character", f_type)
        return String
    elseif f_type == "string"
        return String
    elseif f_type == "integer"
        return Int
    elseif f_type == "logical"
        return Bool
    elseif occursin(".D", f_type)
        return replace(f_type, "D" => "e")
    else
        return Nothing
    end
end
export fort2julia

"""
It's like filter()[1].
"""
function getfirst(f::Function, A)
    for el in A
        if f(el)
            return el
        end
    end
end
export getfirst

"""
    parse_block(f, types...; to_strip=',')

Takes the specified types and parses each line into the types.
When it finds a line where it cannot match all the types, it stops and returns  the parsed values.
The split and strip keywords let the user specify how to first split the line, then strip the splits from the strip char.
"""
function parse_block(f, types...; to_strip = ',')
    output = []
    len_typ = length(types)
    while !eof(f)
        line = strip.(split(readline(f)), to_strip)
        len_lin = length(line)
        if isempty(line)
            continue
        end
        i, j = 1, 1
        tmp = []
        while i <= len_typ && j <= len_lin
            typ = types[i]
            l   = line[j]
            try
                t = Meta.parse(l)
                if typeof(t) == typ
                    push!(tmp, t)
                    i += 1
                    j += 1
                else
                    j += 1
                end
            catch
                j += 1
            end
        end
        if length(tmp) < length(types)
            return output
        end
        push!(output, Tuple{types...}(tmp))
    end
    return output
end
export parse_block

yesterday() = today() - Day(1)
lastweek()  = today() - Week(1)
lastmonth() = today() - Month(1)

username() = read(`whoami`, String)
export yesterday, lastweek, lastmonth, username

function loaded_modules_string()
    all = string.(values(Base.loaded_modules))
    valid = String[]
    ks = keys(Pkg.project().dependencies)
    for p in all
        if p ∈ ks
            push!(valid, p)
        end
    end
    return """using $(join(valid, ", "))\n"""
end
export loaded_modules_string

end