text:
import Test: @test, @testset
import FFTW
import Luna: Grid, Processing, Maths, Fields, settings, PhysData
import Luna.PhysData: wlfreq
import NumericalIntegration: integrate

@testset "normalisation" begin
λ0 = 800e-9
τfwhm = 30e-15
energy = 1e-3
ω0 = wlfreq(λ0)

rg = Grid.RealGrid(1, λ0, (200e-9, 3000e-9), 0.5e-12)
rFT = FFTW.plan_rfft(similar(rg.t), flags=settings["fftw_flag"])
eg = Grid.EnvGrid(1, λ0, (200e-9, 3000e-9), 0.5e-12)
eFT = FFTW.plan_fft(similar(eg.t), flags=settings["fftw_flag"])

itr = ((rg, rFT), (eg, eFT))

for ii = 1:2
    grid, FT = itr[ii]
    input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy)
    energy_t, energy_ω = Fields.energyfuncs(grid)
    Eω = input(grid, FT)
    @test energy_ω(Eω) ≈ energy

    ω, Iω = Processing.getIω(Processing.getEω(grid, Eω)..., :ω)
    f, If = Processing.getIω(Processing.getEω(grid, Eω)..., :f)
    λ, Iλ = Processing.getIω(Processing.getEω(grid, Eω)..., :λ)

    @test integrate(ω, Iω) ≈ energy
    @test integrate(f, If) ≈ energy
    @test isapprox(integrate(λ, Iλ), energy, rtol=1e-4)

    t, Et = Processing.getEt(grid, Eω)
    Pp = 0.94*energy/τfwhm # approximate peak power of gaussian pulse
    @test isapprox(maximum(abs2.(Et)), Pp, rtol=1e-3)
end
end

@testset "arrivaltime" begin
λ0 = 800e-9
ω0 = wlfreq(λ0)

# field
grid = Grid.RealGrid(1, λ0, (200e-9, 3000e-9), 0.5e-12)
Et = zero(grid.t)
It = Maths.gauss.(grid.t, fwhm=20e-15, x0=5e-15)
@. Et += sqrt(It)*cos(grid.t*ω0)
It = Maths.gauss.(grid.t, fwhm=15e-15, x0=-5e-15)
@. Et += sqrt(It)*cos(grid.t*ω0*2)
Eω = FFTW.rfft(Et)
@test isapprox(Processing.arrivaltime(grid, Eω, bandpass=(600e-9, 1000e-9)), 5e-15, rtol=1e-8)
@test isapprox(Processing.arrivaltime(grid, Eω, bandpass=(300e-9, 500e-9)), -5e-15, rtol=1e-8)
δt = grid.t[2] - grid.t[1]
@test abs(Processing.arrivaltime(grid, Eω, bandpass=(600e-9, 1000e-9); method=:peak) - 5e-15) < δt
@test abs(Processing.arrivaltime(grid, Eω, bandpass=(300e-9, 500e-9); method=:peak) + 5e-15) < δt

Et = hcat(zero(grid.t), zero(grid.t))
It = Maths.gauss.(grid.t, fwhm=20e-15, x0=5e-15)
Et[:, 1] .+= sqrt.(It).*cos.(grid.t*ω0)
It = Maths.gauss.(grid.t, fwhm=15e-15, x0=-5e-15)
Et[:, 1] .+= sqrt.(It).*cos.(grid.t*ω0*2)

It = Maths.gauss.(grid.t, fwhm=20e-15, x0=-5e-15)
Et[:, 2] .+= sqrt.(It).*cos.(grid.t*ω0)
It = Maths.gauss.(grid.t, fwhm=15e-15, x0=5e-15)
Et[:, 2] .+= sqrt.(It).*cos.(grid.t*ω0*2)
Eω = FFTW.rfft(Et, 1)
@test isapprox(Processing.arrivaltime(grid, Eω, bandpass=(600e-9, 1000e-9))[1], 5e-15, rtol=1e-8)
@test isapprox(Processing.arrivaltime(grid, Eω, bandpass=(600e-9, 1000e-9))[2], -5e-15, rtol=1e-8)
@test isapprox(Processing.arrivaltime(grid, Eω, bandpass=(300e-9, 500e-9))[1], -5e-15, rtol=1e-8)
@test isapprox(Processing.arrivaltime(grid, Eω, bandpass=(300e-9, 500e-9))[2], 5e-15, rtol=1e-8)

# envelope
grid = Grid.EnvGrid(1, λ0, (200e-9, 3000e-9), 0.5e-12)
Et = complex(zero(grid.t))
It = Maths.gauss.(grid.t, fwhm=20e-15, x0=5e-15)
@. Et += sqrt(It)
It = Maths.gauss.(grid.t, fwhm=15e-15, x0=-5e-15)
@. Et += sqrt(It)*exp(1im*grid.t*ω0)
Eω = FFTW.fft(Et)
@test isapprox(Processing.arrivaltime(grid, Eω, bandpass=(600e-9, 1000e-9)), 5e-15, rtol=1e-8)
@test isapprox(Processing.arrivaltime(grid, Eω, bandpass=(300e-9, 500e-9)), -5e-15, rtol=1e-8)
δt = grid.t[2] - grid.t[1]
@test abs(Processing.arrivaltime(grid, Eω, bandpass=(600e-9, 1000e-9); method=:peak) - 5e-15) < δt
@test abs(Processing.arrivaltime(grid, Eω, bandpass=(300e-9, 500e-9); method=:peak) + 5e-15) < δt

Et = hcat(complex(zero(grid.t)), complex(zero(grid.t)))
It = Maths.gauss.(grid.t, fwhm=20e-15, x0=5e-15)
Et[:, 1] .+= sqrt.(It)
It = Maths.gauss.(grid.t, fwhm=15e-15, x0=-5e-15)
Et[:, 1] .+= sqrt.(It).*exp.(1im*grid.t*ω0)

It = Maths.gauss.(grid.t, fwhm=20e-15, x0=-5e-15)
Et[:, 2] .+= sqrt.(It)
It = Maths.gauss.(grid.t, fwhm=15e-15, x0=5e-15)
Et[:, 2] .+= sqrt.(It).*exp.(1im*grid.t*ω0)
Eω = FFTW.fft(Et, 1)
@test isapprox(Processing.arrivaltime(grid, Eω, bandpass=(600e-9, 1000e-9))[1], 5e-15, rtol=1e-8)
@test isapprox(Processing.arrivaltime(grid, Eω, bandpass=(600e-9, 1000e-9))[2], -5e-15, rtol=1e-8)
@test isapprox(Processing.arrivaltime(grid, Eω, bandpass=(300e-9, 500e-9))[1], -5e-15, rtol=1e-8)
@test isapprox(Processing.arrivaltime(grid, Eω, bandpass=(300e-9, 500e-9))[2], 5e-15, rtol=1e-8)
end

@testset "specres" begin
# field grid
grid = Grid.RealGrid(1.0, 800e-9, (160e-9, 3000e-9), 30e-12)
Eω = (Maths.gauss.(grid.ω, fwhm=2π*PhysData.c/800e-9^2*2e-9, x0=PhysData.wlfreq(grid.referenceλ))
      .+ Maths.gauss.(grid.ω, fwhm=2π*PhysData.c/367e-9^2*20e-9, x0=PhysData.wlfreq(367e-9))
      .+ Maths.gauss.(grid.ω, fwhm=2π*PhysData.c/2000e-9^2*3e-9, x0=PhysData.wlfreq(2000e-9)))
# test energy integral is correct
ωp, Eωp = Processing.getEω(grid, Eω)
λg, Pλ = Processing.getIω(ωp, Eωp, :λ, resolution=10e-9, specrange=(170e-9, 2900e-9))
@test isapprox(Fields.energyfuncs(grid)[2](Eω) / integrate(λg, Pλ), 1.0, rtol=1e-10)
Fg, Pf = Processing.getIω(ωp, Eωp, :f, resolution=5e12)
@test isapprox(Fields.energyfuncs(grid)[2](Eω) / integrate(Fg, Pf), 1.0, rtol=1e-10)
# test resolution
Eω = (Maths.gauss.(grid.ω, fwhm=2π*PhysData.c/800e-9^2*0.1e-9, x0=PhysData.wlfreq(grid.referenceλ)))
ωp, Eωp = Processing.getEω(grid, Eω)
for res in (2e-9, 10e-9, 20e-9)
    λg, Pλ = Processing.getIω(ωp, Eωp, :λ, resolution=res, specrange=(170e-9, 2900e-9))
    @test isapprox(Maths.fwhm(λg, Pλ), res, rtol=1e-2)
end
for res in (1e12, 5e12, 10e12)
    Fg, Pf = Processing.getIω(ωp, Eωp, :f, resolution=res)
    @test isapprox(Maths.fwhm(Fg, Pf), res, rtol=1e-2)
end
# envelope grid
grid = Grid.EnvGrid(1.0, 800e-9, (160e-9, 3000e-9), 30e-12)
Eω = (Maths.gauss.(grid.ω, fwhm=2π*PhysData.c/800e-9^2*2e-9, x0=PhysData.wlfreq(grid.referenceλ))
      .+ Maths.gauss.(grid.ω, fwhm=2π*PhysData.c/367e-9^2*20e-9, x0=PhysData.wlfreq(367e-9))
      .+ Maths.gauss.(grid.ω, fwhm=2π*PhysData.c/2000e-9^2*3e-9, x0=PhysData.wlfreq(2000e-9)))
# test energy integral is correct
ωp, Eωp = Processing.getEω(grid, Eω)
λg, Pλ = Processing.getIω(ωp, Eωp, :λ, resolution=10e-9, specrange=(170e-9, 2900e-9))
@test isapprox(Fields.energyfuncs(grid)[2](Eω) / integrate(λg, Pλ), 1.0, rtol=1e-11)
Fg, Pf = Processing.getIω(ωp, Eωp, :f, resolution=5e12)
@test isapprox(Fields.energyfuncs(grid)[2](Eω) / integrate(Fg, Pf), 1.0, rtol=1e-11)
# test resolution
Eω = (Maths.gauss.(grid.ω, fwhm=2π*PhysData.c/800e-9^2*0.1e-9, x0=PhysData.wlfreq(grid.referenceλ)))
ωp, Eωp = Processing.getEω(grid, Eω)
for res in (2e-9, 10e-9, 20e-9)
    λg, Pλ = Processing.getIω(ωp, Eωp, :λ, resolution=res, specrange=(170e-9, 2900e-9))
    @test isapprox(Maths.fwhm(λg, Pλ), res, rtol=1e-2)
end
for res in (1e12, 5e12, 10e12)
    Fg, Pf = Processing.getIω(ωp, Eωp, :f, resolution=res)
    @test isapprox(Maths.fwhm(Fg, Pf), res, rtol=1e-2)
end
end

@testset "intensity autocorrelation" begin
    grid = Grid.RealGrid(1.0, 800e-9, (160e-9, 3000e-9), 1e-12)
    x = Array{Float64}(undef, length(grid.t))
    FT = FFTW.plan_rfft(x, 1)
    input = Fields.GaussField(λ0=800e-9, τfwhm=30e-15, energy=1e-6)
    Eω = input(grid, FT)
    Et = FT \ Eω
    IAC = Processing.intensity_autocorrelation(Et, grid)
    # Gaussian IAC width is sqrt(2)*FWHM
    @test isapprox(Maths.fwhm(grid.t, IAC), sqrt(2)*30e-15, rtol=3e-6)
    input = Fields.SechField(λ0=800e-9, τfwhm=30e-15, energy=1e-6)
    Eω = input(grid, FT)
    Et = FT \ Eω
    IAC = Processing.intensity_autocorrelation(Et, grid)
    # Numerical AC width from Mathematics: 1.542811*FWHM
    @test isapprox(Maths.fwhm(grid.t, IAC), 1.542811*30e-15, rtol=1e-6)
end

@testset "field autocorrelation" begin
    grid = Grid.EnvGrid(1.0, 800e-9, (160e-9, 3000e-9), 1e-12)
    Δω = 2π*PhysData.c/800e-9^2*10e-9
    Eω = sqrt.(Maths.gauss.(grid.ω, Δω, x0=PhysData.wlfreq(grid.referenceλ)))
    x = Array{Float64}(undef, length(grid.t))
    FT = FFTW.plan_fft(x)
    Et = FFTW.fftshift(FT \ Eω)
    τc = Processing.coherence_time(grid, Et)
    # Gaussian spectrum with natural width Δω has coherence time of 2*sqrt(log(2))/Δω (analytic)
    @test isapprox(τc, 2*sqrt(log(2))/Δω, rtol=6e-6)
    grid = Grid.RealGrid(1.0, 800e-9, (160e-9, 3000e-9), 1e-12)
    Δω = 2π*PhysData.c/800e-9^2*10e-9
    Eω = sqrt.(Maths.gauss.(grid.ω, Δω, x0=PhysData.wlfreq(grid.referenceλ)))
    x = Array{Float64}(undef, length(grid.t))
    FT = FFTW.plan_rfft(x, 1)
    Et = FFTW.fftshift(FT \ Eω)
    τc = Processing.coherence_time(grid, Et)
    # Gaussian spectrum with natural width Δω has coherence time of 2*sqrt(log(2))/Δω (analytic)
    @test isapprox(τc, 2*sqrt(log(2))/Δω, rtol=5e-8)
end

@testset "AutoWindow" begin
    λ0 = 800e-9
    grid = Grid.RealGrid(1, λ0, (100e-9, 3000e-9), 0.5e-12)
    FT = FFTW.plan_rfft(similar(grid.t), flags=settings["fftw_flag"])
    Eω1 = Fields.GaussField(λ0=λ0, τfwhm=10e-15, energy=100e-6)(grid, FT)
    Eω2 = Fields.GaussField(λ0=500e-9, τfwhm=30e-15, energy=100e-6)(grid, FT)
    Eω3 = Fields.GaussField(λ0=200e-9, τfwhm=3e-15, energy=10e-6)(grid, FT)

    pw = Processing.PeakWindow(0.2, 100e-9, 1200e-9; relative=true)
    Eωw = pw(grid.ω, Eω1)
    @test grid.ω[argmax(abs2.(Eωw))] == grid.ω[argmax(abs2.(Eω1))]

    pw = Processing.PeakWindow(0.2, 100e-9, 700e-9; relative=true)
    Eωw = pw(grid.ω, Eω1+Eω2+Eω3)
    @test grid.ω[argmax(abs2.(Eωw))] == grid.ω[argmax(abs2.(Eω2))]

    pw = Processing.PeakWindow(0.2, 100e-9, 400e-9; relative=true)
    Eωw = pw(grid.ω, Eω1+Eω2+Eω3)
    @test grid.ω[argmax(abs2.(Eωw))] == grid.ω[argmax(abs2.(Eω3))]

    pw = Processing.PeakWindow(0.2, 100e-9, 1200e-9; relative=true)
    Eωw = pw(grid.ω, hcat(Eω1, Eω2, Eω3))
    for (ii, Eω) in enumerate((Eω1, Eω2, Eω3))
        λ0 = wlfreq(grid.ω[argmax(abs2.(Eω))])
        l = [0.9, 1.1] .* λ0
        @test l == pw.lims[:, ii]
    end
end

##
@testset "scanproc" begin
    a = 125e-6
    flength = 0.1
    gas = :HeJ

    λ0 = 800e-9
    τfwhm = 10e-15

    λlims = (100e-9, 4e-6)
    trange = 400e-15

    # Scan dimensions:
    energies = collect(range(50e-9, 200e-9; length=8))
    pressures = collect(0.6:0.4:1.4)
    scan = Scan("pressure_energy_example"; energy=energies, pressure=pressures)
    mktempdir() do td
        runscan(scan) do scanidx, energy, pressure
            prop_capillary(a, flength, gas, pressure; λ0, τfwhm, energy,
                           λlims, trange, plasma=false, loss=false, scan, scanidx, filepath=td)
        end
        # multi-output scanproc
        λ, Iλ, zstat, max_peakpower = Processing.scanproc(td) do output
            λ, Iλ = Processing.getIω(output, :λ)
            zstat = Processing.VarLength(output["stats"]["z"])
            max_peakpower = maximum(output["stats"]["peakpower"])
            Processing.Common(λ), Iλ[:, end], zstat, max_peakpower
        end
        @test ndims(λ) == 1
        @test size(Iλ) == (length(λ), length(energies), length(pressures))
        @test size(zstat) == (length(energies), length(pressures))
        @test size(max_peakpower) == (length(energies), length(pressures))
        @test eltype(zstat) == Vector{Float64}

        #single-output scanproc
        energy_in = Processing.scanproc(td) do output
            output["stats"]["energy"][1]
        end
        @test size(energy_in) == (length(energies), length(pressures))
        @test all(isapprox.(energy_in[:, 1], energies, rtol=1e-4))

        # single-array-output scanproc
        Iλ2 = Processing.scanproc(td) do output
            local _, Iλ = Processing.getIω(output, :λ, flength)
            Iλ[:, 1]
        end
        @test Iλ == Iλ2

        #string output scanproc
        strings = Processing.scanproc(td) do output
            "a string"
        end
        @test size(strings) == (length(energies), length(pressures))
    end
end


text:
@testset "simple.jl" begin

    # Load all required packages and define likelihood functions
    ℓ(x) = logpdf(Normal(x, 0.5), 1.0)
    ℓvec(x) = logpdf(MvNormal(x, 0.25 * I), [1.0])
    @everywhere begin
        ℓ(x) = logpdf(Normal(x, 0.5), 1.0)
        ℓvec(x) = logpdf(MvNormal(x, 0.25 * I), [1.0])
    end

    @testset "Scalar model" begin
        Random.seed!(1)

        # model
        prior = Normal(0, 1)

        # true posterior
        μ = 0.8
        σ² = 0.2

        # regular sampling
        samples = sample(ESSModel(prior, ℓ), ESS(), 2_000; progress=false)
        @test samples isa Vector{Float64}
        @test length(samples) == 2_000
        @test mean(samples) ≈ μ atol = 0.05
        @test var(samples) ≈ σ² atol = 0.05

        # parallel sampling
        for alg in (MCMCThreads(), MCMCDistributed(), MCMCSerial())
            samples = sample(ESSModel(prior, ℓ), ESS(), alg, 2_000, 5; progress=false)
            @test samples isa Vector{Vector{Float64}}
            @test length(samples) == 5
            @test all(length(x) == 2_000 for x in samples)
            @test mean(mean, samples) ≈ μ atol = 0.05
            @test mean(var, samples) ≈ σ² atol = 0.05

            # initial parameter
            init_x = randn(5)
            samples = sample(
                ESSModel(prior, ℓ), ESS(), alg, 10, 5; progress=false, init_params=init_x
            )
            @test map(first, samples) == init_x
        end

        # initial parameter
        init_x = randn()
        samples = sample(ESSModel(prior, ℓ), ESS(), 10; progress=false, init_params=init_x)
        @test first(samples) == init_x
    end

    @testset "Scalar model with nonzero mean" begin
        Random.seed!(1)

        # model
        prior = Normal(0.5, 1)

        # true posterior
        μ = 0.9
        σ² = 0.2

        # regular sampling
        samples = sample(ESSModel(prior, ℓ), ESS(), 2_000; progress=false)
        @test samples isa Vector{Float64}
        @test length(samples) == 2_000
        @test mean(samples) ≈ μ atol = 0.05
        @test var(samples) ≈ σ² atol = 0.05

        # parallel sampling
        for alg in (MCMCThreads(), MCMCDistributed(), MCMCSerial())
            samples = sample(ESSModel(prior, ℓ), ESS(), alg, 2_000, 5; progress=false)
            @test samples isa Vector{Vector{Float64}}
            @test length(samples) == 5
            @test all(length(x) == 2_000 for x in samples)
            @test mean(mean, samples) ≈ μ atol = 0.05
            @test mean(var, samples) ≈ σ² atol = 0.05

            # initial parameter
            init_x = randn(5)
            samples = sample(
                ESSModel(prior, ℓ), ESS(), alg, 10, 5; progress=false, init_params=init_x
            )
            @test map(first, samples) == init_x
        end

        # initial parameter
        init_x = randn()
        samples = sample(ESSModel(prior, ℓ), ESS(), 10; progress=false, init_params=init_x)
        @test first(samples) == init_x
    end

    @testset "Scalar model (vectorized)" begin
        Random.seed!(1)

        # model
        prior = MvNormal([0.0], I)

        # true posterior
        μ = [0.8]
        σ² = [0.2]

        # regular sampling
        samples = sample(ESSModel(prior, ℓvec), ESS(), 2_000; progress=false)
        @test samples isa Vector{Vector{Float64}}
        @test length(samples) == 2_000
        @test all(length(x) == 1 for x in samples)
        @test mean(samples) ≈ μ atol = 0.05
        @test var(samples) ≈ σ² atol = 0.05

        # parallel sampling
        for alg in (MCMCThreads(), MCMCDistributed(), MCMCSerial())
            samples = sample(ESSModel(prior, ℓvec), ESS(), alg, 2_000, 5; progress=false)
            @test samples isa Vector{Vector{Vector{Float64}}}
            @test length(samples) == 5
            @test all(length(x) == 2_000 for x in samples)
            @test mean(mean, samples) ≈ μ atol = 0.05
            @test mean(var, samples) ≈ σ² atol = 0.05

            # initial parameter
            init_x = [randn(1) for _ in 1:5]
            samples = sample(
                ESSModel(prior, ℓvec), ESS(), alg, 10, 5; progress=false, init_params=init_x
            )
            @test map(first, samples) == init_x
        end

        # initial parameter
        init_x = randn(1)
        samples = sample(
            ESSModel(prior, ℓvec), ESS(), 10; progress=false, init_params=init_x
        )
        @test first(samples) == init_x
    end

    @testset "Scalar model with nonzero mean (vectorized)" begin
        Random.seed!(1)

        # model
        prior = MvNormal([0.5], I)

        # true posterior
        μ = [0.9]
        σ² = [0.2]

        # regular sampling
        samples = sample(ESSModel(prior, ℓvec), ESS(), 2_000; progress=false)
        @test samples isa Vector{Vector{Float64}}
        @test length(samples) == 2_000
        @test all(length(x) == 1 for x in samples)
        @test mean(samples) ≈ μ atol = 0.05
        @test var(samples) ≈ σ² atol = 0.05

        # parallel sampling
        for alg in (MCMCThreads(), MCMCDistributed(), MCMCSerial())
            samples = sample(ESSModel(prior, ℓvec), ESS(), alg, 2_000, 5; progress=false)
            @test samples isa Vector{Vector{Vector{Float64}}}
            @test length(samples) == 5
            @test all(length(x) == 2_000 for x in samples)
            @test mean(mean, samples) ≈ μ atol = 0.05
            @test mean(var, samples) ≈ σ² atol = 0.05

            # initial parameter
            init_x = [randn(1) for _ in 1:5]
            samples = sample(
                ESSModel(prior, ℓvec), ESS(), alg, 10, 5; progress=false, init_params=init_x
            )
            @test map(first, samples) == init_x
        end

        # initial parameter
        init_x = randn(1)
        samples = sample(
            ESSModel(prior, ℓvec), ESS(), 10; progress=false, init_params=init_x
        )
        @test first(samples) == init_x
    end
end

ext: jl


text:
maxlabel(code) = maximum(l isa LabelNode ? l.label : -1 for l in code)

isgoto(::GotoNode) = true
isgoto(x) = isexpr(x, :gotoifnot)
isjump(x) =
  isgoto(x) || isexpr(x, :return) ||
  (isexpr(x, :call) && x.args[1] == Return())

function blocks(code)
  l = maxlabel(code)
  bs = Dict{Int,Vector{Any}}()
  label = 1
  current = bs[1] = []
  for x in code
    if x isa LabelNode
      label = x.label
      isjump(current[end]) || push!(current, GotoNode(label))
      bs[label] = current = []
    elseif isexpr(x, :gotoifnot)
      label = (l += 1)
      push!(current, x, GotoNode(label))
      bs[label] = current = []
    else
      push!(current, x)
    end
  end
  return bs
end

function branches(code, ignore = ())
  a = Int[]
  if code[end] isa GotoNode
    push!(a, code[end].label)
    length(code) > 1 &&
      isexpr(code[end-1], :gotoifnot) && push!(a, code[end-1].args[2])
  end
  return setdiff(a, ignore)
end

function splitcondition(code)
  code[end] isa GotoNode || return (code, nothing)
  (length(code) > 1 && isexpr(code[end-1], :gotoifnot)) || return (code[1:end-1], nothing)
  code[1:end-2], code[end-1].args[1]
end

function accessible(blocks, label, ignore = ())
  nu = branches(blocks[label], ignore)
  seen = Int[]
  while !isempty(nu)
    nu[1] in seen && (shift!(nu); continue)
    push!(seen, shift!(nu))
    append!(nu, branches(blocks[seen[end]], ignore))
  end
  return seen
end

accessible(blocks, label, labels...) = union(accessible(blocks, label), accessible(blocks, labels...))

function loop_next(bs, label)
  branches(bs[label])[2]
end

# NOTE: this will break if arbitrary gotos are used.
function restructure(bs::Associative, entries = 1, loops = [], cond = nothing)
  if isempty(entries)
    []
  elseif length(entries) == 1
    l = entries[1]
    if !isempty(loops)
      l == loops[end] && return [:(continue)]
      l in loops && return [:(break)]
      l == loop_next(bs, loops[end]) && return [:(break)]
    end
    code, cond = splitcondition(bs[l])
    if l ∉ accessible(bs, l, loops)
      [code..., restructure(bs, branches(bs[l]), loops, cond)...]
    else
      [Expr(:while, true, Expr(:block, code..., restructure(bs, branches(bs[l]), [loops..., l], cond)...)),
       restructure(bs, loop_next(bs, l), loops)...]
    end
  elseif length(entries) == 2 && cond != nothing
    as = accessible(bs, entries...)
    [Expr(:if, cond,
      Expr(:block, restructure(bs, entries[1], loops)...),
      Expr(:block, restructure(bs, entries[2], loops)...))]
  else
    error("Rebuild error")
  end
end

restructure(code) = Expr(:block, restructure(blocks(code))...)


text:
export isethereumaddress

function isethereumaddress(str::AbstractString)::Bool
    ethReg = r"^(0x)[0-9a-f]{40}$"i

    return contains(str, ethReg)
end


text:
export Recall, Precision, MAP, AUC, ReciprocalRank, MPR, NDCG

"""
    Recall

Recall-at-``N`` (Recall@``N``) indicates coverage of truth samples as a result of top-``N`` recommendation. The value is computed by the following equation:
```math
\\mathrm{Recall@}N = \\frac{|\\mathcal{I}^+_u \\cap I_N(u)|}{|\\mathcal{I}^+_u|}.
```
Here, ``|\\mathcal{I}^+_u \\cap I_N(u)|`` is the number of *true positives*.

    measure(
        metric::Recall,
        truth::Array{T},
        pred::Array{T},
        k::Int
    )
"""
struct Recall <: RankingMetric end
function measure(metric::Recall, truth::Array{T}, pred::Array{T}, k::Int) where T
    count_true_positive(truth, pred[1:k]) / length(truth)
end

"""
    Precision

Precision-at-``N`` (Precision@``N``) evaluates correctness of a top-``N`` recommendation list ``I_N(u)`` according to the portion of true positives in the list as:
```math
\\mathrm{Precision@}N = \\frac{|\\mathcal{I}^+_u \\cap I_N(u)|}{|I_N(u)|}.
```

    measure(
        metric::Precision,
        truth::Array{T},
        pred::Array{T},
        k::Int
    )
"""
struct Precision <: RankingMetric end
function measure(metric::Precision, truth::Array{T}, pred::Array{T}, k::Int) where T
    count_true_positive(truth, pred[1:k]) / k
end

"""
    MAP

While the original Precision@``N`` provides a score for a fixed-length recommendation list ``I_N(u)``, mean average precision (MAP) computes an average of the scores over all recommendation sizes from 1 to ``|\\mathcal{I}|``. MAP is formulated with an indicator function for ``i_n``, the ``n``-th item of ``I(u)``, as:
```math
\\mathrm{MAP} = \\frac{1}{|\\mathcal{I}^+_u|} \\sum_{n = 1}^{|\\mathcal{I}|} \\mathrm{Precision@}n \\cdot  \\mathbb{1}_{\\mathcal{I}^+_u}(i_n).
```

It should be noticed that, MAP is not a simple mean of sum of Precision@``1``, Precision@``2``, ``\\dots``, Precision@``|\\mathcal{I}|``, and higher-ranked true positives lead better MAP.

    measure(
        metric::MAP,
        truth::Array{T},
        pred::Array{T}
    )
"""
struct MAP <: RankingMetric end
function measure(metric::MAP, truth::Array{T}, pred::Array{T}, k::Int=0) where T
    tp = accum = 0
    n_pred = length(pred)

    for n = 1:n_pred
        if findfirst(isequal(pred[n]), truth) != nothing
            tp += 1
            accum += tp / n
        end
    end

    accum / length(truth)
end

"""
    AUC

ROC curve and area under the ROC curve (AUC) are generally used in evaluation of the classification problems, but these concepts can also be interpreted in a context of ranking problem. Basically, the AUC metric for ranking considers all possible pairs of truth and other items which are respectively denoted by ``i^+ \\in \\mathcal{I}^+_u`` and ``i^- \\in \\mathcal{I}^-_u``, and it expects that the ``best'' recommender completely ranks ``i^+`` higher than ``i^-``, as follows:

![auc](./assets/images/auc.png)

AUC calculation keeps track the number of true positives at different rank in ``\\mathcal{I}``. At line 8, the function adds the number of true positives which were ranked higher than the current non-truth sample to the accumulated count of correct pairs. Ultimately, an AUC score is computed as portion of the correct ordered ``(i^+, i^-)`` pairs in the all possible combinations determined by ``|\\mathcal{I}^+_u| \\times |\\mathcal{I}^-_u|`` in set notation.

    measure(
        metric::AUC,
        truth::Array{T},
        pred::Array{T}
    )
"""
struct AUC <: RankingMetric end
function measure(metric::AUC, truth::Array{T}, pred::Array{T}, k::Int=0) where T
    tp = correct = 0
    for r in pred
        if findfirst(isequal(r), truth) != nothing
            # keep track number of tp placed before
            tp += 1
        else
            correct += tp
        end
    end
    # number of all possible tp-fp pairs
    pairs = tp * (length(pred) - tp)
    correct / pairs
end

"""
    ReciprocalRank

If we are only interested in the first true positive, reciprocal rank (RR) could be a reasonable choice to quantitatively assess the recommendation lists. For ``n_{\\mathrm{tp}} \\in \\left[ 1, |\\mathcal{I}| \\right]``, a position of the first true positive in ``I(u)``, RR simply returns its inverse:
```math
  \\mathrm{RR} = \\frac{1}{n_{\\mathrm{tp}}}.
```

RR can be zero if and only if ``\\mathcal{I}^+_u`` is empty.

    measure(
        metric::ReciprocalRank,
        truth::Array{T},
        pred::Array{T}
    )
"""
struct ReciprocalRank <: RankingMetric end
function measure(metric::ReciprocalRank, truth::Array{T}, pred::Array{T}, k::Int=0) where T
    n_pred = length(pred)
    for n = 1:n_pred
        if findfirst(isequal(pred[n]), truth) != nothing
            return 1 / n
        end
    end
    return 0
end

"""
    MPR

Mean percentile rank (MPR) is a ranking metric based on ``r_{i} \\in [0, 100]``, the percentile-ranking of an item ``i`` within the sorted list of all items for a user ``u``. It can be formulated as:
```math
\\mathrm{MPR} = \\frac{1}{|\\mathcal{I}^+_u|} \\sum_{i \\in \\mathcal{I}^+_u} r_{i}.
```
``r_{i} = 0\\%`` is the best value that means the truth item ``i`` is ranked at the highest position in a recommendation list. On the other hand, ``r_{i} = 100\\%`` is the worst case that the item ``i`` is at the lowest rank.

MPR internally considers not only top-``N`` recommended items also all of the non-recommended items, and it accumulates the percentile ranks for all true positives unlike MRR. So, the measure is suitable to estimate users' overall satisfaction for a recommender. Intuitively, ``\\mathrm{MPR} > 50\\%`` should be worse than random ranking from a users' point of view.

    measure(
        metric::MPR,
        truth::Array{T},
        pred::Array{T}
    )
"""
struct MPR <: RankingMetric end
function measure(metric::MPR, truth::Array{T}, pred::Array{T}, k::Int=0) where T
    accum = 0
    n_pred = length(pred)
    for t in truth
        r = (coalesce(findfirst(isequal(t), pred), 0) - 1) / n_pred
        accum += r
    end
    accum * 100 / length(truth)
end

"""
    NDCG

Like MPR, normalized discounted cumulative gain (NDCG) computes a score for ``I(u)`` which places emphasis on higher-ranked true positives. In addition to being a more well-formulated measure, the difference between NDCG and MPR is that NDCG allows us to specify an expected ranking within ``\\mathcal{I}^+_u``; that is, the metric can incorporate ``\\mathrm{rel}_n``, a relevance score which suggests how likely the ``n``-th sample is to be ranked at the top of a recommendation list, and it directly corresponds to an expected ranking of the truth samples.

    measure(
        metric::NDCG,
        truth::Array{T},
        pred::Array{T},
        k::Int
    )
"""
struct NDCG <: RankingMetric end
function measure(metric::NDCG, truth::Array{T}, pred::Array{T}, k::Int) where T
    dcg = idcg = 0
    for n = 1:k
        d = 1 / log2(n + 1)
        if findfirst(isequal(pred[n]), truth) != nothing
            dcg += d
        end
        idcg += d
    end
    dcg / idcg
end


text:
"""
    eachchunk(a)

Returns an iterator with `CartesianIndices` elements that mark the index range of each chunk within an array.
"""
function eachchunk end

struct GridChunks{N}
    parentsize::NTuple{N,Int}
    chunksize::NTuple{N,Int}
    chunkgridsize::NTuple{N,Int}
    offset::NTuple{N,Int}
end
GridChunks(a, chunksize; offset = (_->0).(size(a))) = GridChunks(Int.(size(a)), Int.(chunksize), Int.(getgridsize(size(a),chunksize,offset)),Int.(offset))
GridChunks(a::Tuple, chunksize; offset = (_->0).(a)) = GridChunks(Int.(a), Int.(chunksize), Int.(getgridsize(a,chunksize,offset)),Int.(offset))
function getgridsize(a,chunksize,offset)
  map(a,chunksize,offset) do s,cs,of
    fld1(s+of,cs)
  end
end
function Base.show(io::IO, g::GridChunks)
  print(io,"Regular ",join(g.chunksize,"x")," chunks over a ", join(g.parentsize,"x"), " array.")
end
Base.size(g::GridChunks) = g.chunkgridsize
Base.size(g::GridChunks, dim) = g.chunkgridsize[dim]
Base.IteratorSize(::Type{GridChunks{N}}) where N = Base.HasShape{N}()
Base.eltype(::Type{GridChunks{N}}) where N = CartesianIndices{N,NTuple{N,UnitRange{Int64}}}
Base.length(c::GridChunks) = prod(size(c))
@inline function _iterate(g,r)
    if r === nothing
        return nothing
    else
        ichunk, state = r
        outinds = map(ichunk.I, g.chunksize, g.parentsize,g.offset) do ic, cs, ps, of
            max((ic-1)*cs+1-of,1):min(ic*cs-of, ps)
        end |> CartesianIndices
        outinds, state
    end
end
function Base.iterate(g::GridChunks)
    r = iterate(CartesianIndices(g.chunkgridsize))
    _iterate(g,r)
end
function Base.iterate(g::GridChunks, state)
    r = iterate(CartesianIndices(g.chunkgridsize), state)
    _iterate(g,r)
end

#Define the approx default maximum chunk size (in MB)
"The target chunk size for processing for unchunked arrays in MB, defaults to 100MB"
const default_chunk_size = Ref(100)

"""
A fallback element size for arrays to determine a where elements have unknown
size like strings. Defaults to 100MB
"""
const fallback_element_size = Ref(100)

#Here we implement a fallback chunking for a DiskArray although this should normally
#be over-ridden by the package that implements the interface

function eachchunk(a::AbstractArray)
  cs = estimate_chunksize(a)
  GridChunks(a,cs)
end

struct Chunked end
struct Unchunked end
function haschunks end
haschunks(x) = Unchunked()

"""
    element_size(a::AbstractArray)

Returns the approximate size of an element of a in bytes. This falls back to calling `sizeof` on
the element type or to the value stored in `DiskArrays.fallback_element_size`. Methods can be added for
custom containers.
"""
function element_size(a::AbstractArray)
  if isbitstype(eltype(a))
    return sizeof(eltype(a))
  elseif isbitstype(Base.nonmissingtype(eltype(a)))
    return sizeof(Base.nonmissingtype(eltype(a)))
  else
    @warn "Can not determine size of element type. Using DiskArrays.fallback_element_size[] = $(fallback_element_size[]) bytes"
    return fallback_element_size[]
  end
end

estimate_chunksize(a::AbstractArray) = estimate_chunksize(size(a), element_size(a))
function estimate_chunksize(s, si)
  ii = searchsortedfirst(cumprod(collect(s)),default_chunk_size[]*1e6/si)
  ntuple(length(s)) do idim
    if idim<ii
      return s[idim]
    elseif idim>ii
      return 1
    else
      sbefore = idim == 1 ? 1 : prod(s[1:idim-1])
      return floor(Int,default_chunk_size[]*1e6/si/sbefore)
    end
  end
end


text:
""" Cálculo de volume da caixa
    Leia o comprimento a largura e a altura de uma caixa d'água retangular
    e calcule e informe o volume da caixa
"""

println("informe o comprimento")
volume =  parse(Float16, readline())

println("informe a largura")
volume *=  parse(Float16, readline())

println("informe a altura")
volume *=  parse(Float16, readline())

println("Volume = ", volume)

max_stars_repo_name: codinginbrazil/functional-programming

max_issues_repo_name: codinginbrazil/functional-programming

max_forks_repo_name: codinginbrazil/functional-programming


text:
function char_to_magmaInt(option::AbstractChar)
        if option == 'A'
                return MagmaAllVec
        elseif option == 'S'
                return MagmaSomeVec
        elseif option == 'O'
                return MagmaOverwriteVec
        elseif option == 'N'
                return MagmaNoVec
        end
end

for (fname, elty, relty) in    ((:sgesvd, :Float32, :Float32),
                                (:dgesvd, :Float64, :Float64),
                                (:cgesvd, :ComplexF32, :Float32),
                                (:zgesvd, :ComplexF64, :Float64))
        @eval begin
        # magma_int_t magma_cgesvd 	( 	magma_vec_t  	jobu,
        # 	magma_vec_t  	jobvt,
        # 	magma_int_t  	m,
        # 	magma_int_t  	n,
        # 	magmaFloatComplex *  	A,
        # 	magma_int_t  	lda,
        # 	float *  	s,
        # 	magmaFloatComplex *  	U,
        # 	magma_int_t  	ldu,
        # 	magmaFloatComplex *  	VT,
        # 	magma_int_t  	ldvt,
        # 	magmaFloatComplex *  	work,
        # 	magma_int_t  	lwork,
        # 	float *  	rwork,
        # 	magma_int_t *  	info
        # )
                function gesvd!(jobu::AbstractChar, jobvt::AbstractChar, A::AbstractMatrix{$elty})
                        m, n    = size(A)
                        minmn   = min(m, n)
                        lda     = max(1, stride(A, 2))

                        S       = similar(A, $relty, minmn)

                        U       = similar(A, $elty, jobu == 'A' ? (m, m) : (jobu == 'S' ? (m, minmn) : (m, 0)))
                        ldu     = max(1, stride(U, 2))

                        VT      = similar(A, $elty, jobvt == 'A' ? (n, n) : (jobvt == 'S' ? (minmn, n) : (n, 0)))
                        ldvt    = max(1, stride(VT, 2))

                        work    = Vector{$elty}(undef, 1)
                        cmplx   = eltype(A) <: Complex
                        if cmplx
                                rwork = Vector{$relty}(undef, 5minmn)
                        end
                        lwork   = -1
                        info    = Ref{Cint}()

                        jobu_magma      = char_to_magmaInt(jobu)
                        jobvt_magma     = char_to_magmaInt(jobvt)

                        for i in 1:2
                                if cmplx
                                        ccall((@magmafunc($fname), libmagma), Cint,
                                                (Cint, Cint,
                                                Cint, Cint,
                                                PtrOrCuPtr{$elty}, Cint,
                                                PtrOrCuPtr{$relty},
                                                PtrOrCuPtr{$elty}, Cint,
                                                PtrOrCuPtr{$elty}, Cint,
                                                PtrOrCuPtr{$elty}, Cint, PtrOrCuPtr{$relty},
                                                PtrOrCuPtr{Cint}),
                                                  jobu_magma, jobvt_magma,
                                                  m, n,
                                                  A, lda,
                                                  S,
                                                  U, ldu,
                                                  VT, ldvt,
                                                  work, lwork, rwork,
                                                  info)

                                else
                                        ccall((@magmafunc($fname), libmagma), Cint,
                                                  (Cint, Cint,
                                                  Cint, Cint,
                                                  Ptr{$elty}, Cint,
                                                  Ptr{$relty},
                                                  Ptr{$elty}, Cint,
                                                  Ptr{$elty}, Cint,
                                                  Ptr{$elty}, Cint,
                                                  Ptr{Cint}),
                                                    jobu_magma, jobvt_magma,
                                                    m, n,
                                                    A, lda,
                                                    S,
                                                    U, ldu,
                                                    VT, ldvt,
                                                    work, lwork,
                                                    info)
                                end
                                if i == 1
                                        lwork = ceil(Int, (work[1]))
                                        resize!(work, lwork)
                                end
                        end

                        if jobu == 'O'
                                return (A, S, VT)
                        elseif jobvt == 'O'
                                return (U, S, A)
                        else
                                return (U, S, VT)
                        end
                end
        end
end


text:
using Bridge
import Base.valtype

abstract type DiffusionDomain end

struct UnboundedDomain <: DiffusionDomain end

bound_satisfied(::UnboundedDomain, x) = true

struct LowerBoundedDomain{T,N} <: DiffusionDomain
    bounds::NTuple{N,T}
    coords::NTuple{N,Integer}

    function LowerBoundedDomain(bounds::NTuple{N,T},
                                coords::NTuple{N,Integer}) where {N,T}
        new{T,N}(bounds, coords)
    end

    function LowerBoundedDomain(bounds::NTuple{N,T}, coords) where {N,T}
        @assert length(coords) == N
        @assert all([typeof(c) <: Integer for c in coords])
        new{T,N}(bounds, Tuple(coords))
    end

    function LowerBoundedDomain(bounds::Vector{T}, coords) where T
        N = length(bounds)
        @assert length(coords) == N
        @assert all([typeof(c) <: Integer for c in coords])
        new{T,N}(Tuple(bounds), Tuple(coords))
    end
end

function bound_satisfied(d::LowerBoundedDomain{T,N}, x) where {T,N}
    for i in 1:N
        (x[d.coords[i]] < d.bounds[i]) && return false
    end
    true
end

struct UpperBoundedDomain{T,N} <: DiffusionDomain
    bounds::NTuple{N,T}
    coords::NTuple{N,Integer}

    function UpperBoundedDomain(bounds::NTuple{N,T},
                                coords::NTuple{N,Integer})  where {T,N}
        new{T,N}(bounds, coords)
    end

    function UpperBoundedDomain(bounds::NTuple{N,T}, coords) where {T,N}
        @assert length(coords) == N
        @assert all([typeof(c) <: Integer for c in coords])
        new{T,N}(bounds, Tuple(coords))
    end

    function UpperBoundedDomain(bounds::Vector{T}, coords) where T
        N = length(bounds)
        @assert length(coords) == N
        @assert all([typeof(c) <: Integer for c in coords])
        new{T,N}(Tuple(bounds), Tuple(coords))
    end
end

function bound_satisfied(d::UpperBoundedDomain{T,N}, x) where {T,N}
    for i in 1:N
        (x[d.coords[i]] > d.bounds[i]) && return false
    end
    true
end

struct BoundedDomain{T,N1,N2} <: DiffusionDomain
    lowBds::LowerBoundedDomain{T,N1}
    upBds::UpperBoundedDomain{T,N2}

    function BoundedDomain(lowBds::LowerBoundedDomain{T,N1},
                           upBds::UpperBoundedDomain{T,N2}) where {T,N1,N2}
        new{T,N1,N2}(lowBds, upBds)
    end

    function BoundedDomain(lowBds, lowBdsCoords, upBds, upBdsCoords)
        lowBdsObj = lowerBoundedDomain(lowBds, lowBdsCoords)
        upBdsObj = upperBoundedDomain(upBds, upBdsCoords)
        T,N1 = valtype(lowBdsObj)
        S,N2 = valtype(upBdsObj)
        @assert T == S
        new{T,N1,N2}(lowBdsObj, upBdsObj)
    end
end

bound_satisfied(d, x) = bound_satisfied(d.lowBds, x) && bound_satisfied(d.upBds, x)

valtype(::LowerBoundedDomain{T,N}) where {T,N} = T,N
valtype(::UpperBoundedDomain{T,N}) where {T,N} = T,N
valtype(::BoundedDomain{T,N1,N2}) where {T,N1,N2} = T,N1,N2
DomainSomehowBounded = Union{LowerBoundedDomain,UpperBoundedDomain,BoundedDomain}

domain(::Any) = UnboundedDomain() # by default no restrictions

"""
    check_domain_adherence(P::ContinuousTimeProcess, XX::SamplePath)

Check domain calling `domain`
"""
check_domain_adherence(P::ContinuousTimeProcess, XX::SamplePath) =
    check_domain_adherence(P, XX , domain(P.Target))

"""
    check_domain_adherence(P::Vector{ContinuousTimeProcess},
                         XX::Vector{SamplePath}, iRange)

Verify whether all paths in the range `iRange`, i.e. `XX[i].yy`, i in `iRange`
fall on the interior of the domain of diffusions `P[i]`, i in `iRange`
"""
function check_domain_adherence(P::Vector{S}, XX::Vector{T}, iRange
                              ) where {S<:ContinuousTimeProcess, T<:SamplePath}
    for i in iRange
        !check_domain_adherence(P[i], XX[i]) && return false
    end
    true
end

"""
    check_domain_adherence(P::ContinuousTimeProcess, XX::SamplePath,
                         d::UnboundedDomain)

For unrestricted domains there is nothing to check
"""
function check_domain_adherence(P::ContinuousTimeProcess, XX::SamplePath,
                              d::UnboundedDomain)
    true
end

"""
    check_domain_adherence(P::ContinuousTimeProcess, XX::SamplePath,
                         d::DiffusionDomain)

Verify whether path `XX.yy` falls on the interior of the domain of diffusion `P`
"""
function check_domain_adherence(P::ContinuousTimeProcess, XX::SamplePath,
                              d::DiffusionDomain)
    N = length(XX)
    for i in 1:N
        !bound_satisfied(d, XX.yy[i]) && false
    end
    true
end


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#
# this implements the event handling
#

# Return the next scheduled event.
_nextevent(c::Clock) = peek(c.sc.events)[1]

# Return the internal time (unitless) of next scheduled event.
_nextevtime(c::Clock) = peek(c.sc.events)[2]

# Execute or evaluate the next timed event on a clock c.
@inline function _event!(c::Clock)
    c.time = c.tev
    ev = dequeue!(c.sc.events)
    _evaluate(ev.ex)
    c.evcount += 1
    ev.Δt > 0.0 && event!(c, ev.ex, c.time + ev.Δt, cycle=ev.Δt)
end

# First execute all sampling expressions in a schedule, then evaluate all
# conditional events and if conditions are met, execute them.
function _tick!(c::Clock)
    c.time = c.tn
    foreach(x -> _evaluate(x.ex), c.sc.samples)  # exec sampling
    # then lookup conditional events
    ix = findfirst(x->all(_evaluate(x.cond)), c.sc.cevents)
    while ix !== nothing
        _evaluate(splice!(c.sc.cevents, ix).ex)
        if isempty(c.sc.cevents)
            if isempty(c.sc.samples) && isempty(c.ac) # no sampling and active clocks
                c.Δt = 0.0 # delete sample rate
            end
            break
        end
        ix = findfirst(x->all(_evaluate(x.cond)), c.sc.cevents)
    end
    c.scount +=1
end

# ------------------------------------------------------
# step forward to next tick or scheduled event. At a tick evaluate
# 1) all sampling functions or expressions,
# 2) all conditional events, then
# 3) if an event is encountered, trigger the event.
# ------------------------------------------------------
# If sampling rate Δt==0, c.tn is set to 0
# If no events are present, c.tev is set to c.end_time
# -------------------------------------------------------
@inline function _step!(c::Clock)
    c.state = Busy()
    if (c.tev ≤ c.time) && (length(c.sc.events) > 0)
        c.tev = _nextevtime(c)
    end

    if !isempty(c.sc.events)
        if c.Δt > 0.0
            if c.tn <= c.tev     # if t_next_tick  ≤ t_next_event
                _tick!(c)
                if c.tn == c.tev # an event is scheduled at the same time
                    _event!(c)
                end
                c.tn += c.Δt
            else
                _event!(c)
            end
        else
            _event!(c)
        end
    elseif c.Δt > 0.0
        _tick!(c)
        c.tn += c.Δt
    else
        (c.state == Busy()) && (c.state = Idle())
        return 99
        # error("_step!: nothing to evaluate")
    end
    !isempty(c.processes) && yield() # let processes run
    c.tev = !isempty(c.sc.events) ? _nextevtime(c) : c.end_time
    (c.state == Busy()) && (c.state = Idle())
    return 0
end

# ----------------------------------------------------
# execute all events in a clock cycle, then do the periodic actions
# ----------------------------------------------------
function _cycle!(clk::Clock, Δt::Float64, sync::Bool=false)
    clk.state = Busy()
    tcyc = clk.time + Δt
    clk.tev = length(clk.sc.events) ≥ 1 ? _nextevtime(clk) : clk.time
    while clk.time ≤ tcyc
        if (clk.tev ≤ tcyc) && length(clk.sc.events) ≥ 1
            _event!(clk)
        else
            clk.time = tcyc
            break
        end
        length(clk.processes) == 0 || yield() # let processes run
    end
    if !sync
        clk.tn = tcyc
        _tick!(clk)
    end
    (clk.state == Busy()) && (clk.state = Idle())
end


text:
mutable struct UIO{V<:Real}
    F::Matrix{V}
    T::Matrix{V}
    B::Matrix{V}
    H::Matrix{V}
    C::OutputMatrix{V}
    A1::Matrix{V}
    K1::Matrix{V}
    K::Matrix{V}
    z::Vector{V}

    function UIO(A::U, B::U, C::OutputMatrix{V}, E::U) where {V<:Real, U<:Matrix{V}}
        (ny, feasible, lrE) = checkUIOFeasibility(A,C,E);
        if feasible
            H = lrE*pinv(C*lrE)
            T = I - H*C
            A1 = T*A
            nz = size(A1, 1)
            # K2 = F*H;
            new{V}(
                zeros(V, nz, nz),
                T,
                B,
                H,
                C,
                A1,
                zeros(V, nz, ny),
                zeros(V, nz, ny),
                zeros(V, nz));
        else
            error("unfeasible UIO");
        end
    end
end

function UIO(A::U, B::U, C::OutputMatrix{V}, E::U, K1::U) where {V<:Real, U<:Matrix{V}}
    uio = UIO(A, B, C, E)
    setF(uio, K1)
    uio
end

function checkUIOFeasibility(A::Matrix{V},
        C::OutputMatrix{V},
        E::Matrix{V}) where {V<:Real}
    rE = rank(E)
    nx = size(A,1)
    if rE < nx
        SVD = svd(E);
        lrE = SVD.U[:,1:rE] * sqrt(Diagonal(SVD.S[1:rE]))
    else
        lrE = E
    end
    matC = C == I ? 1.0 * Matrix(I, nx, nx) : C
    size(matC, 1), rank(matC*lrE) == rE, lrE
end

function setF(o::UIO{T}, K1::AbstractMatrix{T}) where {T}
    o.F = o.A1 - K1*o.C
    o.K1 = K1;
    o.K = K1 + o.F*o.H
end

function updateState(o::UIO{T}, u::Vector{T}, y::Vector{T}) where {T}
    o.z = o.F*o.z + o.T*o.B*u + o.K*y
end

updateState(o::UIO{T}, u::T, y::Vector{T}) where {T} = updateState(o, [u], y)

outputMap(o::UIO{T}, y::Vector{T}) where {T} = o.z + o.H*y

function Base.getproperty(sys::UIO, s::Symbol)
    if s === :K2
        return sys.K - sys.K1
    else
        return getfield(sys, s)
    end
end

function Base.show(io::IO, sys::UIO)
    # Compose the name vectors
    #inputs = format_names(s.inputnames, "u", "?")
    #outputs = format_names(s.outputnames, "y", "?")
    println(io, "UIO:")
    # println(io, typeof(sys))
    println(io, "F = \n", _string_mat_with_headers(sys.F))
    println(io, "T = \n", _string_mat_with_headers(sys.T))
    println(io, "B = \n", _string_mat_with_headers(sys.B))
    println(io, "K = \n", _string_mat_with_headers(sys.K))
    println(io, "H = \n", _string_mat_with_headers(sys.H))
end


text:
@testset "tmpdir" begin

    #=
    TAKEN FROM JAVIS
    Leftover files from failed tests or errors can cause testing errors.
    Therefore, we remove any files in the `images` directory not pertinent to testing before executing a test.
    =#

    EXCLUDE_FILES = [".keep"]
    for f in readdir("images")
        !(f in EXCLUDE_FILES) && rm("images/$f")
    end
    @test length(readdir("images")) == length(EXCLUDE_FILES)
end

@testset "util" begin

    testvideo = Video(200, 200)
    Background(1:20, ground())
    Object(draw_node(center = Point(50, 50), radius = 20, action = :fill))
    render(testvideo, tempdirectory = "images", pathname = "")

    @test_reference "refs/test_draw_node.png" load("images/$(lpad(5, 10, "0")).png")

    for image in readdir("images", join = true)
        endswith(image, "png") && rm(image)
    end

    @test GB2Luxor(GeometryBasics.Point(10, 10)) == Point(10, 10)
end


text:
using Documenter, SatSolver

makedocs(;
    modules=[SatSolver],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/vutunganh/SatSolver.jl/blob/{commit}{path}#L{line}",
    sitename="SatSolver.jl",
    authors="Tung Anh Vu <vu.tunganh96@gmail.com>",
    assets=String[],
)

deploydocs(;
    repo="github.com/vutunganh/SatSolver.jl",
)


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

JLLWrappers.@generate_main_file_header("scopehal")
JLLWrappers.@generate_main_file("scopehal", UUID("c69b092c-cabd-51ee-983a-b3d58f6c0ce7"))
end  # module scopehal_jll


text:
module SSDPClient

export ssdpquery

using Sockets

"""
    ssdpquery(matchtext=""; timeoutsecs = 10 * 60)

Query local network for devices responding to query with the
Simple Service Discovery Protocol (SSDP).

The `matchtext` argument is an optional text to match, such as "Roku" or "Arduino".
The matching is case-insensitive.

The optional `timeoutsecs` is the number of seconds to monitor the network
before returning from the function. This defaults to 600 seconds or 10 minutes.

If you specify a capture in the match text, that capture will be returned,
otherwise the entry that had a match will be returned in its entirety.

NOTE: If you do not specify an argument all replies will be printed and
the function will loop until timeout (default 10 minutes) or until Control-C
interrupt. This is useful when searching a network.
"""
function ssdpquery(matchtxt=""; timeoutsecs = 24 * 60 * 60)
    MULTICAST = ip"239.255.255.250"
    SERVERPORT = 1900
    SEARCHHEADER = """M-SEARCH * HTTP/1.1\r\nHOST: 239.255.255.250:1900\r\nMAN: "ssdp:discover"\r\nST: ssdp:all\r\nMX: 10\r\n\r\n"""
    endtime = time() + timeoutsecs
    try
        while time() < endtime
            udpsock = UDPSocket()
            send(udpsock, MULTICAST,SERVERPORT, SEARCHHEADER)
            reply = join([Char(ch) for ch in recv(udpsock)], "")
            if matchtxt == ""
                println("Received reply:\r\n$reply")
            else
                regex = Regex(matchtxt, "is")
                if (mat = match(regex, reply)) != nothing
                    if length(mat.captures) > 0
                        return mat.captures
                    end
                    return reply
                end
            end
            close(udpsock)
            sleep(3)
        end
    catch y
        println("Caught error: $y")
        return ""
    end
end

end # module


text:
# -----------------------------------------------------------------------
# Visualization
# -----------------------------------------------------------------------
include("cylinder_structure.jl")
using Plots
gr()

"""
    circleShape(xc, yc, r)

## Description
Create a circunference
## Arguments
- xc and yc: center
- r: radius
"""
function circleShape(xc, yc, r)
    θ = LinRange(0,2*π, 500)
    @. xc + r*sin(θ), yc + r*cos(θ)
end

"""
    function plotCyl(pos::Array{OneCylType,1}; cylAlpha=0.5, B=(-50,20), L=(-10,100), r=1.2, dmfac=50.0, δc=(4,0), ptitle="")

## Description
Plot one cylinder (multiple position)
## Arguments
- pos: array of OneCylType
- cylAlpha: transparency paramenter
- B: tuple with x limits
- L: tuple with y limits
- r: cylinder radius
- dmfac: dimensionless factor
- δc: origin distance from the step
- ptitle: plot title
"""
function plotCyl(pos::Array{OneCylType,1}; cylAlpha=0.5, B=(-50,20), L=(-10,100), r=1.2,
    dmfac=50.0, δc=(4,0), ptitle="")

    n = length(pos)
    xb = [L[1],0.0,0.0] ./ dmfac
    yb = [0.0,0.0,B[1]] ./ dmfac
    r = r / dmfac
    δc = δc ./ dmfac
    L = L ./ dmfac
    B = B ./ dmfac

    if length(cylAlpha) == 1
        cylAlpha = repeat([cylAlpha], n)
    end

    p = plot(xb,yb, lw=2, xlims=L, ylims=B, c=:gray, legend=false, aspect_ratio=1)
    xlabel!("x/B")
    ylabel!("y/B")
    title!(ptitle)

    for i in 1:n
        plot!(circleShape(pos[i].x/dmfac + δc[1], pos[i].y/dmfac + δc[2], r), seriestype=[:shape], lw=0.1,
            c=:blue, linecolor=false, fillalpha=cylAlpha[i])
    end

    display(p)
end

"""
    function plotCyl(pos::Array{TwoCylType,1}; cylAlpha=[0.5], B=(-50,20), L=(-10,120), r=1.2, dmfac=50.0, δc=(4,0,10,0), ptitle="")

## Description
Plot two cylinder (multiple position)
## Arguments
- pos: array of TwoCylType
- cylAlpha: transparency paramenter
- B: tuple with x limits
- L: tuple with y limits
- r: cylinder radius
- dmfac: dimensionless factor
- δc: array with the origin distance from the step and the distance between cylinders at the origin
- ptitle: plot title
"""
function plotCyl(pos::Array{TwoCylType,1}; cylAlpha=[0.5], B=(-50,20), L=(-10,120), r=1.2,
    dmfac=50.0, δc=(4,0,10,0), ptitle="")

    n = length(pos)
    xb = [L[1],0.0,0.0] ./ dmfac
    yb = [0.0,0.0,B[1]] ./ dmfac
    r = r / dmfac
    δc = δc ./ dmfac
    B = B ./ dmfac
    L = L ./ dmfac

    p = plot(xb,yb, lw=2, xlims=L, ylims=B, c=:gray, legend=false, aspect_ratio=1)
    xlabel!("x/B")
    ylabel!("y/B")
    title!(ptitle)

    if length(cylAlpha) == 1
        cylAlpha = repeat([cylAlpha], n)
    end

    for i in 1:n
        Px = [pos[i].x/dmfac + δc[1], pos[i].x/dmfac + pos[i].dx/dmfac + δc[3]]
        Py = [pos[i].y/dmfac + δc[2], pos[i].y/dmfac + pos[i].dy/dmfac + δc[4]]
        plot!(circleShape(Px[1], Py[1], r), seriestype=[:shape], lw=0.1,
            c=:blue, linecolor=false, fillalpha=cylAlpha[i])
        plot!(circleShape(Px[2], Py[2], r), seriestype=[:shape], lw=0.1,
            c=:blue, linecolor=false, fillalpha=cylAlpha[i])
        plot!(Px, Py, c=:red)
    end

    display(p)
end

"""
    function plotCyl(pos::Array{ThreeCylType,1}; cylAlpha=[0.5], B=(-50,20), L=(-10,130), r=1.2, dmfac= 50.0, δc=(4,0,10,0,16,0), ptitle="")

## Description
Plot three cylinder (multiple position)
## Arguments
- pos: array of ThreeCylType
- cylAlpha: transparency paramenter
- B: tuple with x limits
- L: tuple with y limits
- r: cylinder radius
- dmfac: dimensionless factor
- δc: array with the origin distance from the step and the distance between cylinders at the origin
- ptitle: plot title
"""
function plotCyl(pos::Array{ThreeCylType,1}; cylAlpha=[0.5], B=(-50,20), L=(-10,130), r=1.2,
    dmfac= 50.0, δc=(4,0,10,0,16,0), ptitle="")

    n = length(pos)
    xb = [L[1],0.0,0.0] ./ dmfac
    yb = [0.0,0.0,B[1]] ./ dmfac
    r = r / dmfac
    δc = δc ./ dmfac
    B = B ./ dmfac
    L = L ./ dmfac

    p = plot(xb,yb, lw=2, xlims=L, ylims=B, c=:gray, legend=false, aspect_ratio=1)
    xlabel!("x/B")
    ylabel!("y/B")
    title!(ptitle)

    if length(cylAlpha) == 1
        cylAlpha = repeat([cylAlpha], n)
    end

    for i in 1:n
        Px = [pos[i].x/dmfac + δc[1], pos[i].x/dmfac + pos[i].dx1/dmfac + δc[3], pos[i].x/dmfac + pos[i].dx1/dmfac + pos[i].dx2/dmfac + δc[5]]
        Py = [pos[i].y/dmfac + δc[2], pos[i].y/dmfac + pos[i].dy1/dmfac + δc[4], pos[i].y/dmfac + pos[i].dy1/dmfac + pos[i].dy2/dmfac + δc[6]]
        plot!(circleShape(Px[1], Py[1], r), seriestype=[:shape], lw=0.1,
            c=:blue, linecolor=false, fillalpha=cylAlpha[i])
        plot!(circleShape(Px[2], Py[2], r), seriestype=[:shape], lw=0.1,
            c=:blue, linecolor=false, fillalpha=cylAlpha[i])
        plot!(circleShape(Px[3], Py[3], r), seriestype=[:shape], lw=0.1,
            c=:blue, linecolor=false, fillalpha=cylAlpha[i])
        plot!(Px, Py, c=:red)
    end

    display(p)
end
# Methods for one position
function plotCyl(pos::CylType; cylAlpha=0.5, B=(-50,20), L=(-10,130), r=1.2)
    plotCyl([pos], cylAlpha=cylAlpha, B=B, L=L, r=r)
end

"""
    plotGeneration(df::DataFrame, rmcol=missing)

## Description
Plot cylinders for dataframe data format
## Arguments
- df: dataframe with x, y, dx, dy parameters
- rmcol: columns to be removed
"""
function plotGeneration(df::DataFrame; rmcol=missing, dmfac=50.0, ptitle="")

    if !ismissing(rmcol)
        df = df[!,Not(rmcol)]
    end

    nc = ncol(df)
    if nc == 2
        el = convert(OneCylType, df)
    elseif nc == 4
        el = convert(TwoCylType, df)
    elseif nc == 6
        el = convert(ThreeCylType, df)
    else
        throw(DimensionMismatch)
    end

    plotCyl(el, dmfac=dmfac, ptitle=ptitle)
end


text:
import Base.truncate

################################################################
# Verify that support is not empty after applying constraints
################################################################
"""
    verify_nonempty_support(uv::AbstractUncertainScalarKDE,
        constraint::SamplingConstraint)

Check if the support of the distribution furnishing the uncertain value
is empty after applying the sampling constraint.
"""
verify_nonempty_support(uv::AbstractUncertainScalarKDE,
    constraint::SamplingConstraint)

function verify_nonempty_support(uv::AbstractUncertainScalarKDE,
        constraint::TruncateMaximum)
    if constraint.max < minimum(uv.range)
        minrange = minimum(uv.range)
        truncmax = constraint.max

        e = "constraint.max = $truncmax < minimum(uv.range) = $minrange"
        throw(ArgumentError(e))
    end
end

function verify_nonempty_support(uv::AbstractUncertainScalarKDE,
        constraint::TruncateMinimum)

    if constraint.min > maximum(uv.range)
        truncmin = constraint.min
        maxrange = maximum(uv.range)

        e = "maximum(uv.range) = $maxrange < constraint.min = $truncmin"
        throw(ArgumentError(e))
    end
end

function verify_nonempty_support(uv::AbstractUncertainScalarKDE,
        constraint::TruncateRange)
    truncmin = constraint.min
    truncmax = constraint.max
    minrange = minimum(uv.range)
    maxrange = maximum(uv.range)

    if maxrange < truncmin && minrange > truncmax
        e1 = "maximum(uv.range) = $maxrange < TruncateMinimum.min = $truncmin and "
        e2 = "TruncateMaximum.max = $truncmax < minimum(uv.range) = $minrange"
        throw(ArgumentError(string(e1, e2)))
    elseif maxrange < truncmin
        e = "maximum(uv.range) = $maxrange < TruncateMinimum.min = $truncmin"
        throw(ArgumentError(e))
    elseif minrange > truncmax
        e = "TruncateMaximum.max = $truncmax < minimum(uv.range) = $minrange"
        throw(ArgumentError(e))
    end
end

export verify_nonempty_support

################################################################
# Truncating uncertain values based on kernel density estimates
################################################################

"""
    truncate(uv::AbstractUncertainScalarKDE, constraint::SamplingConstraint)

Truncate an uncertain value `uv` furnished by a kernel density estimated
distribution using the supplied `constraint`.
"""
truncate(uv::AbstractUncertainScalarKDE, constraint::SamplingConstraint)

truncate(uv::AbstractUncertainScalarKDE, constraint::NoConstraint) = uv

"""
    truncate(uv::AbstractUncertainScalarKDE, constraint::TruncateLowerQuantile)

Truncate the kernel density estimate to `uv`s distribution using a
`TruncateLowerQuantile` sampling constraint.
"""
function truncate(uv::AbstractUncertainScalarKDE, constraint::TruncateLowerQuantile)

    idx_lower_quantile = getquantileindex(uv, constraint.lower_quantile)

    # Subset the values and weights (values of the pdf at those values)
    idx_min = idx_lower_quantile
    idx_max = length(uv.pdf)
    range_subset = uv.range[idx_min:idx_max]
    pdf_subset = uv.pdf[idx_min:idx_max]

    # Return truncated KDE and the indices used to subset. We have to normalise the pdf
    # here, so that we're still dealing with a probability distribution.
    range_subset, pdf_subset ./ sum(pdf_subset), idx_min, idx_max
end

"""
    truncate(uv::AbstractUncertainScalarKDE, constraint::TruncateUpperQuantile)

Truncate the kernel density estimate to `uv`s distribution using a
`TruncateUpperQuantile` sampling constraint.
"""
function truncate(uv::AbstractUncertainScalarKDE, constraint::TruncateUpperQuantile)

    idx_upper_quantile = getquantileindex(uv, constraint.upper_quantile)

    # Subset the values and weights (values of the pdf at those values)
    idx_min = 1
    idx_max = idx_upper_quantile
    range_subset = uv.range[idx_min:idx_max]
    pdf_subset = uv.pdf[idx_min:idx_max]

    # Return truncated KDE and the indices used to subset. We have to normalise the pdf
    # here, so that we're still dealing with a probability distribution.
    range_subset, pdf_subset ./ sum(pdf_subset), idx_min, idx_max
end

"""
    truncate(uv::AbstractUncertainScalarKDE, constraint::TruncateQuantiles)

Truncate the kernel density estimate to `uv`s distribution using a
`TruncateQuantiles` sampling constraint.
"""
function truncate(uv::AbstractUncertainScalarKDE, constraint::TruncateQuantiles)

    idx_lower_quantile = getquantileindex(uv, constraint.lower_quantile)
    idx_upper_quantile = getquantileindex(uv, constraint.upper_quantile)

    # Subset the values and weights (values of the pdf at those values)
    idx_min = idx_lower_quantile
    idx_max = idx_upper_quantile
    range_subset = uv.range[idx_min:idx_max]
    pdf_subset = uv.pdf[idx_min:idx_max]

    # Return truncated KDE and the indices used to subset. We have to normalise the pdf
    # here, so that we're still dealing with a probability distribution.
    range_subset, pdf_subset ./ sum(pdf_subset), idx_min, idx_max
end

"""
    truncate(uv::AbstractUncertainScalarKDE, constraint::TruncateMinimum)

Truncate the kernel density estimate to `uv`s distribution using a
`TruncateMinimum` sampling constraint.
"""
function truncate(uv::AbstractUncertainScalarKDE, constraint::TruncateMinimum;
        test_support = true)

    # Is the support empty after applying the constraint? If so, throw error.
    test_support ? verify_nonempty_support(uv, constraint) : nothing

    # Subset the values and weights (values of the pdf at those values)
    idx_min = findfirst(uv.range .>= constraint.min)
    idx_max = length(uv.pdf)
    range_subset = uv.range[idx_min:idx_max]
    pdf_subset = uv.pdf[idx_min:idx_max]

    # Return truncated KDE and the indices used to subset. We have to normalise the pdf
    # here, so that we're still dealing with a probability distribution.
    range_subset, pdf_subset ./ sum(pdf_subset), idx_min, idx_max
end

"""
    truncate(uv::AbstractUncertainScalarKDE, constraint::TruncateMaximum)

Truncate the kernel density estimate to `uv`s distribution using a
`TruncateMaximum` sampling constraint.
"""
function truncate(uv::AbstractUncertainScalarKDE, constraint::TruncateMaximum;
        test_support = true)

    # Is the support empty after applying the constraint? If so, throw error.
    test_support ? verify_nonempty_support(uv, constraint) : nothing

    # Subset the values and weights (values of the pdf at those values)
    idx_min = 1
    idx_max = findlast(uv.range .<= constraint.max)
    range_subset = uv.range[idx_min:idx_max]
    pdf_subset = uv.pdf[idx_min:idx_max]

    # Return truncated KDE and the indices used to subset. We have to normalise the pdf
    # here, so that we're still dealing with a probability distribution.
    range_subset, pdf_subset ./ sum(pdf_subset), idx_min, idx_max
end

"""
    truncate(uv::AbstractUncertainScalarKDE, constraint::TruncateRange)

Truncate the kernel density estimate to `uv`s distribution using a
`TruncateRange` sampling constraint.
"""
function truncate(uv::AbstractUncertainScalarKDE, constraint::TruncateRange;
        test_support = true)

    # Is the support empty after applying the constraint? If so, throw error.
    test_support ? verify_nonempty_support(uv, constraint) : nothing

    # Subset the values and weights (values of the pdf at those values)
    idx_min = findfirst(uv.range .>= constraint.min)
    idx_max = findlast(uv.range .<= constraint.max)
    range_subset = uv.range[idx_min:idx_max]
    pdf_subset = uv.pdf[idx_min:idx_max]

    # Return truncated KDE and the indices used to subset. We have to normalise the pdf
    # here, so that we're still dealing with a probability distribution.
    range_subset, pdf_subset ./ sum(pdf_subset), idx_min, idx_max
end

truncate(uv::AbstractUncertainScalarKDE, constraint::TruncateStd) =
    truncate(uv, fallback(uv, constraint))

export truncate


text:
function supportPoint(obj::Composition.Object3D, e::SVector{3,Float64})
    shapeKind = obj.shapeKind
    solid::Modia3D.Solid = obj.feature
    collisionSmoothingRadius = solid.collisionSmoothingRadius

    if shapeKind == Modia3D.SphereKind
        #sphere::Modia3D.Sphere = obj.shape
        return Modia3D.supportPoint_Sphere(obj.shape, obj.r_abs, obj.R_abs, e)
    elseif shapeKind == Modia3D.EllipsoidKind
        #ellipsoid::Modia3D.Ellipsoid = obj.shape
        return Modia3D.supportPoint_Ellipsoid(obj.shape, obj.r_abs, obj.R_abs, e)
    elseif shapeKind == Modia3D.BoxKind
        #box::Modia3D.Box = obj.shape
        return Modia3D.supportPoint_Box(obj.shape, obj.r_abs, obj.R_abs, e, collisionSmoothingRadius)
    elseif shapeKind == Modia3D.CylinderKind
        #cylinder::Modia3D.Cylinder = obj.shape
        return Modia3D.supportPoint_Cylinder(obj.shape, obj.r_abs, obj.R_abs, e, collisionSmoothingRadius)
    elseif shapeKind == Modia3D.ConeKind
        #cone::Modia3D.Cone = obj.shape
        return Modia3D.supportPoint_Cone(obj.shape, obj.r_abs, obj.R_abs, e, collisionSmoothingRadius)
    elseif shapeKind == Modia3D.CapsuleKind
        #capsule::Modia3D.Capsule = obj.shape
        return Modia3D.supportPoint_Capsule(obj.shape, obj.r_abs, obj.R_abs, e, collisionSmoothingRadius)
    elseif shapeKind == Modia3D.BeamKind
        #beam::Modia3D.Beam = obj.shape
        return Modia3D.supportPoint_Beam(obj.shape, obj.r_abs, obj.R_abs, e, collisionSmoothingRadius)
    elseif shapeKind == Modia3D.FileMeshKind
        #fileMesh::Modia3D.FileMesh = obj.shape
        return Modia3D.supportPoint_FileMesh(obj.shape, obj.r_abs, obj.R_abs, e)
    else
        error("not supported shape for support points")
    end
end

function boundingBox!(obj::Composition.Object3D, AABB::Basics.BoundingBox; tight::Bool=true, scaleFactor::Float64=0.01)
    shapeKind = obj.shapeKind
    solid::Modia3D.Solid = obj.feature
    collisionSmoothingRadius = solid.collisionSmoothingRadius

    if shapeKind == Modia3D.SphereKind
        sphere::Modia3D.Sphere = obj.shape
        @inbounds begin
            r = sphere.diameter/2
            if tight    # best fitting AABB which is possible
                AABB.x_min = obj.r_abs[1] - r
                AABB.x_max = obj.r_abs[1] + r
                AABB.y_min = obj.r_abs[2] - r
                AABB.y_max = obj.r_abs[2] + r
                AABB.z_min = obj.r_abs[3] - r
                AABB.z_max = obj.r_abs[3] + r
            else
                # scales each axis with scaleFactor 0.1 means 10 percent bigger than best fitting AABB and so on
                @assert(scaleFactor > 0.0)
                scale = r*scaleFactor/2
                AABB.x_min = obj.r_abs[1] - r - scale
                AABB.x_max = obj.r_abs[1] + r + scale
                AABB.y_min = obj.r_abs[2] - r - scale
                AABB.y_max = obj.r_abs[2] + r + scale
                AABB.z_min = obj.r_abs[3] - r - scale
                AABB.z_max = obj.r_abs[3] + r + scale
            end
            return AABB
        end
    elseif shapeKind == Modia3D.EllipsoidKind
        ellipsoid::Modia3D.Ellipsoid = obj.shape
        xmin = Modia3D.supportPoint_i_Ellipsoid(ellipsoid, obj.r_abs[1], SVector(obj.R_abs[:,1]), -1)
        xmax = Modia3D.supportPoint_i_Ellipsoid(ellipsoid, obj.r_abs[1], SVector(obj.R_abs[:,1]), +1)
        ymin = Modia3D.supportPoint_i_Ellipsoid(ellipsoid, obj.r_abs[2], SVector(obj.R_abs[:,2]), -1)
        ymax = Modia3D.supportPoint_i_Ellipsoid(ellipsoid, obj.r_abs[2], SVector(obj.R_abs[:,2]), +1)
        zmin = Modia3D.supportPoint_i_Ellipsoid(ellipsoid, obj.r_abs[3], SVector(obj.R_abs[:,3]), -1)
        zmax = Modia3D.supportPoint_i_Ellipsoid(ellipsoid, obj.r_abs[3], SVector(obj.R_abs[:,3]), +1)
    elseif shapeKind == Modia3D.BoxKind
        box::Modia3D.Box = obj.shape
        xmin = Modia3D.supportPoint_i_Box(box, obj.r_abs[1], SVector(obj.R_abs[:,1]), -1, collisionSmoothingRadius)
        xmax = Modia3D.supportPoint_i_Box(box, obj.r_abs[1], SVector(obj.R_abs[:,1]), +1, collisionSmoothingRadius)
        ymin = Modia3D.supportPoint_i_Box(box, obj.r_abs[2], SVector(obj.R_abs[:,2]), -1, collisionSmoothingRadius)
        ymax = Modia3D.supportPoint_i_Box(box, obj.r_abs[2], SVector(obj.R_abs[:,2]), +1, collisionSmoothingRadius)
        zmin = Modia3D.supportPoint_i_Box(box, obj.r_abs[3], SVector(obj.R_abs[:,3]), -1, collisionSmoothingRadius)
        zmax = Modia3D.supportPoint_i_Box(box, obj.r_abs[3], SVector(obj.R_abs[:,3]), +1, collisionSmoothingRadius)
    elseif shapeKind == Modia3D.CylinderKind
        cylinder::Modia3D.Cylinder = obj.shape
        xmin = Modia3D.supportPoint_i_Cylinder(cylinder, obj.r_abs[1], SVector(obj.R_abs[:,1]), -1, collisionSmoothingRadius)
        xmax = Modia3D.supportPoint_i_Cylinder(cylinder, obj.r_abs[1], SVector(obj.R_abs[:,1]), +1, collisionSmoothingRadius)
        ymin = Modia3D.supportPoint_i_Cylinder(cylinder, obj.r_abs[2], SVector(obj.R_abs[:,2]), -1, collisionSmoothingRadius)
        ymax = Modia3D.supportPoint_i_Cylinder(cylinder, obj.r_abs[2], SVector(obj.R_abs[:,2]), +1, collisionSmoothingRadius)
        zmin = Modia3D.supportPoint_i_Cylinder(cylinder, obj.r_abs[3], SVector(obj.R_abs[:,3]), -1, collisionSmoothingRadius)
        zmax = Modia3D.supportPoint_i_Cylinder(cylinder, obj.r_abs[3], SVector(obj.R_abs[:,3]), +1, collisionSmoothingRadius)
    elseif shapeKind == Modia3D.ConeKind
        cone::Modia3D.Cone = obj.shape
        xmin = Modia3D.supportPoint_i_Cone(cone, obj.r_abs[1], SVector(obj.R_abs[:,1]), -1, collisionSmoothingRadius)
        xmax = Modia3D.supportPoint_i_Cone(cone, obj.r_abs[1], SVector(obj.R_abs[:,1]), +1, collisionSmoothingRadius)
        ymin = Modia3D.supportPoint_i_Cone(cone, obj.r_abs[2], SVector(obj.R_abs[:,2]), -1, collisionSmoothingRadius)
        ymax = Modia3D.supportPoint_i_Cone(cone, obj.r_abs[2], SVector(obj.R_abs[:,2]), +1, collisionSmoothingRadius)
        zmin = Modia3D.supportPoint_i_Cone(cone, obj.r_abs[3], SVector(obj.R_abs[:,3]), -1, collisionSmoothingRadius)
        zmax = Modia3D.supportPoint_i_Cone(cone, obj.r_abs[3], SVector(obj.R_abs[:,3]), +1, collisionSmoothingRadius)
    elseif shapeKind == Modia3D.CapsuleKind
        capsule::Modia3D.Capsule = obj.shape
        xmin = Modia3D.supportPoint_i_Capsule(capsule, obj.r_abs[1], SVector(obj.R_abs[:,1]), -1, collisionSmoothingRadius)
        xmax = Modia3D.supportPoint_i_Capsule(capsule, obj.r_abs[1], SVector(obj.R_abs[:,1]), +1, collisionSmoothingRadius)
        ymin = Modia3D.supportPoint_i_Capsule(capsule, obj.r_abs[2], SVector(obj.R_abs[:,2]), -1, collisionSmoothingRadius)
        ymax = Modia3D.supportPoint_i_Capsule(capsule, obj.r_abs[2], SVector(obj.R_abs[:,2]), +1, collisionSmoothingRadius)
        zmin = Modia3D.supportPoint_i_Capsule(capsule, obj.r_abs[3], SVector(obj.R_abs[:,3]), -1, collisionSmoothingRadius)
        zmax = Modia3D.supportPoint_i_Capsule(capsule, obj.r_abs[3], SVector(obj.R_abs[:,3]), +1, collisionSmoothingRadius)
    elseif shapeKind == Modia3D.BeamKind
        beam::Modia3D.Beam = obj.shape
        xmin = Modia3D.supportPoint_i_Beam(beam, obj.r_abs[1], SVector(obj.R_abs[:,1]), -1, collisionSmoothingRadius)
        xmax = Modia3D.supportPoint_i_Beam(beam, obj.r_abs[1], SVector(obj.R_abs[:,1]), +1, collisionSmoothingRadius)
        ymin = Modia3D.supportPoint_i_Beam(beam, obj.r_abs[2], SVector(obj.R_abs[:,2]), -1, collisionSmoothingRadius)
        ymax = Modia3D.supportPoint_i_Beam(beam, obj.r_abs[2], SVector(obj.R_abs[:,2]), +1, collisionSmoothingRadius)
        zmin = Modia3D.supportPoint_i_Beam(beam, obj.r_abs[3], SVector(obj.R_abs[:,3]), -1, collisionSmoothingRadius)
        zmax = Modia3D.supportPoint_i_Beam(beam, obj.r_abs[3], SVector(obj.R_abs[:,3]), +1, collisionSmoothingRadius)
    elseif shapeKind == Modia3D.FileMeshKind
        fileMesh::Modia3D.FileMesh = obj.shape
        xmin = Modia3D.supportPoint_i_FileMesh(fileMesh, obj.r_abs[1], SVector(obj.R_abs[:,1]), -1)
        xmax = Modia3D.supportPoint_i_FileMesh(fileMesh, obj.r_abs[1], SVector(obj.R_abs[:,1]), +1)
        ymin = Modia3D.supportPoint_i_FileMesh(fileMesh, obj.r_abs[2], SVector(obj.R_abs[:,2]), -1)
        ymax = Modia3D.supportPoint_i_FileMesh(fileMesh, obj.r_abs[2], SVector(obj.R_abs[:,2]), +1)
        zmin = Modia3D.supportPoint_i_FileMesh(fileMesh, obj.r_abs[3], SVector(obj.R_abs[:,3]), -1)
        zmax = Modia3D.supportPoint_i_FileMesh(fileMesh, obj.r_abs[3], SVector(obj.R_abs[:,3]), +1)
    else
        error("boundingBox! not supported for shape")
    end

    if tight    # best fitting AABB which is possible
        AABB.x_min = xmin
        AABB.x_max = xmax
        AABB.y_min = ymin
        AABB.y_max = ymax
        AABB.z_min = zmin
        AABB.z_max = zmax
    else
        # scales each axis with scaleFactor 0.1 means 10 percent bigger than best fitting AABB and so on
        @assert(scaleFactor > 0.0)
        scale = max(abs(xmax - xmin)*scaleFactor/2,abs(ymax - ymin)*scaleFactor/2,abs(zmax - zmin)*scaleFactor/2)
        AABB.x_min = xmin - scale
        AABB.x_max = xmax + scale
        AABB.y_min = ymin - scale
        AABB.y_max = ymax + scale
        AABB.z_min = zmin - scale
        AABB.z_max = zmax + scale
    end
    return AABB
end


text:
using Documenter, Queryverse

makedocs(
	modules = [Queryverse],
	sitename = "Queryverse.jl",
	analytics="UA-132838790-1",
	pages = [
        "Introduction" => "index.md"
    ]
)

deploydocs(
    repo = "github.com/queryverse/Queryverse.jl.git"
)


text:
struct CountingMeasure{X} <: AbstractDiscreteMeasure{X}
end


text:
@testset "Feature transform/Anisotropic images" begin
    function ind2cart(F)
        s = CartesianIndices(axes(F))
        return map(i -> CartesianIndex(s[i]), F)
    end

    A, w = [false false; false true], (3, 1)
    F = feature_transform(A, w)
    @test F == ind2cart([4 4; 4 4])
    D = distance_transform(F, w)
    @test D ≈ [sqrt(10) 3.0; 1.0 0.0]

    A, w = [false false; false true], (1, 3)
    F = feature_transform(A, w)
    @test F == ind2cart([4 4; 4 4])
    D = distance_transform(F, w)
    @test D ≈ [sqrt(10) 1.0; 3.0 0.0]

    A, w = [true false; false true], (3, 1)
    F = feature_transform(A, w)
    @test F == ind2cart([1 1; 4 4])
    D = distance_transform(F, w)
    @test D == [0 1; 1 0]

    A, w = [true false; false true], (1, 3)
    F = feature_transform(A, w)
    @test F == ind2cart([1 4; 1 4])
    D = distance_transform(F, w)
    @test D == [0 1; 1 0]
end