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

# Convert # of Rata Die days to proleptic Gregorian calendar y,m,d,w
# Reference: http://mysite.verizon.net/aesir_research/date/date0.htm
function yearmonthday(days)
    z = days + 306; h = 100z - 25; a = fld(h, 3652425); b = a - fld(a, 4)
    y = fld(100b + h, 36525); c = b + z - 365y - fld(y, 4); m = div(5c + 456, 153)
    d = c - div(153m - 457, 5); return m > 12 ? (y + 1, m - 12, d) : (y, m, d)
end
function year(days)
   z = days + 306; h = 100z - 25; a = fld(h, 3652425); b = a - fld(a, 4)
   y = fld(100b + h, 36525); c = b + z - 365y - fld(y, 4); m = div(5c + 456, 153)
   return m > 12 ? y + 1 : y
end
function yearmonth(days)
    z = days + 306; h = 100z - 25; a = fld(h,3652425); b = a - fld(a,4)
    y = fld(100b + h, 36525); c = b + z - 365y - fld(y, 4); m = div(5c + 456, 153)
    return m > 12 ? (y + 1, m - 12) : (y, m)
end
function month(days)
    z = days + 306; h = 100z - 25; a = fld(h,3652425); b = a - fld(a,4)
    y = fld(100b + h, 36525); c = b + z - 365y - fld(y, 4); m = div(5c + 456, 153)
    return m > 12 ? m - 12 : m
end
function monthday(days)
    z = days + 306; h = 100z - 25; a = fld(h,3652425); b = a - fld(a,4)
    y = fld(100b + h, 36525); c = b + z - 365y - fld(y, 4); m = div(5c + 456, 153)
    d = c - div(153m - 457, 5); return m > 12 ? (m - 12, d) : (m, d)
end
function day(days)
    z = days + 306; h = 100z - 25; a = fld(h,3652425); b = a - fld(a,4)
    y = fld(100b + h, 36525); c = b + z - 365y - fld(y, 4); m = div(5c + 456, 153)
    return c - div(153m - 457, 5)
end
# https://en.wikipedia.org/wiki/Talk:ISO_week_date#Algorithms
const WEEK_INDEX = (15, 23, 3, 11)
function week(days)
    w = div(abs(days - 1), 7) % 20871
    c, w = divrem((w + (w >= 10435)), 5218)
    w = (w * 28 + WEEK_INDEX[c + 1]) % 1461
    return div(w, 28) + 1
end

# Accessor functions
value(dt::TimeType) = dt.instant.periods.value
value(t::Time) = t.instant.value
days(dt::Date) = value(dt)
days(dt::DateTime) = fld(value(dt), 86400000)
year(dt::TimeType) = year(days(dt))
month(dt::TimeType) = month(days(dt))
week(dt::TimeType) = week(days(dt))
day(dt::TimeType) = day(days(dt))
hour(dt::DateTime)   = mod(fld(value(dt), 3600000), 24)
minute(dt::DateTime) = mod(fld(value(dt), 60000), 60)
second(dt::DateTime) = mod(fld(value(dt), 1000), 60)
millisecond(dt::DateTime) = mod(value(dt), 1000)
hour(t::Time)   = mod(fld(value(t), 3600000000000), Int64(24))
minute(t::Time) = mod(fld(value(t), 60000000000), Int64(60))
second(t::Time) = mod(fld(value(t), 1000000000), Int64(60))
millisecond(t::Time) = mod(fld(value(t), Int64(1000000)), Int64(1000))
microsecond(t::Time) = mod(fld(value(t), Int64(1000)), Int64(1000))
nanosecond(t::Time) = mod(value(t), Int64(1000))

dayofmonth(dt::TimeType) = day(dt)

yearmonth(dt::TimeType) = yearmonth(days(dt))
monthday(dt::TimeType) = monthday(days(dt))
yearmonthday(dt::TimeType) = yearmonthday(days(dt))

# Documentation for exported accessors
for func in (:year, :month)
    name = string(func)
    @eval begin
        @doc """
            $($name)(dt::TimeType) -> Int64

        The $($name) of a `Date` or `DateTime` as an `Int64`.
        """ $func(dt::TimeType)
    end
end

"""
    week(dt::TimeType) -> Int64

Return the [ISO week date](https://en.wikipedia.org/wiki/ISO_week_date) of a `Date` or
`DateTime` as an `Int64`. Note that the first week of a year is the week that contains the
first Thursday of the year which can result in dates prior to January 4th being in the last
week of the previous year. For example `week(Date(2005,1,1))` is the 53rd week of 2004.
"""
week(dt::TimeType)

for func in (:day, :dayofmonth)
    name = string(func)
    @eval begin
        @doc """
            $($name)(dt::TimeType) -> Int64

        The day of month of a `Date` or `DateTime` as an `Int64`.
        """ $func(dt::TimeType)
    end
end

"""
    hour(dt::DateTime) -> Int64

The hour of day of a `DateTime` as an `Int64`.
"""
hour(dt::DateTime)

for func in (:minute, :second, :millisecond)
    name = string(func)
    @eval begin
        @doc """
            $($name)(dt::DateTime) -> Int64

        The $($name) of a `DateTime` as an `Int64`.
        """ $func(dt::DateTime)
    end
end

for parts in (["year", "month"], ["month", "day"], ["year", "month", "day"])
    name = join(parts)
    func = Symbol(name)
    @eval begin
        @doc """
            $($name)(dt::TimeType) -> ($(join(repeated(Int64, length($parts)), ", ")))

        Simultaneously return the $(join($parts, ", ", " and ")) parts of a `Date` or
        `DateTime`.
        """ $func(dt::TimeType)
    end
end

for func in (:hour, :minute, :second, :millisecond, :microsecond, :nanosecond)
    name = string(func)
    @eval begin
        @doc """
            $($name)(t::Time) -> Int64

        The $($name) of a `Time` as an `Int64`.
        """ $func(t::Time)
    end
end


text:
module VectorMathOptInterface

import MathOptInterface
const MOI = MathOptInterface

include("indextypes.jl")
include("attributes.jl")
include("objectives.jl")

########################## MOI "Import & Export" #################################
# get the names defined in MOI and construct and evaluate an expression `_moi_expr`
# such that for every `_moi_name` (not in `IMPORT_FILTER`) we have an assignement
# `_moi_name = VectorMathOptInterface.MOI._moi_name`.

IMPORT_FILTER = [
	:ObjectiveFunction,
	:ObjectiveSense,
	:ObjectiveFunctionType,
	:ObjectiveValue,
	:DualObjectiveValue,
	:ObjectiveBound ]	# TODO RelativeGap?

_moi_names = Symbol[]
for _moi_name in Base.names(MOI; all = true, imported = false)
	_moi_name_str = string(_moi_name)
	if startswith(_moi_name_str, "#")
		_moi_name_str = _moi_name_str[2:end]
	end
	_moi_name = Symbol(_moi_name_str)

	if !(_moi_name in IMPORT_FILTER) && !occursin("#",_moi_name_str) &&
		!isdefined(Base, _moi_name)
		push!(_moi_names, _moi_name)
	end
end

_vmoi_moi_subexpr = [
	Expr(:(=),
		_moi_name,
		Expr(:(.),
			Expr(:(.),
				:VectorMathOptInterface, QuoteNode(:MOI)),
			QuoteNode(_moi_name))
	) for _moi_name in _moi_names
]
_moi_expr = Expr(:toplevel, _vmoi_moi_subexpr...)
eval(_moi_expr)
########################## MOI "Import & Export" #################################

include("utils/ExampleMOP.jl")
using .ExampleMOP

end

ext: jl


text:
# This file is based on LinearAlgebra/test/special.jl a part of Julia. License is MIT: https://julialang.org/license

using Test, SparseArrays, Random, LazyBandedMatrices
import LazyBandedMatrices: Bidiagonal, Tridiagonal, SymTridiagonal
import LinearAlgebra
import LinearAlgebra: UniformScaling

@testset "Tri/Bidiagonal special" begin
    n= 10 #Size of matrix to test
    Random.seed!(1)

    @testset "Interconversion between special matrix types" begin
        a = [1.0:n;]
        A = Diagonal(a)
        @testset for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, Matrix]
            @test Matrix(convert(newtype, A)) == Matrix(A)
            @test Matrix(convert(newtype, Diagonal(GenericArray(a)))) == Matrix(A)
        end

        @testset for isupper in (true, false)
            A = Bidiagonal(a, [1.0:n-1;], ifelse(isupper, :U, :L))
            for newtype in [Bidiagonal, Tridiagonal, Matrix]
            @test Matrix(convert(newtype, A)) == Matrix(A)
            @test Matrix(newtype(A)) == Matrix(A)
            end
            @test_throws ArgumentError convert(SymTridiagonal, A)
            tritype = isupper ? UpperTriangular : LowerTriangular
            @test Matrix(tritype(A)) == Matrix(A)

            A = Bidiagonal(a, zeros(n-1), ifelse(isupper, :U, :L)) #morally Diagonal
            for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, Matrix]
            @test Matrix(convert(newtype, A)) == Matrix(A)
            @test Matrix(newtype(A)) == Matrix(A)
            end
            @test Matrix(tritype(A)) == Matrix(A)
        end

        A = SymTridiagonal(a, [1.0:n-1;])
        for newtype in [Tridiagonal, Matrix]
            @test Matrix(convert(newtype, A)) == Matrix(A)
        end
        for newtype in [Diagonal, Bidiagonal]
            @test_throws ArgumentError convert(newtype,A)
        end
        A = SymTridiagonal(a, zeros(n-1))
        @test Matrix(convert(Bidiagonal,A)) == Matrix(A)

        A = Tridiagonal(zeros(n-1), [1.0:n;], zeros(n-1)) #morally Diagonal
        for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Matrix]
        @test Matrix(convert(newtype, A)) == Matrix(A)
        end
        A = Tridiagonal(fill(1., n-1), [1.0:n;], fill(1., n-1)) #not morally Diagonal
        for newtype in [SymTridiagonal, Matrix]
        @test Matrix(convert(newtype, A)) == Matrix(A)
        end
        for newtype in [Diagonal, Bidiagonal]
            @test_throws ArgumentError convert(newtype,A)
        end
        A = Tridiagonal(zeros(n-1), [1.0:n;], fill(1., n-1)) #not morally Diagonal
        @test Matrix(convert(Bidiagonal, A)) == Matrix(A)
        A = UpperTriangular(Tridiagonal(zeros(n-1), [1.0:n;], fill(1., n-1)))
        @test Matrix(convert(Bidiagonal, A)) == Matrix(A)
        A = Tridiagonal(fill(1., n-1), [1.0:n;], zeros(n-1)) #not morally Diagonal
        @test Matrix(convert(Bidiagonal, A)) == Matrix(A)
        A = LowerTriangular(Tridiagonal(fill(1., n-1), [1.0:n;], zeros(n-1)))
        @test Matrix(convert(Bidiagonal, A)) == Matrix(A)
        @test_throws ArgumentError convert(SymTridiagonal,A)

        A = LowerTriangular(Matrix(Diagonal(a))) #morally Diagonal
        for newtype in [Diagonal, Bidiagonal, SymTridiagonal, LowerTriangular, Matrix]
            @test Matrix(convert(newtype, A)) == Matrix(A)
        end
        A = UpperTriangular(Matrix(Diagonal(a))) #morally Diagonal
        for newtype in [Diagonal, Bidiagonal, SymTridiagonal, UpperTriangular, Matrix]
            @test Matrix(convert(newtype, A)) == Matrix(A)
        end
        A = UpperTriangular(triu(rand(n,n)))
        for newtype in [Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal]
            @test_throws ArgumentError convert(newtype,A)
        end

        # test operations/constructors (not conversions) permitted in the docs
        dl = [1., 1.]
        d = [-2., -2., -2.]
        T = Tridiagonal(dl, d, -dl)
        S = SymTridiagonal(d, dl)
        Bu = Bidiagonal(d, dl, :U)
        Bl = Bidiagonal(d, dl, :L)
        D = Diagonal(d)
        M = [-2. 0. 0.; 1. -2. 0.; -1. 1. -2.]
        U = UpperTriangular(M)
        L = LowerTriangular(Matrix(M'))

        for A in (T, S, Bu, Bl, D, U, L, M)
            Adense = Matrix(A)
            B = Symmetric(A)
            Bdense = Matrix(B)
            for (C,Cdense) in ((A,Adense), (B,Bdense))
                @test Diagonal(C) == Diagonal(Cdense)
                @test Bidiagonal(C, :U) == Bidiagonal(Cdense, :U)
                @test Bidiagonal(C, :L) == Bidiagonal(Cdense, :L)
                @test Tridiagonal(C) == Tridiagonal(Cdense)
                @test UpperTriangular(C) == UpperTriangular(Cdense)
                @test LowerTriangular(C) == LowerTriangular(Cdense)
            end
        end
    end

    @testset "Binary ops among special types" begin
        a=[1.0:n;]
        A=Diagonal(a)
        Spectypes = [Diagonal, Bidiagonal, Tridiagonal, Matrix]
        for (idx, type1) in enumerate(Spectypes)
            for type2 in Spectypes
            B = convert(type1,A)
            C = convert(type2,A)
            @test Matrix(B + C) ≈ Matrix(A + A)
            @test Matrix(B - C) ≈ Matrix(A - A)
        end
        end
        B = SymTridiagonal(a, fill(1., n-1))
        for Spectype in [Diagonal, Bidiagonal, Tridiagonal, Matrix]
            @test Matrix(B + convert(Spectype,A)) ≈ Matrix(B + A)
            @test Matrix(convert(Spectype,A) + B) ≈ Matrix(B + A)
            @test Matrix(B - convert(Spectype,A)) ≈ Matrix(B - A)
            @test Matrix(convert(Spectype,A) - B) ≈ Matrix(A - B)
        end

        C = rand(n,n)
        for TriType in [LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, UpperTriangular, LowerTriangular]
            D = TriType(C)
            for Spectype in [Diagonal, Bidiagonal, Tridiagonal, Matrix]
                @test Matrix(D + convert(Spectype,A)) ≈ Matrix(D + A)
                @test Matrix(convert(Spectype,A) + D) ≈ Matrix(A + D)
                @test Matrix(D - convert(Spectype,A)) ≈ Matrix(D - A)
                @test Matrix(convert(Spectype,A) - D) ≈ Matrix(A - D)
            end
        end

        UpTri = UpperTriangular(rand(20,20))
        LoTri = LowerTriangular(rand(20,20))
        Diag = Diagonal(rand(20,20))
        Tridiag = Tridiagonal(rand(20, 20))
        UpBi = Bidiagonal(rand(20,20), :U)
        LoBi = Bidiagonal(rand(20,20), :L)
        Sym = SymTridiagonal(rand(20), rand(19))
        Dense = rand(20, 20)
        mats = [UpTri, LoTri, Diag, Tridiag, UpBi, LoBi, Sym, Dense]

        for op in (+,-,*)
            for A in mats
                for B in mats
                    @test (op)(A, B) ≈ (op)(Matrix(A), Matrix(B)) ≈ Matrix((op)(A, B))
                end
            end
        end
    end

    @testset "+ and - among structured matrices with different container types" begin
        diag = 1:5
        offdiag = 1:4
        uniformscalingmats = [UniformScaling(3), UniformScaling(1.0), UniformScaling(3//5), UniformScaling(ComplexF64(1.3, 3.5))]
        mats = [Diagonal(diag), Bidiagonal(diag, offdiag, 'U'), Bidiagonal(diag, offdiag, 'L'), Tridiagonal(offdiag, diag, offdiag), SymTridiagonal(diag, offdiag)]
        for T in [ComplexF64, Int64, Rational{Int64}, Float64]
            push!(mats, Diagonal(Vector{T}(diag)))
            push!(mats, Bidiagonal(Vector{T}(diag), Vector{T}(offdiag), 'U'))
            push!(mats, Bidiagonal(Vector{T}(diag), Vector{T}(offdiag), 'L'))
            push!(mats, Tridiagonal(Vector{T}(offdiag), Vector{T}(diag), Vector{T}(offdiag)))
            push!(mats, SymTridiagonal(Vector{T}(diag), Vector{T}(offdiag)))
        end

        for op in (+,*) # to do: fix when operation is - and the matrix has a range as the underlying representation and we get a step size of 0.
            for A in mats
                for B in mats
                    @test (op)(A, B) ≈ (op)(Matrix(A), Matrix(B)) ≈ Matrix((op)(A, B))
                end
            end
        end
        for op in (+,-)
            for A in mats
                for B in uniformscalingmats
                    @test (op)(A, B) ≈ (op)(Matrix(A), B) ≈ Matrix((op)(A, B))
                    @test (op)(B, A) ≈ (op)(B, Matrix(A)) ≈ Matrix((op)(B, A))
                end
            end
        end
    end

    @testset "zero and one for structured matrices" begin
        for elty in (Int64, Float64, ComplexF64)
            D = Diagonal(rand(elty, 10))
            Bu = Bidiagonal(rand(elty, 10), rand(elty, 9), 'U')
            Bl = Bidiagonal(rand(elty, 10), rand(elty, 9), 'L')
            T = Tridiagonal(rand(elty, 9),rand(elty, 10), rand(elty, 9))
            S = SymTridiagonal(rand(elty, 10), rand(elty, 9))
            mats = [D, Bu, Bl, T, S]
            for A in mats
                @test iszero(zero(A))
                @test isone(one(A))
                @test zero(A) == zero(Matrix(A))
                @test one(A) == one(Matrix(A))
            end

            @test zero(D) isa Diagonal
            @test one(D) isa Diagonal

            @test zero(Bu) isa Bidiagonal
            @test one(Bu) isa Bidiagonal
            @test zero(Bl) isa Bidiagonal
            @test one(Bl) isa Bidiagonal
            @test zero(Bu).uplo == one(Bu).uplo == Bu.uplo
            @test zero(Bl).uplo == one(Bl).uplo == Bl.uplo

            @test zero(T) isa Tridiagonal
            @test one(T) isa Tridiagonal
            @test zero(S) isa SymTridiagonal
            @test one(S) isa SymTridiagonal
        end

        # ranges
        D = Diagonal(1:10)
        Bu = Bidiagonal(1:10, 1:9, 'U')
        Bl = Bidiagonal(1:10, 1:9, 'L')
        T = Tridiagonal(1:9, 1:10, 1:9)
        S = SymTridiagonal(1:10, 1:9)
        mats = [D, Bu, Bl, T, S]
        for A in mats
            @test iszero(zero(A))
            @test isone(one(A))
            @test zero(A) == zero(Matrix(A))
            @test one(A) == one(Matrix(A))
        end

        @test zero(D) isa Diagonal
        @test one(D) isa Diagonal

        @test zero(Bu) isa Bidiagonal
        @test one(Bu) isa Bidiagonal
        @test zero(Bl) isa Bidiagonal
        @test one(Bl) isa Bidiagonal
        @test zero(Bu).uplo == one(Bu).uplo == Bu.uplo
        @test zero(Bl).uplo == one(Bl).uplo == Bl.uplo

        @test zero(T) isa Tridiagonal
        @test one(T) isa Tridiagonal
        @test zero(S) isa SymTridiagonal
        @test one(S) isa SymTridiagonal
    end

    @testset "== for structured matrices" begin
        diag = rand(10)
        offdiag = rand(9)
        D = Diagonal(rand(10))
        Bup = Bidiagonal(diag, offdiag, 'U')
        Blo = Bidiagonal(diag, offdiag, 'L')
        Bupd = Bidiagonal(diag, zeros(9), 'U')
        Blod = Bidiagonal(diag, zeros(9), 'L')
        T = Tridiagonal(offdiag, diag, offdiag)
        Td = Tridiagonal(zeros(9), diag, zeros(9))
        Tu = Tridiagonal(zeros(9), diag, offdiag)
        Tl = Tridiagonal(offdiag, diag, zeros(9))
        S = SymTridiagonal(diag, offdiag)
        Sd = SymTridiagonal(diag, zeros(9))

        mats = [D, Bup, Blo, Bupd, Blod, T, Td, Tu, Tl, S, Sd]

        for a in mats
            for b in mats
                @test (a == b) == (Matrix(a) == Matrix(b)) == (b == a) == (Matrix(b) == Matrix(a))
            end
        end
    end

end # module TestSpecial


text:
@testset "module validation" begin
    using Atom: toplevelitems

    path = subjunkspath
    text = read(path, String)

    # basic -- finds every toplevel items with default arguments
    @test filter(toplevelitems(text)) do item
        item isa Atom.ToplevelBinding &&
        item.bind.name == "imwithdoc"
    end |> length === 3

    # don't enter non-target modules, e.g.: submodules
    @test filter(toplevelitems(text; mod = "Junk", inmod = true)) do item
        item isa Atom.ToplevelBinding &&
        item.bind.name == "imwithdoc"
    end |> length === 1 # should only find the `imwithdoc` in Junk module

    # don't include items outside of a module
    @test filter(toplevelitems(text; mod = "SubJunk", inmod = false)) do item
        item isa Atom.ToplevelBinding &&
        item.bind.name == "imwithdoc"
    end |> length === 1 # should only find the `imwithdoc` in SubJunk module
end


text:
# TODO, NOT UPDATED
using DifferentiableStateSpaceModels, BenchmarkTools, Zygote, Random

const RBC_LIKELIHOODS = BenchmarkGroup()
const RBC_LIKELIHOODS["kalman"] = BenchmarkGroup()
const RBC_LIKELIHOODS["joint"] = BenchmarkGroup()
const RBC_LIKELIHOODS["second_order_joint"] = BenchmarkGroup()
const RBC_LIKELIHOODS["perturbation"] = BenchmarkGroup()
const RBC_LIKELIHOODS["solve"] = BenchmarkGroup()

const m_rbc = @include_example_module(Examples.rbc_observables_benchmark)
const m_second_rbc = @include_example_module(Examples.rbc_observables_benchmark, 2)

function generate_benchmark_data(m::FirstOrderPerturbationModel, p, p_f, T, x0,
                                 ϵ = [randn(m.n_ϵ) for _ in 1:T])
    return (solve(generate_perturbation(m, p; p_f), x0, (0, T), LTI(); noise = ϵ).z[2:end],
            ϵ)
end

function generate_benchmark_data(m::SecondOrderPerturbationModel, p, p_f, T, x0,
                                 ϵ = [randn(m.n_ϵ) for _ in 1:T])
    return (solve(generate_perturbation(m, p; p_f), x0, (0, T), QTI(); noise = ϵ).z[2:end],
            ϵ)
end
const p_rbc = [0.4, 0.96]
const p_f_rbc = [0.7, 0.1, 0.01, 0.0001]
const T = 100
const x0_rbc = zeros(m_rbc.n_x)
const z_rbc, ϵ_rbc = generate_benchmark_data(m_rbc, p_rbc, p_f_rbc, T, x0_rbc)
const x0_second_rbc = zeros(m_rbc.n_x)
const z_second_rbc, ϵ_second_rbc = generate_benchmark_data(m_second_rbc, p_rbc, p_f_rbc, T,
                                                           x0_second_rbc)
const cache_rbc = allocate_cache(m_rbc)
const cache_second_rbc = allocate_cache(m_second_rbc)
const settings_skip_hash = PerturbationSolverSettings(; use_solution_cache = false)  #reuses cache storage but doesn't check cache

# conditional likelihood
function benchmark_likelihood_kalman(p, p_f, z, m, cache)
    sol = generate_perturbation(m, p; p_f, cache)
    return solve(sol, sol.x_ergodic, (0, length(z)); observables = z).logpdf
end

# Just the perturbation DifferentiableStateSpaceModels
# WE do not expect the reuse of the memory in the cache to save much time for small models like this
RBC_LIKELIHOODS["perturbation"]["first_order_perturbation"] = @benchmarkable generate_perturbation($m_rbc,
                                                                                $p_rbc;
                                                                                p_f = $p_f_rbc,
                                                                                cache) setup = (cache = allocate_cache(m_rbc))
RBC_LIKELIHOODS["perturbation"]["first_order_perturbation_reused_cache"] = @benchmarkable generate_perturbation($m_rbc,
                                                                                $p_rbc;
                                                                                p_f = $p_f_rbc,
                                                                                cache = $cache_rbc)  # shouldn't do many calculations  since we didn't change the p or p_f
RBC_LIKELIHOODS["perturbation"]["first_order_perturbation_reused_cache_no_hash"] = @benchmarkable generate_perturbation($m_rbc,
                                                                                $p_rbc;
                                                                                p_f = $p_f_rbc,
                                                                                cache = $cache_rbc,
                                                                                settings = $settings_skip_hash)  #reuse memory but no skipping perturbation
RBC_LIKELIHOODS["perturbation"]["second_order_perturbation"] = @benchmarkable generate_perturbation($m_second_rbc,
                                                                                $p_rbc;
                                                                                p_f = $p_f_rbc,
                                                                                cache) setup = (cache = allocate_cache(m_second_rbc))
RBC_LIKELIHOODS["perturbation"]["second_order_perturbation_reused_cache"] = @benchmarkable generate_perturbation($m_second_rbc,
                                                                                $p_rbc;
                                                                                p_f = $p_f_rbc,
                                                                                cache = $cache_second_rbc)  # shouldn't do any significant calculations  since we didn't change the p or p_f
RBC_LIKELIHOODS["perturbation"]["second_order_perturbation_reused_cache_no_hash"] = @benchmarkable generate_perturbation($m_second_rbc,
                                                                                $p_rbc;
                                                                                p_f = $p_f_rbc,
                                                                                cache = $cache_second_rbc,
                                                                                settings = $settings_skip_hash)  #reuse memory but no skipping perturbation

# Testing the solve isolated sequential
const sol_rbc = generate_perturbation(m_rbc, p_rbc; p_f = p_f_rbc)
const sol_second_rbc = generate_perturbation(m_second_rbc, p_rbc; p_f = p_f_rbc)
const tspan = (0, T)
RBC_LIKELIHOODS["solve"]["kalman"] = @benchmarkable solve($sol_rbc, $sol_rbc.x_ergodic,
                                                          $tspan;
                                                          observables = $z_rbc).logpdf
# RBC_LIKELIHOODS["solve"]["kalman_gradient"] = @benchmarkable gradient((A, B) -> solve(A, B, $sol_rbc.x_ergodic, $tspan; h = $sol_rbc.C, D = $sol_rbc.D, observables = $z_rbc).logpdf, $sol_rbc.A, $sol_rbc.B)

# for testing the solve on its own
test_evolution(x, p, t) = p.A * x
test_volatility(x, p, t) = p.B
test_observation(x, p, t) = p.C * x
const p_test = (; sol_rbc.A, sol_rbc.B, sol_rbc.C)
RBC_LIKELIHOODS["solve"]["first_order_joint"] = @benchmarkable solve(test_evolution,
                                                                     test_volatility,
                                                                     $x0_rbc, $tspan,
                                                                     $p_test;
                                                                     h = test_observation,
                                                                     D = $sol_rbc.D,
                                                                     observables = $z_rbc,
                                                                     noise = $ϵ_rbc).logpdf
RBC_LIKELIHOODS["solve"]["first_order_joint_gradient"] = @benchmarkable gradient(p_test -> solve(test_evolution,
                                                                                test_volatility,
                                                                                $x0_rbc,
                                                                                $tspan,
                                                                                p_test;
                                                                                h = test_observation,
                                                                                D = $sol_rbc.D,
                                                                                observables = $z_rbc,
                                                                                noise = $ϵ_rbc).logpdf,
                                                                                $p_test)
RBC_LIKELIHOODS["solve"]["first_order_joint_gradient_ϵ"] = @benchmarkable gradient((p_test, ϵ) -> solve(test_evolution,
                                                                                test_volatility,
                                                                                $x0_rbc,
                                                                                $tspan,
                                                                                p_test;
                                                                                h = test_observation,
                                                                                D = $sol_rbc.D,
                                                                                observables = $z_rbc,
                                                                                noise = ϵ).logpdf,
                                                                                $p_test,
                                                                                $ϵ_rbc)
RBC_LIKELIHOODS["solve"]["second_order_joint"] = @benchmarkable solve($sol_second_rbc,
                                                                      $x0_second_rbc,
                                                                      $tspan;
                                                                      observables = $z_second_rbc,
                                                                      noise = $ϵ_second_rbc).logpdf
RBC_LIKELIHOODS["solve"]["second_order_joint_gradient_ϵ"] = @benchmarkable gradient(ϵ -> solve($sol_second_rbc,
                                                                                $x0_second_rbc,
                                                                                $tspan;
                                                                                observables = $z_second_rbc,
                                                                                noise = ϵ).logpdf,
                                                                                $ϵ_second_rbc)

# likelihood and gradient with cache reallocated each time
RBC_LIKELIHOODS["kalman"]["likelihood"] = @benchmarkable benchmark_likelihood_kalman($p_rbc,
                                                                                $p_f_rbc,
                                                                                $z_rbc,
                                                                                $m_rbc,
                                                                                cache) setup = (cache = allocate_cache(m_rbc))
RBC_LIKELIHOODS["kalman"]["likelihood_gradient"] = @benchmarkable gradient(p -> benchmark_likelihood_kalman(p,
                                                                                $p_f_rbc,
                                                                                $z_rbc,
                                                                                $m_rbc,
                                                                                cache),
                                                                           $p_rbc) setup = (cache = allocate_cache(m_rbc))

# filled cache
generate_perturbation(m_rbc, p_rbc; p_f = p_f_rbc, cache = cache_rbc)  # solve to fill it
RBC_LIKELIHOODS["kalman"]["likelihood_reused_cache"] = @benchmarkable benchmark_likelihood_kalman($p_rbc,
                                                                                $p_f_rbc,
                                                                                $z_rbc,
                                                                                $m_rbc,
                                                                                $cache_rbc)

# gradient of kalman with a cache already filled
const cache_grad_rbc = allocate_cache(m_rbc)
generate_perturbation(m_rbc, p_rbc; p_f = p_f_rbc, cache = cache_grad_rbc)
RBC_LIKELIHOODS["kalman"]["likelihood_gradient_reused_cache"] = @benchmarkable gradient(p -> benchmark_likelihood_kalman(p,
                                                                                $p_f_rbc,
                                                                                $z_rbc,
                                                                                $m_rbc,
                                                                                $cache_grad_rbc),
                                                                                $p_rbc)

# joint-likelihood
function benchmark_joint_likelihood(p, ϵ, x0, p_f, z, m, cache)
    sol = generate_perturbation(m, p; p_f, cache)
    # (sol.retcode == :Success) || return 0.0
    return solve(sol, x0, (0, length(z)); noise = ϵ, observables = z).logpdf
end

# likelihood and gradient with cache reallocated each time
RBC_LIKELIHOODS["joint"]["likelihood"] = @benchmarkable benchmark_joint_likelihood($p_rbc,
                                                                                $ϵ_rbc,
                                                                                $x0_rbc,
                                                                                $p_f_rbc,
                                                                                $z_rbc,
                                                                                $m_rbc,
                                                                                cache) setup = (cache = allocate_cache(m_rbc))
RBC_LIKELIHOODS["joint"]["likelihood_gradient_p"] = @benchmarkable gradient(p -> benchmark_joint_likelihood(p,
                                                                                $ϵ_rbc,
                                                                                $x0_rbc,
                                                                                $p_f_rbc,
                                                                                $z_rbc,
                                                                                $m_rbc,
                                                                                cache),
                                                                            $p_rbc) setup = (cache = allocate_cache(m_rbc))
RBC_LIKELIHOODS["joint"]["likelihood_gradient_p_ϵ"] = @benchmarkable gradient((p, ϵ) -> benchmark_joint_likelihood(p,
                                                                                ϵ,
                                                                                $x0_rbc,
                                                                                $p_f_rbc,
                                                                                $z_rbc,
                                                                                $m_rbc,
                                                                                cache),
                                                                              $p_rbc,
                                                                              $ϵ_rbc) setup = (cache = allocate_cache(m_rbc))
RBC_LIKELIHOODS["joint"]["likelihood_gradient_ϵ"] = @benchmarkable gradient(ϵ -> benchmark_joint_likelihood($p_rbc,
                                                                                ϵ,
                                                                                $x0_rbc,
                                                                                $p_f_rbc,
                                                                                $z_rbc,
                                                                                $m_rbc,
                                                                                cache),
                                                                            $ϵ_rbc) setup = (cache = allocate_cache(m_rbc))

# gradient with a cache already filled
const cache_rbc_joint_grad = allocate_cache(m_rbc)
generate_perturbation(m_rbc, p_rbc; p_f = p_f_rbc, cache = cache_rbc_joint_grad)
RBC_LIKELIHOODS["joint"]["likelihood_reused"] = @benchmarkable benchmark_joint_likelihood($p_rbc,
                                                                                $ϵ_rbc,
                                                                                $x0_rbc,
                                                                                $p_f_rbc,
                                                                                $z_rbc,
                                                                                $m_rbc,
                                                                                $cache_grad_rbc)
RBC_LIKELIHOODS["joint"]["likelihood_gradient_p_reused"] = @benchmarkable gradient(p -> benchmark_joint_likelihood(p,
                                                                                $ϵ_rbc,
                                                                                $x0_rbc,
                                                                                $p_f_rbc,
                                                                                $z_rbc,
                                                                                $m_rbc,
                                                                                $cache_grad_rbc),
                                                                                $p_rbc)
RBC_LIKELIHOODS["joint"]["likelihood_gradient_p_ϵ_reused"] = @benchmarkable gradient((p, ϵ) -> benchmark_joint_likelihood(p,
                                                                                ϵ,
                                                                                $x0_rbc,
                                                                                $p_f_rbc,
                                                                                $z_rbc,
                                                                                $m_rbc,
                                                                                $cache_grad_rbc),
                                                                                $p_rbc,
                                                                                $ϵ_rbc)
RBC_LIKELIHOODS["joint"]["likelihood_gradient_ϵ_reused"] = @benchmarkable gradient(ϵ -> benchmark_joint_likelihood($p_rbc,
                                                                                ϵ,
                                                                                $x0_rbc,
                                                                                $p_f_rbc,
                                                                                $z_rbc,
                                                                                $m_rbc,
                                                                                $cache_grad_rbc),
                                                                                $ϵ_rbc)

# Second order
# likelihood and gradient with cache reallocated each time
RBC_LIKELIHOODS["second_order_joint"]["likelihood"] = @benchmarkable benchmark_joint_likelihood($p_rbc,
                                                                                $ϵ_second_rbc,
                                                                                $x0_second_rbc,
                                                                                $p_f_rbc,
                                                                                $z_second_rbc,
                                                                                $m_second_rbc,
                                                                                cache) setup = (cache = allocate_cache(m_second_rbc))
RBC_LIKELIHOODS["second_order_joint"]["likelihood_gradient_p"] = @benchmarkable gradient(p -> benchmark_joint_likelihood(p,
                                                                                $ϵ_second_rbc,
                                                                                $x0_second_rbc,
                                                                                $p_f_rbc,
                                                                                $z_second_rbc,
                                                                                $m_second_rbc,
                                                                                cache),
                                                                                $p_rbc) setup = (cache = allocate_cache(m_second_rbc))
RBC_LIKELIHOODS["second_order_joint"]["likelihood_gradient_p_ϵ"] = @benchmarkable gradient((p, ϵ) -> benchmark_joint_likelihood(p,
                                                                                ϵ,
                                                                                $x0_second_rbc,
                                                                                $p_f_rbc,
                                                                                $z_second_rbc,
                                                                                $m_second_rbc,
                                                                                cache),
                                                                                $p_rbc,
                                                                                $ϵ_second_rbc) setup = (cache = allocate_cache(m_second_rbc))
RBC_LIKELIHOODS["second_order_joint"]["likelihood_gradient_ϵ"] = @benchmarkable gradient(ϵ -> benchmark_joint_likelihood($p_rbc,
                                                                                ϵ,
                                                                                $x0_second_rbc,
                                                                                $p_f_rbc,
                                                                                $z_second_rbc,
                                                                                $m_second_rbc,
                                                                                cache),
                                                                                $ϵ_second_rbc) setup = (cache = allocate_cache(m_second_rbc))

# gradient with a cache already filled
const cache_second_grad_rbc = allocate_cache(m_second_rbc)
generate_perturbation(m_second_rbc, p_rbc; p_f = p_f_rbc, cache = cache_second_grad_rbc)
RBC_LIKELIHOODS["second_order_joint"]["likelihood_reused"] = @benchmarkable benchmark_joint_likelihood($p_rbc,
                                                                                $ϵ_second_rbc,
                                                                                $x0_second_rbc,
                                                                                $p_f_rbc,
                                                                                $z_second_rbc,
                                                                                $m_second_rbc,
                                                                                $cache_second_grad_rbc)
RBC_LIKELIHOODS["second_order_joint"]["likelihood_gradient_p_reused"] = @benchmarkable gradient(p -> benchmark_joint_likelihood(p,
                                                                                $ϵ_second_rbc,
                                                                                $x0_second_rbc,
                                                                                $p_f_rbc,
                                                                                $z_second_rbc,
                                                                                $m_second_rbc,
                                                                                $cache_second_grad_rbc),
                                                                                $p_rbc)
RBC_LIKELIHOODS["second_order_joint"]["likelihood_gradient_p_ϵ_reused"] = @benchmarkable gradient((p, ϵ) -> benchmark_joint_likelihood(p,
                                                                                ϵ,
                                                                                $x0_second_rbc,
                                                                                $p_f_rbc,
                                                                                $z_second_rbc,
                                                                                $m_second_rbc,
                                                                                $cache_second_grad_rbc),
                                                                                $p_rbc,
                                                                                $ϵ_second_rbc)
RBC_LIKELIHOODS["second_order_joint"]["likelihood_gradient_ϵ_reused"] = @benchmarkable gradient(ϵ -> benchmark_joint_likelihood($p_rbc,
                                                                                ϵ,
                                                                                $x0_second_rbc,
                                                                                $p_f_rbc,
                                                                                $z_second_rbc,
                                                                                $m_second_rbc,
                                                                                $cache_second_grad_rbc),
                                                                                $ϵ_second_rbc)

# return the RBC_LIKELIHOODS benchmarks to the parent
RBC_LIKELIHOODS


text:
# """
# $(DocStringExtensions.TYPEDSIGNATURES)

# Description

# ```jldoctest
# julia> 1 + 1
# 2
# ```
# """
function graph_to_gfa(graph, outfile)
    open(outfile, "w") do io
        println(io, "H\tVN:Z:1.0")
        for vertex in Graphs.vertices(graph)
            if haskey(graph.vprops[vertex], :kmer)
                sequence = graph.vprops[vertex][:kmer]
            else
                sequence = graph.vprops[vertex][:sequence]
            end
#             depth = graph.vprops[vertex][:weight]
            depth = graph.vprops[vertex][:count]
            fields = ["S", "$vertex", sequence, "DP:f:$(depth)"]
            line = join(fields, '\t')
            println(io, line)
        end
        for edge in Graphs.edges(graph)
            overlap = graph.gprops[:k] - 1
            for o in graph.eprops[edge][:orientations]
#                 if !(!o.source_orientation && !o.destination_orientation)
                link = ["L",
                            edge.src,
                            o.source_orientation ? '+' : '-',
                            edge.dst,
                            o.destination_orientation ? '+' : '-',
                            "$(overlap)M"]
                line = join(link, '\t')
                println(io, line)
#                 end
            end
        end
    end
    return outfile
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function open_fastx(path::String)
    if isfile(path)
        io = open(path)
    elseif occursin(r"^ftp", path) || occursin(r"^http", path)
        path = replace(path, r"^ftp:" => "http:")
        io = IOBuffer(HTTP.get(path).body)
    else
        error("unable to locate file $path")
    end
    path_base = basename(path)
    if occursin(r"\.gz$", path_base)
        io = CodecZlib.GzipDecompressorStream(io)
        path_base = replace(path_base, ".gz" => "")
    end
    if occursin(r"\.(fasta|fna|faa|fa)$", path_base)
        fastx_io = FASTX.FASTA.Reader(io)
    elseif occursin(r"\.(fastq|fq)$", path_base)
        fastx_io = FASTX.FASTQ.Reader(io)
    else
        @show path_base
        error()
    end
    return fastx_io
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function list_databases(;address, username, password)
    cmd = "show databases"
    database = "system"
    cmd = cypher(;address, username, password, database, cmd)
    return DataFrames.DataFrame(uCSV.read(open(cmd), header=1, quotes='"', encodings=Dict("FALSE" => false, "TRUE" => true))...)
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function create_database(;database, address, username, password)
    current_databases = list_databases(;address, username, password)
    if database in current_databases[!, "name"]
        return
    else
        f = run
        cmd = "create database $(database)"
        # switch database to system, so that we can create the user-specific database in the system
        database = "system"
        run(cypher(;address, username, password, database, cmd, f))
    end
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function cypher(;address, username, password, database, cmd)
    return `cypher-shell --address $address --username $username --password $password --database $(database) --format auto $(cmd)`
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function open_gff(path::String)
    if isfile(path)
        io = open(path)
    elseif occursin(r"^ftp", path) || occursin(r"^http", path)
        path = replace(path, r"^ftp:" => "http:")
        io = IOBuffer(HTTP.get(path).body)
    else
        error("unable to locate file $path")
    end
    path_base = basename(path)
    if occursin(r"\.gz$", path_base)
        io = CodecZlib.GzipDecompressorStream(io)
#         path_base = replace(path_base, ".gz" => "")
    end
    gff_io = GFF3.Reader(io)
    return gff_io
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function plot_kmer_frequency_spectra(counts; log_scale = log2, kwargs...)
    kmer_counts_hist = StatsBase.countmap(c for c in counts)
    xs = collect(keys(kmer_counts_hist))
    ys = collect(values(kmer_counts_hist))
    if isa(log_scale, Function)
        xs = log_scale.(xs)
        ys = log_scale.(ys)
    end

    p = StatsPlots.plot(
        xs,
        ys,
        xlims = (0, maximum(xs) + max(1, ceil(0.1 * maximum(xs)))),
        ylims = (0, maximum(ys) + max(1, ceil(0.1 * maximum(ys)))),
        seriestype = :scatter,
        legend = false,
        xlabel = isa(log_scale, Function) ? "$(log_scale)(observed frequency)" : "observed frequency",
        ylabel = isa(log_scale, Function) ? "$(log_scale)(# of kmers)" : "observed frequency",
        ;kwargs...
    )
    return p
end

# """
# $(DocStringExtensions.TYPEDSIGNATURES)

# Description

# ```jldoctest
# julia> 1 + 1
# 2
# ```
# """
function convert(args)
    @show args
    in_type = missing
    out_type = missing
    if occursin(r"\.jld2$", args["in"])
        in_type = :jld2
    elseif occursin(r"\.gfa$", args["in"])
        in_type = :gfa
    elseif occursin(r"\.neo4j$", args["in"])
        in_type = :neo4j
    end

    if occursin(r"\.jld2$", args["out"])
        out_type = :jld2
    elseif occursin(r"\.gfa$", args["out"])
        out_type = :gfa
    elseif occursin(r"\.neo4j$", args["out"])
        out_type = :neo4j
    end

    if ismissing(in_type) || ismissing(out_type)
        error("unable to determine in and out types")
    end

    if (in_type == :jld2) && (out_type == :jld2)
        # done
    elseif (in_type == :jld2) && (out_type != :jld2)
        # convert out of jld2
        if out_type == :gfa
            loaded = FileIO.load(args["in"])
            @assert haskey(loaded, "graph")
            graph = loaded["graph"]
            graph_to_gfa(graph, args["out"])
        end
    elseif (in_type != :jld2) && (out_type == :jld2)
        # convert into jld2
    else
        # need to convert from input to jld2 as an intermediate first
    end
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function plot_graph(graph)

#     kmer_counts = MetaGraphs.get_prop(graph, :kmer_counts)
    kmers = [MetaGraphs.get_prop(graph, v, :kmer) for v in Graphs.vertices(graph)]
    counts = [MetaGraphs.get_prop(graph, v, :count) for v in Graphs.vertices(graph)]
    scale = 150

    n = Graphs.nv(graph)
    p = GraphRecipes.graphplot(
        graph,
#         markersize = 1/log2(n),
        markersize = 1/2^2,
        size = (2 * scale * log(n), scale * log(n)),
        node_weights = counts,
        names = kmers)
end

# """
# $(DocStringExtensions.TYPEDSIGNATURES)

# A short description of the function

# ```jldoctest
# julia> 1 + 1
# 2
# ```
# """
# function my_plot(graph::Mycelia.KmerGraph)
#     graph_hash = hash(sort(graph.graph.fadjlist), hash(graph.graph.ne))
#     filename = "/assets/images/$(graph_hash).svg"
#     p = Mycelia.plot_graph(graph)
#     Plots.savefig(p, dirname(pwd()) * filename)
#     display(p)
#     display("text/markdown", "![]($filename)")
# end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function save_graph(graph::Graphs.AbstractGraph, outfile::String)
    if !occursin(r"\.jld2$", outfile)
        outfile *= ".jld2"
    end
    FileIO.save(outfile, Dict("graph" => graph))
    return outfile
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function load_graph(file::String)
    return FileIO.load(file)["graph"]
end

max_stars_repo_name: cjprybol/Eisenia

max_issues_repo_name: cjprybol/Eisenia

max_forks_repo_name: cjprybol/Eisenia


text:
export
    PDPageObject,
    PDPageElement,
    PDPageObjectGroup,
    PDPageTextObject,
    PDPageTextRun,
    PDPageMarkedContent,
    PDPageInlineImage,
    PDPage_BeginGroup,
    PDPage_EndGroup,
    GState

import Base: get, show, getindex, setindex!, delete!, >
import ..Common: CDRect
using ..Cos: CosComment

using LinearAlgebra

abstract type PDPage end
abstract type PDXObject end

"""
```
    PDPageObject
```

The content streams associated with PDF pages contain the objects that can be rendered.
These objects are represented by `PDPageObject`. These objects can contain
a postfix notation based operator prefixed by its operands like:
```
(Draw this text) Tj
```
As can be seen above, the string object is a [`CosString`](@ref) which is a parameter to the
operand `Tj` or draw text. These class of objects are represented by
[`PDPageElement`](@ref).

However, there are certain objects which only provide grouping information or begin and end
markers for grouping information. For example, a text object:
```
BT
    /F1 11 Tf  %selectfont
    (Draw this text) Tj
ET
```
These kind of objects are represented by [`PDPageObjectGroup`](@ref). In this case, the
[`PDPageObjectGroup`](@ref) contains four [`PDPageElement`](@ref). Namely, represented as operators `BT`,
`Tf`, `Tj`, `ET`.

`PDPageElement` and [`PDPageObjectGroup`](@ref) can be extended by composition. Hence, there are
more specialized objects that can be seen as well.

"""
abstract type PDPageObject end

"""
```
    PDPageElement
```
A representation of a content object with operator and operand. See [`PDPageObject`](@ref)
for more details.
"""
mutable struct PDPageElement{S} <: PDPageObject
    t::Symbol
    version::Tuple{Int,Int}
    noperand::Int
    operands::Vector{CosObject}
end

PDPageElement(ts::AbstractString,ver::Tuple{Int,Int},nop::Int=0)=
  PDPageElement{Symbol(ts)}(Symbol(ts),ver,nop,Vector{CosObject}())

function show(io::IO, e::PDPageElement)
    for op in e.operands
        show(io, op)
        print(io, ' ')
    end
    print(io, String(e.t))
end

"""
```
    PDPageObjectGroup
```
A representation of a content object that encloses other content objects. See
[`PDPageObject`](@ref) for more details.
"""
mutable struct PDPageObjectGroup <: PDPageObject
    isEOG::Bool
    objs::Vector{Union{PDPageObject, CosObject}}
    PDPageObjectGroup(isEOG::Bool=false) =
        new(isEOG, Vector{Union{PDPageObject, CosObject}}())
end

Base.isempty(grp::PDPageObjectGroup) = isempty(grp.objs)

function load_objects(grp::PDPageObjectGroup, bis::IO)
    while(!grp.isEOG && !eof(bis))
        obj = parse_value(bis, get_pdfcontentops)
        collect_object(grp, obj, bis)
    end
end

# Ignore comments in content stream
collect_object(grp::PDPageObjectGroup, obj::CosComment, bis::IO) = nothing

collect_object(grp::PDPageObjectGroup, obj::CosObject, bis::IO) =
    push!(grp.objs, obj)

function populate_element(grp::PDPageObjectGroup, elem::PDPageElement)
    #Find operands for the Operator
    if (elem.noperand >= 0)
        for i=1:elem.noperand
            operand=pop!(grp.objs)
            pushfirst!(elem.operands,operand)
        end
    else
        while(isa(grp.objs[end],CosObject))
            operand=pop!(grp.objs)
            pushfirst!(elem.operands,operand)
        end
    end
end

function collect_object(grp::PDPageObjectGroup, elem::PDPageElement,
                        bis::IO)
    populate_element(grp, elem)
    push!(grp.objs, elem)
    return elem
end

"""
```
    PDPageTextObject
```
A [`PDPageObjectGroup`](@ref) object that represents a block of text. See [`PDPageObject`](@ref)
for more details.
"""
mutable struct PDPageTextObject <: PDPageObject
    group::PDPageObjectGroup
    PDPageTextObject()=new(PDPageObjectGroup())
end

"""
```
    PDPageMarkedContent
```
A [`PDPageObjectGroup`](@ref) object that represents a group of a object that is logically
grouped together in case of a structured PDF document.
"""
mutable struct PDPageMarkedContent <: PDPageObject
    group::PDPageObjectGroup
    PDPageMarkedContent()=new(PDPageObjectGroup())
end

"""
```
    PDPageInlineImage
```
Most images in PDF documents are defined in the PDF document and referenced from the page
content stream. [`PDPageInlineImage`](@ref) objects are directly defined in the page
content stream.
"""
mutable struct PDPageInlineImage <: PDPageObject
    params::CosDict
    data::Vector{UInt8}
    isRead::Bool
    PDPageInlineImage()=new(CosDict(),Vector{UInt8}(),false)
end

function show(io::IO, im::PDPageInlineImage)
    println(io, "BI")
    for (k, v) in get(im.params)
        println("$k $v")
    end
    println("ID")
    println("< $(length(im.data)) bytes... >")
    println(io, "EI")
end

#=
```
    PDPage_BeginInlineImage
```
A [`PDPageElement`](@ref) that represents the beginning of an inline image.
=#
mutable struct PDPage_BeginInlineImage <: PDPageObject
    elem::PDPageElement
    PDPage_BeginInlineImage(ts::AbstractString, ver::Tuple{Int,Int}, nop, ::Type)=
        new(PDPageElement(ts,ver,nop))
end

function collect_object(grp::PDPageObjectGroup, beg::PDPage_BeginInlineImage,
                        bis::IO)
    newobj=PDPageInlineImage()

    while(!newobj.isRead)
        value = parse_value(bis, get_pdfcontentops)
        collect_inline_image(newobj, value, bis)
    end
    push!(grp.objs, newobj)
    return newobj
end

"""
```
    PDPageTextRun
```
In PDF text may not be contiguous as there may be chnge of font, style, graphics rendering
parameters. `PDPageTextRun` is a unit of text which can be rendered without any change to
the graphical parameters. There is no guarantee that a text run will represent a meaningful
word or sentence.

`PDPageTextRun` is a composition implementation of [`PDPageElement`](@ref).
"""
mutable struct PDPageTextRun <: PDPageObject
    ss::Vector{Union{CosXString, CosLiteralString, CosInt, CosFloat}}
    elem::PDPageElement
    PDPageTextRun(ts::AbstractString,ver::Tuple{Int,Int},nop::Int=0) =
        new(Vector{Union{CosXString, CosLiteralString, CosInt, CosFloat}}(),
            PDPageElement(ts, ver, nop))
end

function show(io::IO, tr::PDPageTextRun)
    print(io, "[")
    for s in tr.ss
        print(io, s)
    end
    print(io, "]")
end

function collect_object(grp::PDPageObjectGroup, tr::PDPageTextRun,
                        bis::IO)
    elem = collect_object(grp, tr.elem, bis)
    for operand in elem.operands
        if isa(operand, CosString)
            push!(tr.ss, operand)
        elseif isa(operand, CosArray)
            for td in get(operand)
                push!(tr.ss, td)
            end
        end
    end
    val = pop!(grp.objs)
    push!(grp.objs, tr)
    return tr
end

function collect_inline_image(img::PDPageInlineImage, name::CosName, bis::IO)
    value = parse_value(bis, get_pdfcontentops)
    set!(img.params, name, value)
end

function collect_inline_image(img::PDPageInlineImage, elem::PDPageElement,
                              bis::IO)
    elem.t !== :ID && return img
    while(!img.isRead && !eof(bis))
        b1 = _peekb(bis)
        if (b1 == LATIN_UPPER_E)
            mark(bis)
            skip(bis, 1)
            b2 = _peekb(bis)
            if (b2 == LATIN_UPPER_I)
                skip(bis, 1); b3 = _peekb(bis)
                if (ispdfspace(b3))
                    skip(bis, 1)
                    img.isRead = true
                    unmark(bis)
                    break
                else
                    reset(bis)
                end
            else
                reset(bis)
            end
        end
        push!(img.data, b1)
        skip(bis, 1)
    end
    return img
end

"""
```
    PDPage_BeginGroup
```
A [`PDPageElement`](@ref) that represents the beginning of a group object.
"""
mutable struct PDPage_BeginGroup <: PDPageObject
    elem::PDPageElement
    objT::Type
    PDPage_BeginGroup(ts::AbstractString,ver::Tuple{Int,Int},nop,t::Type) =
        new(PDPageElement(ts,ver,nop),t)
end

"""
```
    PDPage_EndGroup
```
A [`PDPageElement`](@ref) that represents the end of a group object.
"""
mutable struct PDPage_EndGroup
    elem::PDPageElement
    PDPage_EndGroup(ts::AbstractString,ver::Tuple{Int,Int},nop) =
        new(PDPageElement(ts,ver,nop))
end

show(io::IO, e::PDPage_BeginGroup) = show(io, e.elem)

show(io::IO, e::PDPage_EndGroup) = show(io, e.elem)

function collect_object(grp::PDPageObjectGroup, beg::PDPage_BeginGroup,
                        bis::IO)
    populate_element(grp,beg.elem)
    newobj=beg.objT()
    push!(newobj.group.objs,beg.elem)
    load_objects(newobj.group,bis)
    push!(grp.objs, newobj)
    return newobj
end

function collect_object(grp::PDPageObjectGroup, elem::PDPage_EndGroup,
                        bis::IO)
    collect_object(grp,elem.elem,bis)
    grp.isEOG = true
    return grp
end

#=
|Operator|PostScript Equivalent|Description|Table|
|:--------|:----------------------|:------------|------:|
|BX||(PDF 1.1) Begin compatibility section|32|
|EX||(PDF 1.1) End compatibility section|32|
|cm|concat|Concatenate matrix to current transformation matrix|57|
|d|setdash|Set line dash pattern|57|
|gs||(PDF 1.2) Set parameters from graphics state parameter dictionary|57|
|i|setflat|Set flatness tolerance|57|
|j|setlinejoin|Set line join style|57|
|J|setlinecap|Set line cap style|57|
|M|setmiterlimit|Set miter limit|57|
|q|gsave|Save graphics state|57|
|Q|grestore|Restore graphics state|57|
|ri||Set color rendering intent|57|
|w|setlinewidth|Set line width|57|
|c|curveto|Append curved segment to path (three control points)|59|
|h|closepath|Close subpath|59|
|l|lineto|Append straight line segment to path|59|
|m|moveto|Begin new subpath|59|
|re||Append rectangle to path|59|
|v|curveto|Append curved segment to path (initial point replicated)|59|
|y|curveto|Append curved segment to path (final point replicated)|59|
|b|closepath, fill, stroke|Close, fill, and stroke path using nonzero winding number rule|60|
|B|fill, stroke|Fill and stroke path using nonzero winding number rule|60|
|b*|closepath, eofill, stroke|Close, fill, and stroke path using even-odd rule|60|
|B*|eofill, stroke|Fill and stroke path using even-odd rule|60|
|f|fill|Fill path using nonzero winding number rule|60|
|F|fill|Fill path using nonzero winding number rule (obsolete)|60|
|f*|eofill|Fill path using even-odd rule|60|
|n||End path without filling or stroking|60|
|s|closepath, stroke|Close and stroke path|60|
|S|stroke|Stroke path|60|
|W|clip|Set clipping path using nonzero winding number rule|61|
|W*|eoclip|Set clipping path using even-odd rule|61|
|CS|setcolorspace|(PDF 1.1) Set color space for stroking operations|74|
|cs|setcolorspace|(PDF 1.1) Set color space for nonstroking operations|74|
|G|setgray|Set gray level for stroking operations|74|
|g|setgray|Set gray level for nonstroking operations|74|
|K|setcmykcolor|Set CMYK color for stroking operations|74|
|k|setcmykcolor|Set CMYK color for nonstroking operations|74|
|RG|setrgbcolor|Set RGB color for stroking operations|74|
|rg|setrgbcolor|Set RGB color for nonstroking operations|74|
|SC|setcolor|(PDF 1.1) Set color for stroking operations|74|
|sc|setcolor|(PDF 1.1) Set color for nonstroking operations|74|
|SCN|setcolor|(PDF 1.2) Set color for stroking operations\n|||(ICCBased and special colour spaces)|74|
|scn|setcolor|(PDF 1.2) Set color for nonstroking operations\n|||(ICCBased and special colour spaces)|74|
|sh|shfill|(PDF 1.3) Paint area defined by shading pattern|77|
|Do||Invoke named XObject|87|
|BI||Begin inline image object|92|
|EI||End inline image object|92|
|ID||Begin inline image data|92|
|BT||Begin text object|107|
|ET||End text object|107|
|T*||Move to start of next text line|108|
|Td||Move text position|108|
|TD||Move text position and set leading|108|
|Tj|show|Show text|109|
|TJ||Show text, allowing individual glyph positioning|109|
|\'||Move to next line and show text|109|
|\"||Set word and character spacing, move to next line, and show text|109|
|Tc||Set character spacing|105|
|Tf|selectfont|Set text font and size|105|
|TL||Set text leading|105|
|Tr||Set text rendering mode|105|
|Ts||Set text rise|105|
|Tw||Set word spacing|105|
|Tz||Set horizontal text scaling|105|
|d0|setcharwidth|Set glyph width in Type 3 font|113|
|d1|setcachedevice|Set glyph width and bounding box in Type 3 font|113|
|BDC||(PDF 1.2) Begin marked-content sequence with property list|320|
|BMC||(PDF 1.2) Begin marked-content sequence|320|
|DP||(PDF 1.2) Define marked-content point with property list|320|
|EMC||(PDF 1.2) End marked-content sequence|320|
|MP||(PDF 1.2) Define marked-content point|320|
=#
const PD_CONTENT_OPERATORS = Dict(
"\'"=>[PDPageTextRun,"\'",(1,0),1],
"\""=>[PDPageTextRun,"\"",(1,0),3],
"b"=>[PDPageElement,"b",(1,0),0],
"b*"=>[PDPageElement,"b*",(1,0),0],
"B"=>[PDPageElement,"B",(1,0),0],
"B*"=>[PDPageElement,"B*",(1,0),0],
"BDC"=>[PDPage_BeginGroup,"BDC",(1,2),2,PDPageMarkedContent],
"BI"=>[PDPage_BeginInlineImage,"BI",(1,0),0,PDPageInlineImage],
"BMC"=>[PDPage_BeginGroup,"BMC",(1,2),1,PDPageMarkedContent],
"BT"=>[PDPage_BeginGroup,"BT",(1,0),0,PDPageTextObject],
"BX"=>[PDPageElement,"BX",(1,1),0],
"c"=>[PDPageElement,"c",(1,0),6],
"cm"=>[PDPageElement,"cm",(1,0),6],
"cs"=>[PDPageElement,"cs",(1,1),1],
"CS"=>[PDPageElement,"CS",(1,1),1],
"d"=>[PDPageElement,"d",(1,0),2],
"d0"=>[PDPageElement,"d0",(1,0),2],
"d1"=>[PDPageElement,"d1",(1,0),6],
"Do"=>[PDPageElement,"Do",(1,0),1],
"DP"=>[PDPageElement,"DP",(1,2),0],
"EI"=>[PDPageElement,"EI",(1,0),0],
"EMC"=>[PDPage_EndGroup,"EMC",(1,2),0],
"ET"=>[PDPage_EndGroup,"ET",(1,0),0],
"EX"=>[PDPageElement,"EX",(1,1),0],
"f"=>[PDPageElement,"f",(1,0),0],
"f*"=>[PDPageElement,"f*",(1,0),0],
"F"=>[PDPageElement,"F",(1,0),0],
"g"=>[PDPageElement,"g",(1,0),1],
"G"=>[PDPageElement,"G",(1,0),1],
"gs"=>[PDPageElement,"gs",(1,2),1],
"h"=>[PDPageElement,"h",(1,0),0],
"i"=>[PDPageElement,"i",(1,0),1],
"ID"=>[PDPageElement,"ID",(1,0),0],
"j"=>[PDPageElement,"j",(1,0),1],
"J"=>[PDPageElement,"J",(1,0),1],
"k"=>[PDPageElement,"k",(1,0),4],
"K"=>[PDPageElement,"K",(1,0),4],
"l"=>[PDPageElement,"l",(1,0),2],
"m"=>[PDPageElement,"m",(1,0),2],
"M"=>[PDPageElement,"M",(1,0),1],
"MP"=>[PDPageElement,"MP",(1,2),0],
"n"=>[PDPageElement,"n",(1,0),0],
"q"=>[PDPageElement,"q",(1,0),0],
"Q"=>[PDPageElement,"Q",(1,0),0],
"re"=>[PDPageElement,"re",(1,0),4],
"rg"=>[PDPageElement,"rg",(1,0),3],
"RG"=>[PDPageElement,"RG",(1,0),3],
"ri"=>[PDPageElement,"ri",(1,0),1],
"s"=>[PDPageElement,"s",(1,0),0],
"S"=>[PDPageElement,"S",(1,0),0],
"sc"=>[PDPageElement,"sc",(1,1),-1],
"SC"=>[PDPageElement,"SC",(1,1),-1],
"scn"=>[PDPageElement,"scn",(1,2),-1],
"SCN"=>[PDPageElement,"SCN",(1,2),-1],
"sh"=>[PDPageElement,"sh",(1,3),1],
"T*"=>[PDPageElement,"T*",(1,0),0],
"Tc"=>[PDPageElement,"Tc",(1,0),1],
"Td"=>[PDPageElement,"Td",(1,0),2],
"TD"=>[PDPageElement,"TD",(1,0),2],
"Tf"=>[PDPageElement,"Tf",(1,0),2],
"Tj"=>[PDPageTextRun,"Tj",(1,0),1],
"TJ"=>[PDPageTextRun,"TJ",(1,0),1],
"TL"=>[PDPageElement,"TL",(1,0),1],
"Tm"=>[PDPageElement,"Tm",(1,0),6],
"Tr"=>[PDPageElement,"Tr",(1,0),1],
"Ts"=>[PDPageElement,"Ts",(1,0),1],
"Tw"=>[PDPageElement,"Tw",(1,0),1],
"Tz"=>[PDPageElement,"Tz",(1,0),1],
"v"=>[PDPageElement,"v",(1,0),4],
"w"=>[PDPageElement,"w",(1,0),1],
"W"=>[PDPageElement,"W",(1,0),0],
"W*"=>[PDPageElement,"W*",(1,0),0],
"y"=>[PDPageElement,"y",(1,0),4]
)

function get_pdfcontentops(b::Vector{UInt8})
    # PDF content operators are never longer than 3 bytes and may not be
    # delimited. Hence, search for the longest 3 byte keyword, then 2 bytes
    # and lastly 1
    arr, l, sb = nothing, length(b), b
    if l > 3
        sb, l = b[1:3], 3
    end
    s = l
    while arr == nothing && s > 0
        arr = get(PD_CONTENT_OPERATORS, String(sb[1:s]), nothing)
        s -= 1
    end
    if arr !== nothing
        return s+1, eval(Expr(:call, arr...))
    end
    error("Invalid content operator: $(String(b))")
end

struct TextLayout
    lbx::Float32
    lby::Float32
    rbx::Float32
    rby::Float32
    rtx::Float32
    rty::Float32
    ltx::Float32
    lty::Float32
    text::String
    fontname::CosName
    fontflags::UInt32
end

CDRect(t::TextLayout) = CDRect(min(t.lbx, t.rbx, t.rtx, t.ltx),
                               min(t.lby, t.rby, t.rty, t.lty),
                               max(t.lbx, t.rbx, t.rtx, t.ltx),
                               max(t.lby, t.rby, t.rty, t.lty))

@inline function width(tl)
    dx = tl.rbx - tl.lbx; dy = tl.rby - tl.lby
    return sqrt(dx*dx + dy*dy)
end

@inline function height(tl)
    dx = tl.ltx - tl.lbx; dy = tl.lty - tl.lby
    return sqrt(dx*dx + dy*dy)
end

Base.:>(tl1::TextLayout, tl2::TextLayout) = Base.isless(tl2, tl1)

@inline function Base.isless(tl1::TextLayout, tl2::TextLayout)
    y2 = max(tl2.lby, tl2.rby, tl2.rty, tl2.lty)
    x2 = min(tl2.lbx, tl2.rbx, tl2.rtx, tl2.ltx)

    y1 = max(tl1.lby, tl1.rby, tl1.rty, tl1.lty)
    x1 = min(tl1.lbx, tl1.rbx, tl1.rtx, tl1.ltx)

    dy = y1 - y2
    dx = x1 - x2
    # This will ensure superscripts with smaller fonts where baseline
    # shift is less than half the height of the baseline fonts is aligned
    # to current line and not an additional line above
    ytol1 = (tl1.lty - tl1.lby)/2
    ytol2 = (tl2.lty - tl2.lby)/2
    ytol = abs(ytol1) > abs(ytol2) ? ytol1 : ytol2
    dy < -ytol && return true
    dy >  ytol && return false
    return dx > 0
end

# using DataStructures

mutable struct GState{T}
    state::Vector{Dict{Symbol, Any}}
    GState{T}() where T = new(init_graphics_state())
end

new_gstate(state::GState{T}) where T = GState{T}()
Base.setindex!(gs::GState, v::V, k::Symbol) where V = (gs.state[end][k] = v)
Base.get(gs::GState, k::Symbol, R::Type) = gs.state[end][k]::R
Base.get(gs::GState, k::Symbol, defval::S, R::Type=S) where S  =
    get(gs.state[end], k, defval)::Union{R, S}
Base.get!(gs::GState, k::Symbol, defval::S, R::Type=S) where S =
    get!(gs.state[end], k, defval)::Union{R, S}
Base.delete!(gs::GState, k::Symbol) = delete!(gs.state[end], k)
save!(gs::GState) = (push!(gs.state, copy(gs.state[end])); gs)
restore!(gs::GState) = (pop!(gs.state); gs)

@inline function init_graphics_state()
    state = Vector{Dict{Symbol, Any}}(undef, 0)
    push!(state, Dict{Symbol, Any}())

    state[end][:text_layout] = Vector{TextLayout}()

    # Histogram along the y-axis. Not used currently.
    state[end][:h_profile] = Dict{Int,Int}()

    # Graphics state
    state[end][:CTM] = Matrix{Float32}(I, 3, 3)

    # Text states
    state[end][:Tc] = 0f0
    state[end][:Tw] = 0f0
    state[end][:Tz] = 100f0
    state[end][:TL] = 0f0
    state[end][:Tr] = 0
    state[end][:Ts] = 0f0
    return state
end

function show_text_layout!(io::IO, state::GState)
    #Make sure to deepcopy. Otherwise the data structures will be lost
    heap = get(state, :text_layout, Vector{TextLayout})
    sort!(heap, lt= > )
    szdict = get(state, :h_profile, Dict{Int, Int})

    x = 0f0
    y = -1f0

    pairs = sort(collect(szdict), lt=(x,y)->(x[2] > y[2]))
    length(pairs) == 0 && return io
    iht, num = pairs[1]
    ht = iht*0.1f0

    # Courier font X width is 600
    xwr = 0.6f0
    ph = 0f0
    npc = 0
    for i = 1:lastindex(heap)
        tlayout = heap[i]
        h = height(tlayout)
        if h > 7f0*ht
            ht = h
        end
        xw = xwr*ht
        nc = length(tlayout.text)
        w = width(tlayout)/nc
        (w > 0.1f0 && h > 0.1f0) || continue
        if (ht > h)
            while (y > tlayout.lty)
                print(io, '\n')
                y -= ht
                x = 0f0
            end
        else
            while (y > tlayout.lby + ht)
                print(io, '\n')
                y -= ht
                x = 0f0
            end
            y = tlayout.lby
        end
        y = tlayout.lby
        # For subscripts and superscripts insert a space on both sides
        # However, ignore the same for dropcaps.
        (x > tlayout.lbx - xw) && (ph < h || (ph > h && npc > 1)) && print(io, ' ')
        while x < tlayout.lbx - xw
            print(io, ' ')
            x += xw
        end
        x = tlayout.lbx
        print(io, tlayout.text)
        x += width(tlayout)
        while x < tlayout.rbx - xw
            print(io, ' ')
            x += xw
        end
        x = tlayout.rbx
        ph = h
        npc = nc
    end
    return io
end

@inline function evalContent!(grp::PDPageObjectGroup, state::GState)
    for obj in grp.objs
        evalContent!(obj, state)
    end
    return state
end

function eval_unicode_mapping(tr::PDPageTextRun, state::GState)
    fontname, font = get(state, :font, (cn"", CosNull), Tuple{CosName, PDFont})
    if font.fum === nothing
        src = get(state, :source, nothing, Union{PDPage, PDXObject})
        src === nothing && error("Graphics state :source is not configured")
        warn_no_unicode_mapping(src, font, fontname)
        @warn "Text run $tr may be decoded as ASCII"
    end
    return fontname, font
end

function warn_no_unicode_mapping(page::PDPage, font, fontname)
    pageno = pdPageGetPageNumber(page)
    @warn "No unicode mapping for font $fontname at Page $pageno"
end

function warn_no_unicode_mapping(xo::PDXObject, font, fontname)
    ref = CosIndirectObjectRef(xo.cosXObj.num, xo.cosXObj.gen)
    @warn "No unicode mapping for font $fontname at XObject at $ref"
end

@inline function evalContent!(tr::PDPageTextRun, state::GState)
    evalContent!(tr.elem, state)
    tfs = get(state, :fontsize, 0f0)
    th  = get(state, :Tz, Float32)/100f0
    ts  = get(state, :Ts, Float32)
    tc  = get(state, :Tc, Float32)
    tw  = get(state, :Tw, Float32)
    tm  = get(state, :Tm, Matrix{Float32})
    ctm = get(state, :CTM, Matrix{Float32})
    trm = tm*ctm

    fontname, font = eval_unicode_mapping(tr, state)

    heap = get(state, :text_layout, Vector{TextLayout})
    text, w, h = get_TextBox(tr.ss, font, tfs, tc, tw, th)

    d = get(state, :h_profile, Dict{Int, Int})
    ih = round(Int, h*10)
    d[ih] = get(d, ih, 0) + length(text)

    tb = [0f0 0f0 1f0; w 0f0 1f0; w h 1f0; 0f0 h 1f0]*trm
    if !get(state, :in_artifact, false)
        tl = TextLayout(tb[1,1], tb[1,2], tb[2,1], tb[2,2],
                        tb[3,1], tb[3,2], tb[4,1], tb[4,2],
                        text, fontname, font.flags)
        push!(heap, tl)
    end
    offset_text_pos!(w, 0f0, state)
    return state
end

@inline function evalContent!(pdo::PDPageTextObject,
                              state::GState)
    state[:Tm]  = Matrix{Float32}(I, 3, 3)
    state[:Tlm] = Matrix{Float32}(I, 3, 3)
    state[:Trm] = Matrix{Float32}(I, 3, 3)
    evalContent!(pdo.group, state)
    delete!(state, :Tm)
    delete!(state, :Tlm)
    delete!(state, :Trm)
    return state
end

@inline function evalContent!(pdo::PDPageMarkedContent, state::GState)
    tag = pdo.group.objs[1].operands[1] # can be used for XML tagging.
    if tag == cn"Artifact"
        state[:in_artifact] = true
        evalContent!(pdo.group, state)
        delete!(state, :in_artifact)
        return state
    end
    return evalContent!(pdo.group, state)
end

evalContent!(pdo::PDPageElement{S}, state::GState) where S = state

evalContent!(pdo::PDPageElement{:q}, state::GState) = save!(state)

evalContent!(pdo::PDPageElement{:Q}, state::GState) = restore!(state)

@inline function evalContent!(pdo::PDPageElement{:cm}, state::GState)
    a = get(pdo.operands[1])
    b = get(pdo.operands[2])
    c = get(pdo.operands[3])
    d = get(pdo.operands[4])
    e = get(pdo.operands[5])
    f = get(pdo.operands[6])
    cm  = [a b 0f0; c d 0f0; e f 1f0]
    ctm = get(state, :CTM, Matrix{Float32})
    ctm = cm*ctm
    state[:CTM] = ctm
    return state
end

@inline function evalContent!(pdo::PDPageElement{:Tm}, state::GState)
    a = get(pdo.operands[1])
    b = get(pdo.operands[2])
    c = get(pdo.operands[3])
    d = get(pdo.operands[4])
    e = get(pdo.operands[5])
    f = get(pdo.operands[6])
    tm  = [a b 0f0; c d 0f0; e f 1f0]
    tlm = copy(tm)
    state[:Tm]  = tm
    state[:Tlm] = tlm
    return state
end

@inline function evalContent!(pdo::PDPageElement{:Tf}, state::GState)
    src = get(state, :source, Union{PDPage, PDXObject})
    fontname = pdo.operands[1]
    font = get_font(src, fontname)
    font === CosNull && return state
    state[:font] = (fontname, font)
    fontsize = get(pdo.operands[2])
    # PDF Spec expects any number so better to standardize to Float32
    state[:fontsize] = Float32(fontsize)
    return state
end

# PDF Spec expects any number so better to standardize to Float32
for op in ["Tc", "Tw", "Tz", "TL", "Tr", "Ts"]
    @eval evalContent!(pdo::PDPageElement{Symbol($op)}, state::GState) =
        (state[Symbol($op)] = Float32(get(pdo.operands[1])); state)
end

@inline function set_text_pos!(tx, ty, state::GState)
    tmul = [1f0 0f0 0f0; 0f0 1f0 0f0; tx ty 1f0]
    #:TL may be called outside of BT...ET block
    tlm = get(state, :Tlm, Matrix{Float32}(I, 3, 3))
    tlm = tmul*tlm
    tm = copy(tlm)

    state[:Tm]  = tm
    state[:Tlm] = tlm
    return state
end

# Affects Tm leaves the Tlm intact
@inline function offset_text_pos!(tx, ty, state::GState)
    tmul = [1f0 0f0 0f0; 0f0 1f0 0f0; tx ty 1f0]
    #:TL may be called outside of BT...ET block
    tm = get(state, :Tm, Matrix{Float32}(I, 3, 3))
    tm = tmul*tm
    state[:Tm]  = tm
    return state
end

@inline function offset_text_leading!(state::GState)
    tl = get(state, :TL, Float32)
    return set_text_pos!(0f0, -tl, state)
end

@inline function evalContent!(pdo::PDPageElement{:TD}, state::GState)
    tx = Float32(get(pdo.operands[1]))
    ty = Float32(get(pdo.operands[2]))

    state[:TL] = -ty
    set_text_pos!(tx, ty, state)
end

@inline function evalContent!(pdo::PDPageElement{:Td}, state::GState)
    tx = Float32(get(pdo.operands[1]))
    ty = Float32(get(pdo.operands[2]))

    set_text_pos!(tx, ty, state)
end

evalContent!(pdo::PDPageElement{Symbol("T*")}, state::GState) =
    offset_text_leading!(state)

evalContent!(pdo::PDPageElement{Symbol("\'")}, state::GState) =
    offset_text_leading!(state)

@inline function evalContent!(pdo::PDPageElement{Symbol("\"")}, state::GState)#"
    state[:Tw] = Float32(get(pdo.operands[1]))
    state[:Tc] = Float32(get(pdo.operands[2]))
    offset_text_leading!(state)
end

function evalContent!(pdo::PDPageElement{:Do}, state::GState)
    xobjname = pdo.operands[1]
    src = get(state, :source, Union{PDPage, PDXObject})
    xobj = get_xobject(src, xobjname)
    xobj === CosNull && return state
    return Do(xobj, state)
end

evalContent!(pdo::PDPageInlineImage, state::GState) = state

evalContent!(pdo::CosObject, state::GState) = state


text:
module CommunalHelperPlayerBubbleRegion

using ..Ahorn, Maple

@mapdef Entity "CommunalHelper/PlayerBubbleRegion" PlayerBubbleRegion(
    x::Integer,
    y::Integer,
    width::Integer=24,
    height::Integer=24,
    nodes::Array{Tuple{Integer,Integer},1}=Tuple{Integer,Integer}[],
)

const placements = Ahorn.PlacementDict(
    "Player Bubble Region (Communal Helper)" => Ahorn.EntityPlacement(
        PlayerBubbleRegion,
        "rectangle",
        Dict{String,Any}(),
        function (entity)
            cx = Int(entity.data["x"]) + Int(entity.data["width"]) / 2
            cy = Int(entity.data["y"]) + Int(entity.data["height"]) / 2
            entity.data["nodes"] = [(cx + 32, cy), (cx + 64, cy)]
        end,
    ),
)

Ahorn.nodeLimits(entity::PlayerBubbleRegion) = 2, 2
Ahorn.resizable(entity::PlayerBubbleRegion) = true, true
Ahorn.minimumSize(entity::PlayerBubbleRegion) = 8, 8

const sprite = "characters/player/bubble"
const color = (45, 103, 111, 200) ./ 255

function Ahorn.selection(entity::PlayerBubbleRegion)
    nodes = entity.data["nodes"]
    controllX, controllY = Int.(nodes[1])
    endX, endY = Int.(nodes[2])

    return [
        Ahorn.getEntityRectangle(entity),
        Ahorn.getSpriteRectangle(sprite, controllX, controllY),
        Ahorn.getSpriteRectangle(sprite, endX, endY),
    ]
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::PlayerBubbleRegion, room::Maple.Room)
    width = Int(entity.data["width"])
    height = Int(entity.data["height"])

    Ahorn.drawRectangle(ctx, 0, 0, width, height, color)

    Ahorn.drawSprite(ctx, sprite, width / 2, height / 2)
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::PlayerBubbleRegion)
    px, py = Ahorn.position(entity)
    nodes = entity.data["nodes"]

    width = Int(entity.data["width"])
    height = Int(entity.data["height"])
    px += width / 2
    py += height / 2

    cx, cy = Int.(nodes[1])
    ex, ey = Int.(nodes[2])

    Ahorn.drawArrow(ctx, px, py, cx, cy, Ahorn.colors.selection_selected_fc, headLength=6)
    Ahorn.drawArrow(ctx, cx, cy, ex, ey, Ahorn.colors.selection_selected_fc, headLength=6)
    Ahorn.drawSprite(ctx, sprite, ex, ey)
end

end


text:
# Use like this:
#
# using LayeredLayouts
# layout(g) = layered_layout(solve_positions(Zarate(), g))
#
function layered_layout(pos)
  xs, ys, _ = pos
  return Point.(zip(xs, ys))
end


text:
module Minimal

using Plots

function test()
    plotly()
    p = scatter(randn(100),rand(100))
    display(p)
end

end # module


text:
using DemoPackageSSD
using Test

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


text:
module StokesDGTests

using Test
using Gridap
import Gridap: ∇
import LinearAlgebra: tr, ⋅

const T = VectorValue{2,Float64}

u(x) = VectorValue(x[1]*x[1], x[2])
∇u(x) = TensorValue(2*x[1],0.0,0.0,1.0)
Δu(x) = VectorValue(2.0,0.0)

p(x) = x[1] - x[2]
∇p(x) = VectorValue(1.0,-1.0)

f(x) = - Δu(x) + ∇p(x)
g(x) = tr(∇u(x))

∇(::typeof(u)) = ∇u
∇(::typeof(p)) = ∇p

L = 1.0
domain = (0.0, L, 0.0, L)
ncellx = 6
partition = (ncellx,ncellx)
model = CartesianDiscreteModel(domain,partition)
model = simplexify(model)

order = 2
const h = L / ncellx
const γ = order*(order+1)
const γ0 = 1.0/10.0

V = TestFESpace(
 model=model,
 order=order,
 reffe=:Lagrangian,
 valuetype=VectorValue{2,Float64},
 conformity=:L2)

Q = TestFESpace(
 model=model,
 order=order,
 reffe=:Lagrangian,
 valuetype=Float64,
 conformity=:L2,
 constraint=:zeromean)

U = TrialFESpace(V)
P = TrialFESpace(Q)

Y = MultiFieldFESpace([V,Q])
X = MultiFieldFESpace([U,P])

trian = get_triangulation(model)
degree = 2*order
quad = CellQuadrature(trian,degree)

btrian = BoundaryTriangulation(model)
bdegree = 2*order
bquad = CellQuadrature(btrian,bdegree)
const nb = get_normal_vector(btrian)

strian = SkeletonTriangulation(model)
sdegree = 2*order
squad = CellQuadrature(strian,sdegree)
const ns = get_normal_vector(strian)

function A_Ω(x,y)
  u, p = x
  v, q = y
  ∇(v)⊙∇(u) - ∇(q)⋅u + v⋅∇(p)
end

function B_Ω(y)
  v, q = y
  v⋅f + q*g
end

function A_∂Ω(x,y)
  u, p = x
  v, q = y
  (γ/h)*v⋅u - v⋅(nb⋅∇(u)) - (nb⋅∇(v))⋅u + 2*(q*nb)⋅u
end

function B_∂Ω(y)
  v, q = y
  (γ/h)*v⋅u - (nb⋅∇(v))⋅u + (q*nb)⋅u
end

function A_Γ(x,y)
  u, p = x
  v, q = y
  (γ/h)*inner( jump(outer(v,ns)), jump(outer(u,ns))) -
    inner( jump(outer(v,ns)), mean(∇(u)) ) -
    inner( mean(∇(v)), jump(outer(u,ns)) ) +
    (γ0*h)*jump(q*ns)⋅jump(p*ns) +
    jump(q*ns)⋅mean(u) -
    mean(v)⋅jump(p*ns)
end

t_Ω = AffineFETerm(A_Ω,B_Ω,trian,quad)
t_∂Ω = AffineFETerm(A_∂Ω,B_∂Ω,btrian,bquad)
t_Γ = LinearFETerm(A_Γ,strian,squad)

op = AffineFEOperator(X,Y,t_Ω,t_∂Ω,t_Γ)

uh, ph = solve(op)

eu = u - uh
ep = p - ph

l2(v) = v⋅v
h1(v) = v⋅v + inner(∇(v),∇(v))

eu_l2 = sqrt(sum(integrate(l2(eu),trian,quad)))
eu_h1 = sqrt(sum(integrate(h1(eu),trian,quad)))
ep_l2 = sqrt(sum(integrate(l2(ep),trian,quad)))

tol = 1.0e-9
@test eu_l2 < tol
@test eu_h1 < tol
@test ep_l2 < tol

end # module


text:
using MCMC_LR_Tests
using Base.Test
using Distributions
using HypothesisTests

# consistent testing
srand(UInt32[0x7f712df6, 0xc9c492e6, 0x4935b672, 0x946d94f9])

@testset "IID normal mean" begin
    N = 10000
    dist = MvNormal(Float64[1, -1, 2], [1,2,3])
    ps = [mean_LR_pvalue(mean(dist), rand(dist, N)') for _ in 1:100]
    t = ExactOneSampleKSTest(ps, Uniform(0,1))
    @test pvalue(t) > 0.05
end

@testset "IID normal cov" begin
    N = 10000
    A = [0.1 0.1 0.0;
         0.0 0.5 0.1;
         0.0 0.0 0.7]
    Σ = A'*A
    dist = MvNormal(Float64[1, -1, 2], Σ)
    @test quantile([cov_LR_pvalue(Σ, rand(dist, N)') for _ in 1:100], 0.1) ≥ 0.1
end

@testset "IID normal cov 2" begin
    Σ = [6.366588528827408 3.25842202827805 0.9729856844746686 3.858746597584623;
         3.25842202827805 3.938997699450515 -2.809133292920283 1.5566006096285196;
         0.9729856844746686 -2.809133292920283 7.344095307429978 -1.2388040483400584;
         3.858746597584623 1.5566006096285196 -1.2388040483400584 9.227552651298026]
    dist = MvNormal(zeros(4), Σ)
    @test quantile([cov_LR_pvalue(Σ, rand(dist, N)') for _ in 1:100], 0.1) ≥ 0.1
end

"Simulate `N` draws of a vector AR(1) process `x = A x + ϵ`, where `ϵ ∼ MvNormal(0, Σ)`."
function simulate_VAR1(A, Σ, N)
    K = Base.LinAlg.checksquare(Σ)
    @assert K == Base.LinAlg.checksquare(A)
    @assert maximum(abs.(eigvals(A))) < 1
    ϵ = MvNormal(zeros(K), Σ)
    x = (I - A) \ rand(ϵ)
    X = similar(x, (N, K))
    for i in 1:N
        x .= A*x .+ rand(ϵ)
        X[i, :] .= x
    end
    X
end

@testset "VAR(1)" begin
    A = [0.1 0.1 0.0;
         0.0 0.5 0.1;
         0.0 0.0 0.7]
    Σ = A'*A
    μ = zeros(3)
    N = 10000
    p10 = quantile([mean_LR_pvalue(μ, simulate_VAR1(A, Σ, N)) for _ in 1:100], 0.1)
    @test p10 ≥ 0.10
end


text:
module Aircraft

using LinearAlgebra
using StaticArrays, ComponentArrays
using UnPack

using Flight.Modeling
import Flight.Modeling: init_x0, init_y0, init_u0, init_d0,f_cont!, f_disc!

using Flight.Plotting
import Flight.Plotting: plots

using Flight.Terrain
using Flight.Atmosphere
using Flight.Airdata

using Flight.Kinematics
import Flight.Kinematics: init!
using Flight.Dynamics

using Flight.Components: get_mass_properties, get_wr_b, get_hr_b

export AbstractAircraftID, GenericAircraftID, AircraftBase
export assign_joystick_inputs!

######################### AircraftBase ########################

abstract type AbstractAircraftID end

struct GenericAircraftID <: AbstractAircraftID end

Base.@kwdef struct AircraftBase{I <: AbstractAircraftID,
                    K <: AbstractKinematics,
                    F <: SystemDescriptor,
                    C <: SystemDescriptor} <: SystemDescriptor

    identifier::I = GenericAircraftID()
    kinematics::K = KinLTF()
    airframe::F = RigidBody()
    controls::C = NullSystemDescriptor()
end

init_x0(ac::AircraftBase) = ComponentVector(
    kin = init_x0(ac.kinematics),
    afr = init_x0(ac.airframe),
    ctl = init_x0(ac.controls),
    )

init_u0(ac::AircraftBase) = init_u0(ac.controls)

init_y0(ac::AircraftBase) = (
    kin = KinData(),
    dyn = DynData(),
    air = AirData(),
    afr = init_y0(ac.airframe),
    ctl = init_y0(ac.controls),
    )

init_d0(ac::AircraftBase) = (
    afr = init_d0(ac.airframe),
    ctl = init_d0(ac.controls),
    )

const AircraftBaseSys{I,K,F,C} = System{AircraftBase{I,K,F,C}} where {I,K,F,C}

function System(ac::AircraftBase, ẋ = init_x0(ac), x = init_x0(ac),
                    y = init_y0(ac), u = init_u0(ac), d = init_d0(ac), t = Ref(0.0))

    params = ()
    subsystems = (
        afr = System(ac.airframe, ẋ.afr, x.afr, y.afr, init_u0(ac.airframe), d.afr, t),
        ctl = System(ac.controls, ẋ.ctl, x.ctl, y.ctl, u, d.ctl, t),)

    System{map(typeof, (ac, x, y, u, d, params, subsystems))...}(
                         ẋ, x, y, u, d, t, params, subsystems)
end

init!(ac::System{<:AircraftBase}, init::KinInit) = init!(ac.x.kin, init)

function f_cont!(sys::AircraftBaseSys, trn::AbstractTerrain, atm::AtmosphericSystem)

    @unpack ẋ, x, subsystems = sys
    @unpack afr, ctl = subsystems

    #update kinematics
    kin = f_kin!(ẋ.kin.pos, x.kin)

    air = AirData(kin, atm)

    #update controls
    f_cont!(ctl, afr, kin, air, trn)
    #update airframe components
    f_cont!(afr, ctl, kin, air, trn)

    mass = get_mass_properties(afr)
    wr_ext_b = get_wr_b(afr)
    hr_b = get_hr_b(afr)

    # update dynamics
    dyn = f_dyn!(ẋ.kin.vel, kin, mass, wr_ext_b, hr_b)

    sys.y = (kin = kin, dyn = dyn, air = air, afr = afr.y, ctl = ctl.y)
    return nothing

end

function f_disc!(sys::AircraftBaseSys)
    @unpack afr, ctl = sys.subsystems

    x_mod = renormalize!(sys.x.kin, 1e-8) |
            f_disc!(afr, ctl) |
            f_disc!(ctl, afr)

    return x_mod
end

assign_joystick_inputs!(args...) = throw(MethodError(assign_joystick_inputs!, args))

#no custom plots function required, since all outputs are namedtuples

end #module


text:
import Random
import Base.getindex, Base.size, Base.eltype

mutable struct SumTree{T}
    capacity::Int64
    tree_struct::Array{Float64, 1}
    data::Array{T, 1}
    write_pos::Int64
    size::Int64
    tree_len::Int64
    SumTree{T}(capacity) where T = new{T}(
        capacity,
        zeros(2*capacity - 1),
        Array{T, 1}(undef, capacity),
        0,
        0,
        2*capacity - 1)
end

total(tree::SumTree) = tree.tree_struct[1]
size(tree::SumTree) = tree.size
capacity(tree::SumTree) = tree.capacity
eltype(tree::SumTree) = typeof(tree.data[1])
randInRange(a, b, size=1; rng=Random.GLOBAL_RNG) = (rand(rng, length(a), size).*(b-a)).+a

function update!(tree::SumTree, idx, weight)
    change = weight - tree.tree_struct[idx + 1]
    tree.tree_struct[idx + 1] = weight
    propagate!(tree, idx, change)
end

function propagate!(tree::SumTree, idx, change)
    parent = div((idx - 1), 2)
    tree.tree_struct[parent + 1] += change
    if parent != 0
        propagate!(tree, parent, change)
    end
end

function add!(tree::SumTree, weight, data)
    idx = tree.write_pos + tree.capacity - 1
    tree.data[tree.write_pos + 1] = data
    update!(tree, idx, weight)

    tree.write_pos += 1
    if tree.write_pos >= tree.capacity
        tree.write_pos = 0
    end

    if tree.size < tree.capacity
        tree.size += 1
    end
end

function get(tree::SumTree, s)
    idx = retrieve(tree, 0, s)
    data_idx = idx - tree.capacity + 1

    return idx + 1, tree.tree_struct[idx + 1], tree.data[data_idx + 1]
end

function getindex(tree::SumTree, idx::Int64)
    tree_idx = (idx - 1) + tree.capacity - 1
    return tree.tree_struct[tree_idx + 1], tree.data[idx]
end

function retrieve(tree::SumTree, idx, s)
    left = 2 * idx + 1
    if left >= tree.tree_len
        return idx
    end
    if s <= tree.tree_struct[left + 1]
        return retrieve(tree, left, s)
    else
        right = left + 1
        return retrieve(tree, right, s-tree.tree_struct[left + 1])
    end
end

function sample(tree::SumTree, batch_size; rng=Random.GLOBAL_RNG)
    batch_idx = zeros(Int64, batch_size)
    batch_priorities = zeros(Float64, batch_size)
    batch = Array{eltype(tree), 1}(undef, batch_size)
    segment = total(tree) / batch_size

    a = convert(Array, 0:(batch_size-1))*segment
    b = convert(Array, 1:batch_size)*segment
    s = randInRange(a, b; rng=rng)
    for i in 1:batch_size
        (x,y,z) = get(tree, s[i])
        batch_idx[i] = x
        batch_priorities[i] = y
        batch[i] = z
    end

    return batch_idx, batch_priorities, batch

end


text:
using Bosonic
using Test
@testset "Bosonic.jl" begin
    @test round(st(State_fixed(ones(10),Op_fixed(4,2)))'*st(State_fixed(ones(10),Op_fixed(4,2))),digits=4) == 1
    @test nume(State_fixed(ones(10),Op_fixed(4,2))) == 2
    @test basis(State_fixed(ones(10),Op_fixed(4,2)))[1,:][4] == 2
    @test rhosp(State_fixed(ones(10),Op_fixed(4,2)))[3,1] == 0.4
end


text:
client = connect(8080)

while true
	print("> ")
	input = readline(STDIN)
	input = replace(input, "\n", "")
	println(client, input)
	if input == ":exit" || input == ":quit"
		close(client)
		exit(1)
	else
		write(STDOUT,readline(client))
	end
end

max_stars_repo_name: openwonk/echo

max_issues_repo_name: openwonk/echo

max_forks_repo_name: openwonk/echo


text:
# Example of logistic map from Keith Briggs's 2006 paper

using ExactReals

"Iterate the function `f` for `n` steps, starting at `x0`."
function iterates(f, n, x0)

    x = x0

    for i in 1:n
        x = f(x)
    end

    return x
end

"Logistic map"
logistic(a, x) = a * x * (1 - x)
logistic(a) = x -> logistic(a, x)

a = rationalize(3.999)
x0 = rationalize(0.9)
n = 54

y = iterates(logistic(ExactReal(a)), n, ExactReal(x0))

@show y > ExactReal(1//2)

# Floating point:

z = iterates(logistic(float(a)), n, float(x0))

@show z > 0.5

nothing

#
# using IntervalArithmetic
#
# x = interval(big"0.9")
#
# iterates(logistic(@biginterval(3.999)), 53, @biginterval(0.9)) |> showfull
#
# iterates(logistic(I"3.999"), 53, I"0.9")
#
# showfull(I"3.999")
#
# # example where double precision is on wrong side of 1/2:
# iterates(logistic(3.999), 54, 0.9)
#
# iterates(logistic(@biginterval(3.999)), 54, @biginterval(0.9)) |> showfull
#
#
# @time y = iterates(logistic(ExactReal(3999//1000)), 54, ExactReal(9//10))
#
# y > ExactReal(1//2)


text:
using Test
using Dates
using ExtensibleScheduler
using TimeFrames

@testset "CustomTrigger" begin
    f = (dt_previous_fire, dt_now) -> dt_now + Dates.Minute(5)
    trigger = Trigger(f)

    n_executed = 0

    function sample_action()
        n_executed += 1
    end

    dt_start = DateTime(2010, 1, 4, 0, 0)
    dt = dt_start
    clock = SimClock(dt)
    sched = BlockingScheduler(clock=clock, delayfunc=NoSleep)
    action = Action(sample_action)

    add(sched, action, trigger)
    run_pending(sched)
    @test n_executed == 0

    for i in 1:13
        dt += Dates.Minute(1)
        #println("Clock set to $dt")
        set(clock, dt)
        run_pending(sched)
    end

    @test n_executed == 2

end


text:
using MemorylessNonlinearities
using Test

x = [-3,-2,-1,0,1,2,3]

@testset "MemorylessNonlinearities.jl" begin

    @info "Test MemorylessNonlinearities"

    @test filt(Arctangent(100000), x) ≈ [-1,-1,-1,0,1,1,1] atol=0.0001
    @test filt(Arctangent(10), x) ≈ [-0.9788,-0.9682,-0.9365,0,0.9365,0.9682,0.9788] atol=0.0001

    @test filt(Blanking(1.0), x) == [0,0,-1,0,1,0,0]

    @test filt(CauchyNL(1.0), x) == [-0.6,-0.8,-1.0,0.0,1.0,0.8,0.6]

    @test filt(Clipping(1.0), x) == [-1,-1,-1,0,1,1,1]

    @test filt(HampelThreePart(1.0,2.0,3.0), x) == [0.0,-1.0,-1.0,0.0,1.0,1.0,0.0]

    @test filt(SαSNL(2.0, 1.0, 0.0, false), x) ≈ x ./ 2 atol=0.0001
    @test filt(SαSNL(1.0, 1.0, 0.0, false), x) ≈ [-0.6,-0.8,-1.0,0.0,1.0,0.8,0.6] atol=0.0001
    @test filt(SαSNL(2.0, 1.0, 0.0, true), x) ≈ x ./ 2 atol=0.0001
    @test filt(SαSNL(1.0, 1.0, 0.0, true), x) ≈ [-0.6,-0.8,-1.0,0.0,1.0,0.8,0.6] atol=0.0001
    @test filt(SαSNL(2.0, 1.0, 0.0), x) ≈ x ./ 2 atol=0.0001
    @test filt(SαSNL(1.0, 1.0, 0.0), x) ≈ [-0.6,-0.8,-1.0,0.0,1.0,0.8,0.6] atol=0.0001
    αs = [1.000001, 1.001000, 1.002000, 1.123457, 1.234567, 1.345678, 1.456789,
          1.567899, 1.678999, 1.789999, 1.899999, 1.998000, 1.999000, 1.999999]
    xtmp = [-999.987654, -500.505050, -100.101010, -30.303030, -10.101010, -5.505050, -3.303030, -1.101010, -0.010101,
            0.0, 0.010101, 1.101010, 3.303030, 5.505050, 10.101010, 30.303030, 100.101010, 500.505050, 999.987654]
    for α in αs
        @test filt(SαSNL(α, 1.0, 0.0, true), xtmp) ≈ filt(SαSNL(α, 1.0, 0.0, false), xtmp) rtol=0.01
        @test filt(SαSNL(α, 2.0, 0.0), xtmp) ≈ filt(SαSNL(α, 2.0, 0.0, false), xtmp) rtol=0.01
    end

    @test filt(SoftClipping(2.0), x) ≈ [-2.0,-2.0,-1.375,0.0,1.375,2.0,2.0]

    @test filt(TurkeyBiweight(3.0), x) ≈ [-0.0,-0.6173,-0.7901,0.0,0.7901,0.6173,0.0] atol=0.0001

    @test filt(InverseHyperbolicSine(1.0), x) ≈ [-1.8184,-1.4436,-0.8814,0.0,0.8814,1.4436,1.8184] atol=0.0001
end

@testset "utils.jl" begin

    @info "Test utils"

    x = randn(1000)
    xrescale = minmaxrescale(x, -0.1, 0.3)
    @test minimum(xrescale) ≈ -0.1 atol=1e-6
    @test maximum(xrescale) ≈ 0.3 atol=1e-6

end


text:
export Connect

function Connect(model;pre=[],post=[],condition="true",expr="",parallel=false)
  var=match(r"^.*(?==)",expr).match
  value=match(r"(?<==).*",expr).match

  if !contains(model.connection,var)
    model.connection*="@parallel for p in workers()  @eval $var=zeros(N,N_tot) end\n"
  end

  for preTemp in [pre;],postTemp in [post;]
    j=model.groups[preTemp][:localindex]
    processor=model.groups[postTemp][:processor]
    condition=replace(condition,"j","(j-$(j[1]-1))")
    value=replace(value,"j","(j-$(j[1]-1))")
    model.connection*="""
    @spawnat $(processor+1) begin
      for j=$j,i=1:N
         if $condition
           $var[i,j]=$value
         end
      end
   end
   """
 end

end # end of function Connect