text:
using Gen

@dist labeled_cat(labels, probs) = labels[categorical(probs)]

@gen function ve_backwards_sampler(latents, fg, ve_result)
    N = length(latents)
    values = Vector{Any}(undef, N)
    for addr in reverse(ve_result.elimination_order)
        var_idx = addr_to_idx(fg, addr)
        fg = ve_result.intermediate_fgs[var_idx]
        dist = conditional_dist(fg, values, addr)
        value = ({addr} ~ labeled_cat(latents[addr].domain, dist))
        values[var_idx] = value
    end
end

# NOTE: we could stage the computation better, so less is done within the
# generative function at runtime
# NOTE: we could in principle emit static IR instead of DML code

"""
    sampler::GenerativeFunction = generate_backwards_sampler_fixed_trace(trace, addrs)

Generate a generative function that takes no arguments that samples from the conditional distribution on the given addresses.

Only applies to traces of generative functions constructing purely using SML and combinators (currently Map and Unfold).
The addresses must be discrete random choices within finite support and must not affect the control flow in the trace.
The sampler is generated using variable elimination followed by backwards sampling, where the order of the provided addresses defines the elimination order.
The sampler takes no arguments and is specialized to the conditional distribution for the given trace.
"""
function generate_backwards_sampler_fixed_trace(
        trace::Union{Gen.StaticIRTrace,Gen.VectorTrace{Gen.MapType},Gen.VectorTrace{Gen.UnfoldType}}, addrs)
    (_, latents, observations) = factor_graph_analysis(trace, addrs)
    fg = compile_trace_to_factor_graph(trace, latents, observations)
    ve_result = variable_elimination(fg, addrs)
    @gen function sampler()
        {*} ~ ve_backwards_sampler(latents, fg, ve_result)
    end
    return sampler
end

"""
    sampler::GenerativeFunction = generate_backwards_sampler_fixed_trace(
        trace, addrs, latents::Dict{Any,Latent}, observations::Dict{Any,Observation})

Generate a generative function that takes no arguments that samples from the conditional distribution on the given addresses.

The addresses must be discrete random choices within finite support and must not affect the control flow in the trace.
The sampler is generated using variable elimination followed by backwards sampling, where the order of the provided addresses defines the elimination order.
The sampler takes no arguments and is specialized to the conditional distribution for the given trace.
"""
function generate_backwards_sampler_fixed_trace(
        trace, addrs, latents::Dict{Any,Latent}, observations::Dict{Any,Observation})
    fg = compile_trace_to_factor_graph(trace, latents, observations)
    ve_result = variable_elimination(fg, addrs)
    @gen function sampler()
        {*} ~ ve_backwards_sampler(latents, fg, ve_result)
    end
    return sampler
end

"""
    sampler::GenerativeFunction = generate_backwards_sampler(trace, addrs)

Generate a generative function that takes a trace argument that samples from the conditional distribution on the given addresses.

Only applies to traces of generative functions constructing purely using SML and combinators (currently Map and Unfold).
The addresses must be discrete random choices within finite support and must not affect the control flow in the trace.
The sampler is generated using variable elimination followed by backwards sampling, where the order of the provided addresses defines the elimination order.
The sampler is specialized to traces with the same control flow path, but not necessarily the same values, as the trace provided to the generation function.
"""
function generate_backwards_sampler_fixed_structure(
        trace::Union{Gen.StaticIRTrace,Gen.VectorTrace{Gen.MapType},Gen.VectorTrace{Gen.UnfoldType}}, addrs)
    (_, latents, observations) = factor_graph_analysis(trace, addrs)
    @gen function sampler(trace)
        fg = compile_trace_to_factor_graph(trace, latents, observations)
        ve_result = variable_elimination(fg, addrs)
        {*} ~ ve_backwards_sampler(latents, fg, ve_result)
    end
    return sampler
end

@gen function backwards_sampler_dml(trace, addrs, latents::Dict{Any,Latent}, observations::Dict{Any,Observation})
    fg = compile_trace_to_factor_graph(trace, latents, observations)
    ve_result = variable_elimination(fg, addrs)
    {*} ~ ve_backwards_sampler(latents, fg, ve_result)
end

@gen function backwards_sampler_sml(trace, addrs)
    (_, latents, observations) = factor_graph_analysis(trace, addrs)
    fg = compile_trace_to_factor_graph(trace, latents, observations)
    ve_result = variable_elimination(fg, addrs)
    {*} ~ ve_backwards_sampler(latents, fg, ve_result)
end

export generate_backwards_sampler_fixed_structure
export generate_backwards_sampler_fixed_trace
export backwards_sampler_dml, backwards_sampler_sml


text:
"""
  A test of the test system itself
"""
function test_empty()
  @testset "Testing empty..." begin

    @testset "Testing if 1 == 1" begin
      @test ( 1 )== 1
    end

  end
end

#test_empty()
add_func1!(EulerTests, test_empty, [TAG_SHORTTEST])

ext: jl


text:
using Statistics
using StatsBase
using Elefridge
using PyCall
xr = pyimport("xarray")
using JLD

path = "/Users/milan/cams/unstructured/"
filelist = filter(x->endswith(x,".grib"),readdir(path))

nvars = length(filelist)
nmethods = 4        # Lin24,Lin16,Log16lin,Log16log
nstats = 7          # min,10%,25%,50%,75%,90%,max

meanerror = fill(0f0,nmethods,nvars)
linerror = fill(0f0,nmethods,nvars,nstats)
logerror = fill(0f0,nmethods,nvars,nstats)
allvars = Array{String,1}(undef,nvars)

for i in 1:nvars

    file = filelist[i]
    allvars[i] = split(split(file,"_")[end],".")[1]

    println("---")
    println("Compressing $(allvars[i])")

    D = xr.open_dataarray(joinpath(path,file),engine="cfgrib").data

    if any(D .<= 0f0)
        println("Negative/zero entries found.")
        D = D[D.>0f0]      # remove non-positive entries
    end

    if any(isnan.(D))
        println("NaNs found.")
        D = D[.~isnan(D)]   # remove NaNs
    end

    Dlin16 = Array(LinQuant16Array(D))
    Dlin24 = Array(LinQuant24Array(D))
    Dlog16log = Array(LogQuant16Array(D))
    Dlog16lin = Array(LogQuant24Array(D))

    # normalised mean error
    Dmean = mean(abs.(D))
    meanerror[1,i] = abs(mean(D - Dlin16)) / Dmean
    meanerror[2,i] = abs(mean(D - Dlin24)) / Dmean
    meanerror[3,i] = abs(mean(D - Dlog16log)) / Dmean
    meanerror[4,i] = abs(mean(D - Dlog16lin)) / Dmean

    # absolute error
    for (j,Dr) in enumerate([Dlin16,Dlin24,Dlog16log,Dlog16lin])
        E = sort(vec(abs.(D-Dr) / Dmean))
        linerror[j,i,1],linerror[j,i,7] = E[1],E[end]
        linerror[j,i,2] = quantile(E,0.1,sorted=true)
        linerror[j,i,3] = quantile(E,0.25,sorted=true)
        linerror[j,i,4] = quantile(E,0.5,sorted=true)
        linerror[j,i,5] = quantile(E,0.75,sorted=true)
        linerror[j,i,6] = quantile(E,0.9,sorted=true)
    end

    # decimal error
    for (j,Dr) in enumerate([Dlin16,Dlin24,Dlog16log,Dlog16lin])
        E = sort(vec(abs.(log10.(abs.(D./Dr)))))    # avoid non-finite entries with abs
        logerror[j,i,1],logerror[j,i,7] = E[1],E[end]
        logerror[j,i,2] = quantile(E,0.1,sorted=true)
        logerror[j,i,3] = quantile(E,0.25,sorted=true)
        logerror[j,i,4] = quantile(E,0.5,sorted=true)
        logerror[j,i,5] = quantile(E,0.75,sorted=true)
        logerror[j,i,6] = quantile(E,0.9,sorted=true)
    end
end

save("/Users/milan/cams/error/linvslog_unstructured.jld",
        "allvars",allvars,
        "meanerror",meanerror,
        "linerror",linerror,
        "logerror",logerror)


text:
using DataFrames
using Distributions
using MixedModels
using PooledArrays
using StableRNGs
using Tables
using Test

using GLM: Link
using MixedModels: dataset

include("modelcache.jl")

@testset "GLMM from MixedModel" begin
    f = first(gfms[:contra])
    d = dataset(:contra)
    @test MixedModel(f, d, Bernoulli()) isa GeneralizedLinearMixedModel
    @test MixedModel(f, d, Bernoulli(), ProbitLink()) isa GeneralizedLinearMixedModel
end

@testset "contra" begin
    contra = dataset(:contra)
    thin = 5
    gm0 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); fast=true, progress=false, thin)
    fitlog = gm0.optsum.fitlog
    @test length(fitlog) == (div(gm0.optsum.feval, thin) + 1) # for the initial value
    @test first(fitlog) == (gm0.optsum.initial, gm0.optsum.finitial)
    @test gm0.lowerbd == zeros(1)
    @test isapprox(gm0.θ, [0.5720734451352923], atol=0.001)
    @test !issingular(gm0)
    @test issingular(gm0, [0])
    @test isapprox(deviance(gm0), 2361.657188518064, atol=0.001)
    # the first 9 BLUPs -- I don't think there's much point in testing all 102
    blups = [-0.5853637711570235, -0.9546542393824562, -0.034754249031292345, # values are the same but in different order
              0.2894692928724314, 0.6381376605845264, -0.2513134928312374,
              0.031321447845204374, 0.10836110432794945, 0.24632286640099466]
    @test only(ranef(gm0))[1:9] ≈ blups atol=1e-4
    retbl = raneftables(gm0)
    @test isone(length(retbl))
    @test isa(retbl, NamedTuple)
    @test Tables.istable(only(retbl))
    @test !dispersion_parameter(gm0)
    @test dispersion(gm0, false) == 1
    @test dispersion(gm0, true) == 1
    @test sdest(gm0) === missing
    @test varest(gm0) === missing
    @test gm0.σ === missing
    @test Distribution(gm0) == Distribution(gm0.resp)
    @test Link(gm0) == Link(gm0.resp)

    gm1 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); progress=false);
    @test isapprox(gm1.θ, [0.573054], atol=0.005)
    @test lowerbd(gm1) == vcat(fill(-Inf, 7), 0.)
    @test isapprox(deviance(gm1), 2361.54575, rtol=0.00001)
    @test isapprox(loglikelihood(gm1), -1180.77288, rtol=0.00001)

    @test dof(gm0) == length(gm0.β) + length(gm0.θ)
    @test nobs(gm0) == 1934
    refit!(gm0; fast=true, nAGQ=7, progress=false)
    @test isapprox(deviance(gm0), 2360.9838, atol=0.001)
    gm1 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); nAGQ=7, progress=false)
    @test isapprox(deviance(gm1), 2360.8760, atol=0.001)
    @test gm1.β == gm1.beta
    @test gm1.θ == gm1.theta
    gm1y = gm1.y
    @test length(gm1y) == size(gm1.X, 1)
    @test eltype(gm1y) == eltype(gm1.X)
    @test gm1y == (MixedModels.dataset(:contra).use .== "Y")
    @test response(gm1) == gm1y
    @test !islinear(gm1)
    @test :θ in propertynames(gm0)

    @testset "GLMM rePCA" begin
        @test length(MixedModels.PCA(gm0)) == 1
        @test length(MixedModels.rePCA(gm0)) == 1
        @test length(gm0.rePCA) == 1
    end
    # the next three values are not well defined in the optimization
    #@test isapprox(logdet(gm1), 75.7217, atol=0.1)
    #@test isapprox(sum(abs2, gm1.u[1]), 48.4747, atol=0.1)
    #@test isapprox(sum(gm1.resp.devresid), 2237.349, atol=0.1)
    show(IOBuffer(), gm1)
    show(IOBuffer(), BlockDescription(gm0))

    gm_slope = fit(MixedModel, gfms[:contra][2], contra, Bernoulli(); progress=false);
    @test !issingular(gm_slope)
    @test issingular(gm_slope, zeros(5))

end

@testset "cbpp" begin
    cbpp = dataset(:cbpp)
    gm2 = fit(MixedModel, first(gfms[:cbpp]), cbpp, Binomial(); wts=float(cbpp.hsz), progress=false)
    @test weights(gm2) == cbpp.hsz
    @test deviance(gm2,true) ≈ 100.09585619892968 atol=0.0001
    @test sum(abs2, gm2.u[1]) ≈ 9.723054788538546 atol=0.0001
    @test logdet(gm2) ≈ 16.90105378801136 atol=0.0001
    @test isapprox(sum(gm2.resp.devresid), 73.47174762237978, atol=0.001)
    @test isapprox(loglikelihood(gm2), -92.02628186840045, atol=0.001)
    @test !dispersion_parameter(gm2)
    @test dispersion(gm2, false) == 1
    @test dispersion(gm2, true) == 1
    @test sdest(gm2) === missing
    @test varest(gm2) === missing
    @test gm2.σ === missing

    @testset "GLMM refit" begin
        gm2r = deepcopy(gm2)
        @test_throws ArgumentError fit!(gm2r; progress=false)

        refit!(gm2r; fast=true, progress=false)
        @test length(gm2r.optsum.final) == 1
        @test gm2r.θ ≈ gm2.θ atol=1e-3

        # swapping successes and failures to give us the same model
        # but with opposite signs. healthy ≈ 1 - response(gm2r)
        # but defining it in terms of the original values avoids some
        # nasty floating point issues
        healthy = @. (cbpp.hsz - cbpp.incid) / cbpp.hsz
        refit!(gm2r, healthy; fast=false, progress=false)
        @test length(gm2r.optsum.final) == 5
        @test gm2r.β ≈ -gm2.β atol=1e-3
        @test gm2r.θ ≈ gm2.θ atol=1e-3
    end
end

@testset "verbagg" begin
    gm3 = fit(MixedModel, only(gfms[:verbagg]), dataset(:verbagg), Bernoulli(); progress=false)
    @test deviance(gm3) ≈ 8151.40 rtol=1e-5
    @test lowerbd(gm3) == vcat(fill(-Inf, 6), zeros(2))
    @test fitted(gm3) == predict(gm3)
    # these two values are not well defined at the optimum
    @test isapprox(sum(x -> sum(abs2, x), gm3.u), 273.29646346940785, rtol=1e-3)
    @test sum(gm3.resp.devresid) ≈ 7156.550941446312 rtol=1e-4
end

@testset "grouseticks" begin
    center(v::AbstractVector) = v .- (sum(v) / length(v))
    grouseticks = DataFrame(dataset(:grouseticks))
    grouseticks.ch = center(grouseticks.height)
    gm4 = fit(MixedModel, only(gfms[:grouseticks]), grouseticks, Poisson(); fast=true, progress=false)  # fails in pirls! with fast=false
    @test isapprox(deviance(gm4), 851.4046, atol=0.001)
    # these two values are not well defined at the optimum
    #@test isapprox(sum(x -> sum(abs2, x), gm4.u), 196.8695297987013, atol=0.1)
    #@test isapprox(sum(gm4.resp.devresid), 220.92685781326136, atol=0.1)
    @test !dispersion_parameter(gm4)
    @test dispersion(gm4, false) == 1
    @test dispersion(gm4, true) == 1
    @test sdest(gm4) === missing
    @test varest(gm4) === missing
    @test gm4.σ === missing
end

@testset "goldstein" begin # from a 2020-04-22 msg by Ben Goldstein to R-SIG-Mixed-Models
    goldstein = (
        group = PooledArray(repeat(string.('A':'J'), outer=10)),
        y = [
            83, 3, 8, 78, 901, 21, 4, 1, 1, 39,
            82, 3, 2, 82, 874, 18, 5, 1, 3, 50,
            87, 7, 3, 67, 914, 18, 0, 1, 1, 38,
            86, 13, 5, 65, 913, 13, 2, 0, 0, 48,
            90, 5, 5, 71, 886, 19, 3, 0, 2, 32,
            96, 1, 1, 87, 860, 21, 3, 0, 1, 54,
            83, 2, 4, 70, 874, 19, 5, 0, 4, 36,
            100, 11, 3, 71, 950, 21, 6, 0, 1, 40,
            89, 5, 5, 73, 859, 29, 3, 0, 2, 38,
            78, 13, 6, 100, 852, 24, 5, 0, 1, 39
            ],
        )
    gform = @formula(y ~ 1 + (1|group))
    m1 = GeneralizedLinearMixedModel(gform, goldstein, Poisson())
    @test !isfitted(m1)
    fit!(m1; progress=false)
    @test isfitted(m1)
    @test deviance(m1) ≈ 193.5587302384811 rtol=1.e-5
    @test only(m1.β) ≈ 4.192196439077657 atol=1.e-5
    @test only(m1.θ) ≈ 1.838245201739852 atol=1.e-5
    m11 = fit(MixedModel, gform, goldstein, Poisson(); nAGQ=11, progress=false)
    @test deviance(m11) ≈ 193.51028088736842 rtol=1.e-5
    @test only(m11.β) ≈ 4.192196439077657 atol=1.e-5
    @test only(m11.θ) ≈ 1.838245201739852 atol=1.e-5
end

@testset "dispersion" begin

    form = @formula(reaction ~ 1 + days + (1+days|subj))
    dat = dataset(:sleepstudy)

    @test_logs (:warn, r"dispersion parameter") GeneralizedLinearMixedModel(form, dat, Gamma())
    @test_logs (:warn, r"dispersion parameter") GeneralizedLinearMixedModel(form, dat, InverseGaussian())
    @test_logs (:warn, r"dispersion parameter") GeneralizedLinearMixedModel(form, dat, Normal(), SqrtLink())

    # notes for future tests when GLMM with dispersion works
    # @test dispersion_parameter(gm)
    # @test dispersion(gm, false) == val
    # @test dispersion(gm, true) == val
    # @test sdest(gm) == dispersion(gm, false) == gm.σ
    # @test varest(gm) == dispersion(gm, true)

end

@testset "mmec" begin
    # Data on "Malignant melanoma in the European community" from the mlmRev package for R
    # The offset of log.(expected) is to examine the ratio of observed to expected, based on population
    mmec = dataset(:mmec)
    mmform = @formula(deaths ~ 1 + uvb + (1|region))
    gm5 = fit(MixedModel, mmform, mmec, Poisson(); offset=log.(mmec.expected), nAGQ=11, progress=false)
    @test isapprox(deviance(gm5), 655.2533533016059, atol=5.e-3)
    @test isapprox(first(gm5.θ), 0.4121684550775567, atol=1.e-3)
    @test isapprox(first(gm5.β), -0.13860166843315044, atol=1.e-3)
    @test isapprox(last(gm5.β), -0.034414458364713504, atol=1.e-3)
end


text:
#=

 Unit tests from libieeep1788 for interval constructors
 (Original author: Marco Nehmeier)
 converted into portable ITL format by Oliver Heimlich.

 Copyright 2013-2015 Marco Nehmeier (nehmeier@informatik.uni-wuerzburg.de)
 Copyright 2015-2017 Oliver Heimlich (oheim@posteo.de)

 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.

=#
#Language imports

#Test library imports
using Test

#Arithmetic library imports
using IntervalArithmetic

#Preamble
setprecision(53)
setprecision(Interval, Float64)
setrounding(Interval, :tight)
# Set full format, and show decorations
@format full

@testset "minimal_nums_to_interval_test" begin
    @test Interval(-1.0, 1.0) == Interval(-1.0, 1.0)
    @test Interval(-Inf, 1.0) == Interval(-Inf, 1.0)
    @test Interval(-1.0, Inf) == Interval(-1.0, Inf)
    @test Interval(-Inf, Inf) == Interval(-Inf, Inf)
    @test Interval(NaN, NaN) == ∅
    @test Interval(1.0, -1.0) == ∅
    @test Interval(-Inf, -Inf) == ∅
    @test Interval(Inf, Inf) == ∅
end

@testset "minimal_nums_to_decorated_interval_test" begin
    @test DecoratedInterval(-1.0, 1.0) == DecoratedInterval(Interval(-1.0, 1.0), com)
    @test DecoratedInterval(-Inf, 1.0) == DecoratedInterval(Interval(-Inf, 1.0), dac)
    @test DecoratedInterval(-1.0, Inf) == DecoratedInterval(Interval(-1.0, Inf), dac)
    @test DecoratedInterval(-Inf, Inf) == DecoratedInterval(Interval(-Inf, Inf), dac)
    @test DecoratedInterval(NaN, NaN) == nai()
    @test DecoratedInterval(1.0, -1.0) == nai()
    @test DecoratedInterval(-Inf, -Inf) == nai()
    @test DecoratedInterval(Inf, Inf) == nai()
end

@testset "minimal_text_to_interval_test" begin
    @test @interval("[ Empty  ]") == ∅
    @test @interval("[ Empty  ]_trv") == ∅
    @test @interval("[  ]") == ∅
    @test @interval("[  ]_trv") == ∅
    @test @interval("[,]") == Interval(-Inf, Inf)
    @test @interval("[,]_trv") == ∅
    @test @interval("[ entire  ]") == Interval(-Inf, Inf)
    @test @interval("[ ENTIRE ]_dac") == ∅
    @test @interval("[ ENTIRE ]") == Interval(-Inf, Inf)
    @test @interval("[ -inf , INF  ]") == Interval(-Inf, Inf)
    @test @interval("[ -inf, INF ]_def") == ∅
    @test @interval("[-1.0,1.0]") == Interval(-1.0, 1.0)
    @test @interval("[  -1.0  ,  1.0  ]") == Interval(-1.0, 1.0)
    @test @interval("[  -1.0  , 1.0]") == Interval(-1.0, 1.0)
    @test @interval("[-1,]") == Interval(-1.0, Inf)
    @test @interval("[-1.0, +inf]") == Interval(-1.0, Inf)
    @test @interval("[-1.0, +infinity]") == Interval(-1.0, Inf)
    @test @interval("[-Inf, 1.000 ]") == Interval(-Inf, 1.0)
    @test @interval("[-Infinity, 1.000 ]") == Interval(-Inf, 1.0)
    @test @interval("[1.0E+400 ]") == Interval(0x1.fffffffffffffp+1023, Inf)
    @test @interval("[ -4/2, 10/5 ]") == Interval(-2.0, 2.0)
    @test @interval("[ -1/10, 1/10 ]") == Interval(-0.1, 0.1)
    @test @interval("0.0?") == Interval(-0.05, 0.05)
    @test @interval("0.0?u") == Interval(0.0, 0.05)
    @test @interval("0.0?d") == Interval(-0.05, 0.0)
    @test @interval("2.5?") == Interval(0x1.3999999999999p+1, 0x1.4666666666667p+1)
    @test @interval("2.5?u") == Interval(2.5, 0x1.4666666666667p+1)
    @test @interval("2.5?d") == Interval(0x1.3999999999999p+1, 2.5)
    @test @interval("0.000?5") == Interval(-0.005, 0.005)
    @test @interval("0.000?5u") == Interval(0.0, 0.005)
    @test @interval("0.000?5d") == Interval(-0.005, 0.0)
    @test @interval("2.500?5") == Interval(0x1.3f5c28f5c28f5p+1, 0x1.40a3d70a3d70bp+1)
    @test @interval("2.500?5u") == Interval(2.5, 0x1.40a3d70a3d70bp+1)
    @test @interval("2.500?5d") == Interval(0x1.3f5c28f5c28f5p+1, 2.5)
    @test @interval("0.0??") == Interval(-Inf, Inf)
    @test @interval("0.0??u") == Interval(0.0, Inf)
    @test @interval("0.0??d") == Interval(-Inf, 0.0)
    @test @interval("2.5??") == Interval(-Inf, Inf)
    @test @interval("2.5??u") == Interval(2.5, Inf)
    @test @interval("2.5??d") == Interval(-Inf, 2.5)
    @test @interval("2.500?5e+27") == Interval(0x1.01fa19a08fe7fp+91, 0x1.0302cc4352683p+91)
    @test @interval("2.500?5ue4") == Interval(0x1.86ap+14, 0x1.8768p+14)
    @test @interval("2.500?5de-5") == Interval(0x1.a2976f1cee4d5p-16, 0x1.a36e2eb1c432dp-16)
    @test @interval("10?3") == Interval(7.0, 13.0)
    @test @interval("10?3e380") == Interval(0x1.fffffffffffffp+1023, Inf)
    @test @interval("1.0000000000000001?1") == Interval(1.0, 0x1.0000000000001p+0)
    # 10?18 + 308 zeros
    @test @interval("10?1800000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000") == Interval(-Inf, Inf)
    @test @interval("[ Nai  ]") == ∅
    @test @interval("[ Nai  ]_ill") == ∅
    @test @interval("[ Nai  ]_trv") == ∅
    @test @interval("[ Empty  ]_ill") == ∅
    @test @interval("[  ]_com") == ∅
    @test @interval("[,]_com") == ∅
    @test @interval("[   Entire ]_com") == ∅
    @test @interval("[ -inf ,  INF ]_com") == ∅
    @test @interval("[  -1.0  , 1.0]_ill") == ∅
    @test @interval("[  -1.0  , 1.0]_fooo") == ∅
    @test @interval("[  -1.0  , 1.0]_da") == ∅
    @test @interval("[-1.0,]_com") == ∅
    @test @interval("[-Inf, 1.000 ]_ill") == ∅
    @test @interval("[-I  nf, 1.000 ]") == ∅
    @test @interval("[-Inf, 1.0  00 ]") == ∅
    @test @interval("[-Inf ]") == ∅
    @test @interval("[Inf , INF]") == ∅
    @test @interval("[ foo ]") == ∅
    # b-textToInterval "0.0??_com" = [empty] signal UndefinedOperation;
    # b-textToInterval "0.0??u_ill" = [empty] signal UndefinedOperation;
    # b-textToInterval "0.0??d_com" = [empty] signal UndefinedOperation;
    @test @interval("[1.0000000000000002,1.0000000000000001]") == Interval(1.0, 0x1.0000000000001p+0)
    @test @interval("[10000000000000001/10000000000000000,10000000000000002/10000000000000001]") == Interval(1.0, 0x1.0000000000001p+0)
    @test @interval("[0x1.00000000000002p0,0x1.00000000000001p0]") == Interval(1.0, 0x1.0000000000001p+0)
end

@testset "minimal_text_to_decorated_interval_test" begin
    @test @decorated("[ Empty  ]") == DecoratedInterval(∅, trv)
    @test @decorated("[ Empty  ]_trv") == DecoratedInterval(∅, trv)
    @test @decorated("[  ]") == DecoratedInterval(∅, trv)
    @test @decorated("[  ]_trv") == DecoratedInterval(∅, trv)
    @test @decorated("[,]") == DecoratedInterval(entireinterval(Float64), dac)
    @test @decorated("[,]_trv") == DecoratedInterval(entireinterval(Float64), trv)
    @test @decorated("[ entire  ]") == DecoratedInterval(entireinterval(Float64), dac)
    @test @decorated("[ ENTIRE ]_dac") == DecoratedInterval(entireinterval(Float64), dac)
    @test @decorated("[ -inf , INF  ]") == DecoratedInterval(entireinterval(Float64), dac)
    @test @decorated("[ -inf, INF ]_def") == DecoratedInterval(entireinterval(Float64), def)
    @test @decorated("[-1.0,1.0]") == DecoratedInterval(Interval(-1.0, 1.0), com)
    @test @decorated("[  -1.0  ,  1.0  ]_com") == DecoratedInterval(Interval(-1.0, 1.0), com)
    @test @decorated("[  -1.0  , 1.0]_trv") == DecoratedInterval(Interval(-1.0, 1.0), trv)
    @test @decorated("[-1,]") == DecoratedInterval(Interval(-1.0, Inf), dac)
    @test @decorated("[-1.0, +inf]_def") == DecoratedInterval(Interval(-1.0, Inf), def)
    @test @decorated("[-1.0, +infinity]_def") == DecoratedInterval(Interval(-1.0, Inf), def)
    @test @decorated("[-Inf, 1.000 ]") == DecoratedInterval(Interval(-Inf, 1.0), dac)
    @test @decorated("[-Infinity, 1.000 ]_trv") == DecoratedInterval(Interval(-Inf, 1.0), trv)
    @test @decorated("[1.0E+400 ]_com") == DecoratedInterval(Interval(0x1.fffffffffffffp+1023, Inf), dac)
    @test @decorated("[ -4/2, 10/5 ]_com") == DecoratedInterval(Interval(-2.0, 2.0), com)
    @test @decorated("[ -1/10, 1/10 ]_com") == DecoratedInterval(Interval(-0.1, 0.1), com)
    @test @decorated("0.0?") == DecoratedInterval(Interval(-0.05, 0.05), com)
    @test @decorated("0.0?u_trv") == DecoratedInterval(Interval(0.0, 0.05), trv)
    @test @decorated("0.0?d_dac") == DecoratedInterval(Interval(-0.05, 0.0), dac)
    @test @decorated("2.5?") == DecoratedInterval(Interval(0x1.3999999999999p+1, 0x1.4666666666667p+1), com)
    @test @decorated("2.5?u") == DecoratedInterval(Interval(2.5, 0x1.4666666666667p+1), com)
    @test @decorated("2.5?d_trv") == DecoratedInterval(Interval(0x1.3999999999999p+1, 2.5), trv)
    @test @decorated("0.000?5") == DecoratedInterval(Interval(-0.005, 0.005), com)
    @test @decorated("0.000?5u_def") == DecoratedInterval(Interval(0.0, 0.005), def)
    @test @decorated("0.000?5d") == DecoratedInterval(Interval(-0.005, 0.0), com)
    @test @decorated("2.500?5") == DecoratedInterval(Interval(0x1.3f5c28f5c28f5p+1, 0x1.40a3d70a3d70bp+1), com)
    @test @decorated("2.500?5u") == DecoratedInterval(Interval(2.5, 0x1.40a3d70a3d70bp+1), com)
    @test @decorated("2.500?5d") == DecoratedInterval(Interval(0x1.3f5c28f5c28f5p+1, 2.5), com)
    @test @decorated("0.0??_dac") == DecoratedInterval(Interval(-Inf, Inf), dac)
    @test @decorated("0.0??u_trv") == DecoratedInterval(Interval(0.0, Inf), trv)
    @test @decorated("0.0??d") == DecoratedInterval(Interval(-Inf, 0.0), dac)
    @test @decorated("2.5??") == DecoratedInterval(Interval(-Inf, Inf), dac)
    @test @decorated("2.5??u_def") == DecoratedInterval(Interval(2.5, Inf), def)
    @test @decorated("2.5??d_dac") == DecoratedInterval(Interval(-Inf, 2.5), dac)
    @test @decorated("2.500?5e+27") == DecoratedInterval(Interval(0x1.01fa19a08fe7fp+91, 0x1.0302cc4352683p+91), com)
    @test @decorated("2.500?5ue4_def") == DecoratedInterval(Interval(0x1.86ap+14, 0x1.8768p+14), def)
    @test @decorated("2.500?5de-5") == DecoratedInterval(Interval(0x1.a2976f1cee4d5p-16, 0x1.a36e2eb1c432dp-16), com)
    @test @decorated("[ Nai  ]") == nai()
    # 10?18 + 308 zeros + _com
    @test @decorated("10?1800000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000_com") == DecoratedInterval(Interval(-Inf, Inf), dac)
    @test @decorated("10?3_com") == DecoratedInterval(Interval(7.0, 13.0), com)
    @test @decorated("10?3e380_com") == DecoratedInterval(Interval(0x1.fffffffffffffp+1023, Inf), dac)
    @test @decorated("[ Nai  ]_ill") == nai()
    @test @decorated("[ Nai  ]_trv") == nai()
    @test @decorated("[ Empty  ]_ill") == nai()
    @test @decorated("[  ]_com") == nai()
    @test @decorated("[,]_com") == nai()
    @test @decorated("[   Entire ]_com") == nai()
    @test @decorated("[ -inf ,  INF ]_com") == nai()
    @test @decorated("[  -1.0  , 1.0]_ill") == nai()
    @test @decorated("[  -1.0  , 1.0]_fooo") == nai()
    @test @decorated("[  -1.0  , 1.0]_da") == nai()
    @test @decorated("[-1.0,]_com") == nai()
    @test @decorated("[-Inf, 1.000 ]_ill") == nai()
    @test @decorated("[-I  nf, 1.000 ]") == nai()
    @test @decorated("[-Inf, 1.0  00 ]") == nai()
    @test @decorated("[-Inf ]") == nai()
    @test @decorated("[Inf , INF]") == nai()
    @test @decorated("[ foo ]") == nai()
    @test @decorated("0.0??_com") == nai()
    @test @decorated("0.0??u_ill") == nai()
    @test @decorated("0.0??d_com") == nai()
    @test @decorated("0.0??_com") == nai()
    @test @decorated("[1.0,2.0") == nai()
    @test @decorated("[1.0000000000000002,1.0000000000000001]") == DecoratedInterval(Interval(1.0, 0x1.0000000000001p+0), com)
    @test @decorated("[10000000000000001/10000000000000000,10000000000000002/10000000000000001]") == DecoratedInterval(Interval(1.0, 0x1.0000000000001p+0), com)
    @test @decorated("[0x1.00000000000002p0,0x1.00000000000001p0]") == DecoratedInterval(Interval(1.0, 0x1.0000000000001p+0), com)
end

@testset "minimal_interval_part_test" begin

end

@testset "minimal_new_dec_test" begin

end

@testset "minimal_set_dec_test" begin

end

@testset "minimal_decoration_part_test" begin
    @test decoration(nai()) == ill
    @test decoration(DecoratedInterval(∅, trv)) == trv
    @test decoration(DecoratedInterval(Interval(-1.0, 3.0), trv)) == trv
    @test decoration(DecoratedInterval(Interval(-1.0, 3.0), def)) == def
    @test decoration(DecoratedInterval(Interval(-1.0, 3.0), dac)) == dac
    @test decoration(DecoratedInterval(Interval(-1.0, 3.0), com)) == com
end
# FactCheck.exitstatus()


text:
# coding: utf-8

# In[1]:

#coin change is a straight forward dynamic programming problem
#you are given one note and a set of coins of different values
#assuming you have infinite supply of coins of different values
#we need to compute different ways of making change
#the order of the coins doesnt matter in this case
#more details can be found in the following link
# https://www.geeksforgeeks.org/coin-change-dp-7/
#the script can also be done in dynamic programming
# https://github.com/je-suis-tm/recursion-and-dynamic-programming/blob/master/coin%20change%20dynamic%20programming.jl

# In[2]:

#to solve this problem via recursion
#we divide one big problem into two sub problems
#the case where coin of η value is excluded in the solutions
#and the case where at least one coin of η value is included
function coin_change(num,choice)

    #base case
    if num==0

        return 1

    end

    #prevent stack overflow
    if num<0

        return 0

    end

    output=0

    #iterate through cases of exclusion
    for i in 1:length(choice)

        #case of inclusion
        include=coin_change(num-choice[i],choice)
        exclude=0

        #case of exclusion
        if i>=2

            exclude=coin_change(num,choice[1:(i-1)])

        end

        #two sub problems merge into a big one
        output=include+exclude

    end

    return output

end

# In[3]:

coin_change(10,[1,2,3])

max_stars_repo_path: coin change recursion.jl

max_issues_repo_path: coin change recursion.jl

max_forks_repo_path: coin change recursion.jl


text:
"""
ZeroKernel([tr=IdentityTransform()])

Create a kernel always returning zero
"""
struct ZeroKernel{Tr} <: Kernel{Tr}
    transform::Tr
end

ZeroKernel() = ZeroKernel(IdentityTransform())

@inline kappa(κ::ZeroKernel, d::T) where {T<:Real} = zero(T)

metric(::ZeroKernel) = Delta()

"""
`WhiteKernel([tr=IdentityTransform()])`

```
    κ(x,y) = δ(x,y)
```
Kernel function working as an equivalent to add white noise.
"""
struct WhiteKernel{Tr} <: Kernel{Tr}
    transform::Tr
end

WhiteKernel() = WhiteKernel(IdentityTransform())

@inline kappa(κ::WhiteKernel,δₓₓ::Real) = δₓₓ

metric(::WhiteKernel) = Delta()

"""
`ConstantKernel([tr=IdentityTransform(),[c=1.0]])`
```
    κ(x,y) = c
```
Kernel function always returning a constant value `c`
"""
struct ConstantKernel{Tr, Tc<:Real} <: Kernel{Tr}
    transform::Tr
    c::Tc
end

params(k::ConstantKernel) = (params(k.transform),k.c)
opt_params(k::ConstantKernel) = (opt_params(k.transform),k.c)

ConstantKernel(c::Real=1.0) = ConstantKernel(IdentityTransform(),c)

ConstantKernel(t::Tr,c::Tc=1.0) where {Tr<:Transform,Tc<:Real} = ConstantKernel{Tr,Tc}(t,c)

@inline kappa(κ::ConstantKernel,x::Real) = κ.c

metric(::ConstantKernel) = Delta()


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdslider

"""
    'feffery'_antdslider(;kwargs...)

An AntdSlider component.

Keyword arguments:
- `id` (String; optional)
- `className` (String; optional)
- `defaultValue` (Real | Array of Reals; optional)
- `disabled` (Bool; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `marks` (Dict; optional)
- `max` (Real; optional)
- `min` (Real; optional)
- `range` (Bool; optional)
- `step` (Real; optional)
- `style` (Dict; optional)
- `tooltipPrefix` (String; optional)
- `tooltipSuffix` (String; optional)
- `tooltipVisible` (Bool; optional)
- `value` (Real | Array of Reals; optional)
"""
function 'feffery'_antdslider(; kwargs...)
        available_props = Symbol[:id, :className, :defaultValue, :disabled, :loading_state, :marks, :max, :min, :range, :step, :style, :tooltipPrefix, :tooltipSuffix, :tooltipVisible, :value]
        wild_props = Symbol[]
        return Component("'feffery'_antdslider", "AntdSlider", "feffery_antd_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: Beddow/feffery-antd-components

max_issues_repo_name: limmyforget/feffery-antd-components

max_forks_repo_name: limmyforget/feffery-antd-components


text:
# The series, 1^1 + 2^2 + 3^3 + ... + 10^10 = 10405071317.
#
# Find the last ten digits of the series, 1^1 + 2^2 + 3^3 + ... + 1000^1000.

using ProjectEulerSolutions

# Multiply using mods since only the last 10 digits matter.  Faster than the
# bigint solution after about n > 30,000, and uses much less memory.
function p048solution_mod(n::Integer=100)::Integer
    total = 1
    for i in 2:n
        s = i
        for j in 2:i
            s = mod(s * i, 10^10)
        end
        total = mod(total + s, 10^10)
    end
    return total
end

# Trivial with Julia's BigInt support, and fast.
function p048solution_bigint(n::Integer=100)::Integer
    return mod(mapreduce(x->x^x, +, BigInt.(1:n)), 10^10)
end

p048 = Problems.Problem(Dict("Modulus" => p048solution_mod,
                             "Bigint" => p048solution_bigint))

Problems.benchmark(p048, 1000)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
import HierarchicalUtils: nodeshow, nodecommshow

function datasummary(n::AbstractMillNode)
    bytes = Base.format_bytes(Base.summarysize(n))
    string("# Summary: ", nobs(n), " obs, ", bytes, ".")
end

function Base.show(io::IO, ::MIME"text/plain", @nospecialize(n::AbstractMillNode))
    HierarchicalUtils.printtree(io, n; htrunc=3, vtrunc=5, breakline=false)
end

nodeshow(io::IO, ::Missing) = print(io, "∅")
nodeshow(io::IO, n::LazyNode{N, Nothing}) where N = print(io, "LazyNode{$N} ∅")

function nodecommshow(io::IO, @nospecialize(n::AbstractMillNode))
    bytes = Base.format_bytes(Base.summarysize(n) - (isleaf(n) ? 0 : Base.summarysize(data(n))))
    print(io, "\t# ", nobs(n), " obs, ", bytes)
end

function Base.show(io::IO, @nospecialize(n::AbstractMillNode))
    print(io, nameof(typeof(n)))
    if !get(io, :compact, false)
        _show_data(io, n)
    end
end

_show_data(io, n::LazyNode{Name}) where {Name} = print(io, "{", Name, "}")
_show_data(io, _) = print(io)

# params summary from https://github.com/FluxML/Flux.jl/blob/master/src/layers/show.jl
function modelsummary(m::AbstractMillModel)
    ps = params(m)
    npars = Flux.underscorise(sum(length, ps))
    bytes = Base.format_bytes(sum(Base.summarysize, collect(ps.params)))
    string("# Summary: ", length(ps), " arrays, ", npars, " params, ", bytes)
end

function Base.show(io::IO, ::MIME"text/plain", @nospecialize(m::AbstractMillModel))
    HierarchicalUtils.printtree(io, m; htrunc=3, vtrunc=5, breakline=false)
end

_levelparams(m::ArrayModel) = params(m.m)
_levelparams(m::BagModel) = params(m.a, m.bm)
_levelparams(m::ProductModel) = params(m.m)
_levelparams(m::LazyModel) = Params([])

# params summary from https://github.com/FluxML/Flux.jl/blob/master/src/layers/show.jl
function nodecommshow(io::IO, @nospecialize(m::AbstractMillModel))
    ps = _levelparams(m)
    if !isempty(ps)
        npars = Flux.underscorise(sum(length, ps))
        print(io, "\t# ", length(ps), " arrays, ", npars, " params")
        if !isempty(ps) && Flux._all(iszero, ps)
            print(io, " (all zero)")
        elseif Flux._any(isnan, ps)
            print(io, " (some NaN)")
        elseif Flux._any(isinf, ps)
            print(io, " (some Inf)")
        end
        bytes = Base.format_bytes(sum(Base.summarysize, collect(ps.params)))
        print(io, ", ", bytes)
    end
end

function Base.show(io::IO, @nospecialize(n::AbstractMillModel))
    print(io, nameof(typeof(n)))
    if !get(io, :compact, false)
        _show_submodels(io, n)
    end
end

_show_submodels(io, m::ArrayModel) = print(io, "(", m.m, ")")
_show_submodels(io, m::BagModel) = print(io, " ↦ ", m.a, " ↦ ", m.bm)
_show_submodels(io, m::ProductModel) = print(io, " ↦ ", m.m)
_show_submodels(io, m::LazyModel{Name}) where {Name} = print(io, "{", Name, "}")
_show_submodels(io, _) = print(io)

function Base.showarg(io::IO, x::MaybeHotVector, toplevel)
    print(io, "MaybeHotVector")
    toplevel && print(io, " with eltype $(eltype(x))")
    return nothing
end

function Base.showarg(io::IO, X::MaybeHotMatrix, toplevel)
    print(io, "MaybeHotMatrix")
    toplevel && print(io, " with eltype $(eltype(X))")
    return nothing
end

# this is from /LinearAlgebra/src/diagonal.jl, official way to print the dots:
function Base.replace_in_print_matrix(x::MaybeHotArray, i::Integer, j::Integer, s::AbstractString)
    ismissing(x[i, j]) || x[i, j] ? s : Base.replace_with_centered_mark(s)
end

function Base.show(io::IO, X::T) where T <: Union{ImputingMatrix, MaybeHotArray, NGramMatrix}
    if get(io, :compact, false)
        if ndims(X) == 1
            print(io, length(X), "-element ", nameof(T))
        else
            print(io, join(size(X), "×"), " ", nameof(T))
        end
    else
        _show_fields(io, X)
    end
end

Base.show(io::IO, ::MIME"text/plain", a::AbstractAggregation) = _show_fields(io, a)

function _show_fields(io, x::T; context=:compact=>true) where T
    o = join(["$f = $(repr(getfield(x, f); context))" for f in fieldnames(T)],", ")
    print(io, nameof(T), "(", o, ")")
end

function Base.print_array(io::IO, A::ImputingMatrix)
    println(io, "W:")
    Base.print_array(io, A.W)
    println(io, "\n\nψ:")
    _print_params(io, A)
end

_print_params(io::IO, A::PreImputingMatrix) = Base.print_array(io, A.ψ |> permutedims)
_print_params(io::IO, A::PostImputingMatrix) = Base.print_array(io, A.ψ)

Base.print_array(io::IO, A::NGramMatrix) = Base.print_array(io, A.S)

function Base.show(io::IO, l::Dense{F, <:ImputingMatrix}) where F
    print(io, _name(l.weight), "Dense(", size(l.weight, 2), ", ", size(l.weight, 1))
    l.σ == identity || print(io, ", ", l.σ)
    l.bias == Flux.Zeros() && print(io, "; bias=false")
    print(io, ")")
end

_name(::PreImputingMatrix) = "[preimputing]"
_name(::PostImputingMatrix) = "[postimputing]"

function Base.show(io::IO, a::T) where T <: AbstractAggregation
    print(io, nameof(T))
    if !get(io, :compact, false)
        print(io, "(", length(a), ")")
    end
end


text:
using Revise, Test

isdefined(@__MODULE__, :do_test) || include("common.jl")

module BackEdgesTest
using Test
flag = false    # this needs to be defined for the conditional part to work
end

do_test("Backedges") && @testset "Backedges" begin
    src = """
    # issue #249
    flag = false
    if flag
        f() = 1
    else
        f() = 2
    end

    # don't do work in the interpreter that isn't needed for function definitions
    # inspired by #300
    const planetdiameters = Dict("Mercury" => 4_878)
    planetdiameters["Venus"] = 12_104

    function getdiameter(name)
        return planetdiameters[name]
    end
    """
    mexs = Revise.parse_source!(Revise.ModuleExprsSigs(BackEdgesTest), src, "backedges_test.jl", BackEdgesTest)
    Revise.moduledeps[BackEdgesTest] = Revise.DepDict()
    Revise.instantiate_sigs!(mexs)
    @test isempty(methods(BackEdgesTest.getdiameter))
    @test !isdefined(BackEdgesTest, :planetdiameters)
    @test length(Revise.moduledeps[BackEdgesTest]) == 1
    @test Revise.moduledeps[BackEdgesTest][:flag] == Set([(BackEdgesTest, first(Iterators.drop(mexs[BackEdgesTest], 1))[1])])

    # issue #399
    src = """
    for jy in ("j","y"), nu in (0,1)
        jynu = Expr(:quote, Symbol(jy,nu))
        jynuf = Expr(:quote, Symbol(jy,nu,"f"))
        bjynu = Symbol("bessel",jy,nu)
        if jy == "y"
            @eval begin
                \$bjynu(x::Float64) = nan_dom_err(ccall((\$jynu,libm),  Float64, (Float64,), x), x)
                \$bjynu(x::Float32) = nan_dom_err(ccall((\$jynuf,libm), Float32, (Float32,), x), x)
                \$bjynu(x::Float16) = Float16(\$bjynu(Float32(x)))
            end
        else
            @eval begin
                \$bjynu(x::Float64) = ccall((\$jynu,libm),  Float64, (Float64,), x)
                \$bjynu(x::Float32) = ccall((\$jynuf,libm), Float32, (Float32,), x)
                \$bjynu(x::Float16) = Float16(\$bjynu(Float32(x)))
            end
        end
        @eval begin
            \$bjynu(x::Real) = \$bjynu(float(x))
            \$bjynu(x::Complex) = \$(Symbol("bessel",jy))(\$nu,x)
        end
    end
    """
    ex = Meta.parse(src)
    @test Revise.methods_by_execution(BackEdgesTest, ex) isa Tuple

    # Issue #428
    src = """
    @testset for i in (1, 2)
        @test i == i
    end
    """
    ex = Meta.parse(src)
    @test Revise.methods_by_execution(BackEdgesTest, ex) isa Tuple

end


text:
using YaoBlocks, YaoArrayRegister, BenchmarkTools, StaticArrays

r = rand_state(20)
U = @SMatrix rand(ComplexF64, 2, 2)

t1 = @benchmark instruct!($(statevec(r)), $U, 1)
t2 = @benchmark instruct!($(statevec(r)), $(Matrix(U)), 1)
(minimum(t1).time - minimum(t2).time) / minimum(t1).time

v = statevec(r)
mU = Matrix(U)

@profiler for i in 1:100
    instruct!(v, U, (1, ))
end

@profiler for i in 1:100
    instruct!(v, mU, (1, ))
end

function instruct2!(st, U, loc)
    a, c, b, d = U
    step = 1 << (loc - 1)
    step_2 = 1 << loc
    kernel(st, step, step_2, a, b, c, d)
end

@inline function kernel(state, step, step_2, a, b, c, d)
    for j in 0:step_2:size(state, 1)-step
        kernel2(state, U, step, step_2, a, b, c, d, j)
    end
    return state
end

@inline function kernel2(state, step, step_2, a, b, c, d, j)
    @inbounds for i in j+1:j+step
        YaoArrayRegister.u1rows!(state, i, i+step, a, b, c, d)
    end
    return state
end

@benchmark instruct2!($(r.state), $U, 1)

@benchmark instruct2!($(r.state), $(Matrix(U)), 1)

a, c, b, d = U
step = 1 << (1 - 1)
step_2 = 1 << 1

@benchmark kernel($(statevec(r)), $step, $step_2, $a, $b, $c, $d)

@benchmark a, c, b, d = $(Matrix(U))


text:
"""
Computes the maximum flow between the source and target vertexes in a flow
graph using [Dinic\'s Algorithm](https://en.wikipedia.org/wiki/Dinic%27s_algorithm)
Returns the value of the maximum flow as well as the final flow matrix.

Use a default capacity of 1 when the capacity matrix isn\'t specified.

Requires arguments:
residual_graph::DiGraph                # the input graph
source::Int                            # the source vertex
target::Int                            # the target vertex
capacity_matrix::AbstractArray{T,2}    # edge flow capacities
"""

function dinic_impl{T<:Number}(
    residual_graph::DiGraph,               # the input graph
    source::Int,                           # the source vertex
    target::Int,                           # the target vertex
    capacity_matrix::AbstractArray{T,2}    # edge flow capacities
    )
    n = nv(residual_graph)                     # number of vertexes

    flow_matrix = zeros(T, n, n)           # initialize flow matrix
    P = zeros(Int, n)                      # Sharable parent vector

    flow = 0

    while true
        augment = blocking_flow!(residual_graph, source, target, capacity_matrix, flow_matrix, P)
        augment == 0 && break
        flow += augment
    end
    return flow, flow_matrix
end

"""
Uses BFS to identify a blocking flow in
the input graph and then backtracks from the targetto the source, aumenting flow
along all possible paths.

Requires arguments:
residual_graph::DiGraph                # the input graph
source::Int                            # the source vertex
target::Int                            # the target vertex
capacity_matrix::AbstractArray{T,2}    # edge flow capacities
flow_matrix::AbstractArray{T,2}        # the current flow matrix
"""
function blocking_flow!{T<:Number}(
    residual_graph::DiGraph,               # the input graph
    source::Int,                           # the source vertex
    target::Int,                           # the target vertex
    capacity_matrix::AbstractArray{T,2},   # edge flow capacities
    flow_matrix::AbstractArray{T,2},       # the current flow matrix
    )
    P = zeros(T, nv(residual_graph))
    return blocking_flow!(residual_graph,
                          source,
                          target,
                          capacity_matrix,
                          flow_matrix,
                          P)
end

"""blocking_flow!
Preallocated version of blocking_flow.Uses BFS to identify a blocking flow in
the input graph and then backtracks from the target to the source, aumenting flow
along all possible paths.

Requires arguments:
residual_graph::DiGraph                # the input graph
source::Int                            # the source vertex
target::Int                            # the target vertex
capacity_matrix::AbstractArray{T,2}    # edge flow capacities
flow_matrix::AbstractArray{T,2}        # the current flow matrix
P::AbstractArray{Int, 1}               # Parent vector to store Level Graph
"""

function blocking_flow!{T<:Number}(
    residual_graph::DiGraph,               # the input graph
    source::Int,                           # the source vertex
    target::Int,                           # the target vertex
    capacity_matrix::AbstractArray{T,2},   # edge flow capacities
    flow_matrix::AbstractArray{T,2},       # the current flow matrix
    P::AbstractArray{Int, 1}               # Parent vector to store Level Graph
    )
    n = nv(residual_graph)                     # number of vertexes

    fill!(P, -1)
    P[source] = -2

    Q = [source]
    sizehint!(Q, n)

    while length(Q) > 0                   # Construct the Level Graph using BFS
        u = pop!(Q)
        for v in fadj(residual_graph, u)
            if P[v] == -1 && capacity_matrix[u,v] > flow_matrix[u,v]
                P[v] = u
                unshift!(Q, v)
            end
        end
    end

    P[target] == -1 && return 0                    # BFS couldn't reach the target

    total_flow = 0

    for bv in badj(residual_graph, target)    # Trace all possible routes to source
        flow = typemax(T)
        v = target
        u = bv
        while v != source
            if u == -1                    # Vertex unreachable from source
                flow = 0
                break
            else
                flow = min(flow, capacity_matrix[u,v] - flow_matrix[u,v])
                v = u
                u = P[u]
            end
        end

        flow == 0 && continue                      # Flow cannot be augmented along path

        v = target
        u = bv
        while v != source             # Augment flow along path
            flow_matrix[u,v] += flow
            flow_matrix[v,u] -= flow
            v = u
            u = P[u]
        end

        total_flow += flow
    end
    return total_flow
end


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license

export add_atom!, remove_atom!
export has_velocities, add_velocities!, positions, velocities
export set_cell!, set_topology!, set_step!, guess_bonds!
export distance, dihedral, out_of_plane

__ptr(frame::Frame) = __ptr(frame.__handle)
__const_ptr(frame::Frame) = __const_ptr(frame.__handle)

# A small wrapper around Array{Float64, 2}, keeping a reference to the
# corresponding frame to prevent garbage collection (see issue
# https://github.com/chemfiles/Chemfiles.jl/issues/37)
struct ChemfilesArray <: AbstractArray{Float64,2}
    data::Array{Float64,2}
    parent::Frame
end

# Implement the Array interface for ChemfilesArray
@inline Base.size(A::ChemfilesArray) = size(A.data)
@inline Base.getindex(A::ChemfilesArray, I::Int) = getindex(A.data, I)
@inline Base.getindex(A::ChemfilesArray, I::Int...) = getindex(A.data, I...)
@inline Base.setindex!(A::ChemfilesArray, v, I::Int) = setindex!(A.data, v, I)
@inline Base.setindex!(A::ChemfilesArray, v, I::Int...) = setindex!(A.data, v, I...)

"""
Create a new empty `Frame`.
"""
function Frame()
    ptr = @__check_ptr(lib.chfl_frame())
    return Frame(CxxPointer(ptr, is_const=false))
end

"""
Get the number of atoms in the `frame`.
"""
function Base.size(frame::Frame)
    count = Ref{UInt64}(0)
    __check(lib.chfl_frame_atoms_count(__const_ptr(frame), count))
    return Int(count[])
end

"""
Get the number of atoms in the `frame`.
"""
function Base.length(frame::Frame)
    size(frame)
end

"""
Resize the positions and the velocities in the `frame`, to make space for
`natoms` atoms. This function may invalidate any pointer to the positions or
the velocities if the new size is bigger than the old one. In all the cases,
previous data is conserved. This function conserve the presence or absence of
velocities.
"""
function Base.resize!(frame::Frame, natoms::Integer)
    __check(lib.chfl_frame_resize(__ptr(frame), UInt64(natoms)))
end

"""
Get the positions in a `Frame` as an array. The positions are readable and
writable from this array. If the frame is resized (by writing to it, or calling
`resize!`), the array is invalidated.
"""
function positions(frame::Frame)
    ptr = Ref{Ptr{Float64}}()
    natoms = Ref{UInt64}(0)
    __check(lib.chfl_frame_positions(__ptr(frame), ptr, natoms))
    data = unsafe_wrap(Array{Float64,2}, ptr[], (3, Int(natoms[])); own=false)
    return ChemfilesArray(data, frame)
end

"""
Get the velocities in a `Frame` as an array. The velocities are readable and
writable from this array. If the frame is resized (by writing to it, or calling
`resize!`), the array is invalidated.

If the frame do not have velocity, this function will error. You can use
`add_velocities!` to add velocities to a frame before calling this function.
"""
function velocities(frame::Frame)
    ptr = Ref{Ptr{Float64}}()
    natoms = Ref{UInt64}(0)
    __check(lib.chfl_frame_velocities(__ptr(frame), ptr, natoms))
    data = unsafe_wrap(Array{Float64,2}, ptr[], (3, Int(natoms[])); own=false)
    return ChemfilesArray(data, frame)
end

"""
Add velocities to this `frame`. The storage is initialized with the result of
`size(frame)` as the number of atoms. If the frame already has velocities, this
does nothing.
"""
function add_velocities!(frame::Frame)
    __check(lib.chfl_frame_add_velocities(__ptr(frame)))
    return nothing
end

"""
Check if a `frame` contains velocity data or not.
"""
function has_velocities(frame::Frame)
    result = Ref{UInt8}(0)
    __check(lib.chfl_frame_has_velocities(__const_ptr(frame), result))
    return convert(Bool, result[])
end

"""
Set the `cell` associated with a `frame`.
"""
function set_cell!(frame::Frame, cell::UnitCell)
    __check(lib.chfl_frame_set_cell(__ptr(frame), __const_ptr(cell)))
    return nothing
end

"""
Set the `topology` associated with a `frame`.
"""
function set_topology!(frame::Frame, topology::Topology)
    __check(lib.chfl_frame_set_topology(__ptr(frame), __const_ptr(topology)))
    return nothing
end

"""
Get the `frame` step, *i.e.* the frame number in the trajectory.
"""
function Base.step(frame::Frame)
    result = Ref{UInt64}(0)
    __check(lib.chfl_frame_step(__const_ptr(frame), result))
    return result[]
end

"""
Set the `frame` step to `step`.
"""
function set_step!(frame::Frame, step::Integer)
    __check(lib.chfl_frame_set_step(__ptr(frame), UInt64(step)))
    return nothing
end

"""
Guess the bonds, angles, and dihedrals in the `frame` using a distance criteria.
"""
function guess_bonds!(frame::Frame)
    __check(lib.chfl_frame_guess_bonds(__ptr(frame)))
    return nothing
end

"""
Calculate the distance between two atoms.
"""
function distance(frame::Frame, i::Integer, j::Integer)
    result = Ref{Float64}(0)
    __check(lib.chfl_frame_distance(__const_ptr(frame), UInt64(i), UInt64(j), result))
    return result[]
end

"""
Calculate the angle made by three atoms.
"""
function Base.angle(frame::Frame, i::Integer, j::Integer, k::Integer)
    result = Ref{Float64}(0)
    __check(lib.chfl_frame_angle(__const_ptr(frame), UInt64(i), UInt64(j), UInt64(k), result))
    return result[]
end

"""
Calculate the dihedral (torsional) angle made by four unbranched atoms.
"""
function dihedral(frame::Frame, i::Integer, j::Integer, k::Integer, m::Integer)
    result = Ref{Float64}(0)
    __check(lib.chfl_frame_dihedral(__const_ptr(frame), UInt64(i), UInt64(j), UInt64(k), UInt64(m), result))
    return result[]
end

"""
Calculate the out-of-plane (improper) angle made by four atoms.
"""
function out_of_plane(frame::Frame, i::Integer, j::Integer, k::Integer, m::Integer)
    result = Ref{Float64}(0)
    __check(lib.chfl_frame_out_of_plane(__const_ptr(frame), UInt64(i), UInt64(j), UInt64(k), UInt64(m), result))
    return result[]
end

"""
Add an `atom` and the corresponding `position` and `velocity` data to a `frame`.
"""
function add_atom!(frame::Frame, atom::Atom, position::Vector{Float64}, velocity::Vector{Float64}=Float64[0.0,0.0,0.0])
    __check(lib.chfl_frame_add_atom(__ptr(frame), __const_ptr(atom), position, velocity))
    return nothing
end

"""
Remove the `atom` at `index` from the `frame`.

This function modifies all the `atoms` indexes after `index`, and invalidates
any array obtained using `positions` or `velocities`.
"""
function remove_atom!(frame::Frame, index::Integer)
    __check(lib.chfl_frame_remove(__ptr(frame), UInt64(index)))
    return nothing
end

"""
Set a named property for the given `Frame`.
"""
function set_property!(frame::Frame, name::String, value)
    property = Property(value)
    __check(lib.chfl_frame_set_property(
        __ptr(frame), pointer(name), __const_ptr(property)
    ))
    return nothing
end

"""
Get a named property for the given atom.
"""
function property(frame::Frame, name::String)::PropertyValue
    ptr = lib.chfl_frame_get_property(__const_ptr(frame), pointer(name))
    return extract(Property(CxxPointer(ptr, is_const=false)))
end

"""
Get the number of properties associated with a frame.
"""
function properties_count(frame::Frame)
    count = Ref{UInt64}(0)
    __check(lib.chfl_frame_properties_count(__const_ptr(frame), count))
    return Int(count[])
end

"""
Get the names of all properties associated with a frame.
"""
function list_properties(frame::Frame)
    count = UInt64(properties_count(frame))
    names = Array{Ptr{UInt8}}(undef, count)
    __check(lib.chfl_frame_list_properties(__const_ptr(frame), pointer(names), count))
    return map(unsafe_string, names)
end

"""
Add an additional bond to the `Frame`'s `Topology`.
"""
function add_bond!(frame::Frame, i::Integer, j::Integer, order=nothing)
    if order === nothing
        __check(lib.chfl_frame_add_bond(__ptr(frame), UInt64(i), UInt64(j)))
    else
        # Check that the order is a valid BondOrder
        order = BondOrder(Integer(order))
        __check(lib.chfl_frame_bond_with_order(
            __ptr(frame), UInt64(i), UInt64(j), lib.chfl_bond_order(order)
        ))
    end
    return nothing
end

"""
Remove a bond from the `Frame`'s `Topology`.
"""
function remove_bond!(frame::Frame, i::Integer, j::Integer)
    __check(lib.chfl_frame_remove_bond(__ptr(frame), UInt64(i), UInt64(j)))
    return nothing
end

"""
Remove all bonds, angles and dihedral angles from the `Frame`'s `Topology`.
"""
function clear_bonds!(frame::Frame)
    __check(lib.chfl_frame_clear_bonds(__ptr(frame)))
    return nothing
end

"""
Add a residue to the `Frame`'s `Topology`.
"""
function add_residue!(frame::Frame, residue::Residue)
    __check(lib.chfl_frame_add_residue(__ptr(frame), __const_ptr(residue)))
    return nothing
end

"""
Make a deep copy of a `Frame`.
"""
function Base.deepcopy(frame::Frame)
    ptr = lib.chfl_frame_copy(__const_ptr(frame))
    return Frame(CxxPointer(ptr, is_const=false))
end

# Indexing support
"""
Get the `Atom` at the given `index` of the `frame`. By default this creates a
copy so as to be safe. To not create a copy, use `@view frame[index]`.

See also [`Base.view(frame::Frame, index::Integer)`](@ref)
"""
Base.getindex(frame::Frame, index::Integer) = Atom(frame, index)

"""
Get the `Atom` at the given `index` of the `frame` without creating a copy.

!!! warning

    This function can lead to unefined behavior when keeping the returned `Atom`
    around. Whith code like this:

    ```
    frame = Frame()
    resize!(frame, 3)

    atom = @view frame[0]
    resize!(frame, 4)
    ```

    `atom` contains a pointer to memory owned by `frame`, but this
    pointer has been invalidated when resizing the frame. Using `atom` after
    the second call to `resize!` might trigger undefined behavior (segmentation
    fault in the best case, silent data corruption in the worst case).
"""
function Base.view(frame::Frame, index::Integer)
    ptr = @__check_ptr(lib.chfl_atom_from_frame(__ptr(frame), UInt64(index)))
    atom = Atom(CxxPointer(ptr, is_const=false))
    return atom
end

# Iteration support
function Base.iterate(frame::Frame, atom=0)
    if atom >= size(frame)
        return nothing
    else
        return (frame[atom], atom + 1)
    end
end
Base.eltype(::Frame) = Atom


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    EstimationProblem(sdata, sdomain, targetvars)

A spatial estimation problem on a given spatial domain `sdomain`
in which the variables to be estimated are listed in `targetvars`.
The data of the problem is stored in spatial data `sdata`.

## Examples

Create an estimation problem for rainfall precipitation measurements:

```julia
julia> EstimationProblem(sdata, sdomain, :precipitation)
```

Create an estimation problem for precipitation and CO₂:

```julia
julia> EstimationProblem(sdata, sdomain, (:precipitation,:CO₂))
```
"""
struct EstimationProblem{S,D,M} <: AbstractProblem
  # input fields
  sdata::S
  sdomain::D
  mapper::M
  targetvars::Dict{Symbol,DataType}

  # state fields
  mappings::Dict{Symbol,Dict{Int,Int}}

  function EstimationProblem{S,D,M}(sdata, sdomain, targetvars, mapper) where {S,D,M}
    probvnames = Tuple(keys(targetvars))
    datavnames = name.(variables(sdata))

    @assert !isempty(probvnames) && probvnames ⊆ datavnames "target variables must be present in spatial data"
    @assert coordtype(sdata) == coordtype(sdomain) "data and domain must have the same coordinate type"

    mappings = map(sdata, sdomain, probvnames, mapper)

    new(sdata, sdomain, mapper, targetvars, mappings)
  end
end

function EstimationProblem(sdata::S, sdomain::D, targetvarnames::NTuple;
                           mapper::M=NearestMapper()) where {S,D,M}
  # build dictionary of target variables
  vars = filter(v -> name(v) ∈ targetvarnames, variables(sdata))
  targetvars = Dict(name(var) => mactype(var) for var in vars)

  EstimationProblem{S,D,M}(sdata, sdomain, targetvars, mapper)
end

EstimationProblem(sdata, sdomain, targetvarname::Symbol; mapper=NearestMapper()) =
  EstimationProblem(sdata, sdomain, (targetvarname,); mapper=mapper)

"""
    data(problem)

Return the spatial data of the estimation `problem`.
"""
data(problem::EstimationProblem) = problem.sdata

"""
    domain(problem)

Return the spatial domain of the estimation `problem`.
"""
domain(problem::EstimationProblem) = problem.sdomain

"""
    mapper(problem)

Return the mapper of the estimation `problem`.
"""
mapper(problem::EstimationProblem) = problem.mapper

"""
    variables(problem)

Return the variable names of the estimation `problem` and their types.
"""
variables(problem::EstimationProblem) = problem.targetvars

"""
    datamap(problem, targetvar)

Return the mapping from domain locations to data locations for the
`targetvar` of the `problem`.
"""
datamap(problem::EstimationProblem, var::Symbol) = problem.mappings[var]

"""
    datamap(problem)

Return the mappings from domain locations to data locations for all
the variables of the `problem`.
"""
datamap(problem::EstimationProblem) = problem.mappings

# ------------
# IO methods
# ------------
function Base.show(io::IO, problem::EstimationProblem)
  N = ncoords(problem.sdomain)
  print(io, "$(N)D EstimationProblem")
end

function Base.show(io::IO, ::MIME"text/plain", problem::EstimationProblem)
  vars = ["$var ($T)" for (var,T) in problem.targetvars]
  println(io, problem)
  println(io, "  data:      ", problem.sdata)
  println(io, "  domain:    ", problem.sdomain)
  print(  io, "  variables: ", join(vars, ", ", " and "))
end


text:
macro gputime(codeblock)
    quote
        local const query        = GLuint[1]
        local const elapsed_time = GLuint64[1]
        local const done         = GLint[0]
        glGenQueries(1, query)
        glBeginQuery(GL_TIME_ELAPSED, query[1])
        value = $(esc(codeblock))
        glEndQuery(GL_TIME_ELAPSED)

        while (done[1] != 1)
            glGetQueryObjectiv(
                query[1],
                GL_QUERY_RESULT_AVAILABLE,
                done
            )
        end
        glGetQueryObjectui64v(query[1], GL_QUERY_RESULT, elapsed_time)
        println("Time Elapsed: ", elapsed_time[1] / 1000000.0, "ms")
    end
end

immutable IterOrScalar{T}
  val::T
end

minlenght(a::Tuple{Vararg{IterOrScalar}}) = foldl(typemax(Int), a) do len, elem
    isa(elem.val, AbstractArray) && len > length(elem.val) && return length(elem.val)
    len
end
getindex{T<:AbstractArray}(A::IterOrScalar{T}, i::Integer) = A.val[i]
getindex(A::IterOrScalar, i::Integer) = A.val

#Some mapping functions for dictionaries
function mapvalues(func, collection::Dict)
    Dict([(key, func(value)) for (key, value) in collection])
end
function mapkeys(func, collection::Dict)
    Dict([(func(key), value) for (key, value) in collection])
end
Base.get{KT, VT}(a::Dict{KT, VT}, keys::Vector{KT}) = [a[key] for key in keys]

function print_with_lines(out::IO, text::AbstractString)
    io = IOBuffer()
    for (i,line) in enumerate(split(text, "\n"))
        println(io, @sprintf("%-4d: %s", i, line))
    end
    println(out, takebuf_string(io))
end
print_with_lines(text::AbstractString) = print_with_lines(STDOUT, text)

"""
Style Type, which is used to choose different visualization/editing styles via multiple dispatch
Usage pattern:
visualize(::Style{:Default}, ...)           = do something
visualize(::Style{:MyAwesomeNewStyle}, ...) = do something different
"""
immutable Style{StyleValue}
end
Style(x::Symbol) = Style{x}()
Style() = Style{:Default}()
mergedefault!{S}(style::Style{S}, styles, customdata) = merge!(copy(styles[S]), Dict{Symbol, Any}(customdata))
macro style_str(string)
    Style{Symbol(string)}
end
export @style_str

"""
splats keys from a dict into variables
"""
macro materialize(dict_splat)
    keynames, dict = dict_splat.args
    keynames = isa(keynames, Symbol) ? [keynames] : keynames.args
    dict_instance = gensym()
    kd = [:($key = $dict_instance[$(Expr(:quote, key))]) for key in keynames]
    kdblock = Expr(:block, kd...)
    expr = quote
        $dict_instance = $dict # handle if dict is not a variable but an expression
        $kdblock
    end
    esc(expr)
end
"""
splats keys from a dict into variables and removes them
"""
macro materialize!(dict_splat)
    keynames, dict = dict_splat.args
    keynames = isa(keynames, Symbol) ? [keynames] : keynames.args
    dict_instance = gensym()
    kd = [:($key = pop!($dict_instance, $(Expr(:quote, key)))) for key in keynames]
    kdblock = Expr(:block, kd...)
    expr = quote
        $dict_instance = $dict # handle if dict is not a variable but an expression
        $kdblock
    end
    esc(expr)
end

"""
Needed to match the lazy gl_convert exceptions.
    `Target`: targeted OpenGL type
    `x`: the variable that gets matched
"""
matches_target{Target, T}(::Type{Target}, x::T) = applicable(gl_convert, Target, x) || T <: Target  # it can be either converted to Target, or it's already the target
matches_target{Target, T}(::Type{Target}, x::Signal{T}) = applicable(gl_convert, Target, x)  || T <: Target
matches_target(::Function, x) = true
matches_target(::Function, x::Void) = false
export matches_target

"""
Takes a dict and inserts defaults, if not already available.
The variables are made accessible in local scope, so things like this are possible:
gen_defaults! dict begin
    a = 55
    b = a * 2 # variables, like a, will get made visible in local scope
    c::JuliaType = X # `c` needs to be of type JuliaType. `c` will be made available with it's original type and then converted to JuliaType when inserted into `dict`
    d = x => GLType # OpenGL convert target. Get's only applied if `x` is convertible to GLType. Will only be converted when passed to RenderObject
    d = x => \"doc string\"
    d = x => (GLType, \"doc string and gl target\")
end
"""
macro gen_defaults!(dict, args)
    args.head == :block || error("second argument needs to be a block of form
    begin
        a = 55
        b = a * 2 # variables, like a, will get made visible in local scope
        c::JuliaType = X # c needs to be of type JuliaType. c will be made available with it's original type and then converted to JuliaType when inserted into data
        d = x => GLType # OpenGL convert target. Get's only applied if x is convertible to GLType. Will only be converted when passed to RenderObject
    end")
    tuple_list = args.args
    dictsym = gensym()
    return_expression = Expr(:block)
    push!(return_expression.args, :($dictsym = $dict)) # dict could also be an expression, so we need to asign it to a variable at the beginning
    push!(return_expression.args, :(gl_convert_targets = get!($dictsym, :gl_convert_targets, Dict{Symbol, Any}()))) # exceptions for glconvert.
    push!(return_expression.args, :(doc_strings = get!($dictsym, :doc_string, Dict{Symbol, Any}()))) # exceptions for glconvert.
    # @gen_defaults can be used multiple times, so we need to reuse gl_convert_targets if already in here
    for (i,elem) in enumerate(tuple_list)
        elem.head == :line && continue
        opengl_convert_target = :() # is optional, so first is an empty expression
        convert_target        = :() # is optional, so first is an empty expression
        doc_strings           = :()
        if elem.head == :(=)
            key_name, value_expr = elem.args
            if isa(key_name, Expr) && key_name.head == :(::) # we need to convert to a julia type
                key_name, convert_target = key_name.args
                convert_target = :(convert($convert_target, $key_name))
            else
                convert_target = :($key_name)
            end
            key_sym = Expr(:quote, key_name)
            if isa(value_expr, Expr) && value_expr.head == :(=>)  # we might need to insert a convert target

                value_expr, target = value_expr.args
                undecided = []
                if isa(target, Expr)
                    undecided = target.args
                else
                    push!(undecided, target)
                end
                for elem in undecided
                    isa(elem, Expr) && continue #
                    if isa(elem, AbstractString) # only docstring
                        doc_strings = :(doc_strings[$key_sym] = $elem)
                    elseif isa(elem, Symbol)
                        opengl_convert_target = quote
                            if matches_target($elem, $key_name)
                                gl_convert_targets[$key_sym] = $elem
                            end
                        end
                    end
                end
            end
            expr = quote
                $key_name = if haskey($dictsym, $key_sym)
                    $dictsym[$key_sym]
                else
                    $value_expr # in case that evaluating value_expr is expensive, we use a branch instead of get(dict, key, default)
                end
                $dictsym[$key_sym] = $convert_target
                $opengl_convert_target
                $doc_strings
            end
            push!(return_expression.args, expr)
        else
            error("all nodes need to be of form a = b OR a::Type = b OR a = b => Type, where a needs to be a var and b any expression. Found: $elem")
        end
    end
    #push!(return_expression.args, :($dictsym[:gl_convert_targets] = gl_convert_targets)) #just pass the targets via the dict
    push!(return_expression.args, :($dictsym)) #return dict
    esc(return_expression)
end
export @gen_defaults!

value(any) = any # add this, to make it easier to work with a combination of signals and constants

makesignal(s::Signal) = s
makesignal(v)         = Signal(v)

@inline const_lift(f::Union{DataType, Function}, inputs...) = map(f, map(makesignal, inputs)...)
export const_lift

function close_to_square(n::Real)
    # a cannot be greater than the square root of n
    # b cannot be smaller than the square root of n
    # we get the maximum allowed value of a
    amax = floor(Int, sqrt(n));
    if 0 == rem(n, amax)
        # special case where n is a square number
        return (amax, div(n, amax))
    end
    # Get its prime factors of n
    primeFactors  = factor(n);
    # Start with a factor 1 in the list of candidates for a
    candidates = Int[1]
    for (f, _) in primeFactors
        # Add new candidates which are obtained by multiplying
        # existing candidates with the new prime factor f
        # Set union ensures that duplicate candidates are removed
        candidates  = union(candidates, f .* candidates);
        # throw out candidates which are larger than amax
        filter!(x-> x <= amax, candidates)
    end
    # Take the largest factor in the list d
    (candidates[end], div(n, candidates[end]))
end

isnotempty(x) = !isempty(x)
AND(a,b) = a&&b
OR(a,b) = a||b

#Uhm I should remove this. Needed for smooth transition between FixedSizeArrays and Number, though
Base.length{T <: Number}(::Type{T}) = 1

#Meshtype holding native OpenGL data.
immutable NativeMesh{MeshType <: HomogenousMesh}
    data::Dict{Symbol, Any}
end
export NativeMesh

@compat (::Type{NativeMesh}){T <: HomogenousMesh}(m::T) = NativeMesh{T}(m)

@compat function (MT::Type{NativeMesh{T}}){T <: HomogenousMesh}(m::T)
    result = Dict{Symbol, Any}()
    attribs = attributes(m)
    @materialize! vertices, faces = attribs
    result[:vertices] = GLBuffer(vertices)
    result[:faces]    = indexbuffer(faces)
    for (field, val) in attribs
        if field in (:texturecoordinates, :normals, :attribute_id, :color)
            if isa(val, Vector)
                result[field] = GLBuffer(val)
            end
        else
            result[field] = Texture(val)
        end
    end
    MT(result)
end

@compat function (MT::Type{NativeMesh{T}}){T <: HomogenousMesh}(m::Signal{T})
    result = Dict{Symbol, Any}()
    mv = value(m)
    attribs = attributes(mv)
    @materialize! vertices, faces = attribs
    result[:vertices] = GLBuffer(vertices)
    result[:faces]    = indexbuffer(faces)
    for (field, val) in attribs
        if field in (:texturecoordinates, :normals, :attribute_id, :color)
            if isa(val, Vector)
                result[field] = GLBuffer(val)
            end
        else
            result[field] = Texture(val)
        end
    end
    preserve(map(m) do mesh
        for (field, val) in attributes(mesh)
            if haskey(result, field)
                update!(result[field], val)
            end
        end
    end)
    MT(result)
end


text:
@testset "polyagammapoisson.jl" begin
    # TODO add more basic test
    d = PolyaGammaPoisson(2, 3.0, 3.0)
    @test length(d) == 2
end


text:
module Hockey

include("shared.jl")
include("json_pbp.jl")
include("json_shifts.jl")

import .JSON_PBP
import .JSON_Shifts

export scrape_game

function scrape_game(game_id::AbstractString, scrape_shifts::Bool=false)
    pbp = JSON_PBP.scrape_game(game_id)
    if scrape_shifts
        shifts = JSON_Shifts.scrape_shifts(game_id)
        return pbp, shifts
    else
        return pbp
    end
end

end


text:
function calc_submatrices(
  nbset  # tuple with the following elements.
    # nodes: array of Strings containing independent node names.
    # num_nodes: number of nodes.
    # branches: array of struct Branch.
    # num_branches: sum of branch phases.
    # num_dyns: number of dynamic branches.
  )  #=

calculates submatrices to build the sparse-tableau equations of the given circuit. The following variables are returned as the elements of the tuple "smset".
- Ac: branch-node connectivity matrix.
- Z: coefficient matrix of i.
- Y: coefficient matrix of v.
- B1: incidence matrix of alpha.
- B2: incidence matrix of [ i v ].

Author: Taku Noda
Started on Jan. 6, 2021
Update on June 29, 2021

Copyright 2021 Central Research Institute of Electric Power Industry

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

=#

( nodes, num_nodes, branches, num_branches, num_dyns ) = nbset
Ac = build_Ac( num_nodes, branches, num_branches )
Z, Y, B1, B2 = build_ZYB!( branches, num_branches, num_dyns )
smset = ( Ac, Z, Y, B1, B2 )

return smset

end

max_stars_repo_name: takunoda/omegasim

max_issues_repo_name: takunoda/omegasim

max_forks_repo_name: takunoda/omegasim


text:
# (C) Copyright IBM 2017, 2018.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.
#
# Modifications copyright Spencer King, 2020

using Qiskit

function build_bell_circuit()
    # Returns a circuit putting 2 qubits in the Bell state.
    qc = Qiskit.QuantumCircuit(2, 2)
    qc.h(0)
    qc.cx(0, 1)
    qc.measure([0, 1], [0, 1])
    return qc
end

# Initialize Qiskit
Qiskit.init()

# Create the circuit
bell_circuit = build_bell_circuit()

# Use the internal .draw() to print the circuit
print(bell_circuit)


text:
#! /usr/bin/julia

# Rosetta Code, Narcissistic decimal number

function isnarcissist{T<:Integer}(n::T, b::Int=10)
    -1 < n || return false
    d = digits(n, b)
    m = length(d)
    n == mapreduce((x)->x^m, +, d)
end

goal = 25
ncnt = 0
println("Finding the first ", goal, " Narcissistic numbers:")
for i in 0:typemax(1)
    isnarcissist(i) || continue
    ncnt += 1
    println(@sprintf "    %2d %7d" ncnt i)
    ncnt < goal || break
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


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

mutable struct FileStatusResult <: SwaggerModel
    FileStatus::Any # spec type: Union{ Nothing, FileStatusProperties } # spec name: FileStatus

    function FileStatusResult(;FileStatus=nothing)
        o = new()
        validate_property(FileStatusResult, Symbol("FileStatus"), FileStatus)
        setfield!(o, Symbol("FileStatus"), FileStatus)
        o
    end
end # type FileStatusResult

const _property_map_FileStatusResult = Dict{Symbol,Symbol}(Symbol("FileStatus")=>Symbol("FileStatus"))
const _property_types_FileStatusResult = Dict{Symbol,String}(Symbol("FileStatus")=>"FileStatusProperties")
Base.propertynames(::Type{ FileStatusResult }) = collect(keys(_property_map_FileStatusResult))
Swagger.property_type(::Type{ FileStatusResult }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_FileStatusResult[name]))}
Swagger.field_name(::Type{ FileStatusResult }, property_name::Symbol) =  _property_map_FileStatusResult[property_name]

function check_required(o::FileStatusResult)
    true
end

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


text:
    using DifferentialEquations
    using OrdinaryDiffEq
    using DiffEqBase
    using Sundials
    using ODEInterfaceDiffEq
    using Turing
    using MCMCChains

        # du -> Derivative
        # u -> State at time t
        # p-> paramter vector
        # t-> time as tuple (init, end)

        function PLacExampleODE!(du,u,p,t)
            Cit_mrna, Cit_foldedP, Cit_fluo, Cit_AU = u;

            alpha1, Vm1, h1, Km1, d1, alpha2, d2, Kf, sc_molec, IPTG = p;

        du[1] = dCit_mrna=alpha1+Vm1*(IPTG^h1/(Km1^h1+IPTG^h1))-d1*Cit_mrna;
        du[2] = dCit_foldedP=alpha2*Cit_mrna-(d2+Kf)*Cit_foldedP;
        du[3] = dCit_fluo=Kf*Cit_foldedP-d2*Cit_fluo;
        du[4] = dCit_AU = sc_molec*dCit_fluo;

        end

        # p -> parameter vector (plus inducer vector at the end)
        # I -> initial Y0 vector

        function PLacExampleSteadyState(p,I)

            expCit_mrna, expCit_foldedP, expCit_fluo, expCit_AU = I;
        alpha1, Vm1, h1, Km1, d1, alpha2, d2, Kf, sc_molec, IPTG = p;

         alp = zeros(4);

        alp[1] = Cit_mrna = (alpha1 + Vm1*(IPTG^h1/(Km1^h1+IPTG^h1)))/d1;;
        alp[2] = Cit_foldedP = (alpha2*Cit_mrna)/(Kf+d2);;
        alp[3] = Cit_fluo = (Kf*Cit_foldedP)/d2;;
        alp[4] = Cit_AU= sc_molec*Cit_fluo;;
        return(alp)

        end

        # ts -> time vector going from t=0 to t=end every 1
        # p -> parameter vector
        # sp -> vector with switching times (t=0 and t=end needs to be included). Length of vector will be length(inducer steps)+1
        # inputs -> Vecor with the inducers (if more than 1 inducer in the system, the shape will be a1,b1,c1,..., a2,b2,c2,... aend,bend,cend,... where each leter represents a different inducer)
        # ivss -> Initial Y0 vector
        # pre -> Vector with the inducers in the ON. It can be empty

        @model function PLacExample_TuringModel(data,prob1)

            p1 ~ Truncated(Distributions.Normal(0.2475194, 0.1237403), 3.88e-5, 0.495)
            p2 ~ Truncated(Distributions.Normal(0.2669, 0.11405), 0.0388, 0.495)
            p3 ~ Truncated(Distributions.Normal(2.7, 1.1), 0.5, 4.9)
            p4 ~ Truncated(Distributions.Normal(6.0, 2.0), 2.0, 10.0)
            p5 ~ Truncated(Distributions.Normal(0.11885000000000001, 0.055575), 0.0077, 0.23)
            p6 ~ Truncated(Distributions.Normal(3.525, 1.6408500000000001), 0.2433, 6.8067)
            p7 ~ Truncated(Distributions.Normal(0.1224799, 0.06121005), 5.98e-5, 0.2449)
            p8 ~ Truncated(Distributions.Normal(0.01685, 0.002425), 0.012, 0.0217)
            p9 ~ Truncated(Distributions.Normal(5.0005, 2.49975), 0.001, 10.0)

            p = [p1,p2,p3,p4,p5,p6,p7,p8,p9,];

            Neq = 4;

            for exp in 1:2

                ts=collect(0.:1:data["finalTime"][exp])
                sp = data["switchT"][exp];
                inputs = data["uInd"][exp];
                ivss = data["y0"][exp];
                pre=data["preInd"][exp];

                maxtime = length(ts);
                Nsp = length(sp);
                Nevents = length(sp)-1;

                if typeof(p) == Array{Float64, 1}
                    if pre != []
                        pSS = vcat(p, pre);
                    else
                        pSS = p;
                    end
                else
                    if pre != []
                        pSS = vcat([p[i].value for i in 1:9], pre);
                    else
                        pSS = [p[i].value for i in 1:9];
                    end
                end

                final = zeros(maxtime, Neq);
                finalBay2 = Array{Any,2}(undef, 1, maxtime)

                y_al = PLacExampleSteadyState(pSS,ivss) # Calculation of initial guesses for steady state
                y0 = y_al;
                initialV = y0;
                i = 1;

                for q in collect(1:Nevents)

                    lts = length(ts[(sp[q]+1):sp[q+1]+1]);  # General way to define the number of elements in each event series

                    Tevent = ts[(sp[q]+1):sp[q+1]+1];  # General way to extract the times of each event
                    I = inputs[i:(i+(0))];
                    pSte = vcat(p, I);

                    if q == 1
                        if typeof(p) == Array{Float64, 1}
                            prob = ODEProblem(PLacExampleODE!,initialV,(ts[(sp[q]+1)],ts[sp[q+1]+1]),pSte);
                            part1 = DifferentialEquations.solve(prob, DifferentialEquations.Tsit5(),reltol=1.0e-9,abstol=1.0e-9,saveat=1);
                        else
                            prob = remake(prob1, p=pSte)
                            prob = remake(prob, u0=initialV)
                            prob = remake(prob, tspan = convert.(Int, (ts[(sp[q]+1)],ts[sp[q+1]+1])))
                            part1 = DifferentialEquations.solve(prob,DifferentialEquations.Tsit5(),reltol=1.0e-9,abstol=1.0e-9,saveat=1)
                        end
                    else
                        if typeof(p) == Array{Float64, 1}
                            prob = ODEProblem(PLacExampleODE!,initialV,(ts[(sp[q]+1)],ts[sp[q+1]+1]),pSte);
                            part1 = DifferentialEquations.solve(prob, DifferentialEquations.Tsit5(),reltol=1.0e-9,abstol=1.0e-9,saveat=1);
                        else
                            prob = remake(prob1, p=pSte)
                            prob = remake(prob, u0=initialV)
                            prob = remake(prob, tspan = (ts[(sp[q]+1)],ts[sp[q+1]+1]))
                            part1 = DifferentialEquations.solve(prob,DifferentialEquations.Tsit5(),reltol=1.0e-9,abstol=1.0e-9,saveat=1)
                        end
                    end

                    if typeof(p) == Array{Float64, 1}
                        initialV = part1[lts];
                    else
                        initialV = [part1[:,lts][h].value for h in 1:4];
                    end

                    i+=1;

                    if typeof(p) == Array{Float64, 1}
                        for d in collect((sp[q]+1):(sp[q]+lts))
                            final[d,:] = part1[d-sp[q]];
                        end
                    else
                        for d in collect((sp[q]+1):(sp[q]+lts))
                            finalBay2[1,d] = 3 .* part1[4, d-sp[q]];

                        end
                    end

                end

                obs = Array{Any,1}(undef, 1)
                if typeof(p) == Array{Float64, 1}
                    obs[1] = 3 .* final[convert.(Int,round(data["tsamps"][exp])).+1, 4];

                else
                    obs[1] = finalBay2[1, convert.(Int,round(data["tsamps"][exp])).+1];
                end

                for ob in 1:1
                    for dat in 1:length(obs)
                        data["DataMean"][ob][dat] ~ Normal(obs[ob][dat], data["DataError"][ob][1][dat])
                    end
                end

            end

        end


text:
using GPLikelihoods
using AbstractGPs
using Test
using Random
using Functors
using Distributions

@testset "GPLikelihoods.jl" begin

    include("test_utils.jl")

    @testset "likelihoods" begin
        include("likelihoods/bernoulli.jl")
        include("likelihoods/binomial.jl")
        include("likelihoods/categorical.jl")
        include("likelihoods/exponential.jl")
        include("likelihoods/gaussian.jl")
        include("likelihoods/poisson.jl")
        include("likelihoods/studentsT.jl")
    end

end


text:
# Read mesh from a triangle generated file
function read_line(ln, types)
    m2 = collect((m.match for m = eachmatch(r"\b((\d*\.)?\d+)\b", ln, overlap=false)))
    return (parse(types[i],x) for (i,x) in enumerate(m2))
end

function parse_nodes!(nodes,root_file)
    open(root_file*".node") do f
        first_line = true
        for ln in eachline(f)
            m = match(r"^\s*(?:#|$)", ln)
            if m == nothing
                if (first_line)   #skip first line
                    first_line = false
                else
                    #parse nodes
                    pln = collect(read_line(ln, (Int,Float64,Float64,Int)))
                    node = Node((pln[2],pln[3]))
                    push!(nodes,node)
                end
            end
        end
    end
end

function parse_faces!(faces_ref,root_file)
    open(root_file*".edge") do f
        first_line = true
        for ln in eachline(f)
            m = match(r"^\s*(?:#|$)", ln)
            if m == nothing
                if (first_line)   #skip first line
                    first_line = false
                else
                    #parse nodes
                    pln = collect(read_line(ln, (Int,Int,Int,Int)))
                    element = (minimum(pln[2:3]),maximum(pln[2:3]))
                    token = ht_keyindex2!(faces_ref, element)
                    if token < 0
                        Base._setindex!(faces_ref, pln[4], element, -token)
                    end
                end
            end
        end
    end
end

function parse_cells!(cells, faces, faces_ref,facesets,nodesets,nodes, root_file)
    #read cell nodes
    open(root_file*".ele") do f
        first_line = true
        n_el = 0
        n_faces = 0
        boundary_faces = Set{Int}()
        boundary_nodes = Set{Int}()
        facesdict = Dict{NTuple{2,Int},Int}()
        el_faces = [0,0,0]
        for ln in eachline(f)
            m = match(r"^\s*(?:#|$)", ln)
            if m == nothing
                if (first_line)   #skip first line
                    first_line = false
                else
                    fill!(el_faces,0)
                    n_el = n_el + 1
                    #parse nodes
                    pln = collect(read_line(ln, (Int,Int,Int,Int)))
                    el_nodes = _check_node_data(nodes, pln[2:4]...)
                    #build faces
                    for (i,fn_id) in enumerate(((2,3),(3,1),(1,2)))
                        v1 = el_nodes[fn_id[1]]; v2 = el_nodes[fn_id[2]]
                        element = (min(v1,v2),max(v1,v2))
                        token = ht_keyindex2!(facesdict, element)
                        if token > 0
                            el_faces[i] = facesdict.vals[token]
                            if !(n_el == faces[facesdict.vals[token],4])
                                faces[facesdict.vals[token],4] = n_el
                            end
                        else
                            ref = -1
                            token2 = ht_keyindex2!(faces_ref, element)
                            if token2 > 0
                                ref = faces_ref.vals[token2]
                            end
                            n_faces = n_faces + 1
                            el_faces[i] = n_faces
                            faces[n_faces,1] = v1
                            faces[n_faces,2] = v2
                            faces[n_faces,3] = n_el
                            faces[n_faces,4] = 0
                            Base._setindex!(facesdict, n_faces, element, -token)
                            if ref > 0
                                push!(boundary_faces, n_faces)
                                push!(boundary_nodes, v1)
                                push!(boundary_nodes, v2)
                            end
                        end
                    end
                    #save cell
                    cell = Cell{2,3,3}(el_nodes, (el_faces[1],el_faces[2],el_faces[3]))
                    push!(cells, cell)
                end
            end
        end
        push!(facesets, "boundary" => boundary_faces)
        push!(nodesets, "boundary" => boundary_nodes)
    end
end

"""
function parse_mesh_triangle(root_file)
Reads mesh generated by triangle with root file name `root_file`
Ex: `parse_mesh_triangle("figure.1")`
"""
function parse_mesh_triangle(root_file)
    nodes = Vector{Node}()
    faces_ref = Dict{NTuple{2,Int},Int}()
    cells = Vector{TriangleCell}()
    facesets = Dict{String,Set{Int}}()
    nodesets = Dict{String,Set{Int}}()
    parse_nodes!(nodes,root_file)
    parse_faces!(faces_ref, root_file)
    faces = Matrix{Int}(undef,length(faces_ref),4)
    parse_cells!(cells, faces, faces_ref,facesets,nodesets,nodes, root_file)
    PolygonalMesh{2,3,3,2,eltype(nodes[1].x)}(cells, nodes, faces, facesets,nodesets)
end


text:
module Batches

import ..Flux

include("batch.jl")

end


text:
function test_epsilon_dominance_func(epsilon_dominates_func)
  describe("with scalar epsilon") do
    test("simple cases with float vector fitness values of length 1") do
      @check epsilon_dominates_func([0.0], [0.0], 1.0)
      @check epsilon_dominates_func([0.99], [0.0], 1.0)

      @check epsilon_dominates_func([-1.0], [-1.0], 1.0)
      @check epsilon_dominates_func([-1.0], [-1.0], 0.00001)

      @check !epsilon_dominates_func([1.0], [0.0], 1.0)
    end

    test("simple cases with float vector fitness values of length 2") do
      @check epsilon_dominates_func([1.0, 2.0], [1.0, 2.0], 0.10)
      @check epsilon_dominates_func([1.1, 2.0], [1.0, 2.0], 0.10)
      @check epsilon_dominates_func([0.0, 0.99], [1.0, 2.0], 1.0)
      @check epsilon_dominates_func([-1.0, 0.99], [0.0, 2.0], 1.0)
      @check epsilon_dominates_func([0.99, -10.0], [2.0, -9.0], 1.0)

      @check !epsilon_dominates_func([1.1, 2.1], [1.0, 2.0], 0.10)
      @check !epsilon_dominates_func([2.1, 1.1], [2.0, 1.0], 0.10)
      @check !epsilon_dominates_func([0.1, 1.1], [0.0, 1.0], 0.10)
      @check !epsilon_dominates_func([-0.9, 0.1], [-1.0, 0.0], 0.10)
      @check !epsilon_dominates_func([2.0, 3.0], [1.0, 2.0], 0.50)
    end

    test("simple cases with float vector fitness values of length 3") do
      @check epsilon_dominates_func([0.0, 1.0, 2.0], [1.0, 2.0, 3.0], 1.0)

      @check !epsilon_dominates_func([1.0, 2.0, 3.0], [0.0, 1.0, 2.0], 1.0)

      @check epsilon_dominates_func([0.99, 2.0,  3.0],  [0.0, 1.0, 2.0], 1.0)
      @check epsilon_dominates_func([1.0,  1.99, 3.0],  [0.0, 1.0, 2.0], 1.0)
      @check epsilon_dominates_func([1.0,  2.0,  2.99], [0.0, 1.0, 2.0], 1.0)
    end
  end

  describe("with vector epsilon") do
    test("float vector fitness values of length 2") do
      @check epsilon_dominates_func([1.0, 2.0], [1.0, 2.0], [0.1, 0.1])
      @check epsilon_dominates_func([1.09, 2.0], [1.0, 2.0], [0.1, 0.1])
      @check epsilon_dominates_func([1.1, 2.0], [1.0, 2.0], [0.1, 0.1])
      @check epsilon_dominates_func([1.1, 2.19], [1.0, 2.0], [0.1, 0.2])
      @check epsilon_dominates_func([-1.0, 2.0], [-1.0, 2.0], [0.1, 0.1])
      @check epsilon_dominates_func([-1.0, -2.0], [-1.0, -2.0], [0.1, 0.1])

      @check !epsilon_dominates_func([1.1, 2.1], [1.0, 2.0], [0.1, 0.1])
      @check !epsilon_dominates_func([1.1, 2.2], [1.0, 2.0], [0.1, 0.2])
      @check !epsilon_dominates_func([-1.0, -2.0], [-1.1, -2.5], [0.1, 0.5])
    end
  end
end

describe("epsilon_dominance") do
  describe("epsilon_dominance_clear") do
    test_epsilon_dominance_func(epsilon_dominates_clear)
  end

  describe("epsilon_dominance_fast") do
    test_epsilon_dominance_func(epsilon_dominates_fast)
  end
end


text:
using OptEvapCool
using StatsPlots
using Printf: @sprintf
using LaTeXStrings

function under_beam_trap(duration = 1.3, input_dir = "")
    # Physical parameters
    Np = 3e7 # Initial atom count
    T₀ = 15e-6 # Initial temperature
    species = Rb87

    # Numerical parameters
    Nt = ceil(Int64, 1e5) # Test particles
    F = Np / Nt
    Nc = 3 # Target number of test particles per cell

    # Beam parameters
    P₁(t) = 20
    P₂(t) = 10
    P_g(t) = (t < 1) ? 15 : 15 + 9*(t-1)/0.3

    w₀ = 130e-6 # Waist (m)
    wg = 100e-6
    θ = ( 22.5 * π / 180 ) / 2 # Half-angle between beams
    ϕ = 5 * π / 180 # Tilt angle of under-beam

    dir1 = [sin(θ), 0, cos(θ)] # Directions
    dir2 = [-sin(θ), 0, cos(θ)]
    dir_g = [0, sin(ϕ), cos(ϕ)]

    λ₁ = 1064e-9 # Wavelengths
    λ₂ = 1090e-9
    λ_g = 1090e-9

    foc_g(t) = (t < 1) ? [0, -0.00035 + t*0.00019, 0] : [0, -0.00016, 0]

    beam1 = GaussianBeam([0,0,0], dir1, P₁, w₀, λ₁)
    beam2 = GaussianBeam([0,0,0], dir2, P₂, w₀, λ₂)
    beam_g = GaussianBeam(foc_g, dir_g, P_g, wg, λ_g)

    # Trapping frequencies
    m = species.m
    κ = kappa(species)

    Uₜ_coeff = 2 * κ / (π * w₀^2)
    Uₜ(t) = Uₜ_coeff * (P₁(t) + P₂(t)) #Trap depth (without ghost beam)

    ωx_coeff = sqrt(4 * cos(θ)^2 / (m * w₀^2))
    ωz_coeff = sqrt(4 * sin(θ)^2 / (m * w₀^2))
    ωy_coeff = sqrt(4 / (m * w₀^2))

    ωx(t) = ωx_coeff * sqrt(Uₜ(t))
    ωz(t) = ωz_coeff * sqrt(Uₜ(t))
    ωy(t) = ωy_coeff * sqrt(Uₜ(t))

    # Cloud initialisation
    positions = harmonic_boltzmann_positions(Nt, m, T₀, ωx(0), ωy(0), ωz(0))
    velocities = boltzmann_velocities(Nt, m, T₀)

    # Adjust vertical equilibrium position
    function y0(t)
        q = -4*g^2 / (ωy(t)^4 * w₀^2)
        return - w₀/2 * sqrt(- q * exp(q))
    end

    centre(t) = [0, y0(t), 0]
    positions .+= centre(0)

    # Three-body loss & background loss
    K = 4.3e-29 * 1e-12
    τ_bg = 180

    # Maximum timestep
    max_dt(t) = 0.05 * 2π / max(ωx(t), ωy(t), ωz(t))

    # Gravitaty
    acc_grav! = acceleration(gravity)
    pot_grav! = potential(gravity)

    # GAUSSIAN BEAM SIMULATION
    acc_b1! = acceleration(beam1)
    acc_b2! = acceleration(beam2)
    acc_g! = acceleration(beam_g)
    pot_b1! = potential(beam1)
    pot_b2! = potential(beam2)
    pot_g! = potential(beam_g)

    function accel(p, s, t, o)
        fill!(o, 0.0)
        acc_b1!(p, s, t, o)
        acc_b2!(p, s, t, o)
        acc_g!(p, s, t, o)
        acc_grav!(p, s, t, o)
        return nothing
    end

    function poten(p, s, t, o)
        fill!(o, 0.0)
        pot_b1!(p, s, t, o)
        pot_b2!(p, s, t, o)
        pot_g!(p, s, t, o)
        pot_grav!(p, s, t, o)
        return nothing
    end

    sensor = GlobalSensor()
    measure = measurer(sensor, 0.001, ωx, ωy, ωz, centre)

    T(t) = (length(sensor.ke) > 0) ? 2 * last(sensor.ke) / (3 * kB) : T₀
    evap = ellipsoid_evap(ωx, ωy, ωz, T, 1e-6, centre)

    conditions = SimulationConditions(species, F, positions, velocities,
        accel, poten, evap = evap, τ_bg = τ_bg, K = K)

    if input_dir == "" # Simulate new results
        final_cloud = evolve(conditions, duration;
            Nc = Nc, max_dt = max_dt, measure = measure)
    else # Load provided results
        sensor = loadsensor("$input_dir/sensor-data.csv")
    end

    # Save sensor data before plotting
    ft = filetime()
    dir = "./results/$ft-under-beam-trap"
    mkpath(dir)
    savecsv(sensor, "$dir/sensor-data.csv")

    # PLOTTING
    default(fontfamily="Computer Modern",
        linewidth=2, framestyle=:box, label=nothing, grid=false)
    #scalefontsizes(1.3)

    window_time = min(duration/2, 1e-2)
    time = sensor.time
    window_size = OptEvapCool.window_time_size(time, window_time)
    rolling_time = rollmean(time, window_size)

    # Colour palette
    col1 = RGB(0.0588, 0, 0.5882)
    #col2 = RGB(0.9000, 0, 0.3490)
    #col3 = RGB(1, 0.6510, 0)
    #col4 = RGB(0.0824, 0.8392, 0)

    # Temperature
    rolling_temp = 2 / (3 * kB) * rollmean(sensor.ke, window_size)

    temp_order = -6
    temp_yformatter(y) = @sprintf("%.2f",y/(10.0^temp_order))

    temperature_plt = plot(rolling_time, rolling_temp,
        xlabel = "Time (s)",
        ylabel = L"\textrm{Temperature\ \ }({}\times10^{%$temp_order}\mathrm{K})",
        label = false,
        linecolor = col1,
        ls = :solid,
        #yscale = :log10,
        yminorticks = 10,
        yformatter = temp_yformatter,
        ylims = (0, Inf),
        dpi = 300)

    # Number
    Np = sensor.Nt .* sensor.F

    num_order = 6
    num_yformatter(y) = @sprintf("%.2f",y/(10.0^num_order))

    number_plt = plot(time, Np,
        xlabel = "Time (s)",
        ylabel = L"\textrm{Number\ \ }({}\times10^{%$num_order})",
        label = false,
        linecolor = col1,
        ls = :solid,
        #yscale = :log10,
        yminorticks=10,
        yformatter = num_yformatter,
        ylims = (0, Inf),
        dpi = 300)

    # Collrate
    rolling_timesteps = (
        time[window_size:end] - time[1:end - window_size + 1]
    )
    rolling_Nt = rollmean(sensor.Nt, window_size)
    rolling_Γ = (2 * rolling(sum, sensor.coll, window_size) ./
        (rolling_Nt .* rolling_timesteps)
    )

    collrate_plt = plot(rolling_time, rolling_Γ,
        xlabel = "Time (s)",
        ylabel = "Collision rate (Hz)",
        label = false,
        linecolor = col1,
        ls = :solid,
        dpi = 300)

    # Phase space density
    rolling_n0 = rollmean(sensor.n0 .* sensor.F, window_size)
    rolling_psd = rolling_n0 .* (
        (2*π*h̄^2 ./ (m * kB .* rolling_temp)) .^1.5
    )

    psd_plt = plot(rolling_time, rolling_psd,
        xlabel = "Time (s)",
        ylabel = "Phase space density",
        label = false,
        linecolor = col1,
        ylims = (0, Inf),
        ls = :solid,
        dpi = 300)

    # Save plots
    savefig(temperature_plt, "$dir/temp.png")
    savefig(collrate_plt, "$dir/collrate.png")
    savefig(number_plt, "$dir/number.png")
    savefig(psd_plt, "$dir/psd.png")

    return nothing
end