text:
using Reexport

using MLStyle
using Distributions

# In Bijectors.jl,
# logpdf_with_trans(dist, x, true) == logpdf(transformed(dist), link(dist, x))

export xform

function xform(m::JointDistribution{A, B}, _data) where {A,B}
    return _xform(m.model, m.args, _data)
end

@gg function _xform(_m::Model{Asub,B}, _args::A, _data) where {Asub, A,B}
    type2model(_m) |> sourceXform(_data) |> loadvals(_args, _data)
end

# function xform(m::Model{EmptyNTtype, B}) where {B}
#     return xform(m,NamedTuple())
# end

export sourceXform

function sourceXform(_data=NamedTuple())
    function(_m::Model)

        _m = canonical(_m)

        _datakeys = getntkeys(_data)
        proc(_m, st::Assign)        = :($(st.x) = $(st.rhs))
        proc(_m, st::Return)     = nothing
        proc(_m, st::LineNumber) = nothing

        function proc(_m, st::Sample)
            if st.x ∈ _datakeys
                return :($(st.x) = _data.$(st.x))
            else
                return (@q begin
                    $(st.x) = rand($(st.rhs))
                    _t = xform($(st.rhs), _data)

                    _result = merge(_result, ($(st.x)=_t,))
                end)
            end

        end

        wrap(kernel) = @q begin
            _result = NamedTuple()
            $kernel
            as(_result)
        end

        buildSource(_m, proc, wrap) |> flatten

    end
end

function xform(d, _data)
    if hasmethod(support, (typeof(d),))
        return asTransform(support(d))
    end

    error("Not implemented:\nxform($d)")
end

using TransformVariables: ShiftedExp, ScaledShiftedLogistic

function asTransform(supp:: RealInterval)
    (lb, ub) = (supp.lb, supp.ub)

    (lb, ub) == (-Inf, Inf) && (return asℝ)
    isinf(ub) && return ShiftedExp(true,lb)
    isinf(lb) && return ShiftedExp(false,lb)
    return ScaledShiftedLogistic(ub-lb, lb)
end

# export xform
# xform(::Normal)       = asℝ
# xform(::Cauchy)       = asℝ
# xform(::Flat)         = asℝ

# xform(::HalfCauchy)   = asℝ₊
# xform(::HalfNormal)   = asℝ₊
# xform(::HalfFlat)     = asℝ₊
# xform(::InverseGamma) = asℝ₊
# xform(::Gamma)        = asℝ₊
# xform(::Exponential)  = asℝ₊

# xform(::Beta)         = as𝕀
# xform(::Uniform)      = as𝕀

function xform(d::For, _data)
    xf1 = xform(d.f(getindex.(d.θ, 1)...), _data)
    return as(Array, xf1, length.(d.θ)...)

    # TODO: Implement case of unequal supports
end

function xform(d::iid, _data)
    as(Array, xform(d.dist, _data), d.size...)
end

xform(d::MvNormal) =  as(Vector, length(d))


text:
function f!(fvec, x)
    fvec[1] = 1 - x[1]
    fvec[2] = 10(x[2]-x[1]^2)
end
function j!(fjac, x)
    fjac[1,1] = -1
    fjac[1,2] = 0
    fjac[2,1] = -20x[1]
    fjac[2,2] = 10
end

function nlsolvers_f(x, Fx, Jx)
	if Fx !== nothing
		f!(Fx, x)
	end
	if Jx !== nothing
		j!(Jx, x)
	end
	objective_return(Fx, Jx)
end
x0 = [-1.2, 1.0]

nlsolve!(OnceDiffed(nlsolvers_f), copy(x0), TrustRegion(NLSolvers.Newton()))
nlsolve!(OnceDiffed(nlsolvers_f), copy(x0), TrustRegion(NLSolvers.Newton(), Dogleg()))
nlsolve!(OnceDiffed(nlsolvers_f), copy(x0), TrustRegion(NLSolvers.Newton(), NTR()))
nlsolve!(OnceDiffed(nlsolvers_f), copy(x0), TrustRegion(NLSolvers.Newton(), NWI()))
# initial convergence
nlsolve!(OnceDiffed(nlsolvers_f), fill(1.0, 2), TrustRegion(NLSolvers.Newton()))
nlsolve!(OnceDiffed(nlsolvers_f), fill(1.0, 2), TrustRegion(NLSolvers.Newton(), Dogleg()))
nlsolve!(OnceDiffed(nlsolvers_f), fill(1.0, 2), TrustRegion(NLSolvers.Newton(), NTR()))
nlsolve!(OnceDiffed(nlsolvers_f), fill(1.0, 2), TrustRegion(NLSolvers.Newton(), NWI()))

ext: jl


text:
struct Mesh
    vertices::Array{Point{3, Float32}}
    edges::Set{Tuple{Int, Int}}
    faces::Array{Face{3,Any}}
end

"""
    distance(v1::Point, v2::Point)

Calculate the distance between two points

"""
function distance(v1::Point, v2::Point)
    return sqrt(reduce(+, map((x,y)->(y-x)^2, v1, v2)))
end

function select_valid_pairs(mesh::Mesh, threshold::Real)
    rs = Set{Tuple{Int, Int}}()
    for i in 1:length(mesh.vertices), e in mesh.edges
        if i in e
            push!(rs, e)
        end
    end
    for i in 1:length(mesh.vertices), j in 1:length(mesh.vertices)
        if i != j && distance(mesh.vertices[i], mesh.vertices[j]) <= threshold
            if !((j,i) in rs)
                push!(rs, (i, j))
            end
        end
    end
    return rs
end

function calculate_contraction_position(v1::Point, v2::Point)
    return Point{3}((v1.+v2)./2)
end

function calculate_contraction_cost(mesh::HomogenousMesh, v1::Point, v2::Point)
    v_bar = contract_vertices(v1, v2)
    cost = Q(mesh, v1)+Q(mesh, v2)
    return cost
end

function get_edges_from_face(face::Face)
    return (face[1], face[2]), (face[1], face[3]), (face[2], face[3])
end

function get_faces_from_vertex(mesh::HomogenousMesh, v_ind::Int)
    return [f for f in mesh.faces if v_ind in f]
end

function Q(mesh::HomogenousMesh, v_ind::Int)
    neighboring_faces = get_faces_from_vertex(mesh, v_ind)
    plane_reps = [get_plane_representation(mesh, f) for f in neighboring_faces]
    error_quadrics = [fundamental_error_quadric(p) for p in plane_reps]
    return reduce(+, error_quadrics)
end

function normalize(n::Vector)
    sq_sum = reduce(+, map(x->x^2, n))
    sq_sum == 0 && error("Can't normalize the 0-vector")
    return map(x->x/(sqrt(sq_sum)), n)
end

function get_plane_representation(mesh::HomogenousMesh, face::Face)
    points_in_plane = mesh.vertices[face]
    v1 = Vector(points_in_plane[2].-points_in_plane[1])
    v2 = Vector(points_in_plane[3].-points_in_plane[1])
    n = v1 × v2
    n = normalize(n)
    d = -points_in_plane[1][1]*n[1] - points_in_plane[1][2]*n[2] - points_in_plane[1][3]*n[3]
    append!(n,d)
    return SVector{4}(n)
end

function convert_from_homogenous(mesh::HomogenousMesh)
    edges =Set{Tuple{Int, Int}}()
    for f in mesh.faces
        edges = union(edges, get_edges_from_face(f))
    end
    return Mesh(mesh.vertices, edges, mesh.faces)
end

function fundamental_error_quadric(p::SVector{4,T}) where {T<:Real}
    return [p[1]^2 p[1]*p[2] p[1]*p[3] p[1]*p[4];
            p[1]*p[2] p[2]^2 p[2]*p[3] p[2]*p[4];
            p[1]*p[3] p[2]*p[3] p[3]^2 p[3]*p[4];
            p[1]*p[4] p[2]*p[4] p[3]*p[4] p[4]^2]
end


text:
using FFTW,LinearAlgebra,Plots,SpecialFunctions;N=128;P=32N;dt=1/6N;T=1024;
W=400;w=W/P*N;ξ=im*zeros(N);x=rand(P);v=rand([-1.0,1.0],P);r=zeros(N);
E=zeros(N);X=copy(x);V=copy(v);ik=2π*im*vcat(1,1:N/2,-N/2+1:-1);D=zeros(T,4);
f(g,c)=erf((g-c)*N)/2;ff(i,c)=f((i+0.5)/N,c)-f((i-0.5)/N,c);l=1e-8;
d(c)=((mod1(i,N),ff(i,c)) for i∈(-6:6).+Int(round(c*N)));F=copy(E);
ρ(x,y)=(r.*=0;for j∈d.((x.+y)./2);for k∈j;r[k[1]]+=k[2]*w;end;end;r);
@gif for t ∈1:T;X.=x;V.=v;F.*=NaN;for _∈0:9;isapprox(F,E,rtol=l)&&break;F.=E;
  x.=X.+(v.+V)/2*dt;E.=real.(ifft((ξ.=fft(ρ(x,X))./ik;ξ[1]*=0;ξ)));
  for j∈1:P;v[j]=V[j]+sum(k->E[k[1]]*k[2],d((x[j]+X[j])/2))*dt;end;end;x.=mod.(x,1);
  D[t,1:2].=(sum(E.^2)/N,sum(v.^2)*W/P)./2;D[t,3:4].=sum.((D[t,1:2],v/P));
  scatter(x,v,ylims=(-3,3));D[t,1:3].*=2/W;plot!((1:t)./T,D[1:t,:],legend=0==1)
end every 8


text:
"""
    gpfitbayes(..., niter::Int=5000, warmup::Int=2000)

Generate a sample from the GP parameters' posterior distribution.

Data provided must be the exceedances above the threshold, *i.e.* the data above the threshold minus
the threshold.

# Arguments
- `niter::Int = 5000`: The total number of MCMC iterations
- `warmup::Int = 2000`: The number of warmup iterations (burn-in).

# Implementation

The function uses the No-U-Turn Sampler (NUTS; [Hoffman and Gelman, 2014](http://jmlr.org/papers/v15/hoffman14a.html))
implemented in the [Mamba.jl](https://mambajl.readthedocs.io/en/latest/index.html)
package to generate a random sample from the posterior distribution.

Currently, only the improper uniform prior is implemented, *i.e.*
```math
f_{(β₂,β₃)}(β₂,β₃) ∝ 1,
```
where
```math
ϕ = X₂ × β₂,
```
```math
ξ = X₃ × β₃.
```
In the stationary case, this improper prior yields to a proper posterior if the
sample size is larger than 2 ([Northrop and Attalides, 2016](https://www.jstor.org/stable/24721296?seq=1)).

The covariates are standardized before estimating the parameters to help fit the
 model. They are transformed back on their original scales before returning the
 fitted model.

See also [`gpfitbayes`](@ref) for the other methods, [`gpfitpwm`](@ref), [`gpfit`](@ref) and [`ThresholdExceedance`](@ref).

# Reference

Hoffman M. D. and Gelman A. (2014). The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. *Journal of Machine Learning Research*, 15:1593–1623.

Paul J. Northrop P. J. and Attalides N. (2016). Posterior propriety in Bayesian extreme value analyses using reference priors. *Statistica Sinica*, 26:721-743.

"""
function gpfitbayes end

"""
    gpfitbayes(y,
        logscalecov = Vector{Variable}(),
        shapecov = Vector{Variable}(),
        niter::Int=5000,
        warmup::Int=2000
        )

Generate a sample from the GP parameters' posterior distribution.

# Arguments

- `y::Vector{<:Real}`: The vector of exceedances.
- `logscalecov::Vector{<:DataItem} = Vector{Variable}()`: The covariates of the log-scale parameter.
- `shapecov::Vector{<:DataItem} = Vector{Variable}()`: The covariates of the shape parameter.
"""
function gpfitbayes(y::Vector{<:Real};
     logscalecov::Vector{<:DataItem} = Vector{Variable}(),
     shapecov::Vector{<:DataItem} = Vector{Variable}(),
     niter::Int=5000, warmup::Int=2000)::BayesianEVA

     logscalecovstd = standardize.(logscalecov)
     shapecovstd = standardize.(shapecov)

    model = ThresholdExceedance(Variable("y", y), logscalecov = logscalecovstd, shapecov = shapecovstd)

    fittedmodel = fitbayes(model, niter=niter, warmup=warmup)

    return transform(fittedmodel)

end

"""
    gpfitbayes(df::DataFrame,
        datacol::Symbol,
        logscalecovid = Vector{Symbol}(),
        shapecovid = Vector{Symbol}(),
        niter::Int=5000,
        warmup::Int=2000)

Generate a sample from the GP parameters' posterior distribution.

# Arguments
- `df::DataFrame`: The dataframe containing the data.
- `datacol::Symbol`: The symbol of the column of `df` containing the exceedances.
- `logscalecovid::Vector{Symbol} = Vector{Symbol}()`: The symbols of the columns of `df` containing the covariates of the log-scale parameter.
- `shapecovid::Vector{Symbol} = Vector{Symbol}()`: The symbols of the columns of `df` containing the covariates of the shape parameter.
"""
function gpfitbayes(df::DataFrame, datacol::Symbol;
    logscalecovid::Vector{Symbol}=Symbol[],
    shapecovid::Vector{Symbol}=Symbol[],
    niter::Int=5000, warmup::Int=2000)::BayesianEVA

    logscalecovstd = standardize.(buildVariables(df, logscalecovid))
    shapecovstd = standardize.(buildVariables(df, shapecovid))

    model = ThresholdExceedance(Variable(string(datacol), df[:, datacol]), logscalecov = logscalecovstd, shapecov = shapecovstd)

    fittedmodel = fitbayes(model, niter=niter, warmup=warmup)

    return transform(fittedmodel)

    return fm

end

"""
    gpfitbayes(model::ThresholdExceedance;
        niter::Int=5000,
        warmup::Int=2000)

Generate a sample from the GP parameters' posterior distribution.

"""
function gpfitbayes(model::ThresholdExceedance; niter::Int=5000, warmup::Int=2000)::BayesianEVA

    return fitbayes(model, niter=niter, warmup=warmup)

end


text:
macro except_str(expr, err_type)
    return quote
        let err = nothing
            try
                $(esc(expr))
            catch err
            end
            err === nothing && error("expected failure, but no exception thrown")
            @test typeof(err) === $(esc(err_type))
            buf = IOBuffer()
            showerror(buf, err)
            String(take!(buf))
        end
    end
end

@testset "error hints" begin
    @testset "Math" begin
        gray = Gray(0.8)
        rgb = RGB{Float32}(1, 0, 0)
        err_str = @except_str gray + rgb MethodError
        @test occursin("no method matching +(::Gray{Float64}, ::RGB{Float32})", err_str)
        @test occursin("Math on colors is deliberately undefined in ColorTypes, but see the ColorVectorSpace package", err_str)

        err_str = @except_str gray * rgb MethodError
        @test occursin("no method matching *(::Gray{Float64}, ::RGB{Float32})", err_str)
        @test occursin("Math on colors is deliberately undefined in ColorTypes, but see the ColorVectorSpace package", err_str)
        @test occursin("You may also need `⋅`, `⊙`, or `⊗`.", err_str)

        err_str = @except_str rgb * 0.5 MethodError
        @test !occursin("You may also need `⋅`, `⊙`, or `⊗`.", err_str)
    end
end


text:
_call_if_something(f, arg, dflt = nothing) = isnothing(arg) ? dflt : f(arg)

# currbook
bookdir() = _call_if_something(bookdir, currbook(), "")

# currdoc
get_doi() = _call_if_something(get_doi, currdoc(), "")
get_year() = _call_if_something(get_year, currdoc(), "")
get_author() = _call_if_something(get_author, currdoc(), "")
get_bibkey() = _call_if_something(get_bibkey, currdoc(), "")

# currobj
get_title() = _call_if_something(get_title, currobj(), "")
get_tags() = _call_if_something(get_tags, currobj(), nothing)
localpath() = _call_if_something(localpath, currobj(), "")
localrelpath() = _call_if_something(localrelpath, currobj(), "")

# crossrefs
crossrefs(; force = false) = crossrefs(get_doi(); force)


text:
"""
     EmbeddingLags

Stores the information on which embedding lags were used to construct an
embedding.

## Fields
- **`lags`**: The lags which were used to construct the embedding. The i-th
    lag gives the lag for the i-th embedding variable.
"""
struct EmbeddingLags
    lags::Vector{Int}
end

export EmbeddingLags


text:
module CIFAR100_Tests
using Test
using ColorTypes
using ImageCore
using FixedPointNumbers
using MLDatasets
using DataDeps

function readimage(path, m, i)
    CIFAR100.Reader.readdata(path, m, i)[1]
end

@testset "Constants" begin
    @test CIFAR100.Reader.NROW === 32
    @test CIFAR100.Reader.NCOL === 32
    @test CIFAR100.Reader.NCHAN === 3
    @test CIFAR100.Reader.NBYTE === 3074

    @test CIFAR100.TRAINSET_SIZE === 50_000
    @test CIFAR100.TESTSET_SIZE  === 10_000

    @test DataDeps.registry["CIFAR100"] isa DataDeps.DataDep
end

@testset "convert2images" begin
    @test CIFAR100.convert2image === CIFAR10.convert2image
end

# NOT executed on CI. only executed locally.
# This involves dataset download etc.
if parse(Bool, get(ENV, "CI", "false"))
    @info "CI detected: skipping dataset download"
else
    data_dir = withenv("DATADEPS_ALWAY_ACCEPT"=>"true") do
        datadep"CIFAR100"
    end

    @testset "classnames" begin
        @test CIFAR100.classnames_coarse() isa Vector{String}
        @test length(CIFAR100.classnames_coarse()) == 20
        @test length(unique(CIFAR100.classnames_coarse())) == 20
        @test CIFAR100.classnames_fine() isa Vector{String}
        @test length(CIFAR100.classnames_fine()) == 100
        @test length(unique(CIFAR100.classnames_fine())) == 100
    end

    @testset "Images" begin
        # Sanity check that the first trainimage is not the first testimage
        @test CIFAR100.traintensor(1) != CIFAR100.testtensor(1)
        # Make sure other integer types work as indicies
        @test CIFAR100.traintensor(0xBAE) == CIFAR100.traintensor(2990)

        @testset "Test that traintensor are the train images" begin
            path = joinpath(data_dir, "cifar-100-binary", "train.bin")
            @test CIFAR100.traintensor(1) == reinterpret(N0f8, readimage(path, 10_000, 1))
            @test CIFAR100.traintensor(Float64, 1) ≈ readimage(path, 10_000, 1) ./ 255
            @test CIFAR100.traintensor(Int, 1) == Int.(readimage(path, 10_000, 1))
            @test CIFAR100.traintensor(UInt8, 1) == readimage(path, 10_000, 1)
            @test CIFAR100.traintensor(UInt8, 2) == readimage(path, 10_000, 2)
            @test CIFAR100.traintensor(UInt8, 10_000) == readimage(path, 10_000, 10_000)
        end

        @testset "Test that testtensor are the test images" begin
            path = joinpath(data_dir, "cifar-100-binary", "test.bin")
            @test CIFAR100.testtensor(1) == reinterpret(N0f8, readimage(path, 10_000, 1))
            @test CIFAR100.testtensor(Float64, 1) ≈ readimage(path, 10_000, 1) ./ 255
            @test CIFAR100.testtensor(Int, 1) == Int.(readimage(path, 10_000, 1))
            @test CIFAR100.testtensor(UInt8, 1) == readimage(path, 10_000, 1)
            @test CIFAR100.testtensor(UInt8, 2) == readimage(path, 10_000, 2)
            @test CIFAR100.testtensor(UInt8, 10_000) == readimage(path, 10_000, 10_000)
        end

        # These tests check if the functions return internaly
        # consistent results for different parameters (e.g. index
        # as int or as vector). That means no matter how you
        # specify an index, you will always get the same result
        # for a specific index.
        for (image_fun, T, nimages) in (
                (CIFAR100.traintensor, Float32, 50_000),
                (CIFAR100.traintensor, Float64, 50_000),
                (CIFAR100.traintensor, N0f8,    50_000),
                (CIFAR100.traintensor, Int,     50_000),
                (CIFAR100.traintensor, UInt8,   50_000),
                (CIFAR100.testtensor,  Float32, 10_000),
                (CIFAR100.testtensor,  Float64, 10_000),
                (CIFAR100.testtensor,  N0f8,    10_000),
                (CIFAR100.testtensor,  Int,     10_000),
                (CIFAR100.testtensor,  UInt8,   10_000)
            )
            @testset "$image_fun with T=$T" begin
                # whole image set
                A = @inferred image_fun(T)
                @test typeof(A) <: Union{Array{T,4},Base.ReinterpretArray{T,4}}
                @test size(A) == (32,32,3,nimages)

                @test_throws AssertionError image_fun(T,-1)
                @test_throws AssertionError image_fun(T,0)
                @test_throws AssertionError image_fun(T,nimages+1)

                @testset "load single images" begin
                    # Sample a few random images to compare
                    for i = rand(1:nimages, 10)
                        A_i = @inferred image_fun(T,i)
                        @test typeof(A_i) <: Union{Array{T,3},Base.ReinterpretArray{T,3}}
                        @test size(A_i) == (32,32,3)
                        @test A_i == A[:,:,:,i]
                    end
                end

                @testset "load multiple images" begin
                    A_5_10 = @inferred image_fun(T,5:10)
                    @test typeof(A_5_10) <: Union{Array{T,4},Base.ReinterpretArray{T,4}}
                    @test size(A_5_10) == (32,32,3,6)
                    for i = 1:6
                        @test A_5_10[:,:,:,i] == A[:,:,:,i+4]
                    end

                    # also test edge cases `1`, `nimages`
                    indices = [10,3,9,1,nimages]
                    A_vec   = image_fun(T,indices)
                    A_vec_f = image_fun(T,Vector{Int32}(indices))
                    @test typeof(A_vec) <: Union{Array{T,4},Base.ReinterpretArray{T,4}}
                    @test typeof(A_vec_f) <: Union{Array{T,4},Base.ReinterpretArray{T,4}}
                    @test size(A_vec)   == (32,32,3,5)
                    @test size(A_vec_f) == (32,32,3,5)
                    for i in 1:5
                        @test A_vec[:,:,:,i] == A[:,:,:,indices[i]]
                        @test A_vec[:,:,:,i] == A_vec_f[:,:,:,i]
                    end
                end
            end
        end
    end

    @testset "Labels" begin
        # Sanity check that the first trainlabel is not also
        # the first testlabel
        @test CIFAR100.trainlabels(1) != CIFAR100.testlabels(1)

        # Check a few hand picked examples. I looked at both the
        # pictures and the native output to make sure these
        # values are correspond to the image at the same index.
        @test CIFAR100.trainlabels(1) === (11, 19)
        @test CIFAR100.trainlabels(2) === (15, 29)
        @test CIFAR100.trainlabels(1337) === (11, 19)
        @test CIFAR100.trainlabels(0xCAF) === (18, 58)
        @test CIFAR100.trainlabels(50_000) === (1, 73)
        @test CIFAR100.testlabels(1) === (10, 49)
        @test CIFAR100.testlabels(4) === (4, 51)
        @test CIFAR100.testlabels(0xDAD) === (9, 12)
        @test CIFAR100.testlabels(10_000) === (2, 70)

        # These tests check if the functions return internaly
        # consistent results for different parameters (e.g. index
        # as int or as vector). That means no matter how you
        # specify an index, you will always get the same result
        # for a specific index.
        # -- However, technically these tests do not check if
        #    these are the actual CIFAR100 labels of that index!
        for (label_fun, nlabels) in
                    ((CIFAR100.trainlabels, 50_000),
                     (CIFAR100.testlabels,  10_000))
            @testset "$label_fun" begin
                # whole label set
                A, B = @inferred label_fun()
                @test typeof(A) <: Vector{Int64}
                @test size(A) == (nlabels,)
                @test typeof(B) <: Vector{Int64}
                @test size(B) == (nlabels,)

                @testset "load single label" begin
                    # Sample a few random labels to compare
                    for i = rand(1:nlabels, 10)
                        A_i, B_i = @inferred label_fun(i)
                        @test typeof(A_i) <: Int64
                        @test A_i == A[i]
                        @test typeof(B_i) <: Int64
                        @test B_i == B[i]
                    end
                end

                @testset "load multiple labels" begin
                    A_5_10, B_5_10 = @inferred label_fun(5:10)
                    @test typeof(A_5_10) <: Vector{Int64}
                    @test size(A_5_10) == (6,)
                    @test typeof(B_5_10) <: Vector{Int64}
                    @test size(B_5_10) == (6,)
                    for i = 1:6
                        @test A_5_10[i] == A[i+4]
                        @test B_5_10[i] == B[i+4]
                    end

                    # also test edge cases `1`, `nlabels`
                    indices = [10,3,9,1,nlabels]
                    A_vec, B_vec     = @inferred label_fun(indices)
                    A_vec_f, B_vec_f = @inferred label_fun(Vector{Int32}(indices))
                    @test typeof(A_vec)   <: Vector{Int64}
                    @test typeof(A_vec_f) <: Vector{Int64}
                    @test size(A_vec)   == (5,)
                    @test size(A_vec_f) == (5,)
                    @test typeof(B_vec)   <: Vector{Int64}
                    @test typeof(B_vec_f) <: Vector{Int64}
                    @test size(B_vec)   == (5,)
                    @test size(B_vec_f) == (5,)
                    for i in 1:5
                        @test A_vec[i] == A[indices[i]]
                        @test A_vec[i] == A_vec_f[i]
                        @test B_vec[i] == B[indices[i]]
                        @test B_vec[i] == B_vec_f[i]
                    end
                end
            end
        end
    end

    # Check against the already tested tensor and labels functions
    @testset "Data" begin
        for (data_fun, feature_fun, label_fun, nobs) in
                ((CIFAR100.traindata, CIFAR100.traintensor, CIFAR100.trainlabels, 50_000),
                 (CIFAR100.testdata,  CIFAR100.testtensor,  CIFAR100.testlabels,  10_000))
            @testset "check $data_fun against $feature_fun and $label_fun" begin
                data, l1, l2 = @inferred data_fun()
                @test data == @inferred feature_fun()
                @test (l1, l2) == @inferred label_fun()

                for i = rand(1:nobs, 10)
                    d_i, l1_i, l2_i = @inferred data_fun(i)
                    @test d_i == @inferred feature_fun(i)
                    @test (l1_i, l2_i) == @inferred label_fun(i)
                end

                data, l1, l2 = @inferred data_fun(5:10)
                @test data == @inferred feature_fun(5:10)
                @test (l1, l2) == @inferred label_fun(5:10)

                data, l1, l2 = @inferred data_fun(Int, 5:10)
                @test data == @inferred feature_fun(Int, 5:10)
                @test (l1, l2) == @inferred label_fun(5:10)

                indices = [10,3,9,1,nobs]
                data, l1, l2 = @inferred data_fun(indices)
                @test data == @inferred feature_fun(indices)
                @test (l1, l2) == @inferred label_fun(indices)
            end
        end
    end
end

end # module


text:
# This file illustrates the use of memory blocks to break algebraic loops

using Causal
using Plots

# Simualation settings
t0, dt, tf = 0, 1 / 64, 1.

# Construct the model
@defmodel model begin
    @nodes begin
        gen = FunctionGenerator(readout=identity)
        adder = Adder(signs = (+, -))
        mem = Memory(delay = dt)
        writer = Writer(input=Inport(2))
    end
    @branches begin
        gen[1]      =>      adder[1]
        adder[1]    =>      mem[1]
        mem[1]      =>      adder[2]
        gen[1]      =>      writer[1]
        adder[1]    =>      writer[2]
    end
end

# Simulate the model
sim = simulate!(model, t0, dt, tf)

# Diplay model taskmanager
display(model.taskmanager.pairs)

# Read the simulation data
t, x = read(getnode(model, :writer).component)

# Plot the results
p1 = plot(t, x[:, 1], label=:u, marker=(:circle, 1))
    plot!(t, x[:, 2], label=:y, marker=(:circle, 1))
display(p1)


text:
module DualNumbers

"""
Dual numbers with a scalar dual part
"""
struct Dual{T} <: Real
    primal::T
    dual::T
    Dual{T}(x, y) where {T} = new{T}(T(x), T(y))
end
Dual(x::T1, y::T2) where {T1,T2} = Dual{promote_type(T1, T2)}(x, y)
Dual{T}(x) where {T} = Dual{T}(x, 0)
Dual(x) = Dual(x, 0)
export Dual

# Only compare primal parts
Base.:(==)(x::Dual, y::Dual) = x.primal == y.primal
Base.:(<)(x::Dual, y::Dual) = x.primal < y.primal

Base.:(==)(x::Dual, y::Real) = x == Dual(y)
Base.:(<)(x::Dual, y::Real) = x < Dual(y)
Base.:(==)(x::Real, y::Dual) = Dual(x) == y
Base.:(<)(x::Real, y::Dual) = Dual(x) < y

Base.promote_rule(::Type{Dual{T1}}, ::Type{Dual{T2}}) where {T1,T2} = Dual{promote_type(T1, T2)}
Base.promote_rule(::Type{T1}, ::Type{Dual{T2}}) where {T1<:Real,T2} = Dual{promote_type(T1, T2)}
Base.promote_rule(::Type{Bool}, ::Type{Dual{T2}}) where {T2} = Dual{promote_type(Bool, T2)}

Base.zero(x::Dual) = Dual(zero(x.primal))
Base.one(x::Dual) = Dual(one(x.primal))

Base.:+(x::Dual) = Dual(+x.primal, +x.dual)
Base.:-(x::Dual) = Dual(-x.primal, -x.dual)

Base.:+(x::Dual, y::Dual) = Dual(x.primal + y.primal, x.dual + y.dual)
Base.:-(x::Dual, y::Dual) = Dual(x.primal - y.primal, x.dual - y.dual)

Base.inv(x::Dual) = Dual(inv(x.primal), -x.dual / x.primal^2)

Base.:*(x::Dual, y::Dual) = Dual(x.primal * y.primal, x.dual * y.primal + x.primal * y.dual)
Base.:/(x::Dual, y::Dual) = Dual(x.primal / y.primal, x.dual / y.primal - x.primal / y.primal^2 * y.dual)
Base.:\(x::Dual, y::Dual) = Dual(x.primal \ y.primal, x.primal \ y.dual - (x.primal^2 \ y.primal) * x.dual)

Base.literal_pow(::typeof(^), x::Dual, ::Val{0}) = one(x)
function Base.literal_pow(::typeof(^), x::Dual, ::Val{n}) where {n}
    return Dual(Base.literal_pow(^, x.primal, Val(n)), n * Base.literal_pow(^, x.primal, Val(n - 1)) * x.dual)
end
Base.:^(x::Dual, n::Integer) = n == 0 ? one(x) : Dual(x.primal^n, n * x.primal^(n - 1) * x.dual)
Base.:^(x::Dual, y::Real) = (@assert !(y isa Dual); Dual(x.primal^y, y * x.primal^(y - 1) * x.dual))
Base.:^(x::Dual, y::Dual) = error("not implemented")

Base.cos(x::Dual) = Dual(cos(x.primal), -x.dual * sin(x.primal))
Base.log(x::Dual) = Dual(log(x.primal), x.dual / x.primal)
Base.sin(x::Dual) = Dual(sin(x.primal), x.dual * cos(x.primal))
Base.sqrt(x::Dual) = Dual(sqrt(x.primal), x.dual / (2 * sqrt(x.primal)))

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

"""
    derivative(f, x)
    derivative(f, xs, i)
    derivative(f, x, c)
    derivative(f, xs, i, c)

Evaluate the partial derivative of the scalar-argument function `f(x)`
at `x`, or of the array-argument function `f(xs)` at `xs[i]`.

When the argument `x` (or `xs[i]`) is complex, then an additional
argument `c` chooses whether the partial derivative is with respect to
the real (`c==1`) or imaginary (`c==2`) part of the argument.
"""
function derivative end
export derivative

derivative(f, x::T) where {T<:Real} = f(Dual(x, 1)).dual

function derivative(f, xs::AbstractArray{T}, n::Integer) where {T<:Real}
    xs′ = similar(xs, Dual{T})
    for i in 1:length(xs′)
        xs′[i] = Dual(xs[i], i == n)
    end
    return f(xs′).dual
end

function derivative(f, x::Complex{T}, c::Integer) where {T<:Real}
    fx′ = f(Complex(Dual(real(x), c == 1), Dual(imag(x), c == 2)))
    fx′::Complex
    return Complex(real(fx′).dual, imag(fx′).dual)
end

function derivative(f, xs::AbstractArray{Complex{T}}, n::Integer, c::Integer) where {T<:Real}
    xs′ = similar(xs, Complex{Dual{T}})
    for i in 1:length(xs′)
        xs′[i] = Complex(Dual(real(xs[i]), i == n && c == 1), Dual(imag(xs[i]), i == n && c == 2))
    end
    fx′ = f(xs′)
    fx′::Complex
    return Complex(real(fx′).dual, imag(fx′).dual)
end

end


text:
function step_names_params(team::Int64,project::Int64,experiment::Int64,task::Int64,protocol::Int64)
    data = get_steps(team::Int64,project::Int64,experiment::Int64,task::Int64,protocol::Int64)
    step_ids = [parse(Int,data[i]["id"]) for i in 1:length(data)]
    comments = [comment2measurement(html2string(data[i]["attributes"]["description"])) for i in 1:length(data)]
    params = [get_step_table(team,project,experiment,task,protocol,step_ids[i]) for i in 1:length(step_ids)]
    df = DataFrame()
        df.Names  = [data[i]["attributes"]["name"] for i in 1:length(data)]
        df.Scan_Rate = [params[i][2] for i in 1:length(params)]
        df.Pump = [params[i][2] for i in 1:length(params)]
        df.Comment

end
"""
    html2string(html)

Convert a html to a normal String.
"""
function html2string(html_raw)
    html_parsed = parsehtml(html_raw)
    s = text(html_parsed.root)

end

"""
comment2measurement(str)
Extract Measurment and other comments from string.
"""
function comment2measurement(str)

    _idx = findall(x-> isspace(x)==true,str)
    _idx_s =[_idx[i] for i in 2:length(_idx) if _idx[i] !== _idx[i-1]+1]

    regex = r"(\w*|$)(?:[ ]{2,})"

    if isempty(_idx)==true
        measurement = str
        pump = "Pump: n/a"
        rest = "comment n/a"
        return measurement,pump,rest
    else
        if length(_idx_s) == 0
            matches = collect(eachmatch(regex,str))
            measurement = first(matches).captures[1]
            second = last(matches)
            if occursin("pump",lowercase(str))== true
                pump = "Pump: "* second.captures[1]
            else
                pump = second.captures[1]
            end
            rest= "comment n/a"

            return measurement,pump,rest

        elseif length(_idx_s) >= 1
            regex2 = r"(?:[ ]{2,})(.*)"

            matches = collect(eachmatch(regex,str))
            measurement = first(matches).captures[1]
            second = last(matches)
            if occursin("pump",lowercase(str))== true
                pump = "Pump: "* second.captures[1]
            else
                pump = second.captures[1]
            end
            rest =  match(regex2,str[second.offset+length(second.captures[1]):end]).captures[1]

            return measurement,pump,rest
        end
    end

end

max_stars_repo_name: lambcount/EC

max_issues_repo_name: lambcount/EC

max_forks_repo_name: lambcount/EC


text:
module PSIS

using Turing, DynamicPPL

export psislw
export psisloo

# Copyright (c) 2015 Aki Vehtari, Tuomas Sivula
# Original Matlab version by Aki Vehtari. Translation to Python
# by Tuomas Sivula.

# This software is distributed under the GNU General Public
# License (version 3 or later); please refer to the file
# License.txt, included with the software, for details.

"""
    psisloo(log_lik, [wcpp, wtrunc])
PSIS leave-one-out log predictive densities.
Computes the log predictive densities given posterior samples of the log
likelihood terms p(y_i|\theta^s) in input parameter `log_lik`. Returns a
sum of the leave-one-out log predictive densities `loo`, individual
leave-one-out log predictive density terms `loos` and an estimate of Pareto
tail indeces `ks`. If tail index k>0.5, variance of the raw estimate does
not exist and if tail index k>1 the mean of the raw estimate does not exist
and the PSIS estimate is likely to have large variation and some bias.
# Arguments
* `log_lik::Union{AbstractArray, Mamba.Chains}`: Array of size n x m containing n posterior samples of the log likelihood terms p(y_i|\theta^s).
* `wcpp::Real`: Percentage of samples used for GPD fit estimate (default is 20).
* `wtrunc::Float64`: Positive parameter for truncating very large weights to n^wtrunc. Providing False or 0 disables truncation. Default values is 3/4.
# Returns
* `loo::Real`: sum of the leave-one-out log predictive densities.
* `loos::AbstractArray`: Individual leave-one-out log predictive density terms.* `ks::AbstractArray`: Estimated Pareto tail indeces.
"""
# Compute LOO and standard error
function psisloo(model:: Model, chain:: Chains, wcpp::Int64=20, wtrunc::Float64=3/4)

    # log raw weights from log_lik

    lw = pointwise_loglikelihoods(model, chain)

    # compute Pareto smoothed log weights given raw log weights
    lwp, ks = psislw(-lw, wcpp, wtrunc)

    lwp += lw
    loos = logsumexp(lwp, 1)
    loo = sum(loos)

    return loo, loos, ks
end

"""
    psislw(lw, [wcpp, wtrunc])
Compute the Pareto smoothed importance sampling (PSIS).
# Arguments
* `lw::Union{AbstractArray, Mamba.Chains}`: Array of size n x m containing m sets of n log weights. It is also possible to provide one dimensional array of length n.
* `wcpp::Real`: Percentage of samples used for GPD fit estimate (default is 20).
* `wtrunc::Float64`: Positive parameter for truncating very large weights to n^wtrunc. Providing False or 0 disables truncation. Default values is 3/4.
# Returns
* `lw_out::AbstractArray`: Smoothed log weights
* `kss::AbstractArray`: Pareto tail indices
"""
function psislw(lw::Union{AbstractArray, Mamba.Chains},
                wcpp::Int64=20, wtrunc::Float64=3/4)

    if isa(lw, Mamba.Chains)
        lw_out = Mamba.combine(lw)
    else
        lw_out = copy(lw)
    end

    if ~(1 <= ndims(lw_out) <= 2)
        throw(DimensionMismatch("Argument `lw` must be 1 or 2 dimensional."))
    end
    if size(lw_out,1) <= 1
        error("More than one log-weight needed.")
    end
    return _psislw(lw_out, wcpp, wtrunc)
end

function _psislw(lw_out::Array{Float64}, wcpp::Int64, wtrunc::Float64)

    if ndims(lw_out) == 2
        n, m = size(lw_out)
    elseif ndims(lw_out) == 1
        n = length(lw_out)
        m = 1
    end
    kss = zeros(Float64, m)

    # precalculate constants
    cutoffmin = log(realmin(Float64))
    logn = log(n)

    # loop over sets of log weights
    for i in 1:m
        x = view(lw_out, :, i)
        x .-= maximum(x)
        xcutoff = max(quantile(x, 1 - wcpp/100), cutoffmin)
        expxcutoff = exp(xcutoff)
        x2 = x[x .> xcutoff]
        n2 = length(x2)
        if n2 <= 4
            # not enough tail samples for gpdfitnew
            k = Inf
        else
            # order of tail samples
            x2si = sortperm(x2)
            # fit generalized Pareto distribution to the right tail samples
            x2 .= exp.(x2) .- expxcutoff
            k, sigma = gpdfitnew(x2)
            # compute ordered statistic for the fit
            sti = collect(0.5:n2) / n2
            qq = gpinv(sti, k, sigma)
            qq .= log.(qq .+ expxcutoff)
            # place the smoothed tail into the output array
            x[x .> xcutoff] = qq[x2si]
        end
        if wtrunc > 0
            # truncate too large weights
            lwtrunc = wtrunc * logn - logn + logsumexp(x)
            x[x .> lwtrunc] = lwtrunc
        end
        # renormalize weights
        x .-= logsumexp(x)
        # store tail index k
        kss[i] = k
    end
    # If the provided input array is one dimensional, return kss as scalar.
    if ndims(lw_out) == 1
        kss = kss[1]
    end
    return lw_out, kss
end

"""
    gpdfitnew(x)
Estimate the paramaters for the Generalized Pareto Distribution (GPD). Returns empirical Bayes estimate for the parameters of the two-parameter generalized Parato distribution given the data.
# Arguments
* `x::AbstractArray`: One dimensional data array.
# Returns
* `k::Float64, sigma::Float64`: Estimated parameter values.
# Notes
* This function returns a negative of Zhang and Stephens's k, because it is more common parameterization.
"""
function gpdfitnew{T<:Real}(x::AbstractArray{T})
    if ndims(x) != 1 || length(x) <= 1
        throw(DimensionMismatch())
    end
    sort!(x)

    n = length(x)
    m = 80 + trunc(Int, sqrt(n))
    prior = 3

    bs = collect(1:m) - 0.5
    bs .= (1 .- sqrt.(m ./ bs)) ./ prior ./ x[trunc(Int64, n/4 + 0.5)] .+ 1 / x[end]

    ks = mean(log1p.(-bs .* x'), 2)[:,1]
    L = n * (log.(bs ./ -ks) .- ks .- 1)
    w = similar(L)
    w .= (1 ./ sum(exp.(L .- L'), 1)')[:,1]

    # remove negligible weights
    dii = w .>= (10 * eps(Float64))
    if ~all(dii)
        w = w[dii]
        bs = bs[dii]
    end
    # normalise w
    w ./= sum(w)
    # posterior mean for b
    b = sum(bs .* w)
    # Estimate for k, note that we return a negative of Zhang and
    # Stephens's k, because it is more common parameterisation.
    k = mean(log1p.(-b .* x))
    # estimate for sigma
    sigma = -k / b

    return k, sigma

end

"""
    gpinv(p, k, sigma)
Inverse Generalised Pareto distribution function.
"""
function gpinv(p, k, sigma)
    x = similar(p)
    if sigma <= 0
        return x
    end
    ok = (p .> 0) & (p .< 1)
    pv = view(p, ok)
    xv = view(x,ok)
    if abs(k) < eps(Float64)
        xv .= -log1p.(-pv)
    else
        xv .= expm1.(log1p.(-pv) .* (-k)) ./ k
    end
    x .*= sigma
    if ~all(ok)
        x[p .== 0] = 0
        if k >= 0
            x[p .== 1] = Inf
        else
            x[p .== 1] = -sigma / k
        end
    end
    x
end

"""
    logsumexp(x[, d])
Compute `log(sum(exp(x), d))` of `x`.
"""
function logsumexp{T<:Real}(x::AbstractArray{T}, d::Int64=1)
    u = maximum(x)
    abs(u) == Inf && return any(isnan, x) ? T(NaN) : u
    result = u .+ log.(sum(exp.(x-u), d))
    if length(result) == 1
        result = result[1]
    end
    result
end

end


text:
function test_rgb_over_rgb(bm, expected_f64)
    @test all(t->bm(t[1]...) ≈ t[2],
              zip(((c1, c2) for c1 in rgb1_f64, c2 in rgb2_f64), expected_f64))

    @test all(c->bm(c[1], c[2]) ≈ bm(c[1]|>f64, c[2]|>f64)|>f32,
              ((c1, c2) for c1 in rgb1_f32, c2 in rgb2_f32))
    @test all(c->bm(c[1], c[2]) === bm(c[1]|>f64, c[2]|>f64)|>n0f8,
              ((c1, c2) for c1 in rgb1_n0f8, c2 in rgb2_n0f8))
    @test all(c->bm(c[1], c[2]) === bm(c[1]|>f64, c[2]|>f64)|>n0f16,
              ((c1, c2) for c1 in rgb1_n0f16, c2 in rgb2_n0f16))

    @test all(c->bm(c[1], c[2]) ≈ bm(c[1]|>f64, c[2]|>f64)|>f32,
              ((c2, c1) for c1 in rgb1_f32, c2 in rgb2_f32))
    @test all(c->bm(c[1], c[2]) === bm(c[1]|>f64, c[2]|>f64)|>n0f8,
              ((c2, c1) for c1 in rgb1_n0f8, c2 in rgb2_n0f8))
    @test all(c->bm(c[1], c[2]) === bm(c[1]|>f64, c[2]|>f64)|>n0f16,
              ((c2, c1) for c1 in rgb1_n0f16, c2 in rgb2_n0f16))
end

@testset "normal" begin
    @testset "RGB over RGB" begin
        test_rgb_over_rgb(BlendNormal, [c2 for c1 in rgb1_f64, c2 in rgb2_f64])
    end

    @testset "RGBA over RGB: $T" for T in Ts
        @test BlendNormal(RGB{T}(1, 0.75, 0), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGB{T}(0.4, 0.6, 0.6)
    end

    @testset "RGB over RGBA: $T" for T in Ts
        @test BlendNormal(RGBA{T}(1, 0.75, 0, 0.6), RGB{T}(0, 0.5, 1)) ≈ RGBA{T}(0, 0.5, 1, 1)
    end

    @testset "RGBA over RGBA: $T" for T in Ts
        @test BlendNormal(RGBA{T}(1, 0.75, 0, 0.6), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGBA{T}(2/7, 4/7, 5/7, 0.84)
    end
end

@testset "multiply" begin
    @testset "RGB over RGB" begin
        test_rgb_over_rgb(BlendMultiply, RGB{Float64}[
            RGB(0.0, 0.0,                     0.0)                    RGB(0.0,                     0.0,                     0.0)                     RGB(0.0,                     0.0,                    0.0)
            RGB(0.0, 0.0,                     5.562684646268003e-309) RGB(1.1125369292536007e-308, 1.1125369292536007e-308, 1.1125369292536007e-308) RGB(2.2250738585072014e-308, 1.668805393880401e-308, 2.2250738585072014e-308)
            RGB(0.0, 5.562684646268003e-309,  0.0625)                 RGB(0.125,                   0.12500000000000003,     0.12499999999999999)     RGB(0.25,                    0.1875,                 0.24999999999999997)
            RGB(0.0, 1.1125369292536007e-308, 0.12499999999999999)    RGB(0.24999999999999997,     0.25,                    0.24999999999999994)     RGB(0.49999999999999994,     0.37499999999999994,    0.4999999999999999)
            RGB(0.0, 1.1125369292536007e-308, 0.125)                  RGB(0.25,                    0.25000000000000006,     0.24999999999999997)     RGB(0.5,                     0.375,                  0.49999999999999994)
            RGB(0.0, 1.1125369292536007e-308, 0.12500000000000003)    RGB(0.25000000000000006,     0.2500000000000001,      0.25)                    RGB(0.5000000000000001,      0.3750000000000001,     0.5)
            RGB(0.0, 1.668805393880401e-308,  0.1875)                 RGB(0.375,                   0.3750000000000001,      0.37499999999999994)     RGB(0.75,                    0.5625,                 0.7499999999999999)
            RGB(0.0, 2.2250738585072014e-308, 0.24999999999999997)    RGB(0.49999999999999994,     0.5,                     0.4999999999999999)      RGB(0.9999999999999999,      0.7499999999999999,     0.9999999999999998)
            RGB(0.0, 2.2250738585072014e-308, 0.25)                   RGB(0.5,                     0.5000000000000001,      0.49999999999999994)     RGB(1.0,                     0.75,                   0.9999999999999999)
        ])
    end

    @testset "RGBA over RGB: $T" for T in Ts
        @test BlendMultiply(RGB{T}(1, 0.75, 0), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGB{T}(0.4, 0.525, 0.0)
    end

    @testset "RGB over RGBA: $T" for T in Ts
        @test BlendMultiply(RGBA{T}(1, 0.75, 0, 0.6), RGB{T}(0, 0.5, 1)) ≈ RGBA{T}(0, 0.425, 0.4, 1)
    end

    @testset "RGBA over RGBA: $T" for T in Ts
        @test BlendMultiply(RGBA{T}(1, 0.75, 0, 0.6), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGBA{T}(2/7, 29/56, 2/7, 0.84)
    end
end

@testset "screen" begin
    @testset "RGB over RGB" begin
        test_rgb_over_rgb(BlendScreen, RGB{Float64}[
            RGB(0.0,                 0.0,                 0.25)               RGB(0.5,   0.5000000000000001, 0.49999999999999994) RGB(1.0, 0.75,   0.9999999999999999)
            RGB(0.0,                 0.0,                 0.25)               RGB(0.5,   0.5000000000000001, 0.49999999999999994) RGB(1.0, 0.75,   0.9999999999999999)
            RGB(0.25,                0.25,                0.4375)             RGB(0.625, 0.6250000000000001, 0.625)               RGB(1.0, 0.8125, 0.9999999999999999)
            RGB(0.49999999999999994, 0.49999999999999994, 0.625)              RGB(0.75,  0.75,               0.7499999999999999)  RGB(1.0, 0.875,  0.9999999999999999)
            RGB(0.5,                 0.5,                 0.625)              RGB(0.75,  0.75,               0.75)                RGB(1.0, 0.875,  1.0)
            RGB(0.5000000000000001,  0.5000000000000001,  0.6250000000000001) RGB(0.75,  0.7500000000000001, 0.75)                RGB(1.0, 0.875,  1.0)
            RGB(0.75,                0.75,                0.8125)             RGB(0.875, 0.875,              0.875)               RGB(1.0, 0.9375, 1.0)
            RGB(0.9999999999999999,  0.9999999999999999,  0.9999999999999999) RGB(1.0,   1.0,                0.9999999999999999)  RGB(1.0, 1.0,    1.0)
            RGB(1.0,                 1.0,                 1.0)                RGB(1.0,   1.0,                1.0)                 RGB(1.0, 1.0,    1.0)
        ])
    end

    @testset "RGBA over RGB: $T" for T in Ts
        @test BlendScreen(RGB{T}(1, 0.75, 0), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGB{T}(1.0, 0.825, 0.6)
    end

    @testset "RGB over RGBA: $T" for T in Ts
        @test BlendScreen(RGBA{T}(1, 0.75, 0, 0.6), RGB{T}(0, 0.5, 1)) ≈ RGBA{T}(0.6, 0.725, 1.0, 1)
    end

    @testset "RGBA over RGBA: $T" for T in Ts
        @test BlendScreen(RGBA{T}(1, 0.75, 0, 0.6), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGBA{T}(5/7, 41/56, 5/7, 0.84)
    end
end

@testset "overlay" begin
    @testset "RGB over RGB" begin
        test_rgb_over_rgb(BlendOverlay, RGB{Float64}[
            RGB(0.0,                   0.0,                     0.0)                     RGB(0.0,                     0.0,                    0.0)                     RGB(0.0,                    0.0,                    0.0)
            RGB(0.0,                   0.0,                     1.1125369292536007e-308) RGB(2.2250738585072014e-308, 2.225073858507202e-308, 2.2250738585072014e-308) RGB(4.450147717014403e-308, 3.337610787760802e-308, 4.4501477170144023e-308)
            RGB(0.0,                   1.1125369292536007e-308, 0.125)                   RGB(0.25,                    0.25000000000000006,    0.24999999999999997)     RGB(0.5,                    0.375,                  0.49999999999999994)
            RGB(0.0,                   2.2250738585072014e-308, 0.24999999999999997)     RGB(0.49999999999999994,     0.5,                    0.4999999999999999)      RGB(0.9999999999999999,     0.7499999999999999,     0.9999999999999998)
            RGB(0.0,                   2.2250738585072014e-308, 0.25)                    RGB(0.5,                     0.5000000000000001,     0.49999999999999994)     RGB(1.0,                    0.75,                   0.9999999999999999)
            RGB(2.220446049250313e-16, 2.220446049250313e-16,   0.25000000000000017)     RGB(0.5000000000000001,      0.5000000000000002,     0.5000000000000001)      RGB(1.0,                    0.75,                   0.9999999999999999)
            RGB(0.5,                   0.5,                     0.625)                   RGB(0.75,                    0.75,                   0.75)                    RGB(1.0,                    0.875,                  1.0)
            RGB(0.9999999999999998,    0.9999999999999998,      0.9999999999999998)      RGB(0.9999999999999999,      0.9999999999999999,     0.9999999999999999)      RGB(1.0,                    1.0,                    1.0)
            RGB(1.0,                   1.0,                     1.0)                     RGB(1.0,                     1.0,                    1.0)                     RGB(1.0,                    1.0,                    1.0)
        ])
    end

    @testset "RGBA over RGB: $T" for T in Ts
        @test BlendOverlay(RGB{T}(1, 0.75, 0), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGB{T}(1.0, 0.75, 0.0)
    end

    @testset "RGB over RGBA: $T" for T in Ts
        @test BlendOverlay(RGBA{T}(1, 0.75, 0, 0.6), RGB{T}(0, 0.5, 1)) ≈ RGBA{T}(0.6, 0.65, 0.4, 1)
    end

    @testset "RGBA over RGBA: $T" for T in Ts
        @test BlendOverlay(RGBA{T}(1, 0.75, 0, 0.6), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGBA{T}(5/7, 19/28, 2/7, 0.84)
    end
end

@testset "darken" begin
    @testset "RGB over RGB" begin
        test_rgb_over_rgb(BlendDarken, RGB{Float64}[
            RGB(0.0, 0.0,                     0.0)                     RGB(0.0,                     0.0,                     0.0)                     RGB(0.0,                     0.0,                     0.0)
            RGB(0.0, 2.2250738585072014e-308, 2.2250738585072014e-308) RGB(2.2250738585072014e-308, 2.2250738585072014e-308, 2.2250738585072014e-308) RGB(2.2250738585072014e-308, 2.2250738585072014e-308, 2.2250738585072014e-308)
            RGB(0.0, 2.2250738585072014e-308, 0.25)                    RGB(0.25,                    0.25,                    0.25)                    RGB(0.25,                    0.25,                    0.25)
            RGB(0.0, 2.2250738585072014e-308, 0.25)                    RGB(0.49999999999999994,     0.49999999999999994,     0.49999999999999994)     RGB(0.49999999999999994,     0.49999999999999994,     0.49999999999999994)
            RGB(0.0, 2.2250738585072014e-308, 0.25)                    RGB(0.5,                     0.5,                     0.49999999999999994)     RGB(0.5,                     0.5,                     0.5)
            RGB(0.0, 2.2250738585072014e-308, 0.25)                    RGB(0.5,                     0.5000000000000001,      0.49999999999999994)     RGB(0.5000000000000001,      0.5000000000000001,      0.5000000000000001)
            RGB(0.0, 2.2250738585072014e-308, 0.25)                    RGB(0.5,                     0.5000000000000001,      0.49999999999999994)     RGB(0.75,                    0.75,                    0.75)
            RGB(0.0, 2.2250738585072014e-308, 0.25)                    RGB(0.5,                     0.5000000000000001,      0.49999999999999994)     RGB(0.9999999999999999,      0.75,                    0.9999999999999999)
            RGB(0.0, 2.2250738585072014e-308, 0.25)                    RGB(0.5,                     0.5000000000000001,      0.49999999999999994)     RGB(1.0,                     0.75,                    0.9999999999999999)
        ])
    end

    @testset "RGBA over RGB: $T" for T in Ts
        @test BlendDarken(RGB{T}(1, 0.75, 0), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGB{T}(0.4, 0.6, 0.0)
    end

    @testset "RGB over RGBA: $T" for T in Ts
        @test BlendDarken(RGBA{T}(1, 0.75, 0, 0.6), RGB{T}(0, 0.5, 1)) ≈ RGBA{T}(0.0, 0.5, 0.4, 1)
    end

    @testset "RGBA over RGBA: $T" for T in Ts
        @test BlendDarken(RGBA{T}(1, 0.75, 0, 0.6), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGBA{T}(2/7, 4/7, 2/7, 0.84)
    end
end

@testset "lighten" begin
    @testset "RGB over RGB" begin
        test_rgb_over_rgb(BlendLighten, RGB{Float64}[
            RGB(0.0,                     2.2250738585072014e-308, 0.25)                RGB(0.5,                0.5000000000000001, 0.49999999999999994) RGB(1.0, 0.75,               0.9999999999999999)
            RGB(2.2250738585072014e-308, 2.2250738585072014e-308, 0.25)                RGB(0.5,                0.5000000000000001, 0.49999999999999994) RGB(1.0, 0.75,               0.9999999999999999)
            RGB(0.25,                    0.25,                    0.25)                RGB(0.5,                0.5000000000000001, 0.49999999999999994) RGB(1.0, 0.75,               0.9999999999999999)
            RGB(0.49999999999999994,     0.49999999999999994,     0.49999999999999994) RGB(0.5,                0.5000000000000001, 0.49999999999999994) RGB(1.0, 0.75,               0.9999999999999999)
            RGB(0.5,                     0.5,                     0.5)                 RGB(0.5,                0.5000000000000001, 0.5)                 RGB(1.0, 0.75,               0.9999999999999999)
            RGB(0.5000000000000001,      0.5000000000000001,      0.5000000000000001)  RGB(0.5000000000000001, 0.5000000000000001, 0.5000000000000001)  RGB(1.0, 0.75,               0.9999999999999999)
            RGB(0.75,                    0.75,                    0.75)                RGB(0.75,               0.75,               0.75)                RGB(1.0, 0.75,               0.9999999999999999)
            RGB(0.9999999999999999,      0.9999999999999999,      0.9999999999999999)  RGB(0.9999999999999999, 0.9999999999999999, 0.9999999999999999)  RGB(1.0, 0.9999999999999999, 0.9999999999999999)
            RGB(1.0,                     1.0,                     1.0)                 RGB(1.0,                1.0,                1.0)                 RGB(1.0, 1.0,                1.0)
        ])
    end

    @testset "RGBA over RGB: $T" for T in Ts
        @test BlendLighten(RGB{T}(1, 0.75, 0), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGB{T}(1.0, 0.75, 0.6)
    end

    @testset "RGB over RGBA: $T" for T in Ts
        @test BlendLighten(RGBA{T}(1, 0.75, 0, 0.6), RGB{T}(0, 0.5, 1)) ≈ RGBA{T}(0.6, 0.65, 1.0, 1)
    end

    @testset "RGBA over RGBA: $T" for T in Ts
        @test BlendLighten(RGBA{T}(1, 0.75, 0, 0.6), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGBA{T}(5/7, 19/28, 5/7, 0.84)
    end
end

@testset "color-dodge" begin
    @testset "RGB over RGB" begin
        test_rgb_over_rgb(BlendColorDodge, RGB{Float64}[
            RGB(0.0,                     0.0,                     0.0)                     RGB(0.0,                    0.0,                    0.0)                     RGB(0.0, 0.0,                    0.0)
            RGB(2.2250738585072014e-308, 2.2250738585072014e-308, 2.9667651446762683e-308) RGB(4.450147717014403e-308, 4.450147717014404e-308, 4.4501477170144023e-308) RGB(1.0, 8.900295434028806e-308, 2.004168360008973e-292)
            RGB(0.25,                    0.25,                    0.3333333333333333)      RGB(0.5,                    0.5000000000000001,     0.49999999999999994)     RGB(1.0, 1.0,                    1.0)
            RGB(0.49999999999999994,     0.49999999999999994,     0.6666666666666666)      RGB(0.9999999999999999,     1.0,                    0.9999999999999998)      RGB(1.0, 1.0,                    1.0)
            RGB(0.5,                     0.5,                     0.6666666666666666)      RGB(1.0,                    1.0,                    0.9999999999999999)      RGB(1.0, 1.0,                    1.0)
            RGB(0.5000000000000001,      0.5000000000000001,      0.6666666666666669)      RGB(1.0,                    1.0,                    1.0)                     RGB(1.0, 1.0,                    1.0)
            RGB(0.75,                    0.75,                    1.0)                     RGB(1.0,                    1.0,                    1.0)                     RGB(1.0, 1.0,                    1.0)
            RGB(0.9999999999999999,      0.9999999999999999,      1.0)                     RGB(1.0,                    1.0,                    1.0)                     RGB(1.0, 1.0,                    1.0)
            RGB(1.0,                     1.0,                     1.0)                     RGB(1.0,                    1.0,                    1.0)                     RGB(1.0, 1.0,                    1.0)
        ])
    end

    @testset "RGBA over RGB: $T" for T in Ts
        @test BlendColorDodge(RGB{T}(1, 0.75, 0), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGB{T}(1.0, 0.9, 0.0)
    end

    @testset "RGB over RGBA: $T" for T in Ts
        @test BlendColorDodge(RGBA{T}(1, 0.75, 0, 0.6), RGB{T}(0, 0.5, 1)) ≈ RGBA{T}(0.6, 0.8, 0.4, 1)
    end

    @testset "RGBA over RGBA: $T" for T in Ts
        @test BlendColorDodge(RGBA{T}(1, 0.75, 0, 0.6), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGBA{T}(5/7, 11/14, 2/7, 0.84)
    end
end

@testset "color-burn" begin
    @testset "RGB over RGB" begin
        test_rgb_over_rgb(BlendColorBurn, RGB{Float64}[
            RGB(0.0, 0.0, 0.0)                RGB(0.0,                   0.0,                    0.0)                    RGB(0.0,                 0.0,                 0.0)
            RGB(0.0, 0.0, 0.0)                RGB(0.0,                   0.0,                    0.0)                    RGB(0.0,                 0.0,                 0.0)
            RGB(0.0, 0.0, 0.0)                RGB(0.0,                   0.0,                    0.0)                    RGB(0.25,                0.0,                 0.24999999999999992)
            RGB(0.0, 0.0, 0.0)                RGB(0.0,       #=1.11e-16=#2.2204460492503126e-16, 0.0)                    RGB(0.49999999999999994, 0.33333333333333326, 0.4999999999999999)
            RGB(0.0, 0.0, 0.0)                RGB(0.0,                   2.2204460492503126e-16, 0.0)                    RGB(0.5,                 0.3333333333333333,  0.49999999999999994)
            RGB(0.0, 0.0, 0.0)                RGB(2.220446049250313e-16, 4.440892098500625e-16,  1.1102230246251568e-16) RGB(0.5000000000000001,  0.3333333333333335,  0.5000000000000001)
            RGB(0.0, 0.0, 0.0)                RGB(0.5,                   0.5000000000000001,     0.49999999999999994)    RGB(0.75,                0.6666666666666666,  0.75)
            RGB(0.0, 0.0, 0.9999999999999996) RGB(0.9999999999999998,    0.9999999999999998,     0.9999999999999998)     RGB(0.9999999999999999,  0.9999999999999999,  0.9999999999999999)
            RGB(1.0, 1.0, 1.0)                RGB(1.0,                   1.0,                    1.0)                    RGB(1.0,                 1.0,                 1.0)
        ])
    end

    @testset "RGBA over RGB: $T" for T in Ts
        @test BlendColorBurn(RGB{T}(1, 0.75, 0), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGB{T}(1.0, 0.6, 0.0)
    end

    @testset "RGB over RGBA: $T" for T in Ts
        @test BlendColorBurn(RGBA{T}(1, 0.75, 0, 0.6), RGB{T}(0, 0.5, 1)) ≈ RGBA{T}(0.6, 0.5, 0.4, 1)
    end

    @testset "RGBA over RGBA: $T" for T in Ts
        @test BlendColorBurn(RGBA{T}(1, 0.75, 0, 0.6), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGBA{T}(5/7, 4/7, 2/7, 0.84)
    end
end

@testset "hard-light" begin
    @testset "RGB over RGB" begin
        test_rgb_over_rgb(BlendHardLight, RGB{Float64}[
            RGB(0.0, 0.0,                     0.0)                     RGB(0.0,                     2.220446049250313e-16, 0.0)                     RGB(1.0, 0.5,   0.9999999999999998)
            RGB(0.0, 0.0,                     1.1125369292536007e-308) RGB(2.2250738585072014e-308, 2.220446049250313e-16, 2.2250738585072014e-308) RGB(1.0, 0.5,   0.9999999999999998)
            RGB(0.0, 1.1125369292536007e-308, 0.125)                   RGB(0.25,                    0.25000000000000017,   0.24999999999999997)     RGB(1.0, 0.625, 0.9999999999999998)
            RGB(0.0, 2.2250738585072014e-308, 0.24999999999999997)     RGB(0.49999999999999994,     0.5000000000000001,    0.4999999999999999)      RGB(1.0, 0.75,  0.9999999999999999)
            RGB(0.0, 2.2250738585072014e-308, 0.25)                    RGB(0.5,                     0.5000000000000001,    0.49999999999999994)     RGB(1.0, 0.75,  0.9999999999999999)
            RGB(0.0, 2.225073858507202e-308,  0.25000000000000006)     RGB(0.5000000000000001,      0.5000000000000002,    0.5)                     RGB(1.0, 0.75,  0.9999999999999999)
            RGB(0.0, 3.337610787760802e-308,  0.375)                   RGB(0.75,                    0.75,                  0.7499999999999999)      RGB(1.0, 0.875, 1.0)
            RGB(0.0, 4.4501477170144023e-308, 0.49999999999999994)     RGB(0.9999999999999999,      0.9999999999999999,    0.9999999999999998)      RGB(1.0, 1.0,   1.0)
            RGB(0.0, 4.450147717014403e-308,  0.5)                     RGB(1.0,                     1.0,                   0.9999999999999999)      RGB(1.0, 1.0,   1.0)
        ])
    end

    @testset "RGBA over RGB: $T" for T in Ts
        @test BlendHardLight(RGB{T}(1, 0.75, 0), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGB{T}(0.4, 0.75, 0.6)
    end

    @testset "RGB over RGBA: $T" for T in Ts
        @test BlendHardLight(RGBA{T}(1, 0.75, 0, 0.6), RGB{T}(0, 0.5, 1)) ≈ RGBA{T}(0.0, 0.65, 1.0, 1)
    end

    @testset "RGBA over RGBA: $T" for T in Ts
        @test BlendHardLight(RGBA{T}(1, 0.75, 0, 0.6), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGBA{T}(2/7, 19/28, 5/7, 0.84)
    end
end

@testset "soft-light" begin
    @testset "RGB over RGB" begin
        test_rgb_over_rgb(BlendSoftLight, RGB{Float64}[
            RGB(0.0,                 0.0,                 0.0)                     RGB(0.0,                     0.0,                    0.0)                     RGB(0.0,                    0.0,                    0.0)
            RGB(0.0,                 0.0,                 1.1125369292536007e-308) RGB(2.2250738585072014e-308, 2.225073858507203e-308, 2.2250738585072014e-308) RGB(8.900295434028806e-308, 5.562684646268003e-308, 8.900295434028804e-308)
            RGB(0.0625,              0.0625,              0.15625)                 RGB(0.25,                    0.25000000000000006,    0.24999999999999997)     RGB(0.5,                    0.375,                  0.49999999999999994)
            RGB(0.24999999999999994, 0.24999999999999994, 0.37499999999999994)     RGB(0.49999999999999994,     0.5,                    0.49999999999999994)     RGB(0.7071067811865475,     0.6035533905932737,     0.7071067811865475)
            RGB(0.25,                0.25,                0.375)                   RGB(0.5,                     0.5,                    0.5)                     RGB(0.7071067811865476,     0.6035533905932737,     0.7071067811865475)
            RGB(0.2500000000000001,  0.2500000000000001,  0.3750000000000001)      RGB(0.5000000000000001,      0.5000000000000001,     0.5000000000000001)      RGB(0.7071067811865476,     0.6035533905932738,     0.7071067811865476)
            RGB(0.5625,              0.5625,              0.65625)                 RGB(0.75,                    0.75,                   0.75)                    RGB(0.8660254037844386,     0.8080127018922193,     0.8660254037844386)
            RGB(0.9999999999999998,  0.9999999999999998,  0.9999999999999999)      RGB(0.9999999999999999,      0.9999999999999999,     0.9999999999999999)      RGB(0.9999999999999999,     0.9999999999999999,     0.9999999999999999)
            RGB(1.0,                 1.0,                 1.0)                     RGB(1.0,                     1.0,                    1.0)                     RGB(1.0,                    1.0,                    1.0)
        ])
    end

    @testset "RGBA over RGB: $T" for T in Ts
        @test BlendSoftLight(RGB{T}(1, 0.75, 0), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGB{T}(1.0, 0.75, 0.0)
    end

    @testset "RGB over RGBA: $T" for T in Ts
        @test BlendSoftLight(RGBA{T}(1, 0.75, 0, 0.6), RGB{T}(0, 0.5, 1)) ≈ RGBA{T}(0.6, 0.65, 0.4, 1)
    end

    @testset "RGBA over RGBA: $T" for T in Ts
        @test BlendSoftLight(RGBA{T}(1, 0.75, 0, 0.6), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGBA{T}(5/7, 19/28, 2/7, 0.84)
    end
end

@testset "difference" begin
    @testset "RGB over RGB" begin
        test_rgb_over_rgb(BlendDifference, RGB{Float64}[
            RGB(0.0,                     2.2250738585072014e-308, 0.25)                RGB(0.5,                    0.5000000000000001,     0.49999999999999994)    RGB(1.0,                    0.75,                0.9999999999999999)
            RGB(2.2250738585072014e-308, 0.0,                     0.25)                RGB(0.5,                    0.5000000000000001,     0.49999999999999994)    RGB(1.0,                    0.75,                0.9999999999999999)
            RGB(0.25,                    0.25,                    0.0)                 RGB(0.25,                   0.2500000000000001,     0.24999999999999994)    RGB(0.75,                   0.5,                 0.7499999999999999)
            RGB(0.49999999999999994,     0.49999999999999994,     0.24999999999999994) RGB(5.551115123125783e-17,  1.6653345369377348e-16, 0.0)                    RGB(0.5,                    0.25000000000000006, 0.49999999999999994)
            RGB(0.5,                     0.5,                     0.25)                RGB(0.0,                    1.1102230246251565e-16, 5.551115123125783e-17)  RGB(0.5,                    0.25,                0.4999999999999999)
            RGB(0.5000000000000001,      0.5000000000000001,      0.2500000000000001)  RGB(1.1102230246251565e-16, 0.0,                    1.6653345369377348e-16) RGB(0.4999999999999999,     0.2499999999999999,  0.4999999999999998)
            RGB(0.75,                    0.75,                    0.5)                 RGB(0.25,                   0.2499999999999999,     0.25000000000000006)    RGB(0.25,                   0.0,                 0.2499999999999999)
            RGB(0.9999999999999999,      0.9999999999999999,      0.7499999999999999)  RGB(0.4999999999999999,     0.4999999999999998,     0.49999999999999994)    RGB(1.1102230246251565e-16, 0.2499999999999999,  0.0)
            RGB(1.0,                     1.0,                     0.75)                RGB(0.5,                    0.4999999999999999,     0.5)                    RGB(0.0,                    0.25,                1.1102230246251565e-16)
        ])
    end

    @testset "RGBA over RGB: $T" for T in Ts
        @test BlendDifference(RGB{T}(1, 0.75, 0), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGB{T}(1.0, 0.45, 0.6)
    end

    @testset "RGB over RGBA: $T" for T in Ts
        @test BlendDifference(RGBA{T}(1, 0.75, 0, 0.6), RGB{T}(0, 0.5, 1)) ≈ RGBA{T}(0.6, 0.35, 1.0, 1)
    end

    @testset "RGBA over RGBA: $T" for T in Ts
        @test BlendDifference(RGBA{T}(1, 0.75, 0, 0.6), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGBA{T}(5/7, 13/28, 5/7, 0.84)
    end
end

@testset "exclusion" begin
    @testset "RGB over RGB" begin
        test_rgb_over_rgb(BlendExclusion, RGB{Float64}[
            RGB(0.0,                     2.2250738585072014e-308, 0.25)  RGB(0.5, 0.5000000000000001,  0.49999999999999994) RGB(1.0,                    0.75,                0.9999999999999999)
            RGB(2.2250738585072014e-308, 4.450147717014403e-308,  0.25)  RGB(0.5, 0.5000000000000001,  0.49999999999999994) RGB(1.0,                    0.75,                0.9999999999999999)
            RGB(0.25,                    0.25,                    0.375) RGB(0.5, 0.5,                 0.5)                 RGB(0.75,                   0.625,               0.75)
            RGB(0.49999999999999994,     0.49999999999999994,     0.5)   RGB(0.5, 0.5,                 0.5)                 RGB(0.5,                    0.5,                 0.5)
            RGB(0.5,                     0.5,                     0.5)   RGB(0.5, 0.5,                 0.5)                 RGB(0.5,                    0.5,                 0.5)
            RGB(0.5000000000000001,      0.5000000000000001,      0.5)   RGB(0.5, 0.5,                 0.5)                 RGB(0.4999999999999999,     0.49999999999999994, 0.4999999999999999)
            RGB(0.75,                    0.75,                    0.625) RGB(0.5, 0.49999999999999994, 0.5)                 RGB(0.25,                   0.375,               0.25000000000000006)
            RGB(0.9999999999999999,      0.9999999999999999,      0.75)  RGB(0.5, 0.4999999999999999,  0.5)                 RGB(1.1102230246251565e-16, 0.25000000000000006, 2.2204460492503128e-16)
            RGB(1.0,                     1.0,                     0.75)  RGB(0.5, 0.4999999999999999,  0.5)                 RGB(0.0,                    0.25,                1.1102230246251565e-16)
        ])
    end

    @testset "RGBA over RGB: $T" for T in Ts
        @test BlendExclusion(RGB{T}(1, 0.75, 0), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGB{T}(1.0, 0.6, 0.6)
    end

    @testset "RGB over RGBA: $T" for T in Ts
        @test BlendExclusion(RGBA{T}(1, 0.75, 0, 0.6), RGB{T}(0, 0.5, 1)) ≈ RGBA{T}(0.6, 0.5, 1.0, 1)
    end

    @testset "RGBA over RGBA: $T" for T in Ts
        @test BlendExclusion(RGBA{T}(1, 0.75, 0, 0.6), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGBA{T}(5/7, 4/7, 5/7, 0.84)
    end
end

@testset "hue" begin
    @testset "RGB over RGB" begin
        test_rgb_over_rgb(BlendHue, RGB{Float64}[
            RGB(0.0,                     0.0,                     0.0)                     RGB(0.0,                     0.0,                     0.0)                     RGB(0.0,                     0.0,                     0.0)
            RGB(2.2250738585072014e-308, 2.2250738585072014e-308, 2.2250738585072014e-308) RGB(2.2250738585072014e-308, 2.2250738585072014e-308, 2.2250738585072014e-308) RGB(2.2250738585072014e-308, 2.2250738585072014e-308, 2.2250738585072014e-308)
            RGB(0.25,                    0.25,                    0.25)                    RGB(0.25,                    0.25,                    0.25)                    RGB(0.25,                    0.25,                    0.25)
            RGB(0.49999999999999994,     0.49999999999999994,     0.49999999999999994)     RGB(0.49999999999999994,     0.49999999999999994,     0.49999999999999994)     RGB(0.49999999999999994,     0.49999999999999994,     0.49999999999999994)
            RGB(0.5,                     0.5,                     0.5)                     RGB(0.5,                     0.5,                     0.5)                     RGB(0.5,                     0.5,                     0.5)
            RGB(0.5000000000000001,      0.5000000000000001,      0.5000000000000001)      RGB(0.5000000000000001,      0.5000000000000001,      0.5000000000000001)      RGB(0.5000000000000001,      0.5000000000000001,      0.5000000000000001)
            RGB(0.75,                    0.75,                    0.75)                    RGB(0.75,                    0.75,                    0.75)                    RGB(0.75,                    0.75,                    0.75)
            RGB(0.9999999999999999,      0.9999999999999999,      0.9999999999999999)      RGB(0.9999999999999999,      0.9999999999999999,      0.9999999999999999)      RGB(0.9999999999999999,      0.9999999999999999,      0.9999999999999999)
            RGB(1.0,                     1.0,                     1.0)                     RGB(1.0,                     1.0,                     1.0)                     RGB(1.0,                     1.0,                     1.0)
        ])
    end

    @testset "RGBA over RGB: $T" for T in Ts
        @test BlendHue(RGB{T}(1, 0.75, 0), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGB{T}(0.740, 0.770, 0.6) atol=max(0.0005, eps(T))
    end

    @testset "RGB over RGBA: $T" for T in Ts
        @test BlendHue(RGBA{T}(1, 0.75, 0, 0.6), RGB{T}(0, 0.5, 1)) ≈ RGBA{T}(0.340, 0.670, 1.0, 1) atol=max(0.0005, eps(T))
    end

    @testset "RGBA over RGBA: $T" for T in Ts
        @test BlendHue(RGBA{T}(1, 0.75, 0, 0.6), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGBA{T}(0.529, 0.693, 5/7, 0.84) atol=max(0.0005, eps(T))
    end
end

@testset "saturation" begin
    @testset "RGB over RGB" begin
        test_rgb_over_rgb(BlendSaturation, RGB{Float64}[
            RGB(0.0,                     0.0,                     0.0)                     RGB(0.0,                     0.0,                     0.0)                     RGB(0.0,                     0.0,                     0.0)
            RGB(2.2250738585072014e-308, 2.2250738585072014e-308, 2.2250738585072014e-308) RGB(2.2250738585072014e-308, 2.2250738585072014e-308, 2.2250738585072014e-308) RGB(2.2250738585072014e-308, 2.2250738585072014e-308, 2.2250738585072014e-308)
            RGB(0.25,                    0.25,                    0.25)                    RGB(0.25,                    0.25,                    0.25)                    RGB(0.25,                    0.25,                    0.25)
            RGB(0.49999999999999994,     0.49999999999999994,     0.49999999999999994)     RGB(0.49999999999999994,     0.49999999999999994,     0.49999999999999994)     RGB(0.49999999999999994,     0.49999999999999994,     0.49999999999999994)
            RGB(0.5,                     0.5,                     0.5)                     RGB(0.5,                     0.5,                     0.5)                     RGB(0.5,                     0.5,                     0.5)
            RGB(0.5000000000000001,      0.5000000000000001,      0.5000000000000001)      RGB(0.5000000000000001,      0.5000000000000001,      0.5000000000000001)      RGB(0.5000000000000001,      0.5000000000000001,      0.5000000000000001)
            RGB(0.75,                    0.75,                    0.75)                    RGB(0.75,                    0.75,                    0.75)                    RGB(0.75,                    0.75,                    0.75)
            RGB(0.9999999999999999,      0.9999999999999999,      0.9999999999999999)      RGB(0.9999999999999999,      0.9999999999999999,      0.9999999999999999)      RGB(0.9999999999999999,      0.9999999999999999,      0.9999999999999999)
            RGB(1.0,                     1.0,                     1.0)                     RGB(1.0,                     1.0,                     1.0)                     RGB(1.0,                     1.0,                     1.0)
        ])
    end

    @testset "RGBA over RGB: $T" for T in Ts
        @test BlendSaturation(RGB{T}(1, 0.75, 0), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGB{T}(1.0, 0.75, 0.0)
    end

    @testset "RGB over RGBA: $T" for T in Ts
        @test BlendSaturation(RGBA{T}(1, 0.75, 0, 0.6), RGB{T}(0, 0.5, 1)) ≈ RGBA{T}(0.6, 0.65, 0.4, 1)
    end

    @testset "RGBA over RGBA: $T" for T in Ts
        @test BlendSaturation(RGBA{T}(1, 0.75, 0, 0.6), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGBA{T}(5/7, 19/28, 2/7, 0.84)
    end
end

@testset "color" begin
    @testset "RGB over RGB" begin
        test_rgb_over_rgb(BlendColor, RGB{Float64}[
            RGB(0.0,                0.0,                0.0)                     RGB(0.0,                     0.0,                    0.0)                 RGB(0.0,                    0.0,                 0.0)
            RGB(0.0,                0.0,                2.0227944168247287e-307) RGB(1.0749149074914015e-308, 3.224744722474205e-308, 0.0)                 RGB(5.427009410993175e-308, 0.0,                 5.427009410993172e-308)
            RGB(0.2225,             0.2225,             0.4725)                  RGB(0.24999999999999994,     0.25000000000000006,    0.2499999999999999)  RGB(0.3975,                 0.14750000000000002, 0.3974999999999999)
            RGB(0.4724999999999999, 0.4724999999999999, 0.7224999999999999)      RGB(0.4999999999999999,      0.5,                    0.49999999999999983) RGB(0.6475,                 0.39749999999999996, 0.6474999999999999)
            RGB(0.4725,             0.4725,             0.7225)                  RGB(0.49999999999999994,     0.5,                    0.4999999999999999)  RGB(0.6475,                 0.3975,              0.6474999999999999)
            RGB(0.4725000000000001, 0.4725000000000001, 0.7225000000000001)      RGB(0.5,                     0.5000000000000001,     0.5)                 RGB(0.6475000000000001,     0.39750000000000013, 0.6475)
            RGB(0.7225,             0.7225,             0.9725)                  RGB(0.7499999999999999,      0.75,                   0.7499999999999999)  RGB(0.8975,                 0.6475,              0.8974999999999999)
            RGB(0.9999999999999999, 1.0,                1.0)                     RGB(0.9999999999999998,      0.9999999999999999,     0.9999999999999998)  RGB(1.0,                    1.0,                 1.0)
            RGB(1.0,                1.0,                1.0)                     RGB(1.0,                     1.0,                    1.0)                 RGB(1.0,                    1.0,                 1.0)
        ])
    end

    @testset "RGBA over RGB: $T" for T in Ts
        @test BlendColor(RGB{T}(1, 0.75, 0), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGB{T}(0.740, 0.770, 0.6) atol=max(0.0005, eps(T))
    end

    @testset "RGB over RGBA: $T" for T in Ts
        @test BlendColor(RGBA{T}(1, 0.75, 0, 0.6), RGB{T}(0, 0.5, 1)) ≈ RGBA{T}(0.340, 0.670, 1.0, 1) atol=max(0.0005, eps(T))
    end

    @testset "RGBA over RGBA: $T" for T in Ts
        @test BlendColor(RGBA{T}(1, 0.75, 0, 0.6), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGBA{T}(0.529, 0.693, 5/7, 0.84) atol=max(0.0005, eps(T))
    end
end

@testset "luminosity" begin
    @testset "RGB over RGB" begin
        test_rgb_over_rgb(BlendLuminosity, RGB{Float64}[
            RGB(0.0275, 0.0275, 0.0275) RGB(0.5000000000000001, 0.5000000000000001, 0.5000000000000001) RGB(0.8525, 0.8525, 0.8525)
            RGB(0.0275, 0.0275, 0.0275) RGB(0.5000000000000001, 0.5000000000000001, 0.5000000000000001) RGB(0.8525, 0.8525, 0.8525)
            RGB(0.0275, 0.0275, 0.0275) RGB(0.5000000000000001, 0.5000000000000001, 0.5000000000000001) RGB(0.8525, 0.8525, 0.8525)
            RGB(0.0275, 0.0275, 0.0275) RGB(0.5000000000000001, 0.5000000000000001, 0.5000000000000001) RGB(0.8525, 0.8525, 0.8525)
            RGB(0.0275, 0.0275, 0.0275) RGB(0.5000000000000001, 0.5000000000000001, 0.5000000000000001) RGB(0.8525, 0.8525, 0.8525)
            RGB(0.0275, 0.0275, 0.0275) RGB(0.5000000000000001, 0.5000000000000001, 0.5000000000000001) RGB(0.8525, 0.8525, 0.8525)
            RGB(0.0275, 0.0275, 0.0275) RGB(0.5000000000000001, 0.5000000000000001, 0.5000000000000001) RGB(0.8525, 0.8525, 0.8525)
            RGB(0.0275, 0.0275, 0.0275) RGB(0.5000000000000001, 0.5000000000000001, 0.5000000000000001) RGB(0.8525, 0.8525, 0.8525)
            RGB(0.0275, 0.0275, 0.0275) RGB(0.5000000000000001, 0.5000000000000001, 0.5000000000000001) RGB(0.8525, 0.8525, 0.8525)
        ])
    end

    @testset "RGBA over RGB: $T" for T in Ts
        @test BlendLuminosity(RGB{T}(1, 0.75, 0), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGB{T}(0.727, 0.545, 0.0) atol=max(0.0005, eps(T))
    end

    @testset "RGB over RGBA: $T" for T in Ts
        @test BlendLuminosity(RGBA{T}(1, 0.75, 0, 0.6), RGB{T}(0, 0.5, 1)) ≈ RGBA{T}(0.327, 0.445, 0.4, 1) atol=max(0.0005, eps(T))
    end

    @testset "RGBA over RGBA: $T" for T in Ts
        @test BlendLuminosity(RGBA{T}(1, 0.75, 0, 0.6), RGBA{T}(0, 0.5, 1, 0.6)) ≈ RGBA{T}(0.519, 0.532, 2/7, 0.84) atol=max(0.0005, eps(T))
    end
end


text:
# config.jl

"""
Contains configuration data.
#TODO add something

$(TYPEDEF)
"""
struct PnmlCfg

end


text:
### reconstruct procedure for Compressed Stacks

## Look for the first level of information available to reconstruct
function reconstruct!(stack::CompressedStack)
  for i in 1:stack.depth # lvl = h - i; 0 is for the fully compressed part
    lvl = stack.depth - i
    if lvl == 0
      sign = get(stack.compressed)
      return reconstruct!(stack, sign, lvl + 1)
    else
      if !isempty(stack.first_partial[lvl])
        sign = stack.first_partial[lvl][end]
        return reconstruct!(stack, sign, lvl + 1)
      elseif !isempty(stack.second_partial[lvl])
        sign = stack.second_partial[lvl][end]
        return reconstruct!(stack, sign, lvl + 1)
      end
    end
  end
end

## Reconstruct an auxiliary stack based on the signature found above
function reconstruct!(stack::CompressedStack, sign::Signature, lvl::Int)
  pos_reminder = position(stack.input)
  size = stack.space ^ (stack.depth + 1 - lvl)
  aux = CompressedStack(stack, size, sign.context,sign.first-1,sign.pos)

## Call the run! function in run.jl
  run!(aux, sign.last - sign.first)
  δ = stack.depth - aux.depth

  for i in 1:(aux.depth-1)
    stack.second_partial[δ + i] = aux.first_partial[i]
  end
  if lvl == 1
    if !isempty(aux.second_partial)
      compress!(stack, aux.compressed, aux.second_partial[1])
    else
      stack.compressed = aux.compressed
    end
  end
  stack.second_explicit = aux.first_explicit
  seek(stack.input, pos_reminder)

## Clean the memory
  aux = 0
  gc()
end


text:
using PPInterpolation, LaverySpline, Test

@testset "HaganCurveLavery" begin
    #for quintic: first pass compute cubic first derivative, second pass, use cubic b, and optimize second derivative.
    t = [0.0, 0.1, 1.0, 2.0, 3.0, 4.0, 9.0, 20.0, 30.0]
    y = [0.06, 0.06, 0.06, 0.05, 0.05, 0.05, 0.065, 0.06, 0.06]
    t = [0.0, 0.1, 1.0, 4.0, 9.0, 20.0, 30.0]
    y = [0.081, 0.081, 0.07, 0.044, 0.07, 0.04, 0.03]
    t = [0, 0.1, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 2.0, 2.5, 3.0, 3.5]
    y = [0.012, 0.008, 0.0104, 0.00995, 0.010633333333333333, 0.010849999999999999, 0.010179999999999998, 0.009816666666666665, 0.008612499999999999, 0.007789999999999998, 0.007741666666666665, 0.007635714285714285]

    # t = [0.0, 0.1, 1, 2, 3, 4, 5, 6, 7, 30.0]
    # f = [0.01, 0.01, 0.01, 0.01, 0.01, 0.03, 0.03, 0.03, 0.03, 0.02]
    # y = zeros(length(t))
    # y[1] = f[1]
    # y[2:end] = (y[1:end-1] .* t[1:end-1] + (t[2:end]-t[1:end-1]) .* f[2:end]) ./ t[2:end]
    # y[3:end-2] = f[3:end-2] ./ t[3:end-2]
    # println(y)
    yt = zeros(length(t))
    yt[1:end] = y[1:end] .* t[1:end]
    spline = makeQuinticLaverySpline(t, yt, 1e-3, 1e-2, 3)
    nSample = 1000
    for i = 1:nSample
        ti = (t[end]) * (i) / (nSample)
        val = evaluateDerivative(spline, ti)
        zval = evaluate(spline, ti) / ti
        println(ti, " QuinticLavery3 ", val, " ", zval)
    end
end


text:
import SMCExamples.Particles.Float64Particle
using SMCExamples.LinearGaussian
import SMCExamples.LinearGaussian: LGTheta, simulateLGModel, kalman,
  makeLGModel
using RNGPool
using CoupledConditionalSMC

import Statistics: mean, var

function testCouplingTime(n::Int64, N::Int64, maxit::Int64)
  theta = LGTheta(1.0, 1.0, 1.0, 1.0, 0.0, 1.0)
  ys = simulateLGModel(theta, n)

  model = makeLGModel(theta, ys)
  lM = LinearGaussian.makelM(theta)

  ccsmcio = CCSMCIO{model.particle, model.pScratch}(N, model.maxn)

  timeBS = CoupledConditionalSMC.couplingTime(model, lM, ccsmcio, :BS)
  @test timeBS > 0
  timeAS = CoupledConditionalSMC.couplingTime(model, lM, ccsmcio, :AS)
  @test timeAS > 0
  timeAT = CoupledConditionalSMC.couplingTime(model, ccsmcio)
  @test timeAT > 0

  timesBS = CoupledConditionalSMC.couplingTimes(model, lM, N, model.maxn, 100,
    :BS)
  @test all(timesBS .> 0)
end

setRNGs(12345)

@time @testset "Time test" begin
  testCouplingTime(10, 16, 100000)
  testCouplingTime(20, 32, 100000)
end


text:
## Example for using LSTNet

# Make sure all the required packages are available
cd(@__DIR__)
using Pkg; Pkg.activate(".")
Pkg.instantiate()

@info "Loading packages"
using FluxArchitectures
using Plots

# Load some sample data
@info "Loading data"
poollength = 10
horizon = 15
datalength = 500
input, target = get_data(:solar, poollength, datalength, horizon) |> gpu

# Define the network architecture
@info "Creating model and loss"
inputsize = size(input, 1)
encodersize = 10
decodersize = 10

# Define the neural net
model = DARNN(inputsize, encodersize, decodersize, poollength, 1) |> gpu

# MSE loss
function loss(x, y)
    Flux.reset!(model)
    return Flux.mse(model(x), permutedims(y))
end

# Callback for plotting the training
cb = function ()
    Flux.reset!(model)
    pred = model(input) |> permutedims |> cpu
    Flux.reset!(model)
    p1 = plot(pred, label="Predict")
    p1 = plot!(cpu(target), label="Data", title="Loss $(loss(input, target))")
    display(plot(p1))
end

# Training loop
@info "Start loss" loss = loss(input, target)
@info "Starting training"
Flux.train!(loss, Flux.params(model),Iterators.repeated((input, target), 10),
            ADAM(0.007), cb=cb)

@info "Finished"
@info "Final loss" loss = loss(input, target)


text:
module Yota

export
    grad,
    update!,
    @diffrule,
    @diffrule_kw,
    @nodiff,
    best_available_device,
    to_device,
    CPU,
    GPU

include("core.jl")

end


text:
#!/usr/bin/env julia
#using Plots, StatsPlots, Distributions, LaTeXStrings, Turing

"""
    function qsdgwas(nid, nlc, μ, σ, ϵ; maf = .1)
`Q`uick (and dirty) `s`imulation `d`ata for GWAS.
This function uses a uniform distribution of allele frequencies.
It is closer to reality than using just one unique frequency.

Here I used some smart, may be stupid later, one line code to generate genotypes.
- Ca = `rand.(Binomial.(2, p)` will generate genotypes of one individual.
- Cb = `[Ca for i in 1:nid]` then generate genotypes for all ID.
- `reduce(hcat, Cb)` will convert the results to a matrix of dimension `nid × nlc`.
"""
function qsdgwas(nid, nlc, μ, σ, ϵ; maf = .2)
    p = rand(Uniform(maf, 1-maf), nlc)              # allele frequencies,
    training_snp   = reduce(hcat, [rand.(Binomial.(2, p)) for i in 1:nid])'
    validation_snp = reduce(hcat, [rand.(Binomial.(2, p)) for i in 1:nid])'
    T = [ones(nid) training_snp]
    V = [ones(nid) validation_snp]
    b = rand(Normal(μ, σ), nlc + 1)                 # true SNP effects
    y = T * b + rand(Normal(0., ϵ), nid)            # phenotypes
    T, y, V, b, p
end

function naive_mcmc(T, y, σₑ; nit = 10_000)
    np = size(T, 2)             # number of parameters to estimate
    b = zeros(np)
    mcmcSample = zeros(nit, np)
    sp, tt = zeros(np), zeros(np)
    for i in 1:np
        tt[i] = T[:, i]'T[:, i]
        sp[i] = σₑ / sqrt(tt[i])
    end
    for i in 1:nit
        for t in 1:np
            b[t] = .0           # nullify this parameter
            yadj = y - T * b    # adjust y with every one of other parameters
            tcol = view(T, :, t)
            bhat = tcol'yadj / tt[t] # ̂b conditional on other parameters
            b[t] = rand(Normal(bhat, sp[t]))
        end
        mcmcSample[i, :] = b
    end
    mcmcSample
end

#=
function using_turing(T, y; nit = 2_000, nthreads = 8)
    # must using Turing to be valid
    @model linreg(X, y; predictors=size(X, 2)) = begin
        α ~ Normal(mean(y), 2.5std(y))
        β ~ filldist(TDist(3), predictors)
        σ ~ Exponential(1)
        y ~ MvNormal(α .+ X * β, σ)
    end
    model = linreg(T, y)
    chain = sample(model, NUTS(), MCMCThreads(), nit, nthreads)
    rst = summarystats(chain)
    rst[:, :mean][1:end-1]
end
=#

function _tst_gwas()
    nid, nlc, μ, σ, ϵ = 1000, 5, 0., .5, 1.0
    T, y, V, b, p = qsdgwas(nid, nlc, μ, σ, ϵ)
    sample = naive_mcmc(T, y, σ^2)

    # show the estimates side by side
    [b mean(sample, dims=1)' rst[:, :mean][1:1+nlc]]
end

"""
    function qsgp(nlc, nid, μ, σ, ϵ; maf = .2)
`Q`uick (and dirty) `s`imulation of (012) `g`enotypes and `p`henotypes.
The results are to evaluate SNP effects, so genotypes are locus majored.
That is, it is of `n_Loci × n_ID`.
"""
function qsgp(nlc, nid, μ, σ, ϵ; maf = .2)
    ps = rand(Uniform(maf, 1-maf), nlc) # allele frequencies
    T = reduce(hcat, [rand(Binomial(2, p), nid) for p in ps])'
    b = rand(Normal(μ, σ), nlc)
    y = T'b + rand(Normal(0., ϵ), nid)
    Int8.(T), y, b
end


text:
function write_tsp(x::Vector{Float64}, y::Vector{Float64}; dist="EUC_2D")
    n_nodes = length(x)
    @assert length(x) == length(y)

    name = randstring(10)
    filepath = joinpath(pwd(), name * ".tsp")

    open(filepath, "w") do io
        write(io, "NAME: $(name)\n")
        write(io, "TYPE: TSP\n")
        write(io, "COMMENT: $(name)\n")
        write(io, "DIMENSION: $(n_nodes)\n")
        write(io, "EDGE_WEIGHT_TYPE: $(dist)\n")
        write(io, "EDGE_WEIGHT_FORMAT: FUNCTION \n")
        write(io, "NODE_COORD_TYPE: TWOD_COORDS \n")
        write(io, "NODE_COORD_SECTION\n")
        for i in 1:n_nodes
            write(io, "$i $(x[i]) $(y[i])\n")
        end
        write(io, "EOF\n")
    end

    return name, filepath
end

function write_tsp(dist_mtx::Matrix{Int})
    n_nodes = size(dist_mtx, 1)

    name = randstring(10)
    filepath = joinpath(pwd(), name * ".tsp")

    lower_diag_row = Int[]
    for i in 1:n_nodes
        for j in 1:i
            push!(lower_diag_row, dist_mtx[i, j])
        end
    end
    buf = 10
    n_rows = length(lower_diag_row) / buf |> ceil |> Int
    rows = String[]
    for i in 1:n_rows
        s = buf * (i-1) + 1
        t = min(buf * i, length(lower_diag_row))
        push!(rows, join(lower_diag_row[s:t], " "))
    end

    open(filepath, "w") do io
        write(io, "NAME: $(name)\n")
        write(io, "TYPE: TSP\n")
        write(io, "COMMENT: $(name)\n")
        write(io, "DIMENSION: $(n_nodes)\n")
        write(io, "EDGE_WEIGHT_TYPE: EXPLICIT\n")
        write(io, "EDGE_WEIGHT_FORMAT: LOWER_DIAG_ROW \n")
        write(io, "EDGE_WEIGHT_SECTION\n")
        for r in rows
            write(io, "$r\n")
        end
        write(io, "EOF\n")
    end

    return name, filepath
end

function write_atsp(dist_mtx::Matrix{Int})
    n_nodes = size(dist_mtx, 1)

    name = randstring(10)
    filepath = joinpath(pwd(), name * ".tsp")

    rows = String[]
    for i in 1:n_nodes
        push!(rows, join(dist_mtx[i, :], " "))
    end

    open(filepath, "w") do io
        write(io, "NAME: $(name)\n")
        write(io, "TYPE: ATSP\n")
        write(io, "COMMENT: $(name)\n")
        write(io, "DIMENSION: $(n_nodes)\n")
        write(io, "EDGE_WEIGHT_TYPE: EXPLICIT\n")
        write(io, "EDGE_WEIGHT_FORMAT: FULL_MATRIX \n")
        write(io, "EDGE_WEIGHT_SECTION\n")
        for r in rows
            write(io, "$r\n")
        end
        write(io, "EOF\n")
    end

    return name, filepath
end

function tour_length(tour, M)
    n_nodes = length(tour)
    len = 0
    for i in 1:n_nodes
        j = i + 1
        if i == n_nodes
            j = 1
        end

        len += M[tour[i], tour[j]]
    end

    return len
end

function filepath(name, ext)
    return joinpath(pwd(), "$(name).$(ext)")
end

function write_par(name; kwargs...)
    par_filepath = filepath(name, "par")
    out_filepath = filepath(name, "out")

    open(par_filepath, "w") do io
        write(io, "PROBLEM_FILE = $(name).tsp\n")
        write(io, "TOUR_FILE = $(name).out\n")
        if length(kwargs) > 0
            for (key, val) in kwargs
                write(io, key, " = ", string(val), "\n")
            end
        end
    end

    return name
end

function read_output(name)
    output = readlines(filepath(name, "out"))
    idx1 = findfirst(x -> x=="TOUR_SECTION", output) + 1
    idx2 = findfirst(x -> x=="-1", output) - 1
    tour = parse.(Int, output[idx1:idx2])

    len = parse(Int, split(output[2])[end])

    return tour, len
end

function cleanup(name)
    exts = ["par", "tsp", "out"]
    rm.(filepath.(name, exts))
end


text:
module PGQueryTests

using AbstractTrees
using JSON
using PGQuery
using Test

const DEFAULT_SHOW_JSON = false

function try_query(qry; show_json::Bool=DEFAULT_SHOW_JSON)
    println("query:")
    println(qry)
    if show_json
        parse_query_json(qry) do parsed_qry_json_str
            parsed_qry_json = JSON.parse(parsed_qry_json_str)
            println("parsed_qry_json:")
            JSON.print(parsed_qry_json, 2)
        end
    end
    parse_query(qry) do parsed_qry
        AbstractTrees.print_tree(parsed_qry, 999)
    end
    println("-----------------------------------------------")
end

function sql_files(dir)
    [file for file in readdir(dir) if endswith(file, ".sql")]
end

@testset "simple SQL parsing" begin
    try_query("SELECT 1")
    try_query("SELECT 'abc'")
    try_query("SELECT B'101010111'")
    try_query("SELECT 3.14")
    try_query("SELECT NULL")
    try_query("SELECT a, t.b from tbl t")
    for imdb_job_qry in sql_files(joinpath(@__DIR__, "imdb_job"))
        println("query name: $imdb_job_qry")
        try_query(read(joinpath(@__DIR__, "imdb_job", imdb_job_qry), String))
    end
    for tpch_qry in sql_files(joinpath(@__DIR__, "tpch"))
        println("query name: $tpch_qry")
        try_query(read(joinpath(@__DIR__, "tpch", tpch_qry), String))
    end
    @test_throws SQLParserException parse_query("INSERT FROM DOES NOT WORK")
end

end


text:
module AnalyticalSolutions

include("AbstractTypes.jl")
include("prototype.jl")

using Reexport #reexports must be placed after includes
@reexport using .AbstractType #passes thru all exports from specified module

#=
using .AbstractType

struct g  <: GoverningDiffEqn end
struct an <: Analysis end
struct ap <: Approach end
struct d  <: Dimensions end
struct m  <: Material end
struct b  <: BoundaryCondition end
struct s  <: Source end
struct r  <: ResultSpec end

analyticalSolution(g(),an(),ap(),d(),m(),b(),s(),r())
=#

#=
modules:
--------
- Structural mechanics
  - Plate
    - Thin
    - Thick
  - Membrane
    - Plane stress
    - Plane strain
  - Shell
  - Frame
- CFD (?)
- Heat transport (?)
- Electrical (?)
  - Circuit
- Porous media (?)
- etc.

Definitions:
------------
- Abstract types
- Fall thru for not implemented function

=#

end


text:
struct GradientDescent{T1,TP} <: QuasiNewton{T1}
    approx::T1
    P::TP
end
GradientDescent() = GradientDescent(Direct(), nothing)
GradientDescent(m) = GradientDescent(m, nothing)
hasprecon(::GradientDescent{<:Any,<:Nothing}) = NoPrecon()
hasprecon(::GradientDescent{<:Any,<:Any}) = HasPrecon()

summary(gr::GradientDescent) = "Gradient Descent"
update!(scheme::GradientDescent, A, s, y) = A
update(scheme::GradientDescent, A, s, y) = A

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
"""
    console()

Lorem ipsum dolor sit amet.
"""
function console(args::Dict)
  console()
end

function console()

  target_name = bump()

  Base.run(`julia -ie "using $(target_name)"`)

end


text:
#!/usr/bin/env julia

using Pkg
using Dates
using DelimitedFiles
using Distributed
using Statistics
using MIToS.Utils.Scripts

Args = parse_commandline(
    # TO DO ----------------------------------------------------------------------
    ["--format", "-f"],
    Dict(
        :help => "Format of the MSA: Stockholm, Raw or FASTA",
        :arg_type => String,
        :default => "Stockholm"
    ),
    ["--savelist", "-s"],
    Dict(
        :help => "Create and pidlist.csv file with the percentage identity for each pairwise comparison.",
        :action => :store_true
    ),
    # Keywords...
    description="""
    Calculates the percentage identity between all the sequences of an MSA and creates an *.pidstats.csv file with:
    The number of columns and sequences. The mean, standard deviation, median, minimum and maximum values and first and third quantiles of the percentage identity.
    It could also create and pidlist.csv file with the percentage identity for each pairwise comparison.
    """,
    output=".pidstats.csv"
    # ----------------------------------------------------------------------------
    )

set_parallel(Args["parallel"])

@everywhere begin

    const args = remotecall_fetch(()->Args,1)

    import MIToS.Utils.Scripts: script

    # TO DO ----------------------------------------------------------------------
    using MIToS.MSA
    using PairwiseListMatrices
    # ----------------------------------------------------------------------------

    function script(input::Union{Base.LibuvStream,  AbstractString},
                    args,
                    fh_out::Union{Base.LibuvStream, IO})
        # TO DO ------------------------------------------------------------------
        println(fh_out, "# MIToS ", Pkg.installed()["MIToS"], " PercentIdentity.jl ", now())
        println(fh_out, "# used arguments:")
        for (key, value) in args
            println(fh_out, "# \t", key, "\t\t", value)
        end
        form = string(args["format"])
        if form == "Stockholm"
            msa = readorparse(input, Stockholm)
        elseif form == "FASTA"
            msa = readorparse(input, FASTA)
        elseif form == "Raw"
            msa = readorparse(input, Raw)
        else
            throw(ErrorException("--format should be Stockholm, Raw or FASTA."))
        end
        savelist = args["savelist"]
        println(fh_out, "ncol,nseq,mean,std,min,firstq,median,thirdq,max")
        plm = percentidentity(msa, Float16)
        mean_pid = mean(getlist(plm))
        min_pid, max_pid = extrema(getlist(plm))
        first, med, third = quantile(getlist(plm), Float64[0.25,0.5,0.75])
        println(fh_out, size(msa, 2), ",", size(msa, 1), ",", mean_pid, ",", stdm(getlist(plm), mean_pid), ",", min_pid, ",", first, ",", med, ",", third, ",", max_pid)
        flush(fh_out)
        if savelist
            writedlm("pidlist.csv", to_table(plm, diagonal=false), ',')
        end
        # ------------------------------------------------------------------------
    end

end

runscript(args)


text:
module Visualize_TextFonts

using Modia3D

@assembly DisplayFontFamily(fontFamily; cs=0.1) begin
   frame0         = Modia3D.Object3D()
   fontNormal     = Modia3D.Object3D(frame0, Modia3D.TextShape(fontFamily; alignment = Modia3D.Left, offset=[0.0,0.0,10*cs],
                                                            font = Modia3D.Font(fontFamily=fontFamily, bold=false, italic=false, charSize=cs)))
   fontItalic     = Modia3D.Object3D(frame0, Modia3D.TextShape(fontFamily * " italic"; alignment = Modia3D.Left, offset=[0.0,0.0,8.5*cs],
                                                            font = Modia3D.Font(fontFamily=fontFamily, bold=false, italic=true, charSize=cs)))
   fontBold       = Modia3D.Object3D(frame0, Modia3D.TextShape(fontFamily * " bold"; alignment = Modia3D.Left, offset=[0.0,0.0,7*cs],
                                                            font = Modia3D.Font(fontFamily=fontFamily, bold=true, italic=false, charSize=cs)))
   fontBoldItalic = Modia3D.Object3D(frame0, Modia3D.TextShape(fontFamily * " bold + italic"; alignment = Modia3D.Left, offset=[0.0,0.0,5.5*cs],
                                                            font = Modia3D.Font(fontFamily=fontFamily, bold=true, italic=true, charSize=cs)))
   fontNormalRed  = Modia3D.Object3D(frame0, Modia3D.TextShape(fontFamily * " red color"; alignment = Modia3D.Left, offset=[0.0,0.0,4*cs],
                                                            font = Modia3D.Font(fontFamily=fontFamily, color="Red", charSize=cs)))
end

@assembly DisplayFonts begin
   Arial         = DisplayFontFamily("Arial")
   ArialNarrow   = DisplayFontFamily("ArialNarrow")
   TimesNewRoman = DisplayFontFamily("TimesNewRoman")
   Verdana       = DisplayFontFamily("Verdana")
   CourierNew    = DisplayFontFamily("CourierNew")
   FreeSans      = DisplayFontFamily("FreeSans")
   FreeSerif     = DisplayFontFamily("FreeSerif")

   world         = Modia3D.Object3D()

   Modia3D.connect(world, Arial.frame0        ; r=[0.0,-4.5,0.0])
   Modia3D.connect(world, ArialNarrow.frame0  ; r=[0.0,-3.5,0.0])
   Modia3D.connect(world, TimesNewRoman.frame0; r=[0.0,-2.0,0.0])
   Modia3D.connect(world, Verdana.frame0      ; r=[0.0,-0.5,0.0])
   Modia3D.connect(world, CourierNew.frame0   ; r=[0.0, 1.0,0.0])
   Modia3D.connect(world, FreeSans.frame0     ; r=[0.0, 2.5,0.0])
   Modia3D.connect(world, FreeSerif.frame0    ; r=[0.0, 4.0,0.0])
end
Modia3D.visualizeAssembly!( DisplayFonts() )

println("... success of Visualize_TextFonts.jl!")

end


text:
monicderivative(p, i::Int) = derivative(p, degree(p)-i)*factorial(i)/factorial(degree(p))

isapproxreal(x::Real) = true
isapproxreal(x::Complex{T}; atol::Real=1e-7) where {T} = abs(imag(x)) < atol
realroots(poly) = real.(filter(isapproxreal, _ROOTS[](poly)))

# define a function that returns a Plots.Shape
# code from https://discourse.julialang.org/t/how-to-draw-a-rectangular-region-with-plots-jl/1719/4
rectangle(x1, y1, x2, y2) = ([x1,x2,x2,x1], [y1,y1,y2,y2])
function ithrectangles(sortrealroots, i;
                       mininf=-5, inf=5, w=.5)
    xs = [mininf; sortrealroots]
    xs2 = [sortrealroots; inf]
    return [rectangle(x1, i-w, x2, i+w) for (x1,x2) in zip(xs, xs2)]
end

function ithcolors(_iseven::Bool, nboxes, gray, black)
    grays = _iseven ? iseven.(1:nboxes) : isodd.(1:nboxes)
    blacks = (!).(grays)
    return gray*grays + black*blacks
end

function find_bounds(allsortvectors)
    nonempty_vectors = filter((!)∘isempty, allsortvectors)
    mins, maxs = first.(nonempty_vectors), last.(nonempty_vectors)
    return minimum(mins), maximum(maxs)
end