text:
# SOME NOTES TO REVISIT
# ----------------------
#  * Let's try to get away without definint cutoff for the TBModel ?!?!

export AbstractTBModel,
      TBModel,
      isorthogonal,
      TBHamiltonian,
      ChemicalPotential,
      BZQuadratureRule

# =======================  Abstract Hamiltonian Types =========================

const ORTHOGONAL = true
const NONORTHOGONAL = false

"""
abstract Tight Binding Hamiltonian operator; the parameter
`ISORTH` is used to  allow dispatch w.r.t. whether the TB model is orthogonal
   or non-orthogonal
"""
abstract type TBHamiltonian{ISORTH} end

"""
An auxiliary hamiltonian that doesn't do anything but lets us
construct an empty TBModel.
"""
mutable struct NullHamiltonian <: TBHamiltonian{ORTHOGONAL}
end

isorthogonal{ISORTH}(::TBHamiltonian{ISORTH}) = ISORTH == ORTHOGONAL
isorth(H::TBHamiltonian) = isorthogonal(H)

"""
`ndofs(::TBHamiltonian, ::AbstractAtoms)`

returns total number of electron dofs
"""
@protofun ndofs(::TBHamiltonian, ::AbstractAtoms)

@protofun evaluate(::TBHamiltonian, ::AbstractAtoms, ::AbstractVector)

evaluate(H::TBHamiltonian, at::AbstractAtoms) = evaluate(H, at, JVecF([0.0,0.0,0.0]))

"""
`SKHamiltonian{ISORTH, NORB}`

abstract supertype for all Slater-Koster type TB hamiltonians. Type parameters
are

* `ISORTH` ∈ {`true`, `false`} : determine whether TB model is orthogonal or non-orthogonal
* `NORB` ∈ {1,4,9} : integer, number of orbitals.
"""
abstract type SKHamiltonian{ISORTH, NORB} <: TBHamiltonian{ISORTH} end

norbitals{ISORTH,NORB}(::SKHamiltonian{ISORTH, NORB}) = NORB

nbonds{ISORTH}(::SKHamiltonian{ISORTH, 1}) = 1
nbonds{ISORTH}(::SKHamiltonian{ISORTH, 4}) = 4
nbonds{ISORTH}(::SKHamiltonian{ISORTH, 9}) = 10

ndofs(H::SKHamiltonian, at::AbstractAtoms) = norbitals(H) * length(at)

# =======================  Smearing Functions =============================

"""
`ChemicalPotential`: abstract supertype for different chemical potentials.
`TightBinding.jl` implements:
# TODO: make list
"""
abstract type ChemicalPotential end

abstract type FiniteTPotential <: ChemicalPotential end
abstract type ZeroTPotential <: ChemicalPotential end

"""
auxiliary smearing function that doesn't do anything but lets us construct
and empty TBModel.
"""
mutable struct NullPotential <: ChemicalPotential end

# ======================= BZ Quadrature supertype =====================
# see `bzintegration.jl` for implementations

"""
`BZQuadratureRule`: abstract BZ quadrature supertype. Quadrature rules
can be applied either useing `w_and_pts` or through iterators, e.g.,
```
for (w, k) in tbm.bzquad
    ...
end
```
"""
abstract type BZQuadratureRule end

mutable struct NullBZQ <: BZQuadratureRule end

# ===================  Standard TightBinding Calculator =====================

"""
supertype for all TB model type calcualtors
"""
abstract type AbstractTBModel <: AbstractCalculator end

"""
`TBModel`: basic non-self consistent tight binding calculator.
"""
mutable struct TBModel{HT <: TBHamiltonian, ST <: ChemicalPotential} <: AbstractTBModel
   H::HT                 # hamiltonian
   Vrep::SitePotential   # additional MM potential (typically but not necessarily pair)
   potential::ST         # chemical potential / smearing function
   bzquad::BZQuadratureRule      # k-point sampling
   # -------------- internals ------------------
   hfd::Float64          # step used for finite-difference approximations
end

const TightBindingModel = TBModel

TBModel() = TBModel(NullHamiltonian(), ZeroSitePotential(),
                    NullPotential(), NullBZQ(), 0.0)

isorthogonal(tbm::TBModel) = isorthogonal(tbm.H)
isorth(tbm::TBModel) = isorthogonal(tbm)

get_eF(tbm::TBModel) = get_eF(tbm.potential)

norbitals(tbm::TBModel) = norbitals(tbm.H)

"""
`hamiltonian`: computes the hamiltonitan and overlap matrix for a tight
binding model.

#### Parameters:

* `tbm::TBModel`
* `atm::AbstractAtoms`
* `k = [0.;0.;0.]` : k-point at which the hamiltonian is evaluated

### Output: (H, M)

* `H` : hamiltonian in suitable format (typically CSC)
* `M` : overlap matrix in suitable format (typically CSC or I if orthogonal)
"""
hamiltonian(tbm::AbstractTBModel, at::AbstractAtoms, args...) =
   evaluate(tbm.H, at, args...)

# TODO: return full or sparse hamiltonian depending on sparsity: 0.06% seems a good heuristic!

# ==============================================================================


text:
## Code from aamini/FastConv.jl for convolutions of small size

using Base.Cartesian

# direct version (do not check if threshold is satisfied)
function fastconv(E::Array{T,N}, k::Array{T,N}) where {T,N}
    retsize = ntuple(n -> size(E, n) + size(k, n) - 1, Val{N}())
    ret = zeros(T, retsize)
    convn!(ret, E, k)
    return ret
end

# in place helper operation to speedup memory allocations
@generated function convn!(out::Array{T}, E::Array{T,N}, k::Array{T,N}) where {T,N}
    quote
        @inbounds begin
            @nloops $N x E begin
                @nloops $N i k begin
                    (@nref $N out d -> (x_d + i_d - 1)) += (@nref $N E x) * (@nref $N k i)
                end
            end
        end
        return out
    end
end

## Code from Julia 1.4   (https://github.com/JuliaLang/julia/blob/master/base/math.jl#L101 on 4/8/20)
## cf. https://github.com/JuliaLang/julia/pull/32753
## Slight modification when `x` is a matrix
## Need to keep as we use _one and _muladd to work with x a vector or matrix
## e.g. 1 => I
module EvalPoly
using LinearAlgebra
function evalpoly(x::S, p::Tuple) where {S}
    p == () && return zero(S)
    if @generated
        N = length(p.parameters)
        ex = :(p[end]*_one(S))
        for i in N-1:-1:1
            ex = :(_muladd($ex, x, p[$i]))
        end
        ex
    else
        _evalpoly(x, p)
    end
end

evalpoly(x, p::AbstractVector) = _evalpoly(x, p)

function _evalpoly(x::S, p) where {S}
    N = length(p)
    ex = p[end]*_one(x)
    for i in N-1:-1:1
        ex = _muladd(ex, x, p[i])
    end
    ex
end

function evalpoly(z::Complex, p::Tuple)
    if @generated
        N = length(p.parameters)
        a = :(p[end]*_one(z))
        b = :(p[end-1]*_one(z))
        as = []
        for i in N-2:-1:1
            ai = Symbol("a", i)
            push!(as, :($ai = $a))
            a = :(muladd(r, $ai, $b))
            b = :(p[$i] - s * $ai)
        end
        ai = :a0
        push!(as, :($ai = $a))
        C = Expr(:block,
                 :(x = real(z)),
                 :(y = imag(z)),
                 :(r = x + x),
                 :(s = muladd(x, x, y*y)),
                 as...,
                 :(muladd($ai, z, $b)))
    else
        _evalpoly(z, p)
    end
end
evalpoly(z::Complex, p::Tuple{<:Any}) = p[1]*_one(z)

evalpoly(z::Complex, p::AbstractVector) = _evalpoly(z, p)

function _evalpoly(z::Complex, p)
    length(p) == 1 && return p[1]*_one(z)
    N = length(p)
    a = p[end]*_one(z)
    b = p[end-1]*_one(z)

    x = real(z)
    y = imag(z)
    r = 2x
    s = muladd(x, x, y*y)
    for i in N-2:-1:1
        ai = a
        a = muladd(r, ai, b)
        b = p[i] - s * ai
    end
    ai = a
    muladd(ai, z, b)
end

## modify muladd, as needed
_muladd(a,b,c) = muladd(a,b,c)
_muladd(a, b::Vector, c) = a.*b .+ c
_muladd(a, b::Matrix, c) = (a*I)*b + c*I

# try to get y = P(c::T)(x::S) = P{T}(c)(x::S) to
# have y == one(T)*one(S)*x
_one(P::Type{<:Matrix}) = one(eltype(P))*I
_one(x::Matrix) = one(eltype(x))*I
_one(x) = one(x)
end

## get type of parametric composite type without type parameters
## this is needed when the underlying type changes, e.g. with integration
## where T=Int might become T=Float64
##
## trick from [ConstructionBase.jl](https://github.com/JuliaObjects/ConstructionBase.jl/blob/b5686b755bd3bee29b181b3cb18fe2effa0f10a2/src/ConstructionBase.jl#L25)
## as noted in https://discourse.julialang.org/t/get-new-type-with-different-parameter/37253/4
##
#@generated function constructorof(::Type{T}) where T
#    getfield(parentmodule(T), nameof(T))
#end

# https://discourse.julialang.org/t/how-do-a-i-get-a-type-stripped-of-parameters/73465/11
constructorof(::Type{T}) where T = Base.typename(T).wrapper

# Define our own minimal Interval type, inspired by Intervals.jl.
# We vendor it in to avoid adding the heavy Intervals.jl dependency and
# using IntervalSets leads to many needs for type piracy that may interfere
# with uses by its many dependent packages.
# likely this command will be needed to use outside of this package:
# import Polynomials: domain, Interval, Open, Closed, bounds_types

abstract type Bound end
struct Closed <: Bound end
struct Open <: Bound end
struct Unbounded <: Bound end

"""
    Interval{T, L <: Bound, R <: Bound}

Very bare bones Interval following `Intervals.jl` assuming `T<:Real`.
"""
struct Interval{T, L <: Bound, R <: Bound}
    first::T
    last::T
    function Interval{T,L,R}(f::T, l::T) where {T, L <: Bound, R <: Bound}
        f > l && throw(ArgumentError("first not less than last"))
        𝐋 = isinf(f) ? Unbounded : L
        𝐑 = isinf(l) ? Unbounded : R
        return new{T,𝐋,𝐑}(f, l)
    end
    function Interval{L,R}(f, l) where {L <: Bound, R <: Bound}
        𝒇, 𝒍 = promote(f,l)
        T = eltype(𝒇)
        new{T,L,R}(𝒇, 𝒍)
    end
    Interval(f, l) = Interval{Closed, Closed}(f, l)
end

bounds_types(x::Interval{T,L,R}) where {T,L,R} = (L, R)

Base.broadcastable(I::Interval) = Ref(I)

function Base.show(io::IO, I::Interval{T,L,R}) where {T,L,R}
    l,r = extrema(I)
    print(io, L == Closed ? "[" : "(")
    print(io, l, ", ", r)
    print(io, R == Closed ? "]" : ")")
end

Base.first(I::Interval) = I.first
Base.last(I::Interval) = I.last
Base.extrema(I::Interval) = (first(I), last(I))

function Base.in(x, I::Interval{T,L,R}) where {T, L, R}
    a, b = extrema(I)
    (L == Open ? a < x : a <= x) && (R == Open ? x < b : x <= b)
end

Base.isopen(I::Interval{T,L,R}) where {T,L,R} = (L != Closed && R != Closed)

ext: jl


text:
using StatsFuns

@testset "Testing the CIR specific function μπh, t and T" begin
    @test DualOptimalFiltering.logμπh(2, 1.2, 3.1, 4) ≈ -1.9336877133784256 atol=10.0^(-5)
    @test DualOptimalFiltering.logμπh(2, 1.2, 3.1, [4]) ≈ -1.9336877133784256 atol=10.0^(-5)
    @test DualOptimalFiltering.logμπh(2, 1.2, 3.1, [4, 6]) ≈ -4.261806748257223 atol=10.0^(-5)
    @test DualOptimalFiltering.logμπh([2], 1.2, 3.1, [4, 6])[1] ≈ -4.261806748257223 atol=10.0^(-5)
    tmp = DualOptimalFiltering.logμπh([2, 3], 1.2, 3.1, [4, 6])
    @test tmp[1] ≈ -4.261806748257223 atol = 10.0^(-7)
    @test tmp[2] ≈ -4.15718179717835 atol = 10.0^(-7)

    @test DualOptimalFiltering.t_CIR(5, 3) == 8
    @test DualOptimalFiltering.t_CIR([5], 3) == 8
    @test DualOptimalFiltering.t_CIR([5,6], 3) == 14
    tmp = DualOptimalFiltering.t_CIR(5, [3,2,1])
end

@testset "Testing the transition functions" begin

    @test DualOptimalFiltering.Λ_prime_1D([3,2,1]) == 0:3

    tmp = DualOptimalFiltering.next_Λ_from_Λ_prime([3,2,1], 6, DualOptimalFiltering.t_CIR)
    for t in 1:3
        @test tmp[t] == [9, 8, 7][t]
    end
    tmp = DualOptimalFiltering.next_Λ_from_Λ_prime([3,2,1], [6], DualOptimalFiltering.t_CIR)
    for t in 1:3
        @test tmp[t] == [9, 8, 7][t]
    end
    tmp = DualOptimalFiltering.next_Λ_from_Λ_prime([3,2,1], [5,6], DualOptimalFiltering.t_CIR)
    for t in 1:3
        @test tmp[t] == [14, 13, 12][t]
    end

    @test DualOptimalFiltering.θ_prime_from_θ_CIR(1.4, 0.3, 1.2, 2.3) ≈ 0.38310545974007354 atol = 10.0^(-7)

    @test DualOptimalFiltering.θ_from_θ_prime(5, 0.3, DualOptimalFiltering.T_CIR)  == 1.3
    @test DualOptimalFiltering.θ_from_θ_prime([5], 0.3, DualOptimalFiltering.T_CIR)  == 1.3
    @test DualOptimalFiltering.θ_from_θ_prime([5,6], 0.3, DualOptimalFiltering.T_CIR)  == 2.3

    tmp = DualOptimalFiltering.next_logwms_from_log_wms_prime1D(log.([0.5,0.25,0.25]), [1,4,2], 6, 1.3, 0.7)
    for t in 1:3
        @test tmp[t] ≈ [-1.74651, -0.506833, -1.49961][t] atol = 5*10.0^(-5)
    end
    tmp = DualOptimalFiltering.next_logwms_from_log_wms_prime1D(log.([0.5,0.25,0.25]), [1,4,2], [6], 1.3, 0.7)
    for t in 1:3
        @test tmp[t] ≈ [-1.74651, -0.506833, -1.49961][t] atol = 5*10.0^(-5)
    end
    tmp = DualOptimalFiltering.next_logwms_from_log_wms_prime1D(log.([0.5,0.25,0.25]), [1,4,2], [5,6], 1.3, 0.7)
    ref = [-1.956560511980828, -0.42974982819492613, -1.5702921541433446]
    for t in 1:3
        @test tmp[t] ≈ ref[t] atol = 5*10.0^(-5)
    end
    tmp = DualOptimalFiltering.next_logwms_from_log_wms_prime1D(log.([0.5,0.25,0.25]), [1,4,2], [6,5], 1.3, 0.7)
    for t in 1:3
        @test tmp[t] ≈ ref[t] atol = 5*10.0^(-5)
    end
    tmp = DualOptimalFiltering.next_wms_prime_from_wms([0.5,0.25,0.25], [1,4,2], 0.2, 1.3, 1.1, 0.6)
    for t in 1:5
        @test tmp[t] ≈ [0.104581, 0.487464, 0.199507, 0.101034, 0.107413][t] atol = 5*10.0^(-5)
    end
    tmp = DualOptimalFiltering.next_log_wms_prime_from_log_wms1D(log.([0.5,0.25,0.25]), [1,4,2], 0.2, 1.3, 1.1, 0.6)
    for t in 1:5
        @test tmp[t] ≈ [-2.25779, -0.718539, -1.6119, -2.2923, -2.23107][t] atol = 5*10.0^(-5)
    end
    tmp = log.(DualOptimalFiltering.next_wms_prime_from_wms([0.5,0.25,0.25], [1,4,2], 0.2, 1.3, 1.1, 0.6))
    for t in 1:5
        @test tmp[t] ≈ [-2.25779, -0.718539, -1.6119, -2.2923, -2.23107][t] atol = 5*10.0^(-5)
    end

end;

@testset "Testing the likelihood functions" begin

    Random.seed!(1)

    δ = 3.
    γ = 2.5
    σ = 4.
    Nobs = 3
    dt = 0.011
    Nsteps = 10
    λ = 1.

    α = δ/2
    β = γ/σ^2

    function rec_rcCIR(Dts, x, δ, γ, σ)
        x_new = DualOptimalFiltering.rCIR(1, Dts[1], x[end], δ, γ, σ)
        if length(Dts) == 1
            return Float64[x; x_new]
        else
            return Float64[x; rec_rcCIR(Dts[2:end], x_new, δ, γ, σ)]
        end
    end

    function generate_CIR_trajectory2(times, x0, δ, γ, σ)
        θ1 = δ*σ^2
        θ2 = 2*γ
        θ3 = 2*σ
        Dts = diff(times)
        return rec_rcCIR(Dts, [x0], δ, γ, σ)
    end

    time_grid = [k*dt for k in 0:(Nsteps-1)]
    X = generate_CIR_trajectory2(time_grid, 3, δ*1.2, γ/1.2, σ*0.7)
    Y = map(λ -> rand(Poisson(λ), Nobs), X);
    data = zip(time_grid, Y) |> Dict;

    @test_nowarn DualOptimalFiltering.log_likelihood_CIR(δ, γ, σ, λ, data)

end;

@testset "Testing the approximate likelihood functions" begin

    Random.seed!(1)

    δ = 3.
    γ = 2.5
    σ = 4.
    Nobs = 2
    dt = 0.011
    Nsteps = 10
    λ = 1.

    α = δ/2
    β = γ/σ^2

    function rec_rcCIR(Dts, x, δ, γ, σ)
        x_new = DualOptimalFiltering.rCIR(1, Dts[1], x[end], δ, γ, σ)
        if length(Dts) == 1
            return Float64[x; x_new]
        else
            return Float64[x; rec_rcCIR(Dts[2:end], x_new, δ, γ, σ)]
        end
    end

    function generate_CIR_trajectory2(times, x0, δ, γ, σ)
        θ1 = δ*σ^2
        θ2 = 2*γ
        θ3 = 2*σ
        Dts = diff(times)
        return rec_rcCIR(Dts, [x0], δ, γ, σ)
    end

    time_grid = [k*dt for k in 0:(Nsteps-1)]
    X = generate_CIR_trajectory2(time_grid, 3, δ*1.2, γ/1.2, σ*0.7)
    Y = map(λ -> rand(Poisson(λ), Nobs), X);
    data = zip(time_grid, Y) |> Dict;

    @test_nowarn tmp = DualOptimalFiltering.log_likelihood_CIR_fixed_fraction(δ, γ, σ, λ, data, 0.9; silence = false)
    @test_nowarn tmp = DualOptimalFiltering.log_likelihood_CIR_keep_above_threshold(δ, γ, σ, λ, data, 0.001; silence = false)
    @test_nowarn tmp = DualOptimalFiltering.log_likelihood_CIR_keep_fixed_number(δ, γ, σ, λ, data, 10; silence = false)

end;

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
struct JobStatistics{A}
    elapsed_time::A
    bytes_allocated::A
end
JobStatistics() = JobStatistics( Series( Mean(), Variance() ), Series( Mean(), Variance() ))

struct JobStatisticsSample
    elapsed_time::Union{Missing, Float64}
    bytes_allocated::Union{Missing, Float64}
end
JobStatisticsSample() = JobStatisticsSample(missing, missing)

@enum WORKER_STATE begin
    available   = 0
    ready       = 1
    launched    = 2
    hasdata     = 3
    failed      = 4
end

mutable struct WorkerCommunication
    task_stats::MrPhelps.JobStatisticsSample
    last_task::Int
    state::WORKER_STATE
end

#All thunk code was copied with only superficial changes from ChainRulesCore.jl: Permalink below.
#https://github.com/JuliaDiff/ChainRulesCore.jl/blob/47f5354191773d73a5dc372cd049b01556f6145f/src/differentials/thunks.jl#L73

struct Thunk{F}
    f::F
end
(x::Thunk)() = x.f()

"""
    @thunk expression

Define a [`Thunk`](@ref) wrapping the `expression`, to lazily defer its evaluation.
https://github.com/JuliaDiff/ChainRulesCore.jl/blob/47f5354191773d73a5dc372cd049b01556f6145f/src/differentials/thunks.jl#L73
"""
macro thunk(expression)
    func = Expr(:->, Expr(:tuple), Expr(:block, __source__, expression))
    return :(Thunk($(esc(func))))
end

"""
    dispatch_task( fn::Union{Thunk, Function},
                        remote_hook::RemoteChannel{ Channel{ Any } },
                        local_hook::RemoteChannel{ Channel{ WorkerCommunication } },
                        task_ID::Int,
                        src = nothing )

Executes a worker task on a remote machine. Statistics are sent back to the local scheduler,
data is stored locally.

"""
function dispatch_task( fn::Union{Thunk, Function},
                        remote_hook::RemoteChannel{ Channel{ Any } },
                        local_hook::RemoteChannel{ Channel{ WorkerCommunication } },
                        task_ID::Int,
                        src = nothing )
    try
        stats = Base.gc_num()
        elapsedtime = Base.time_ns()
        #if the channel is empty call the base function to add stuff too it
        #Note State_channel is a global that gets created in execute_mission() on every worker!!!!
        if !isnothing( fn )
            if !isready( remote_hook ) #channel is empty! should be our first iteration
                if isnothing(src)
                    put!( remote_hook, fn()() )
                else
                    put!( remote_hook, fn()( src ) )
                end
            else
                #if the channel is full, put the previous result into the next function...
                curval = fn()( take!( remote_hook ) )
                put!( remote_hook, curval)
                println("new thing: $curval")
            end
        else
            put!(remote_hook, src )#is actually a string or metadata?
        end
        elapsedtime = Base.time_ns() - elapsedtime
        diff = Base.GC_Diff( Base.gc_num(), stats )
        jobstats = JobStatisticsSample( elapsedtime * 1e-9, diff.allocd * 1e-6 )
        put!( local_hook, WorkerCommunication( jobstats, task_ID, ready ) )
    catch #uh oh
        @info("Worker failed on task node: $task_ID .")
        put!( local_hook, WorkerCommunication( JobStatisticsSample(), task_ID, failed ) )
    end
end

"""
    reset_worker(  remote_hook::RemoteChannel{ Channel{ Any } },
                    local_hook::RemoteChannel{ Channel{ WorkerCommunication } })

Resets the remote and local hooks for a worker which has completed its tasking.

"""
function reset_worker(  remote_hook::RemoteChannel{ Channel{ Any } },
                        local_hook::RemoteChannel{ Channel{ WorkerCommunication } } )
    try
        #think the memory is already reset but the user should have access to this
        #isready(remote_hook) && put!( remote_hook, nothing ) #Reset memory
        put!( local_hook, WorkerCommunication( JobStatisticsSample(), 0, available ) )
    catch #uh oh
        @info("Failed to reset worker...")
        put!( local_hook, WorkerCommunication( JobStatisticsSample(), 0, failed ) )
    end
end

max_stars_repo_name: caseykneale/MrPhelps

max_issues_repo_name: caseykneale/MrPhelps

max_forks_repo_name: caseykneale/MrPhelps


text:
module Segmentations

using ..Types
using ..Domains
using Colors
using FixedPointNumbers

export seg2aff, singleton2boundary!, relabel_seg
export add_seg_boundary!, seg2rgb, seg_overlay_img, seg_overlay_img!
export seg2sgm, seg2segMST, segid1N!, segid1N!_V3
export downsample

const NORMED8_HALF  = Normed{UInt8,8}(0.5)
const NORMED8_ONE   = Normed{UInt8,8}(1.0)

"""
construct affinity map from segmentation
"""
function seg2aff(seg::Array{T,3}) where T
  aff = zeros(Float32, (size(seg)..., 3))
  aff[2:end, :,:,1] = (seg[2:end, :,:] .== seg[1:end-1, :,:])
  aff[:, 2:end,:,2] = (seg[:, 2:end,:] .== seg[:, 1:end-1,:])
  aff[:,:, 2:end,3] = (seg[:,:, 2:end] .== seg[:,:, 1:end-1])
  aff
end

"""
label all the singletones as boundary
use original segmentation as mask to lable boundary
some singletons could be created by cropping margins,
this protected the marginal singletons.
"""
function singleton2boundary!( seg::Array{T,3}, ref::Array{T,3} ) where T
    Threads.@threads for i in eachindex(seg)
        seg[i] = ref[i]>0x00000000 ? seg[i] : 0x00000000
    end
end

function singleton2boundary!( seg::Array{T,3} ) where T

    # a flag array indicating whether it is segment
    flg = falses(size(seg))
    # size
    X,Y,Z = size(seg)

    # traverse the segmentation
    for z in 1:Z
        for y in 1:Y
            for x in 1:X
                if flg[x,y,z]
                    continue
                end
                if x>1 && seg[x,y,z]==seg[x-1,y,z]
                    continue
                end
                if x<X && seg[x,y,z]==seg[x+1,y,z]
                    flg[x+1,y,z] = true
                    continue
                end
                if y>1 && seg[x,y,z]==seg[x,y-1,z]
                    continue
                end
                if y<Y && seg[x,y,z]==seg[x,y+1,z]
                    flg[x,y+1,z] = true
                    continue
                end
                if z>1 && seg[x,y,z]==seg[x,y,z-1]
                    continue
                end
                if z<Z && seg[x,y,z]==seg[x,y,z+1]
                    flg[x,y,z+1] = true
                    continue
                end
                # it is a singletone
                seg[x,y,z] = T(0)
            end
        end
    end
end

# relabel the segment according to connectivity
# where N is the total number of segments
# Note that this is different from relabel1N in segerror package, which relabeles in 2D and labeled the segment ID to 1-N, where N is the total number of segments.
function relabel_seg( seg::Array{T,3} ) where T
    N = length(seg)
    X,Y,Z = size(seg)
    X = T(X);   Y = T(Y);   Z = T(Z);

    # initialize the disjoint sets
    djs = Tdjsets(N)

    # x affinity
    for x in 0x00000002:X
        for y in 0x00000001:Y
            for z in 0x00000001:Z
                if seg[x,y,z]>0 && seg[x,y,z]==seg[x-1,y,z]
                    # should union these two sets
                    vid1 = x            + (y-0x00000001)*X + (z-0x00000001)*X*Y
                    vid2 = x-0x00000001 + (y-0x00000001)*X + (z-0x00000001)*X*Y
                    # find tree root
                    r1 = find!(djs, vid1)
                    r2 = find!(djs, vid2)
                    # union two sets
                    union!(djs, r1, r2)
                end
            end
        end
    end

    # y affinity
    for x in 0x00000001:T(X)
        for y in 0x00000002:T(Y)
            for z in 0x00000001:T(Z)
                if seg[x,y,z]>0x00000000 && seg[x,y,z]==seg[x,y-0x00000001,z]
                    vid1 = x + (y-0x00000001)*X + (z-0x00000001)*X*Y
                    vid2 = x + (y-0x00000002)*X + (z-0x00000001)*X*Y
                    r1 = find!(djs, vid1)
                    r2 = find!(djs, vid2)
                    union!(djs, r1, r2)
                end
            end
        end
    end

    # z affinity
    for x in 0x00000001:T(X)
        for y in 0x00000001:T(Y)
            for z in 0x00000002:T(Z)
                if seg[x,y,z]>0x00000000 && seg[x,y,z] == seg[x,y,z-0x00000001]
                    vid1 = x + (y-0x00000001)*X + (z-0x00000001)*X*Y
                    vid2 = x + (y-0x00000001)*X + (z-0x00000002)*X*Y
                    r1 = find!(djs, vid1)
                    r2 = find!(djs, vid2)
                    union!(djs, r1, r2)
                end
            end
        end
    end

    # get current segmentation
    setallroot!( djs )
    ret = djs.sets
    ret = reshape(ret, size(seg))
    # mark all the singletons to 0 as boundary
    # singleton2boundary!(ret, seg) # this will remove the singletons created by cropping
    # should use original segmentation as a mask to remove boundary regions
    Threads.@threads for i in eachindex(seg)
        if seg[i] == 0x00000000
            ret[i] = 0x00000000
        end
    end
    return ret
end

# reassign segment ID as 1-N
function segid1N!( seg::Array{T,3} ) where T
    # dictionary of ids
    did = Dict{T, T}(T(0)=>T(0))
    sizehint!(did, div(length(seg),16))

    # number of segments
    N = 0x00000000
    v = 0x00000000
    # assign the map to a new segment
    for z in 1:size(seg, 3)
        for y in 1:size(seg, 2)
            for x in 1:size(seg, 1)
                v = seg[x,y,z]
                if !haskey(did, v)
                    # a new segment ID
                    N += 0x00000001
                    did[v] = N
                    seg[x,y,z] = N
                else
                    seg[x,y,z] = did[ v ]
                end
            end
        end
    end
    return N
end

"""
downsampling segmentation by simply sampling voxels
"""
function downsample( seg::Union{Array{UInt32,3},Array{UInt64,3}};
                    scale::Vector{Int} = [2,2,1])
    start = map(x->cld(x,2), scale)
    seg[start[1]:scale[1]:end, start[2]:scale[2]:end, start[3]:scale[3]:end]
end

#==
"""
multiple threads version runs without stop, used all the cpu!
"""
function segid1N!_V3{T}( seg::Array{T,3} )
    # dictionary of ids
    did = Dict{T, T}(T(0)=>T(0))
    sizehint!(did, div(length(seg),32))

    # number of segments
    N = Atomic{T}(0)
    v = T(0)

    # lock of critical section
    critical = SpinLock()
    # assign the map to a new segment
    gc_enable(false)
    for z in 1:size(seg, 3)
        for y in 1:size(seg, 2)
            @threads for x in 1:size(seg, 1)
                v = seg[x,y,z]
                if !haskey(did, v)
                    # a new segment ID
                    lock(critical)
                    # atomic_add!(N, 1)
                    N += 1
                    did[v] = N
                    seg[x,y,z] = N
                    unlock(critical)
                else
                    seg[x,y,z] = did[ v ]
                end
            end
        end
    end
    gc_enable(true)
    return N
end
==#

# add boundary between contacting segments
function add_lbl_boundary!(seg::Array{T,3}, conn=8) where T
    # neighborhood definition
    @assert conn==8 || conn==4
    sx,sy,sz = size(seg)
    for z = 1:sz
        for y = 1:sy
            for x = 1:sx
                if seg[x,y,z]==T(0)
                    # ignore the existing boundary
                    continue
                end
                # flag of central pixel
                cf = false
                # x direction
                if x<sx && seg[x+1,y,z]>T(0) && seg[x,y,z]!=seg[x+1,y,z]
                    cf = true
                    seg[x+1,y,z] = T(0)
                end
                # y direction
                if y<sy && seg[x,y+1,z]>T(0) && seg[x,y,z]!=seg[x,y+1,z]
                    cf = true
                    seg[x,y+1,z] = T(0)
                end
                if x>1 && seg[x-1,y,z]>T(0) && seg[x,y,z]!=seg[x-1,y,z]
                    cf = true
                    seg[x-1,y,z] = T(0)
                end
                if y>1 && seg[x,y-1,z]>T(0) && seg[x,y,z]!=seg[x,y-1,z]
                    cf = true
                    seg[x,y-1,z] = T(0)
                end

                if conn==8
                    if x<sx && y<sy && seg[x+1,y+1,z]>T(0) && seg[x,y,z]!=seg[x+1,y+1,z]
                        cf = true
                        seg[x+1,y+1,z] = T(0)
                    end

                    if x>1 && y<sy && seg[x-1,y+1,z]>T(0) && seg[x,y,z]!=seg[x-1,y+1,z]
                        cf = true
                        seg[x-1,y+1,z] = T(0)
                    end
                    if x<sx && y>1 && seg[x+1,y-1,z]>T(0) && seg[x,y,z]!=seg[x+1,y-1,z]
                        cf = true
                        seg[x+1,y-1,z] = T(0)
                    end
                    if x>1 && y>1 && seg[x-1,y-1,z]>T(0) && seg[x,y,z]!=seg[x-1,y-1,z]
                        cf = true
                        seg[x-1,y-1,z] = T(0)
                    end
                end
                if cf
                    print("$x,$y, ")
                    seg[x,y,z] = T(0)
                end
            end
        end
    end
end

"""
transform segmentation to domains
Inputs:
seg: a segmentation or label of image volume

Outputs:
dms: domains for fast union-find algorithm defined in "domains.jl"
"""
function segmentation2domains(seg::Array{T,3}; is_merge = true) where T
    # initialize a domain as singletons
    @assert ndims(seg)==2 || ndims(seg)==3
    dms = Tdomains( length(seg) )

    # if do not merge, return directly
    if !is_merge
        return dms
    end

    # volume size
    sx,sy,sz = size(seg)
    sx = T(sx); sy = T(sy); sz = T(sz)

    # union all the voxel with same segment ID
    for z in 0x00000001:sz
        for y in 0x00000001:sy
            for x in 0x00000001:sx
                # voxel id
                vid1 = x + (y-0x00000001)*sx + (z-0x00000001)*sx*sy
                # segmentation ID
                sid1 = seg[x,y,z]

                # x affinity
                if x>0x00000001 && sid1==seg[x-0x00000001,y,z]
                    vid2 = x-0x00000001 + (y-0x00000001)*sx + (z-0x00000001)*sx*sy
                    union!(dms, vid1, vid2)
                end

                # y affinity
                if y>0x00000001 && sid1==seg[x,y-0x00000001,z]
                    vid2 = x + (y-0x00000002)*sx + (z-0x00000001)*sx*sy
                    union!(dms, vid1, vid2)
                end

                # z affinity
                if z>0x00000001 && sid1==seg[x,y,z-0x00000001]
                    vid2 = x + (y-0x00000001)*sx + (z-0x00000002)*sx*sy
                    union!(dms, vid1, vid2)
                end
            end
        end
    end
    return dms
end

"""
transform indexed segmentation image to RGB image with random color label
seg: segmentation, an indexed array

Outputs:
ret: rgb image array with a size of X x Y x Z x 3, the color dim is the last one
"""
function seg2rgb(seg::Segmentation)
    # the color dict, key is segment id, value is color
    dcol = Dict{UInt32, RGB{N0f8}}()
    # set the boundary color to be black
    dcol[0] = RGB{N0f8}(0,0,0)

    # create RGB image
    sx,sy,sz = size(seg)
    ret = Array(RGB{N0f8}, (sx,sy,sz))
    # assign random color
    for z = 1:sz
        for y = 1:sy
            for x = 1:sx
                key = seg[x,y,z]
                if !haskey(dcol, key)
                    dcol[key] = rand(RGB{N0f8})
                end
                ret[x,y,z] = dcol[key]
            end
        end
        ret[i] = dcol[key]
    end
    return ret
end

"""
overlay segmentation to gray image using Alpha compositing
https://en.wikipedia.org/wiki/Alpha_compositing
Inputs:
img: gray image array
seg: segmentation, an indexed array
alpha1: the alpha value of the image
alpha2: the alpha value of the segmentation

Outputs:
ret: composited RGBA image array
"""
function seg_overlay_img(img::Array{UInt8, 3}, seg::Array{Ts,3};
                            alpha1::Normed = NORMED8_HALF,
                            alpha2::Normed = NORMED8_HALF ) where Ts
    @assert size(img)==size(seg)
    @assert alpha1>0.0f0 && alpha1<1.0f0
    @assert alpha2>0.0f0 && alpha2<1.0f0

    # initialize the returned RGBA image
    # ret = Array{RGB{U8},3}(sx,sy,sz)
    ret = similar(img, RGB{U8})

    # colorful segmentation image
    cseg = seg2rgb(Segmentation(seg))

    # transform img to 0-1
    fimg = reinterpret(N0f8, img)
    fimg = ( fimg-minimum(fimg) ) ./ (maximum(fimg) - minimum(fimg))

    # Threads.@threads for z in 1:sz
    # # for z in 1:sz
    #     for y in 1:sy
    #         for x in 1:sx
    Threads.@threads for i in eachindex( img )
        if seg[i]==0x00000000
            ret[i] = RGB{U8}(fimg[i], fimg[i], fimg[i])
        else
            ret[i] = RGB{U8}(   (fimg[i]*alpha1 + cseg[i].r * alpha2*(NORMED8_ONE-alpha1) ) / (alpha1+alpha2*(NORMED8_ONE-alpha1)),
                                (fimg[i]*alpha1 + cseg[i].g * alpha2*(NORMED8_ONE-alpha1) ) / (alpha1+alpha2*(NORMED8_ONE-alpha1)),
                                (fimg[i]*alpha1 + cseg[i].b * alpha2*(NORMED8_ONE-alpha1) ) / (alpha1+alpha2*(NORMED8_ONE-alpha1)) )
        end
    end
    return ret
end

"""
transform segmentation to sgm by making fake mst
"""
function seg2segMST(seg::Segmentation)
    # making fake mst
    segmentPairs = zeros(UInt32, (1,2))
    segmentPairs[1] = seg[1]
    segmentPairs[2] = seg[end]
    segmentPairAffinities = Vector{Float32}( [0.001] )

    return SegMST(seg, segmentPairs, segmentPairAffinities)
end
seg2sgm = seg2segMST

end # end of module


text:
"""
    ReachSet{SN}

Type that wraps a reach set, which is a set representing (an approximation of)
the reachable states for a given time interval.

### Fields

- `X`       -- set
- `t_start` -- time interval lower bound
- `t_end`   -- time interval upper bound
"""
struct ReachSet{SN}
    X::SN
    t_start::Float64
    t_end::Float64
end


text:
arr = [1,2,3,4,5]
@show arr
another = arr
@show another
another[1]=-1
@show another,arr
arr[end]=-100
@show another,arr

arr = [1,2,3,4,5]
@show arr
cparr = copy(arr)
@show cparr
cparr[1]=-1
@show cparr,arr
arr[end]=-100
@show cparr,arr

function copyarr(n)
    s = zeros(n,n)
    for i in 1:n
        sample=copy(ones(n,n))
        s+=(-1)^i*sample
    end
    s
end

ret=copyarr(11)

function copyarr2(n)
    s = zeros(n,n)
    sample = similar(s)
    for i in 1:n
        copy!(sample,ones(n,n))
        s+=(-1)^i*sample
    end
    s
end

ret2=copyarr2(11)

@assert ret2==ret

using BenchmarkTools
copyarr(303)
copyarr2(303)

p = 1+1==2 ? "seyana" : "seyaroka"
print(p)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
export parent, has_parent, root, isroot,
       lca, find_leaf, depth, find_inode,
       Tree, isequal_local, print_tree

#####################
# types and traits
#####################

"A node in a tree (of which it is the root)"
abstract type Tree <: Dag end

#####################
# methods
#####################

import Base: parent

"Get the parent of a given tree node (or nothing if the node is root)"
parent(n::Tree) = n.parent

"Does the node have a parent?"
@inline has_parent(n::Tree)::Bool = issomething(parent(n))

"Get the root of the given tree node"
root(n::Tree) =
    has_parent(n) ? root(parent(n)) : n

"Is the node the root of its tree?"
@inline isroot(n::Tree)::Bool = !has_parent(n)

"Is one node equal to another locally, ignoring children?"
isequal_local(::Tree, ::Tree)::Bool = false #default to only `===`` equality; can be extended to support value equality for certain types of trees

function Base.:(==)(n1::Tree, n2::Tree)::Bool
    (n1 === n2) && return true
    !isequal_local(n1,n2) && return false
    isleaf(n1) && isleaf(n2) && return true
    (num_children(n1) != num_children(n2)) && return false
    return all(cs -> (cs[1] == cs[2]), zip(children(n1), children(n2)))
end

"""
Find the least common ancestor. Assumes the `Tree` has access to a `parent`.
A given `descends_from` function is required to quickly check whether a node is an ancestor.
"""
function lca(v::Tree, w::Tree, descends_from::Function)::Tree
    v === w && return v
    descends_from(w,v) && return v
    candidate::Union{Dag,Nothing} = w
    while issomething(candidate)
        descends_from(v, candidate) && return candidate
        candidate = parent(candidate)
    end
    error("First argument is not contained in the root of second argument. There is no LCA.")
end
lca(::Nothing, v::Tree, ::Function)::Tree = v
lca(v::Tree, ::Nothing, ::Function)::Tree = v
lca(::Nothing, ::Nothing, ::Function)::Nothing = nothing
lca(v::Tree, w::Tree, u::Tree, r::Tree...) = lca(lca(v,w), u, r...)

"Find the leaf in the tree by follwing the branching function"
find_leaf(n::Tree, branch::Function) = find_leaf(n, NodeType(n), branch)
find_leaf(n, ::Leaf, _) = n
@inline function find_leaf(n, ::Inner, branch)::Tree
    for c in children(n)
        branch(c) && return find_leaf(c, branch)
    end
    error("Could not find a branch that contains the desired leaf")
end

"Compute the length of the path from a tree node to the leaf following the branching function"
depth(n::Tree, branch::Function)::Int = depth(n, NodeType(n), branch)
depth(_, ::Leaf, _) = 0

function depth(n::Tree, ::Inner, branch)::Int
    for c in children(n)
        branch(c) && return 1 + depth(c, branch)
    end
    error("Could not find a branch that contains the desired leaf")
end

"""
Find a binary inner node that has `left` in its left subtree and `right` in its right subtree.
Supports `nothing` as a catch-all for either left or right. Returns nothing if no such node exists.
"""
function find_inode(left, right, descends_from::Function)
    isnothing(left) && isnothing(right) && return nothing
    candidate = lca(left, right, descends_from)
    while issomething(candidate)
        if isinner(candidate) && num_children(candidate) == 2 &&
            (isnothing(left) ||
                descends_from(left, children(candidate)[1])) &&
            (isnothing(right) ||
                descends_from(right, children(candidate)[2]))
            return candidate
        end
        candidate = parent(candidate)
    end
    nothing # there is no solution
end

"Print the given tree"
print_tree(root::Tree, io::IO=stdout) =
    print_tree(io, root)

print_tree(io::IO, root::Tree, prefix="", onceprefix="", laterprefix="") =
    print_tree(io, root, NodeType(root), prefix, onceprefix, laterprefix)

function print_tree(io::IO, root::Tree, ::Leaf, prefix, onceprefix, laterprefix)
    println(io, "$prefix$(onceprefix)━$root")
end

function print_tree(io::IO, root::Tree, ::Inner, prefix, onceprefix, laterprefix)
    println(io, "$prefix$(onceprefix)━$root")
    for c in children(root)[1:end-1]
        print_tree(io, c, "$prefix$laterprefix", " ┣"," ┃")
    end
    print_tree(io, children(root)[end], "$prefix$laterprefix", " ┗", "  ")
end

#####################
# traversal
#####################

function foreach(f::Function, node::Tree)
    if isinner(node)
        for c in children(node)
            foreach(f, c)
        end
    end
    f(node)
    nothing
end

function foldup(node::Tree, f_leaf::Function, f_inner::Function, ::Type{T})::T where T
    v = if isinner(node)
        callback(c) = (foldup(c, f_leaf, f_inner, T)::T)
        f_inner(node, callback)::T
    else
        f_leaf(node)::T
    end
    return v
end

function foldup_aggregate(node::Tree, f_leaf::Function, f_inner::Function, ::Type{T})::T where T
    v = if isinner(node)
        child_values = Vector{T}(undef, num_children(node))
        map!(c -> foldup_aggregate(c, f_leaf, f_inner, T)::T, child_values, children(node))
        f_inner(node, child_values)::T
    else
        f_leaf(node)::T
    end
    return v
end


text:
import StaticArrays: SVector

#########################################################
# Helper functions to plot 3D rectangles.
#########################################################
"""
    rectangle3dvertices(x, y, z, ϵx, ϵy, ϵz)

Compute the coordinates of the vertices of a 3D rectangle, given the
coordinates of the origin (`x`, `y`, and `z`) and the corresponding edge lengths
(`ϵx`, `ϵy`, `ϵz`).
"""
function rectangle3dvertices(x, y, z, ϵx, ϵy, ϵz)
    v1b = [x,    y,    z]
    v2b = [x+ϵx, y,    z]
    v3b = [x+ϵx, y+ϵy, z]
    v4b = [x,    y+ϵy, z]
    v1t = [x,    y,    z+ϵz]
    v2t = [x+ϵx, y,    z+ϵz]
    v3t = [x+ϵx, y+ϵy, z+ϵz]
    v4t = [x,    y+ϵy, z+ϵz]
    (v1b, v2b, v3b, v4b, v1t, v2t, v3t, v4t)
end

"""
    rectangle3dvertices(o, ϵ)

Given the origin `o` (3-element vector) and edge lengths `ϵ` (also a
3-element vector) of a 3D rectangle, return the coordinates of its
vertices.
"""
function rectangle3dvertices(o, ϵ)
    x, y, z = (o...,)
    ϵx, ϵy, ϵz = (ϵ...,)
    rectangle3dvertices(x, y, z, ϵx, ϵy, ϵz)
end

"""
    connectvertices3d(o, ϵ)

Given the origin `o` (3-element vector) and edge lengths `ϵ` (also a
3-element vector) of a 3D rectangle, return a vector of line segments
(each a dim-by-n_vertices array) that when plotted together yields
the entire rectangle.
"""
function connectvertices3d(o, ϵ)
    # Get the vertices
    v1b, v2b, v3b, v4b, v1t, v2t, v3t, v4t = rectangle3dvertices(o, ϵ)

    # Connect vertices in top and bottom planes
    connect_top = zeros(3, 5)
    connect_bottom = zeros(3, 5)

    # Connect corners of the top and bottom planes
    corner1 = zeros(3, 2)
    corner2 = zeros(3, 2)
    corner3 = zeros(3, 2)
    corner4 = zeros(3, 2)
    connect_bottom[:, 1:5] = [v1b v2b v3b v4b v1b]
    connect_top[:, 1:5]    = [v1t v2t v3t v4t v1t]
    corner1[:, 1:2] = [v1b v1t]
    corner2[:, 1:2] = [v2b v2t]
    corner3[:, 1:2] = [v3b v3t]
    corner4[:, 1:2] = [v4b v4t]

    return [connect_top, connect_bottom,
            corner1, corner2, corner3, corner4]
end

export connectvertices3d

"""
    splitaxes(x)

Return a vector of the individual components of an array of points
provided as an array where each point is a column.
"""
splitaxes(x) = ([x[k, :] for k = 1:size(x, 1)]...,)

"""
    plot_3D_rect!(p, origin, edgelengths;
             lc = :black, lw = 1.0, ls = :solid)

Append a a 3D rectangle to an existing plot `p`, given the
origin `o` (3-element vector) and edge lengths
`ϵ` (3-element vector) of the rectangle.
"""
function plot_3D_rect!(p, o, ϵ;
             lc = :black, lw = 1.0, ls = :solid, lα = 0.7)

    linesegments = connectvertices3d(o, ϵ)
    for segment in linesegments
        plot!(p, splitaxes(segment),
            lc = lc, lw = lw, ls = ls, lα = lα)
    end
end

"""
Given the `origin` of a hyperrectangle and the `stepsizes` along each
coordinate axis, divide each axis into `N` stepsizes starting
from `origin[i] + (stepsizes[i]/N)/2`
"""
function fill_hyperrectangle_3D(origin, stepsizes, N = 10)
    xs, ys, zs = ([LinRange(origin[i] + stepsizes[i]/N/2, origin[i] + stepsizes[i] - stepsizes[i]/N/2, N) for i = 1:length(origin)]...,)
    pts = zeros(Float64, 3, length(xs)^3)
    i = 0
    for x in xs
        for y in ys
            for z in zs
                i += 1
                pts[:, i] = [x, y, z]
            end
        end
    end
    pts
end

struct Rectangle{D, T}
    origin::SVector{D, T}
    edgelengths::SVector{D, T}
end

rectangle3dvertices(r::Rectangle) = rectangle3dvertices((r.origin...,), (r.edgelengths...,))

function linesegments(r::Rectangle)
    # Get the vertices
    v1b, v2b, v3b, v4b, v1t, v2t, v3t, v4t = rectangle3dvertices(r.origin, r.edgelengths)

    # Connect vertices in top and bottom planes
    connect_top = zeros(3, 5)
    connect_bottom = zeros(3, 5)

    # Connect corners of the top and bottom planes
    corner1 = zeros(3, 2)
    corner2 = zeros(3, 2)
    corner3 = zeros(3, 2)
    corner4 = zeros(3, 2)
    connect_bottom[:, 1:5] = [v1b v2b v3b v4b v1b]
    connect_top[:, 1:5]    = [v1t v2t v3t v4t v1t]
    corner1[:, 1:2] = [v1b v1t]
    corner2[:, 1:2] = [v2b v2t]
    corner3[:, 1:2] = [v3b v3t]
    corner4[:, 1:2] = [v4b v4t]

    linesegments = [connect_top, connect_bottom, corner1, corner2, corner3, corner4]
    return linesegments
end

function connectrect(r::Rectangle)
    v1b, v2b, v3b, v4b, v1t, v2t, v3t, v4t = rectangle3dvertices(r.origin, r.edgelengths)

    x = hcat(v1b, v2b, v2t, v3t, v3b, v4b, v4t, v1t, v1b, v1t, v2t, v2b, v3b, v3t, v4t, v4b, v1b)
    splitaxes(x)
end

export
Rectangle,
connectrect,
rectangle3dvertices,
linesegments,
connectvertices3d


text:
module FrostHelperCustomRisingLavaStartHeightTriggerModule

using ..Ahorn, Maple

@mapdef Trigger "FrostHelper/CustomRisingLavaStartHeightTrigger" FrostHelperCustomRisingLavaStartHeightTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16)

Ahorn.nodeLimits(entity::FrostHelperCustomRisingLavaStartHeightTrigger) = 1, 1

const placements = Ahorn.PlacementDict(
	"Custom Rising Lava Start Height (Frost Helper)" => Ahorn.EntityPlacement(
		FrostHelperCustomRisingLavaStartHeightTrigger,
		"rectangle",
        Dict{String, Any}(),
        function(entity)
            entity.data["nodes"] = [(Int(entity.data["x"]), Int(entity.data["y"]) + 16)]
        end
	)
)

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
# Test automatic differentiation
@testset "autodiff" begin
function f!(fvec, x)
    fvec[1] = (x[1]+3)*(x[2]^3-7)+18
    fvec[2] = sin(x[2]*exp(x[1])-1)
end

r = nlsolve(f!, [ -0.5; 1.4], autodiff = true)
@test converged(r)
@test norm(r.zero - [ 0; 1]) < 1e-8

end


text:
# In this example, we use EKI to tune the closure parameters of a HydrostaticFreeSurfaceModel
# with a CATKEBasedVerticalDiffusivity closure in order to align the predictions of the model
# to those of a high-resolution LES data generated in LESbrary.jl. Here `predictions` refers to the
# 1-D profiles of temperature, velocity, and turbulent kinetic energy horizontally averaged over a
# 3-D physical domain.

pushfirst!(LOAD_PATH, joinpath(@__DIR__, ".."))
pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", "projects", "OceanBoundaryLayerParameterizations", "src"))

using Oceananigans
using OceanTurbulenceParameterEstimation
using OceanBoundaryLayerParameterizations

# CATKE parameters involved in setting field diffusivities
@free_parameters StabilityFnParameters CᴷRiʷ CᴷRiᶜ Cᴷu⁻ Cᴷuʳ Cᴷc⁻ Cᴷcʳ Cᴷe⁻ Cᴷeʳ

parameters = Parameters(
    RelevantParameters = CATKEParametersRiDependent,  # Parameters that are used in CATKE
    ParametersToOptimize = StabilityFnParameters    # Subset of RelevantParameters that we want to optimize
)

two_day_suite_dir = "/home/greg/Projects/OceanTurbulenceParameterEstimation/data/2DaySuite"
two_day_suite = TwoDaySuite(two_day_suite_dir)

# InverseProblem represents the model, data, loss function, and parameters
calibration = InverseProblem(two_day_suite, # "Truth data" for model calibration
                             parameters;   # Model parameters
                             # Loss function parameters
                             relative_weights = relative_weight_options["all_but_e"],
                             # Model (hyper)parameters
                             architecture = GPU(),
                             ensemble_size = 10,
                             Nz = 16,
                             Δt = 1.0)

# Run the model forward and store the solution
output = model_time_series(calibration, calibration.default_parameters)


text:
"""
Class for sparse variational Gaussian Processes

```julia
SVGP(X::AbstractArray{T1},y::AbstractVector{T2},kernel::Kernel,
    likelihood::LikelihoodType,inference::InferenceType, nInducingPoints::Int;
    verbose::Int=0,optimizer::Union{Optimizer,Nothing,Bool}=Adam(α=0.01),atfrequency::Int=1,
    mean::Union{<:Real,AbstractVector{<:Real},PriorMean}=ZeroMean(),
    Zoptimizer::Union{Optimizer,Nothing,Bool}=false,
    ArrayType::UnionAll=Vector)
```

Argument list :

**Mandatory arguments**
 - `X` : input features, should be a matrix N×D where N is the number of observation and D the number of dimension
 - `y` : input labels, can be either a vector of labels for multiclass and single output or a matrix for multi-outputs (note that only one likelihood can be applied)
 - `kernel` : covariance function, can be either a single kernel or a collection of kernels for multiclass and multi-outputs models
 - `likelihood` : likelihood of the model, currently implemented : Gaussian, Student-T, Laplace, Bernoulli (with logistic link), Bayesian SVM, Multiclass (softmax or logistic-softmax) see [`Likelihood`](@ref likelihood_user)
 - `inference` : inference for the model, can be analytic, numerical or by sampling, check the model documentation to know what is available for your likelihood see the [`Compatibility table`](@ref compat_table)
 - `nInducingPoints` : number of inducing points
**Optional arguments**
 - `verbose` : How much does the model print (0:nothing, 1:very basic, 2:medium, 3:everything)
 - `optimizer` : Optimizer for kernel hyperparameters (to be selected from [GradDescent.jl](https://github.com/jacobcvt12/GradDescent.jl)) or set it to `false` to keep hyperparameters fixed
 - `atfrequency` : Choose how many variational parameters iterations are between hyperparameters optimization
 - `mean` : PriorMean object, check the documentation on it [`MeanPrior`](@ref meanprior)
 - `IndependentPriors` : Flag for setting independent or shared parameters among latent GPs
 - `optimizer` : Optimizer for inducing point locations (to be selected from [GradDescent.jl](https://github.com/jacobcvt12/GradDescent.jl))
 - `ArrayType` : Option for using different type of array for storage (allow for GPU usage)
"""
mutable struct SVGP{T<:Real,TLikelihood<:Likelihood{T},TInference<:Inference,TGP<:Abstract_GP{T},N} <: AbstractGP{T,TLikelihood,TInference,TGP,N}
    X::Matrix{T} #Feature vectors
    y::Vector #Output (-1,1 for classification, real for regression, matrix for multiclass)
    nSamples::Int64 # Number of data points
    nDim::Int64 # Number of covariates per data point
    nFeatures::Int64 # Number of features of the GP (equal to number of points)
    nLatent::Int64 # Number pf latent GPs
    f::NTuple{N,TGP}
    likelihood::TLikelihood
    inference::TInference
    verbose::Int64
    atfrequency::Int64
    Trained::Bool
end

function SVGP(X::AbstractArray{T1},y::AbstractVector{T2},kernel::Kernel,
            likelihood::TLikelihood,inference::TInference, nInducingPoints::Int;
            verbose::Int=0,optimizer::Union{Optimizer,Nothing,Bool}=Adam(α=0.01),atfrequency::Int=1,
            mean::Union{<:Real,AbstractVector{<:Real},PriorMean}=ZeroMean(), variance::Real = 1.0,
            Zoptimizer::Union{Optimizer,Nothing,Bool}=false,
            ArrayType::UnionAll=Vector) where {T1<:Real,T2,TLikelihood<:Likelihood,TInference<:Inference}

            X,y,nLatent,likelihood = check_data!(X,y,likelihood)
            @assert check_implementation(:SVGP,likelihood,inference) "The $likelihood is not compatible or implemented with the $inference"

            nSamples = size(X,1); nDim = size(X,2);
            if isa(optimizer,Bool)
                optimizer = optimizer ? Adam(α=0.01) : nothing
            end

            @assert nInducingPoints > 0 && nInducingPoints < nSamples "The number of inducing points is incorrect (negative or bigger than number of samples)"
            Z = KMeansInducingPoints(X,nInducingPoints,nMarkov=10)
            if isa(Zoptimizer,Bool)
                Zoptimizer = Zoptimizer ? Adam(α=0.01) : nothing
            end
            Z = InducingPoints(Z,Zoptimizer)

            nFeatures = nInducingPoints

            if typeof(mean) <: Real
                mean = ConstantMean(mean)
            elseif typeof(mean) <: AbstractVector{<:Real}
                mean = EmpiricalMean(mean)
            end

            nMinibatch = nSamples
            if inference.Stochastic
                @assert inference.nMinibatch > 0 && inference.nMinibatch < nSamples "The size of mini-batch $(inference.nMinibatch) is incorrect (negative or bigger than number of samples), please set nMinibatch correctly in the inference object"
                nMinibatch = inference.nMinibatch
            end

            latentf = ntuple( _ -> _SVGP{T1}(nFeatures,nMinibatch,Z,kernel,mean,variance,optimizer),nLatent)

            likelihood = init_likelihood(likelihood,inference,nLatent,nMinibatch,nFeatures)
            inference = tuple_inference(inference,nLatent,nFeatures,nSamples,nMinibatch)
            inference.xview = view(X,1:nMinibatch,:)
            inference.yview = view_y(likelihood,y,1:nMinibatch)

            model = SVGP{T1,TLikelihood,typeof(inference),_SVGP{T1},nLatent}(X,y,
                    nSamples, nDim, nFeatures, nLatent,
                    latentf,likelihood,inference,
                    verbose,atfrequency,false)
            # if isa(inference.optimizer,ALRSVI)
                # init!(model.inference,model)
            # end
            # return model
end

function Base.show(io::IO,model::SVGP{T,<:Likelihood,<:Inference}) where {T}
    print(io,"Sparse Variational Gaussian Process with a $(model.likelihood) infered by $(model.inference) ")
end

const SVGP1 = SVGP{<:Real,<:Likelihood,<:Inference,<:Abstract_GP,1}

get_y(model::SVGP) = model.inference.yview
get_Z(model::SVGP) = getproperty.(getproperty.(model.f,:Z),:Z)


text:
@testset "Visualization Forecast" begin
    air_passengers = CSV.File(StateSpaceModels.AIR_PASSENGERS) |> DataFrame
    log_air_passengers = log.(air_passengers.passengers)
    model = BasicStructural(log_air_passengers, 12)
    fit!(model)
    # forecasting
    forec = forecast(model, 12)
    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), model, forec)
    @test length(rec) == 4
    # simulating
    scen = simulate_scenarios(model, 12, 100)
    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), model, scen)
    @test length(rec) == 2
end


text:
"""
    Segment{T<:Real}

The intersection of a mesh element or cell and a [`Track`](@ref) yields to a `Segment` with
entry point `p`, exit point `q` and length `ℓ`.
"""
struct Segment{T<:Real}
    p::Point2D{T}
    q::Point2D{T}
    ℓ::T
    τ::Vector{T}
    element::Int32 # TODO: deberia llamarlo cell_id o algo asi
end

function Segment(p::Point2D{T}, q::Point2D{T}, element::Int32=Int32(-1)) where {T}
    return Segment(p, q, norm(p - q), Vector{T}(undef, 0), element)
end

ℓ(segment::Segment) = segment.ℓ

in(x::Point2D, segment::Segment) = point_in_segment(segment.p, segment.q, x)

function point_in_segment(p::Point2D, q::Point2D, x::Point2D)
    lpx = norm(p - x)
    lqx = norm(q - x)
    lpq = norm(p - q)
    return isapprox(lpx + lqx, lpq)
end


text:
include("cg.jl")
using CG;
include("render.jl");

v1 = Vec4f(1, 0, 0, 0);
v2 = Vec4f(0.0, 42.0, 0.0, 0.0);

println(v1 - v2);

m1 = Mat4f(Vec4f(1.0, 0.0, 0.0, 0.0),
           Vec4f(0.0, 1.0, 0.0, 0.0),
           Vec4f(0.0, 0.0, 1.0, 0.0),
           Vec4f(0.0, 0.0, 0.0, 1.0));

println(2 * (m1 * v1));
println(m1 * v2);

v1 = Vec4f(1, 0, 0, 1);
v2 = Vec4f(0, 0, 0, 1);
v3 = Vec4f(0, 0, 1, 1);

triangle = Object(v1, v2, v3);

house = Object(Vec4f(-1, -1, 0, 1),
Vec4f(1, -1, 0, 1),
Vec4f(-1, 1, 0, 1),
Vec4f(0, 2, 0, 1),
Vec4f(1, 1, 0, 1),
Vec4f(1, 1, 0, 1),
Vec4f(-1, -1, 0, 1),
Vec4f(-1, 1, 0, 1),
Vec4f(1, 1, 0, 1),
Vec4f(1, -1, 0, 1))

#render(translation(Vec4f(1.5, 0, 0, 1)) * house, figAxis=[-2, 2, -2, 2]);
#render(scaling(Vec4f(0.5, 0.5, 0.5, 1)) * house, figAxis=[-2, 2, -2, 2]);
render(roty(pi/4) * house, figAxis=[-2, 2, -2, 2]);

max_stars_repo_name: fpoppinga/cg

max_issues_repo_name: fpoppinga/cg

max_forks_repo_name: fpoppinga/cg


text:
@doc """
 COST 231- Cost-Waldrosch-Ikegami Model
 freq: signal frequency
 hb: base station's height
 hm: mobile's height
 ws: average width of the street in meters
 bs: average setback of buildings in meters
 hr: mean height of houses in meters
 areaKind: area type (1-rural, 2-suburban e 3-urban).
 cityKind : cyte type(1-small, 2-medium e 3-large).
"""
type Cost231Model<:Model
  freq::Float64
  txH::Float64
  rxH::Float64
  ws::Float64
  bs::Float64
  hr::Float64
  cityKind::Int
  checkFreqRange::Bool
  Cost231Model() = new(800,50.0,1.5,20,10,35,CityKind.Medium,true)
end

"Calculate Path loss for distance dist (km) with COST 231 Model"
function pathloss_op(model::Cost231Model, dist::AbstractFloat)
  #150 a 2000 MHz
  hm::Float64 = model.rxH
  hb::Float64 = model.txH
  hr::Float64 = model.hr
  d::Float64 = dist
  f::Float64 = model.freq
  ws::Float64 = model.ws
  bs::Float64 = model.bs
  # Checking the frequency range
  if (model.checkFreqRange)
    if (f < 150 || f > 2000)
     error("The frequency range for Ecc-33 Model is 150MHz-2000Mhz")
    end
  end

  deltaH::Float64 = hm/hb #relaction between heighths
  Lbsh::Float64 = 18*log(1+deltaH) # Loss due to difference of heights
  Ka::Float64=54.0  #Coefficient of proximity Buildings
  Kd::Float64=18.0  #Coeficiente of proximidade Edifica??es
  Kf::Float64=4.0  #Coeficient of environment(urban or not)

  #Coeficient's calculate
  if (hr > hb)
    Lbsh = 0
  end
  if (hb <= hr && d >= 0.5)
    ka=Ka-0.8*deltaH
  elseif (hb<=hr && d<0.5)
    ka= Ka-0.8*deltaH*(d/0.5)
  end
  if (hb < hr)
    Kd= Kd-15*(hb-hr)/(hr-hm)
  end
  if (model.cityKind==CityKind.Small)
    Kf=Kf+0.7*(f/925-1)
  else
    Kf=Kf+1.5*(f/925-1)
  end
  #path loss's calculate
  Lo = 32.4+20*log10(d)+20*log10(f)                     #free space path loss
  Lrts = 8.2+10*log(ws) + 10*log10(f) + 10*log(deltaH) # roofTop loss
  Lmsd =Lbsh +Ka +Kd*log10(d)+Kf*log(f)-9*log10(bs)    #Multpath loss
  #final path loss
  PL = Lo + Lrts + Lmsd;
  return PL
end


text:
using Test
using DECAES
using Logging

function compare_epg(
        work₁::DECAES.AbstractEPGWorkspace{T,ETL},
        work₂::DECAES.AbstractEPGWorkspace{T,ETL},
    ) where {T,ETL}
    α, TE, T2, T1, β = T(163.0), T(11e-3), T(39e-3), T(1.1), T(151.0)
    θ   = DECAES.EPGOptions((; α, TE, T2, T1, β), Val(ETL), T)
    dc₁ = DECAES.EPGdecaycurve!(work₁, θ)
    dc₂ = DECAES.EPGdecaycurve!(work₂, θ)
    if !(dc₁ ≈ dc₂)
        @info "Comparing: $((nameof(typeof(work₁)), nameof(typeof(work₂))))"
        @info "    max error:   $(maximum(abs, dc₁ .- dc₂))"
        @info "    diff vector: $(round.(abs.(dc₁ .- dc₂)'; sigdigits = 4))"
    end
    @test dc₁ ≈ dc₂
end

@testset "EPG algorithms" begin
    # NOTE: Generated function approach is extremely slow to compile for large ETL (around 16)
    epg_algs = DECAES.EPG_Algorithms
    for i in 1:length(epg_algs), j in 1:i-1, T in [Float32, Float64], ETL in [4,5,6,7] # test four ETL >= 4 which are unique mod 4
        algᵢ, algⱼ = epg_algs[i], epg_algs[j]
        compare_epg(algᵢ(T, ETL), algⱼ(T, ETL))
    end
end

nothing


text:
# This file was generated, do not modify it. # hide
using MLJ
import RDatasets: dataset
import DataFrames: DataFrame, describe, select, Not
import StatsBase: countmap, cor, var
MLJ.color_off() # hide
using PrettyPrinting

smarket = dataset("ISLR", "Smarket")
@show size(smarket)
@show names(smarket)


text:
using DistributedFactorGraphs
using IncrementalInference
using Interpolations
using Manifolds
using StaticArrays
using Test
import IncrementalInference: LevelSetGridNormal

## define new local variable types for testing

@defVariable TranslationGroup2 TranslationGroup(2) [0.0, 0.0]

# @defVariable SpecialEuclidean2 SpecialEuclidean(2) ProductRepr(@MVector([0.0,0.0]), @MMatrix([1.0 0.0; 0.0 1.0]))
@defVariable SpecialEuclidean2 SpecialEuclidean(2) ProductRepr([0.0,0.0], [1.0 0.0; 0.0 1.0])

##

@testset "Test SpecialEuclidean(2)" begin

##

Base.convert(::Type{<:Tuple}, M::SpecialEuclidean{2}) = (:Euclid, :Euclid, :Circular)
Base.convert(::Type{<:Tuple}, ::IIF.InstanceType{SpecialEuclidean{2}})  = (:Euclid, :Euclid, :Circular)

M = getManifold(SpecialEuclidean2)
@test M == SpecialEuclidean(2)
pT = getPointType(SpecialEuclidean2)
@test pT == ProductRepr{Tuple{Vector{Float64}, Matrix{Float64}}}
# @test pT == ProductRepr{Tuple{MVector{2, Float64}, MMatrix{2, 2, Float64, 4}}}
pϵ = getPointIdentity(SpecialEuclidean2)
# @test_broken pϵ == ProductRepr(@MVector([0.0,0.0]), @MMatrix([1.0 0.0; 0.0 1.0]))
@test all(isapprox.(pϵ,ProductRepr([0.0,0.0], [1.0 0.0; 0.0 1.0])).parts)

@test is_point(getManifold(SpecialEuclidean2), getPointIdentity(SpecialEuclidean2))

##
fg = initfg()

v0 = addVariable!(fg, :x0, SpecialEuclidean2)

# mp = ManifoldPrior(SpecialEuclidean(2), ProductRepr(@MVector([0.0,0.0]), @MMatrix([1.0 0.0; 0.0 1.0])), MvNormal([0.01, 0.01, 0.01]))
mp = ManifoldPrior(SpecialEuclidean(2), ProductRepr(@MVector([0.0,0.0]), @MMatrix([1.0 0.0; 0.0 1.0])), MvNormal(Diagonal(abs2.([0.01, 0.01, 0.01]))))
p = addFactor!(fg, [:x0], mp)

##

doautoinit!(fg, :x0)

##
vnd = getVariableSolverData(fg, :x0)
@test all(isapprox.(mean(vnd.val), ProductRepr([0.0,0.0], [1.0 0.0; 0.0 1.0]), atol=0.1).parts)
@test all(is_point.(Ref(M), vnd.val))

##
v1 = addVariable!(fg, :x1, SpecialEuclidean2)
mf = ManifoldFactor(SpecialEuclidean(2), MvNormal([1,2,pi/4], [0.01,0.01,0.01]))
f = addFactor!(fg, [:x0, :x1], mf)

doautoinit!(fg, :x1)

vnd = getVariableSolverData(fg, :x1)
@test all(isapprox.(mean(vnd.val), ProductRepr([1.0,2.0], [0.7071 -0.7071; 0.7071 0.7071]), atol=0.1).parts)
@test all(is_point.(Ref(M), vnd.val))

##
smtasks = Task[]
solveTree!(fg; smtasks, verbose=true) #, recordcliqs=ls(fg))
# hists = fetchCliqHistoryAll!(smtasks);

vnd = getVariableSolverData(fg, :x0)
@test all(isapprox.(mean(vnd.val), ProductRepr([0.0,0.0], [1.0 0.0; 0.0 1.0]), atol=0.1).parts)
@test all(is_point.(Ref(M), vnd.val))

vnd = getVariableSolverData(fg, :x1)
@test all(isapprox.(mean(vnd.val), ProductRepr([1.0,2.0], [0.7071 -0.7071; 0.7071 0.7071]), atol=0.1).parts)
@test all(is_point.(Ref(M), vnd.val))

v1 = addVariable!(fg, :x2, SpecialEuclidean2)
mf = ManifoldFactor(SpecialEuclidean(2), MvNormal([1,2,pi/4], [0.01,0.01,0.01]))
f = addFactor!(fg, [:x1, :x2], mf)

##

#test new error from solvetree
# smtasks = Task[]
@test solveTree!(fg; smtasks, verbose=true) isa AbstractBayesTree

## test partial prior issue

fg = initfg()

v0 = addVariable!(fg, :x0, SpecialEuclidean2)
mp = PartialPrior(MvNormal([0.01, 0.01]), (1,2))

p = addFactor!(fg, [:x0], mp, graphinit=false)

##

pbel_ = approxConvBelief(fg, :x0f1, :x0)

@test isPartial(pbel_)

@test pbel_._partial == [1;2]
@test length(pbel_.infoPerCoord) == 3

##
end

##
struct ManifoldFactorSE2{T <: SamplableBelief} <: IIF.AbstractManifoldMinimize
    Z::T
end

ManifoldFactorSE2() = ManifoldFactorSE2(MvNormal(Diagonal([1,1,1])))
DFG.getManifold(::ManifoldFactorSE2) = SpecialEuclidean(2)

IIF.selectFactorType(::Type{<:SpecialEuclidean2}, ::Type{<:SpecialEuclidean2}) = ManifoldFactorSE2

function IIF.getSample(cf::CalcFactor{<:ManifoldFactorSE2})
  M = SpecialEuclidean(2)
  ϵ = identity_element(M)
  X = sampleTangent(M, cf.factor.Z, ϵ)
  return X
end

function (cf::CalcFactor{<:ManifoldFactorSE2})(X, p, q)
    M = SpecialEuclidean(2)
    q̂ = Manifolds.compose(M, p, exp(M, identity_element(M, p), X)) #for groups
    Xc = zeros(3)
    vee!(M, Xc, q, log(M, q, q̂))
    return Xc
end

##

@testset "Test Pose2 like hex as SpecialEuclidean2" begin
##

M = getManifold(SpecialEuclidean2)
fg = initfg()
v0 = addVariable!(fg, :x0, SpecialEuclidean2)

mp = ManifoldPrior(SpecialEuclidean(2), ProductRepr(@MVector([10.0,10.0]), @MMatrix([-1.0 0.0; 0.0 -1.0])), MvNormal([0.1, 0.1, 0.01]))
p = addFactor!(fg, [:x0], mp)

##

for i in 0:5
    psym = Symbol("x$i")
    nsym = Symbol("x$(i+1)")
    addVariable!(fg, nsym, SpecialEuclidean2)
    mf = ManifoldFactor(SpecialEuclidean(2), MvNormal([10.0,0,pi/3], [0.5,0.5,0.05]))
    f = addFactor!(fg, [psym;nsym], mf)
end

addVariable!(fg, :l1, SpecialEuclidean2, tags=[:LANDMARK;])
mf = ManifoldFactor(SpecialEuclidean(2), MvNormal([10.0,0,0], [0.1,0.1,0.01]))
addFactor!(fg, [:x0; :l1], mf)

mf = ManifoldFactor(SpecialEuclidean(2), MvNormal([10.0,0,0], [0.1,0.1,0.01]))
addFactor!(fg, [:x6; :l1], mf)

vnd = getVariableSolverData(fg, :x0)
@test isapprox(M, mean(M, vnd.val), ProductRepr([10.0,10.0], [-1.0 0.0; 0.0 -1.0]), atol=0.2)

vnd = getVariableSolverData(fg, :x1)
@test isapprox(M, mean(M, vnd.val), ProductRepr([0.0,10.0], [-0.5 0.866; -0.866 -0.5]), atol=0.4)

vnd = getVariableSolverData(fg, :x6)
@test isapprox(M, mean(M, vnd.val), ProductRepr([10.0,10.0], [-1.0 0.0; 0.0 -1.0]), atol=0.5)

##

smtasks = Task[]
solveTree!(fg; smtasks);

## Special test for manifold based messages

#FIXME this may show some bug in propagateBelief caused by empty factors
fg.solverParams.useMsgLikelihoods = true
@test_broken solveTree!(fg; smtasks) isa AbstractBayesTree

##
end

@testset "test deconv on <:AbstractManifoldMinimize" begin
##

fg = initfg()
getSolverParams(fg).useMsgLikelihoods = true

addVariable!(fg, :x0, SpecialEuclidean2)
addVariable!(fg, :x1, SpecialEuclidean2)

mp = ManifoldPrior(SpecialEuclidean(2), ProductRepr(@MVector([10.0,10.0]), @MMatrix([-1.0 0.0; 0.0 -1.0])), MvNormal([0.1, 0.1, 0.01]))
p = addFactor!(fg, [:x0], mp)

doautoinit!(fg,:x0)

addFactor!(fg, [:x0;:x1], ManifoldFactorSE2(MvNormal([10.0,0,0.1], diagm([0.5,0.5,0.05].^2))))

initAll!(fg)

# now check deconv

pred, meas = approxDeconv(fg, :x0x1f1)

@test mmd(SpecialEuclidean(2), pred, meas) < 1e-1

p_t = map(x->x.parts[1], pred)
m_t = map(x->x.parts[1], meas)
p_θ = map(x->x.parts[2][2], pred)
m_θ = map(x->x.parts[2][2], meas)

@test isapprox(mean(p_θ), 0.1, atol=0.02)
@test isapprox(std(p_θ), 0.05, atol=0.02)

@test isapprox(mean(p_t), [10,0], atol=0.3)
@test isapprox(std(p_t), [0.5,0.5], atol=0.3)

@test isapprox(mean(p_θ), mean(m_θ), atol=0.02)
@test isapprox(std(p_θ), std(m_θ), atol=0.02)

@test isapprox(mean(p_t), mean(m_t), atol=0.3)
@test isapprox(std(p_t), std(m_t), atol=0.3)

##
end

## ======================================================================================
##
## ======================================================================================

struct ManiPose2Point2{T <: SamplableBelief} <: IIF.AbstractManifoldMinimize
    Z::T
end

function IIF.getSample(cf::CalcFactor{<:ManiPose2Point2})
    return rand(cf.factor.Z)
end

DFG.getManifold(::ManiPose2Point2) = TranslationGroup(2)

# define the conditional probability constraint
function (cfo::CalcFactor{<:ManiPose2Point2})(measX, p, q)
  #
    M = SpecialEuclidean(2)
    q_SE = ProductRepr(q, identity_element(SpecialOrthogonal(2), p.parts[2]))

    X_se2 = log(M, identity_element(M, p), Manifolds.compose(M, inv(M, p), q_SE))
    X = X_se2.parts[1]
    # NOTE wrong for what we want X̂ = log(M, p, q_SE)
    return measX - X
end

##
@testset "Test SpecialEuclidean(2)" begin
##

# Base.convert(::Type{<:Tuple}, M::TranslationGroup{Tuple{2},ℝ}) = (:Euclid, :Euclid)
# Base.convert(::Type{<:Tuple}, ::IIF.InstanceType{TranslationGroup{Tuple{2},ℝ}})  = (:Euclid, :Euclid)

##
fg = initfg()

v0 = addVariable!(fg, :x0, SpecialEuclidean2)

mp = ManifoldPrior(SpecialEuclidean(2), ProductRepr(@MVector([0.0,0.0]), @MMatrix([1.0 0.0; 0.0 1.0])), MvNormal([0.01, 0.01, 0.01]))
p = addFactor!(fg, [:x0], mp)

##
v1 = addVariable!(fg, :x1, TranslationGroup2)
mf = ManiPose2Point2(MvNormal([1,2], [0.01,0.01]))
f = addFactor!(fg, [:x0, :x1], mf)

doautoinit!(fg, :x1)

vnd = getVariableSolverData(fg, :x1)
@test all(isapprox.(mean(vnd.val), [1.0,2.0], atol=0.1))

# ##
smtasks = Task[]
solveTree!(fg; smtasks, verbose=true, recordcliqs=ls(fg))
# # hists = fetchCliqHistoryAll!(smtasks);

vnd = getVariableSolverData(fg, :x0)
@test all(isapprox.(mean(vnd.val), ProductRepr([0.0,0.0], [1.0 0.0; 0.0 1.0]), atol=0.1).parts)

vnd = getVariableSolverData(fg, :x1)
@test all(isapprox.(mean(vnd.val), [1.0,2.0], atol=0.1))

##
end

@testset "test propagateBelief w HeatmapSampler and init for PartialPriorPassThrough w Priors" begin
##

fg = initfg()

v0 = addVariable!(fg, :x0, SpecialEuclidean2)

img_ = rand(10,10).+5.0
x_,y_ = ([-9:2.0:9;],[-9:2.0:9;])

hmd = LevelSetGridNormal(img_, (x_,y_), 5.5, 0.1, N=120)
pthru = PartialPriorPassThrough(hmd, (1,2))

## quick

pf = convert( AbstractPackedFactor, pthru )
upf = convert( AbstractFactor, pf )

@test pthru.partial == upf.partial
@test isapprox( pthru.Z.heatmap.data, upf.Z.heatmap.data )
@test isapprox( pthru.Z.heatmap.domain[1], upf.Z.heatmap.domain[1] )
@test isapprox( pthru.Z.heatmap.domain[2], upf.Z.heatmap.domain[2] )
@test isapprox( pthru.Z.level, upf.Z.level )
@test isapprox( pthru.Z.sigma, upf.Z.sigma )
@test isapprox( pthru.Z.sigma_scale, upf.Z.sigma_scale )

## test without nullhyp

f0 = addFactor!(fg, [:x0], pthru, graphinit=false)

## test the inference functions

bel, infd = propagateBelief(fg, v0, [f0])

@test isPartial(bel)
@test length(getPoints(bel)) == 120

## repeat test with nullhypo

fg = initfg()

v0 = addVariable!(fg, :x0, SpecialEuclidean2)
# test with nullhypo
f0 = addFactor!(fg, [:x0], pthru, graphinit=false, nullhypo=0.2)

## test the inference functions

bel, infd = propagateBelief(fg, v0, [f0])
@test isPartial(bel)

##

doautoinit!(fg, :x0)

@test length(getPoints(getBelief(fg, :x0))) == getSolverParams(fg).N # 120
# @info "PassThrough transfers the full point count to the graph, unless a product is calculated during the propagateBelief step."

##

solveGraph!(fg);

@test 120 == length(getPoints(fg, :x0))

@warn "must still check if bandwidths are recalculated on many points (not necessary), or lifted from this case single prior"

##

mp = ManifoldPrior(SpecialEuclidean(2), ProductRepr(@MVector([0.0,0.0]), @MMatrix([1.0 0.0; 0.0 1.0])), MvNormal([0.01, 0.01, 0.01],[1 0 0;0 1 0;0 0 1.]))
f1 = addFactor!(fg, [:x0], mp, graphinit=false)

@test length(ls(fg, :x0)) == 2

##

prp, infd = propagateBelief(fg, v0, [f0;f1])

@test length(getPoints(prp)) == getSolverParams(fg).N

## check that solve corrects the point count on graph variable

@test 120 == length(getPoints(fg, :x0))

solveGraph!(fg);

# this number should drop down to usual, 100 at time of writing
@test getSolverParams(fg).N == length(getPoints(fg, :x0))

## check saveDFG (check consistency of packing converters above)

saveDFG("/tmp/passthru", fg)
fg_ = loadDFG("/tmp/passthru.tar.gz")
Base.rm("/tmp/passthru.tar.gz")

##
end

@testset "test point count on propagate, solve, init for PartialPriorPassThrough w Relative" begin
##

fg = initfg()

v0 = addVariable!(fg, :x0, SpecialEuclidean2)

img_ = rand(10,10).+5.0
x_,y_ = ([-9:2.0:9;],[-9:2.0:9;])

hmd = LevelSetGridNormal(img_, (x_,y_), 5.5, 0.1, N=120)
pthru = PartialPriorPassThrough(hmd, (1,2))

# test without nullhyp
f0 = addFactor!(fg, [:x0], pthru, graphinit=false)

## test the inference functions
addVariable!(fg, :x1, SpecialEuclidean2)
mp = ManifoldPrior(SpecialEuclidean(2), ProductRepr(@MVector([0.0,0.0]), @MMatrix([1.0 0.0; 0.0 1.0])), MvNormal([0.01, 0.01, 0.01]))
f1 = addFactor!(fg, [:x1], mp, graphinit=false)

doautoinit!(fg, :x1)

## connect with relative and check calculation size on x0

mf = ManifoldFactor(SpecialEuclidean(2), MvNormal([1,2,pi/4], [0.01,0.01,0.01]))
f2 = addFactor!(fg, [:x0, :x1], mf, graphinit=false)

##

bel, infd = propagateBelief(fg, v0, [f0;f2])

@test !isPartial(bel)
@test getSolverParams(fg).N == length(getPoints(bel))

## check other functions

solveTree!(fg);

@test getSolverParams(fg).N == length(getPoints(fg, :x0))
@test getSolverParams(fg).N == length(getPoints(fg, :x1))

## and check that initAll! works the same (different init sequences may change code execution path)

fg = initfg()
v0 = addVariable!(fg, :x0, SpecialEuclidean2)
img_ = rand(10,10).+5.0
x_,y_ = ([-9:2.0:9;],[-9:2.0:9;])
hmd = LevelSetGridNormal(img_, (x_,y_), 5.5, 0.1, N=120)
pthru = PartialPriorPassThrough(hmd, (1,2))
f0 = addFactor!(fg, [:x0], pthru, graphinit=false)
addVariable!(fg, :x1, SpecialEuclidean2)
mp = ManifoldPrior(SpecialEuclidean(2), ProductRepr(@MVector([0.0,0.0]), @MMatrix([1.0 0.0; 0.0 1.0])), MvNormal([0.01, 0.01, 0.01]))
f1 = addFactor!(fg, [:x1], mp, graphinit=false)
mf = ManifoldFactor(SpecialEuclidean(2), MvNormal([1,2,pi/4], [0.01,0.01,0.01]))
f2 = addFactor!(fg, [:x0, :x1], mf, graphinit=false)

##

bel, infd = propagateBelief(fg, v0, [f0;f2])

@test !isPartial(bel)
@test getSolverParams(fg).N == length(getPoints(bel))

##

initAll!(fg)

@test getSolverParams(fg).N == length(getPoints(fg, :x0))
@test getSolverParams(fg).N == length(getPoints(fg, :x1))

##
end

@testset "Test SpecialEuclidean(2) to TranslationGroup(2) multihypo" begin

##
fg = initfg()
# fg.solverParams.attemptGradients=false

v0 = addVariable!(fg, :x0, SpecialEuclidean2)

mp = ManifoldPrior(SpecialEuclidean(2), ProductRepr(@MVector([0.0,0.0]), @MMatrix([1.0 0.0; 0.0 1.0])), MvNormal([0.01, 0.01, 0.01]))
p = addFactor!(fg, [:x0], mp)

##
addVariable!(fg, :x1a, TranslationGroup2)
addVariable!(fg, :x1b, TranslationGroup2)
mf = ManiPose2Point2(MvNormal([1,2], [0.01,0.01]))
f = addFactor!(fg, [:x0, :x1a, :x1b], mf; multihypo=[1,0.5,0.5])

solveTree!(fg)

vnd = getVariableSolverData(fg, :x0)
@test isapprox(SpecialEuclidean(2), mean(SpecialEuclidean(2), vnd.val), ProductRepr([0.0,0.0], [1.0 0; 0 1]), atol=0.1)

#FIXME I would expect close to 50% of particles to land on the correct place
# Currently software works so that 33% should land there so testing 20 for now
pnt = getPoints(fg, :x1a)
@test sum(isapprox.(pnt, Ref([1.0,2.0]), atol=0.1)) > 20

#FIXME I would expect close to 50% of particles to land on the correct place
pnt = getPoints(fg, :x1b)
@test sum(isapprox.(pnt, Ref([1.0,2.0]), atol=0.1)) > 20

## other way around

fg = initfg()
fg.solverParams.attemptGradients=false

addVariable!(fg, :x0, SpecialEuclidean2)
addVariable!(fg, :x1a, TranslationGroup2)
addVariable!(fg, :x1b, TranslationGroup2)

mp = ManifoldPrior(SpecialEuclidean(2), ProductRepr(@MVector([0.0,0.0]), @MMatrix([1.0 0.0; 0.0 1.0])), MvNormal([10, 10, 0.01]))
p = addFactor!(fg, [:x0], mp)
mp = ManifoldPrior(TranslationGroup(2), [1.,1], MvNormal([0.01, 0.01]))
p = addFactor!(fg, [:x1a], mp)
mp = ManifoldPrior(TranslationGroup(2), [-1.,1], MvNormal([0.01, 0.01]))
p = addFactor!(fg, [:x1b], mp)

mf = ManiPose2Point2(MvNormal([0., 1], [0.01,0.01]))
f = addFactor!(fg, [:x0, :x1a, :x1b], mf; multihypo=[1,0.5,0.5])

solveTree!(fg)

pnts = getPoints(fg, :x0)
# c = getCoordinates.(SpecialEuclidean2, pnts)
# @cast p[i,j] := c[i][j]
# scatter(p[:,1], p[:,2])

#FIXME
@test 10 < sum(isapprox.(Ref(SpecialEuclidean(2)), pnts, Ref(ProductRepr([-1.0,0.0], [1.0 0; 0 1])), atol=0.5))
@test 10 < sum(isapprox.(Ref(SpecialEuclidean(2)), pnts, Ref(ProductRepr([1.0,0.0], [1.0 0; 0 1])), atol=0.5))

end

@testset "Test SpecialEuclidean(2) to SpecialEuclidean(2) multihypo" begin
##
fg = initfg()
# fg.solverParams.attemptGradients=false

v0 = addVariable!(fg, :x0, SpecialEuclidean2)

mp = ManifoldPrior(SpecialEuclidean(2), ProductRepr(@MVector([0.0,0.0]), @MMatrix([1.0 0.0; 0.0 1.0])), MvNormal([0.01, 0.01, 0.01]))
p = addFactor!(fg, [:x0], mp)

##
addVariable!(fg, :x1a, SpecialEuclidean2)
addVariable!(fg, :x1b, SpecialEuclidean2)
mf = ManifoldFactor(SpecialEuclidean(2), MvNormal([1,2,pi/4], [0.01,0.01,0.01]))
f = addFactor!(fg, [:x0, :x1a, :x1b], mf; multihypo=[1,0.5,0.5])

solveTree!(fg)

vnd = getVariableSolverData(fg, :x0)
@test isapprox(SpecialEuclidean(2), mean(SpecialEuclidean(2), vnd.val), ProductRepr([0.0,0.0], [1.0 0; 0 1]), atol=0.1)

#FIXME I would expect close to 50% of particles to land on the correct place
# Currently software works so that 33% should land there so testing 20 for now
pnt = getPoints(fg, :x1a)
@test sum(isapprox.(Ref(SpecialEuclidean(2)), pnt, Ref(ProductRepr([1.0,2.0], [0.7071 -0.7071; 0.7071 0.7071])), atol=0.1)) > 20

#FIXME I would expect close to 50% of particles to land on the correct place
pnt = getPoints(fg, :x1b)
@test sum(isapprox.(Ref(SpecialEuclidean(2)), pnt, Ref(ProductRepr([1.0,2.0], [0.7071 -0.7071; 0.7071 0.7071])), atol=0.1)) > 20

end
#


text:
# @testset "GGPlot(::RObject{VecSxp})" begin
#     p = R"""
#     suppressPackageStartupMessages(library("ggplot2"))
#     df <- data.frame(a = rnorm(100), b = rnorm(100))
#     p <- ggplot(df) +
#         geom_point(aes(a, b))

#     pdf(file = NULL)
#     p
#     """

#     p = Escape.GGPlot(p)
#     @test p isa Escape.GGPlot
# end

# @testset "save(::GGPlot)" begin
#     p = R"""
#     suppressPackageStartupMessages(library("ggplot2"))
#     df <- data.frame(a = rnorm(100), b = rnorm(100))
#     p <- ggplot(df) +
#         geom_point(aes(a, b))

#     pdf(file = NULL)
#     p
#     """

#     p = Escape.GGPlot(p)
#     temp_path = tempname() * "test.pdf"

#     @test !isfile(temp_path)
#     Escape.save(p, temp_path)
#     @test isfile(temp_path)
#     rm(temp_path)
# end


text:
module LiterateTest

export @dedent, @evaltest, @evaltest_throw, @testset_error

import Test
using Base.Meta: isexpr

"""
    LiterateTest.config(; overloads...) -> config::Dict{Symbol,Any}

Construct a `config` dictionary that can be passed to
`Literate.markdown` etc.

`LiterateTest.config` post-composes given `preprocess` function
(default: `identity`) with `LiterateTest.preprocess`.
"""
config(; preprocess = identity, kw...) =
    Dict{Symbol,Any}(:preprocess => (@__MODULE__).preprocess ∘ preprocess, kw...)

nitems(itr) = count(_ -> true, itr)

if isdefined(Iterators, :takewhile)
    using Base.Iterators: takewhile
else
    takewhile(pred, xs) = TakeWhile(pred, xs)
    struct TakeWhile{I,P<:Function}
        pred::P
        xs::I
    end
    function Base.iterate(ibl::TakeWhile, itr...)
        y = iterate(ibl.xs, itr...)
        y === nothing && return nothing
        ibl.pred(y[1]) || return nothing
        y
    end
    Base.IteratorSize(::Type{<:TakeWhile}) = Base.SizeUnknown()
    Base.IteratorEltype(::Type{TakeWhile{I,P}}) where {I,P} = Base.IteratorEltype(I)
end

#! format: off
const INTERPOLATION_NOT_SUPPORTED_ERROR = ArgumentError(
    raw"""
    Interpolation with `$` is not supported.
    Use `raw"..."` to avoid parsing `$` as interpolation.
    """
)
#! format: on

"""
    LiterateTest.preprocess(code::AbstractString) -> code′::String

Remove testing related code from `code`.

See [Testing LiterateTest.jl](@ref tests) for examples of how it works.

!!! note

    Currently, `LiterateTest.preprocess` is implemented as a plain
    line-based transformation.  It may break with complex expressions.
"""
function preprocess(original::AbstractString)
    source = Iterators.Stateful(eachline(IOBuffer(original)))
    io = IOBuffer()
    while true
        ln = try
            popfirst!(source)
        catch
            break
        end
        if ln == "ans = begin"
            println(io, "ans = begin # hide")
            print_deindent_until(io, source) do ln
                ln == "end"
            end
            println(io, "end # hide")
        elseif (m = match(r"^( *)@test +begin$", ln)) !== nothing
            spaces = m[1]
            re = Regex("^" * spaces * raw"end\b")
            print_deindent_until(io, source) do ln
                match(re, ln) !== nothing
            end
        elseif (m = match(r"^( *)@test\b", ln)) !== nothing
            # remove `@test ...`
        elseif startswith(ln, "@testset ")
            consume_until_end(source)
        elseif (
            m = match(
                r"""
                ^(?<macroname>@evaltest|@evaltest_throw)\ +
                (?<str>
                    (?:raw)? (?<left>\"\"\"|\") (?<code>.*?) (?<right>\"\"\"|\")
                )\ + begin$
                """x,
                ln,
            )
        ) !== nothing
            ex = Meta.parse(m[:str])
            if ex isa String
                code = ex
            elseif ex.head === :string
                throw(INTERPOLATION_NOT_SUPPORTED_ERROR)
            else
                @assert ex.head == :macrocall
                @assert ex.args[1] == Symbol("@raw_str")
                code = ex.args[end]
            end
            if m[:macroname] == "@evaltest"
                println(io, code)
                consume_until_end(source)
            elseif m[:macroname] == "@evaltest_throw"
                print(io, THROWING_HEADER)
                println(io, code)
                print(io, THROWING_FOOTER)
                consume_until_end(source)
            end
        elseif (m = match(r"^@testset_error .*try$", ln)) !== nothing
            print(io, THROWING_HEADER)
            # TODO: less manual recursive processing
            inner = sprint() do io
                print_deindent_until(io, source) do ln
                    startswith(ln, "catch ")
                end
            end
            print(io, preprocess(inner))
            print(io, THROWING_FOOTER)
            consume_until_end(source)
        #! format: off
        #=
        elseif (m = match(r"^@throwing *(.*)$", ln)) !== nothing
            print(io, THROWING_HEADER)
            println(io, m[1])
            print(io, THROWING_FOOTER)
        elseif (m = match(r"^@throwing\((.*)\)$", ln)) !== nothing
            print(io, THROWING_HEADER)
            println(io, m[1])
            print(io, THROWING_FOOTER)
        elseif ln == "^@throwing begin"
            print(io, THROWING_HEADER)
            print_deindent_until(io, source) do ln
                startswith(ln, "end ")
            end
            print(io, THROWING_FOOTER)
        =#
        #! format: on
        elseif (m = match(r"^( *)@dedent +([^ ].*)$", ln)) !== nothing
            println(io, m[2], " # hide")
            spaces = m[1]
            re = Regex("^" * spaces * raw"(end\b.*)$")
            mend = nothing
            print_deindent_until(io, source) do ln
                mend = match(re, ln)
                mend !== nothing
            end
            println(io, mend[1], " # hide")
        elseif match(r"^ *#! format: (on|off)$", ln) !== nothing
            # Filter out JuliaFormatter directives:
            # https://domluna.github.io/JuliaFormatter.jl/dev/skipping_formatting/
        else
            println(io, ln)
        end
    end
    return String(take!(io))
end

function consume_until_end(source)
    consume_until(source) do ln
        match(r"^end\b$", ln) !== nothing
    end
end

function consume_until(f, source)
    while true
        ln = popfirst!(source)
        f(ln) && return ln
    end
end

function print_deindent_until(f, io, source)
    indent = typemax(Int)
    while true
        ln = popfirst!(source)
        f(ln) && break
        indent = min(indent, nitems(takewhile(isequal(' '), ln)))
        println(io, remove_global(ln[nextind(ln, indent):end]))
    end
end

function remove_global(ln)
    m = match(r"^global +([^ ]+ =(\s.*)?)$", ln)
    m === nothing || return m[1]
    return ln
end

"""
    LiterateTest.preprocess_juliaformatter(code::AbstractString) -> code′::String

Remove `#! format: off` and `#! format: on`.  Note that
[`LiterateTest.preprocess`](@ref) does the same preprocessing.
"""
function preprocess_juliaformatter(original::AbstractString)
    io = IOBuffer()
    for ln in eachline(IOBuffer(original))
        if match(r"^ *#! format: (on|off)$", ln) === nothing
            println(io, ln)
        end
    end
    return String(take!(io))
end

"""
    @evaltest "\$code" begin \$tests end
    @evaltest raw"\$code" begin \$tests end

Evaluate `code`, assign it to the variable `ans`, and then run `tests`.

This macro is meant to be processed by `LiterateTest.preprocess`.
"""
macro evaltest(code::AbstractString, tests::Expr)
    ex = :($Test.@testset $code begin
        ans = $Base.include_string(@__MODULE__, $(QuoteNode(code)))
        $tests
    end)
    @assert ex.head == :macrocall
    @assert ex.args[2] isa LineNumberNode
    @assert ex.args[4].head === :block
    @assert ex.args[4].args[1] isa LineNumberNode
    ex.args[2] = __source__
    ex.args[4].args[1] = __source__
    return esc(ex)
end

macro evaltest(str::Expr, tests::Expr)
    code = macroexpand(__module__, str)
    code::AbstractString
    return esc(:($(@__MODULE__).@evaltest $code $tests))
end

struct Success
    value::Any
end

const ∉ᵢₛₐ = !isa

function check_exception(ans)
    Test.@test ans ∉ᵢₛₐ Success
end

"""
    @evaltest_throw "\$code" begin \$tests end
    @evaltest_throw raw"\$code" begin \$tests end

Evaluate `code`, assign the exception thrown to the variable `ans`,
and then run `tests`.

This macro is meant to be processed by `LiterateTest.preprocess`.
"""
macro evaltest_throw(code::AbstractString, tests::Expr)
    code_catch = """
    import LiterateTest
    try
        LiterateTest.Success($code)
    catch err
        err
    end
    """
    ex = :($Test.@testset $code begin
        ans = $Base.include_string(@__MODULE__, $(QuoteNode(code_catch)))
        $check_exception(ans)
        $tests
    end)
    @assert ex.head == :macrocall
    @assert ex.args[2] isa LineNumberNode
    @assert ex.args[4].head === :block
    @assert ex.args[4].args[1] isa LineNumberNode
    ex.args[2] = __source__
    ex.args[4].args[1] = __source__
    return esc(ex)
end

macro evaltest_throw(str::Expr, tests::Expr)
    code = macroexpand(__module__, str)
    code::AbstractString
    return esc(:($(@__MODULE__).@evaltest_throw $code $tests))
end

"""
    @testset_error expr

`expr` must contain a `try`-`catch` block.
"""
macro testset_error(expr)
    return esc(_testset_error(__source__, __module__, nothing, expr))
end

macro testset_error(label, expr)
    return esc(_testset_error(__source__, __module__, label, expr))
end

function _testset_error(__source__, __module__, label, expr)
    if !isexpr(expr, :try, 3)
        error("Expect a `try`-`catch` block. Got:\n", expr)
    end
    tryblock, error_symbol, testblock = expr.args
    @gensym err1 err2
    testset_body = quote
        $err2 =
            try
                $tryblock
                nothing
            catch $err1
                Some($err1)
            end
        $error_symbol = something($err2)
        $testblock
    end
    args = [:macrocall, getfield(Test, Symbol("@testset")), __source__]
    if label !== nothing
        push!(args, label)
    end
    push!(args, Expr(:block, __source__, testset_body))
    return Expr(args...)
end

#=
"""
    @throwing ex

Catch and return the exception thrown in `ex`.
"""
macro throwing(ex)
    quote
        try
            $(esc(ex))
            nothing
        catch err
            err
        end
    end
end
=#

const THROWING_HEADER = """
err = try # hide
"""

const THROWING_FOOTER = """
catch _err; _err; end # hide
print(stdout, "ERROR: ") # hide
showerror(stdout, err) # hide
#-
"""

"""
    @dedent begin
        ex
    end

A no-op macro that just evaluates `ex`.

This macro is meant to be processed by `LiterateTest.preprocess`.
"""
macro dedent(ex)
    esc(ex)
end

"""
    LiterateTest.AssertAsTest

A module that exports `Test.@test` as `@assert`.

An intended usage is to import `@assert` before including example file
in test:

```julia
module MyTestExampels
using LiterateTest.AssertAsTest: @assert
include("../examples/file.jl")
end
```

This way, `examples/file.jl` can still use `@assert` for clean
documentation while a better test failure messages are shown when
running the tests.

# Examples
```jldoctest
julia> using LiterateTest.AssertAsTest: @assert

julia> @assert 1 == 1
Test Passed
```
"""
module AssertAsTest
export @assert
using Test: @test
@eval const $(Symbol("@assert")) = $(Symbol("@test"))
end

# Use README as the docstring of the module:
@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), r"^```julia"m => "```jldoctest README")
end LiterateTest

end


text:
type ComplexSpeciesDefinition <: SpeciesDefinition
  fitness::FitnessScheme
  stages::Vector{SpeciesStageDefinition}

  ComplexSpeciesDefinition() =
    new(fitness.scalar())
  ComplexSpeciesDefinition(st::Vector{SpeciesStageDefinition}) =
    new(fitness.scalar(), st)
  ComplexSpeciesDefinition(f::FitnessScheme, st::Vector{SpeciesStageDefinition}) =
    new(f, st)
end

"""
Composes a complex species from a provided definition.
"""
function compose!(c::ComplexSpeciesDefinition)
  stages = Dict{String, SpeciesStage}()
  for st in c.stages
    stages[st.label] = compose!(st)
  end
  Species(stages, c.fitness)
end

"""
Unlike a simple species, a complex species may be composed of numerous
developmental stages, or different views on a particular stage.

For example, in the case of grammatical evolution, one may wish to represent
use different stages to represent the individual as a bit vector, codon
vector, program string, and compiled program.
"""
complex(fitness::FitnessScheme, stages::Vector{SpeciesStageDefinition}) =
  ComplexSpeciesSpecification(fitness, stages)

complex(stages::Vector{SpeciesStageDefinition}) =
  ComplexSpeciesDefinition(stages)

function complex(spec::Function)
  ss = ComplexSpeciesDefinition()
  spec(ss)
  ss
end