text:
"""
 This contains the results from a function evaluation in the FixedPointAcceleration framework.
 It containing the following fields:
 * `Input_` - The input
 * `Output_` - The output of the  function. May be `missing` if function could not complete without error.
 * `Error_` - A symbol representing what error occured.
"""
struct FunctionEvaluationResult{T<:Real,R}
    Input_::Vector{T}
    Output_::Union{Missing,Vector{Missing},Vector{Union{Missing,R}},Vector{R}}
    Other_Output_::Union{Missing,NamedTuple}
    Error_::Symbol
    function FunctionEvaluationResult(Input::Vector{T}, Output_::Missing, Error_::Symbol, Other_Output_::Union{Missing,NamedTuple} = missing) where T<:Real
        return new{T,T}(Input, Output_, Other_Output_, Error_)
    end
    function FunctionEvaluationResult(Input::Vector{T}, Output::Vector{R}, Error_::Symbol, Other_Output_::Union{Missing,NamedTuple} = missing) where T<:Real where R<:Real
        return new{T,R}(Input, convert(Array{Union{Missing,R},1}, Output), Other_Output_, Error_)
    end
    function FunctionEvaluationResult(Input::Vector{T}, Output::Vector{Union{Missing,R}}, Error_::Symbol, Other_Output_::Union{Missing,NamedTuple} = missing) where T<:Real where R<:Real
        return new{T,R}(Input, Output, Other_Output_, Error_)
    end
    function FunctionEvaluationResult(Input::Vector{T}, Output::Vector{Missing}, Error_::Symbol, Other_Output_::Union{Missing,NamedTuple} = missing) where T<:Real
        return new{T,T}(Input, Output, Other_Output_, Error_)
    end
end

"""
 This contains the results from a the fixed_point function.

 It containing the following fields:
 * `FixedPoint_` - The `Vector` with the fixed point that has been found.
 * `Other_Output_` - The other output of the fixedpoint function.
 * `Convergence_` - A real number showing how close the `FixedPoint_` is to the input (to the function) that created it as output.
  * `TerminationCondition_' - Why did the fixedpoint acceleration stop.
  * `Iterations_' - How many iterations were undertaken
  * `ConvergenceVector_` - What is the convergence value at each iteration.
  * `FailedEvaluation_` - Why did the fixedpoint iteration fail (missing if it did not fail)
  * `Inputs_` - What were all of the inputs tried
  * `Outputs_` - What were all the corresponding outputs.
"""
struct FixedPointResults{R<:Real}
    FixedPoint_::Union{Missing,Vector{R}}
    Other_Output_::Union{Missing,NamedTuple}
    Convergence_::Union{Missing,R}          #
    TerminationCondition_::Symbol
    Iterations_::Integer
    ConvergenceVector_::Union{Missing,Array{R,1}}
    FailedEvaluation_::Union{Missing,FunctionEvaluationResult}
    Inputs_::Array{R,2}
    Outputs_::Array{R,2}
    function FixedPointResults(Inputs_::Array{R,2}, Outputs_::Array{R,2}, TerminationCondition_::Symbol;
                               ConvergenceVector_::Union{Missing,Array{R,1}} = missing,
                               FailedEvaluation_::Union{Missing,FunctionEvaluationResult} = missing,
                               Other_Output::Union{Missing,NamedTuple} = missing) where R<:Real
        Iterations_ = size(Outputs_)[2]
        FixedPoint_ = missing
        Convergence_ = missing
        if (!(ismissing(ConvergenceVector_))) && !(isempty(ConvergenceVector_))
            Convergence_ = ConvergenceVector_[Iterations_]
        end
        if TerminationCondition_ == :ReachedConvergenceThreshold
            FixedPoint_ = Outputs_[:,Iterations_]
        end
        return new{R}(FixedPoint_, Other_Output, Convergence_, TerminationCondition_, Iterations_, ConvergenceVector_, FailedEvaluation_, Inputs_, Outputs_)
    end
end


text:
"""
`Crop <: ImageTransformation`

Description
============

Crops out the area of the specified pixel dimensions starting
at a specified position that denotes the top-left corner
of the crop. A position of `x = 1`, and `y = 1` would mean that
the crop is located in the top-left corner of the original image

Usage
======

    Crop(1, 1, 64, 64)

    Crop(x = 1, y = 1, width = 64, height = 64)

Arguments
==========

- **`x`** : The horizontal offset of the top left corner of the
window that should be cropped out

- **`y`** : The vertical offset of the top left corner of the
window that should be cropped out

- **`width`** : The desired width or the cropped out sub-image
in pixels

- **`height`** : The desired height or the cropped out sub-image
in pixels

Methods
========

- **`transform`** : Applies the transformation to the given Image
or set of images

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

Examples
========

    using Augmentor
    using TestImages

    # load an example image
    img = testimage("lena")

    # Crop out a 64x32 area starting at x=10, y=20
    img_new = transform(Crop(10, 20, 64, 32), img)

see also
=========

`ImageTransformation`, `transform`
"""
immutable Crop <: ImageTransformation
    x::Int
    y::Int
    width::Int
    height::Int

    function Crop(x::Int, y::Int, width::Int, height::Int)
        x > 0 || throw(ArgumentError("x has to be greater than 0"))
        y > 0 || throw(ArgumentError("y has to be greater than 0"))
        width > 0 || throw(ArgumentError("width has to be greater than 0"))
        height > 0 || throw(ArgumentError("height has to be greater than 0"))
        new(x, y, width, height)
    end
end

Crop(; x::Int = 1, y::Int = 1, width::Int = 64, height::Int = 64) = Crop(x, y, width, height)

Base.showcompact(io::IO, op::Crop) = print(io, "Crop at x=$(op.x) y=$(op.y) a $(op.width)x$(op.height) window.")
multiplier(::Crop) = 1

function transform{T<:AbstractImage}(op::Crop, img::T)
    w = width(img)
    h = height(img)
    op.x <= w || throw(ArgumentError("x of crop window is out of bounds"))
    op.y <= h || throw(ArgumentError("y of crop window is out of bounds"))

    x_end = op.x + op.width - 1
    y_end = op.y + op.height - 1
    x_end <= w || throw(ArgumentError("width of crop window is out of bounds"))
    y_end <= h || throw(ArgumentError("height of crop window is out of bounds"))

    result = crop(img, op.x:x_end, op.y:y_end)
    _log_operation!(result, op)
end

ext: jl


text:
using LoopVectorization, OffsetArrays
using Test

@testset "dot" begin
    dotq = :(for i ∈ eachindex(a,b)
             s += a[i]*b[i]
             end)
    lsdot = LoopVectorization.loopset(dotq);
    @test LoopVectorization.choose_order(lsdot) == (Symbol[:i], :i, Symbol("##undefined##"), :i, 4, -1)
    function mydot(a::AbstractVector, b::AbstractVector)
        s = zero(eltype(a))
        za = OffsetArray(a, OffsetArrays.Origin(0))
        zb = OffsetArray(b, OffsetArrays.Origin(0))
        @inbounds @simd for i ∈ LoopVectorization.CloseOpen(min(length(a),length(b)))
            s += za[i]*zb[i]
        end
        s
    end
    function mydotavx(a::AbstractVector, b::AbstractVector)
        s = zero(eltype(a))
        za = OffsetArray(a, OffsetArrays.Origin(0))
        zb = OffsetArray(b, OffsetArrays.Origin(0))
        @turbo for i ∈ LoopVectorization.CloseOpen(min(length(a),length(b)))
            s += za[i]*zb[i]
        end
        s
    end
    @test LoopVectorization.ArrayInterface.static_step(LoopVectorization.CloseOpen(-5,10)) === LoopVectorization.One()
    function mydot(a, b)
        s = zero(eltype(a))
        @inbounds @simd for i ∈ eachindex(a,b)
            s += a[i]*b[i]
        end
        s
    end
    function mydotavx(a, b)
        s = zero(eltype(a))
        @turbo for i ∈ eachindex(a,b)
            s += a[i]*b[i]
        end
        s
    end

    function mydot_avx(a, b)
        s = zero(eltype(a))
        @_avx for i ∈ eachindex(a,b)
            aᵢ = getindex(a, i)
            s += aᵢ*b[i]
        end
        s
    end

    selfdotq = :(for i ∈ indices(a,1)
                 s += a[i]*a[i]
                 end)
    lsselfdot = LoopVectorization.loopset(selfdotq);
    @test LoopVectorization.choose_order(lsselfdot) == (Symbol[:i], :i, Symbol("##undefined##"), :i, 8, -1)

    function myselfdot(a)
        s = zero(eltype(a))
        @inbounds @simd for i ∈ eachindex(a)
            s += getindex(a, i) * a[i]
        end
        s
    end
    function myselfdotavx(a)
        s = zero(eltype(a))
        @turbo for i ∈ eachindex(a)
            s += a[i]*a[i]
        end
        s
    end
    function myselfdotavx_v2(a)
        s = zero(eltype(a))
        @turbo for i ∈ 1:length(a)
            s += a[i]*a[i]
        end
        s
    end
    function myselfdotavx_range(a)
        s = zero(eltype(a))
        rng = axes(a, 1)
        @turbo for i ∈ rng
            s += a[i]*a[i]
        end
        s
    end
    function myselfdot_avx(a)
        s = zero(eltype(a))
        @_avx for i ∈ eachindex(a)
            s += a[i]*a[i]
        end
        s
    end
    function dot_unroll2avx(x::Vector{T}, y::Vector{T}) where {T<:Number}
        z = zero(T)
        @turbo unroll=2 for i ∈ 1:length(x)
            z += x[i]*y[i]
        end
        z
    end
    function dot_unroll3avx(x::Vector{T}, y::Vector{T}) where {T<:Number}
        z = zero(T)
        @turbo unroll=3 for i ∈ 1:length(x)
            z += x[i]*y[i]
        end
        z
    end
    # @macroexpand @turbo inline=false unroll=2 for i ∈ 1:length(x)
    #         z += x[i]*y[i]
    #     end

    function dot_unroll2avx_noinline(x::Vector{T}, y::Vector{T}) where {T<:Number}
        z = zero(T)
        @turbo inline=true unroll=2 for i ∈ 1:length(x)
            z += x[i]*y[i]
        end
        z
    end
    function dot_unroll3avx_inline(x::Vector{T}, y::Vector{T}) where {T<:Number}
        z = zero(T)
        @turbo unroll=3 inline=true check_empty=true for i ∈ 1:length(x)
            z += x[i]*y[i]
        end
        z
    end
    function dot_unroll2_avx(x::Vector{T}, y::Vector{T}) where {T<:Number}
        z = zero(T)
        @_avx unroll=2 for i ∈ 1:length(x)
            z += x[i]*y[i]
        end
        z
    end
    function dot_unroll3_avx(x::Vector{T}, y::Vector{T}) where {T<:Number}
        z = zero(T)
        @_avx unroll=3 for i ∈ 1:length(x)
            z += x[i]*y[i]
        end
        z
    end
    function complex_dot_soa(
        xre::AbstractVector{T}, xim::AbstractVector{T},
        yre::AbstractVector{T}, yim::AbstractVector{T}
    ) where {T}
        zre = zero(T)
        zim = zero(T)
        @turbo for i ∈ 1:length(xre)
            zre += xre[i]*yre[i] - xim[i]*yim[i]
            zim += xre[i]*yim[i] + xim[i]*yre[i]
        end
        Complex{T}(zre,zim)
    end
    qc = :(for i ∈ 1:length(xre)
           zre += xre[i]*yre[i] - xim[i]*yim[i]
           zim += xre[i]*yim[i] + xim[i]*yre[i]
           end);
    lsc = LoopVectorization.loopset(qc);
    function complex_mul_with_index_offset!(c_re, c_im, a_re, a_im, b_re, b_im)
        @inbounds @simd ivdep for i = 1:length(a_re) - 1
            c_re[i] = b_re[i] * a_re[i + 1] - b_im[i] * a_im[i + 1]
            c_im[i] = b_re[i] * a_im[i + 1] + b_im[i] * a_re[i + 1]
        end
    end
    function complex_mul_with_index_offsetavx!(c_re, c_im, a_re, a_im, b_re, b_im)
        @turbo for i = 1:length(a_re) - 1
            c_re[i] = b_re[i] * a_re[i + 1] - b_im[i] * a_im[i + 1]
            c_im[i] = b_re[i] * a_im[i + 1] + b_im[i] * a_re[i + 1]
        end
    end
    function complex_mul_with_index_offset_avx!(c_re, c_im, a_re, a_im, b_re, b_im)
        @_avx for i = 1:length(a_re) - 1
            c_re[i] = b_re[i] * a_re[i + 1] - b_im[i] * a_im[i + 1]
            setindex!(c_im, b_re[i] * a_im[i + 1] + b_im[i] * a_re[i + 1], i)
        end
    end
    # q = :(for i = 1:length(a_re) - 1
    #         c_re[i] = b_re[i] * a_re[i + 1] - b_im[i] * a_im[i + 1]
    #         c_im[i] = b_re[i] * a_im[i + 1] + b_im[i] * a_re[i + 1]
    #       end);
    # ls = LoopVectorization.loopset(q)

    function mcpi(x, y)
        acc = 0
        @inbounds @simd for i ∈ eachindex(x)
            acc += (x[i]*x[i] + y[i]*y[i]) < 1.0
        end
        4acc/length(x)
    end
    function mcpiavx(x, y)
        acc = 0
        @turbo for i ∈ eachindex(x)
            acc += (x[i]*x[i] + y[i]*y[i]) < 1.0
        end
        4acc/length(x)
    end
    function mcpiavx_u4(x, y)
        acc = 0
        @turbo unroll=4 for i ∈ eachindex(x)
            acc += (x[i]*x[i] + y[i]*y[i]) < 1.0
        end
        4acc/length(x)
    end
    function mcpi_avx(x, y)
        acc = 0
        @_avx for i ∈ eachindex(x)
            acc += (x[i]*x[i] + y[i]*y[i]) < 1.0
        end
        4acc/length(x)
    end
    function mcpi_avx_u4(x, y)
        acc = 0
        @_avx unroll=4 for i ∈ 1:length(x)
            acc += (x[i]*x[i] + y[i]*y[i]) < 1.0
        end
        4acc/length(x)
    end

    function dotloopinductvarpow(x::AbstractArray{T}) where {T}
        s = zero(T)
        for i ∈ eachindex(x)
            s += x[i] * T(i)^3
        end
        s
    end
    function dotloopinductvarpowavx(x)
        s = zero(eltype(x))
        @turbo for i ∈ eachindex(x)
            s += x[i] * i^3
        end
        s
    end
    function dot_from_n_to_100(a, b, n)
        s = zero(eltype(a))
        @turbo for i ∈ n:100
            s += a[i] * b[i]
        end
        s
    end
    function dot33(a,b)
        s = zero(eltype(a))
        @turbo for i ∈ 1:33
            s += a[i] * b[i]
        end
        s
    end
    function dot17(a,b)
        s = zero(eltype(a))
        @turbo for i ∈ 1:17
            s += a[i] * b[i]
        end
        s
    end
    # @macroexpand @_avx for i = 1:length(a_re) - 1
    #     c_re[i] = b_re[i] * a_re[i + 1] - b_im[i] * a_im[i + 1]
    #     c_im[i] = b_re[i] * a_im[i + 1] + b_im[i] * a_re[i + 1]
    # end

    # a = rand(400);
    for T ∈ (Float32, Float64, Int32, Int64)
        @show T, @__LINE__
        N = 143
        R = T <: Integer ? (T(-100):T(100)) : T
        a = rand(R, N); b = rand(R, N);
        ao = OffsetArray(a, -60:N-61); bo = OffsetArray(b, -60:N-61);
        s = mydot(a, b)
        @test mydotavx(a,b) ≈ s
        @test mydot_avx(a,b) ≈ s
        @test mydotavx(ao,bo) ≈ s
        @test mydot_avx(ao,bo) ≈ s
        @test dot_unroll2avx(a,b) ≈ s
        @test dot_unroll3avx(a,b) ≈ s
        @test dot_unroll2_avx(a,b) ≈ s
        @test dot_unroll3_avx(a,b) ≈ s
        @test dot_unroll2avx_noinline(a,b) ≈ s
        @test dot_unroll3avx_inline(a,b) ≈ s
        s = myselfdot(a)
        @test myselfdotavx(a) ≈ s
        @test myselfdotavx_v2(a) ≈ s
        @test myselfdotavx_range(a) ≈ s
        @test myselfdot_avx(a) ≈ s
        @test myselfdotavx(a) ≈ s

        A = OffsetArray(rand(37, 61), -5, 10);
        s = myselfdot(A);
        @test myselfdotavx(A) ≈ myselfdotavx(A') ≈ s
        @test myselfdotavx_v2(A) ≈ myselfdotavx_v2(A') ≈ s
        # @test myselfdot_avx(A) ≈ myselfdot_avx(A') ≈ s

        @test dot17(a,b) ≈ @view(a[1:17])' * @view(b[1:17])
        @test dot33(a,b) ≈ @view(a[1:33])' * @view(b[1:33])

        if T <: Union{Float32,Float64}
            πest = mcpi(a, b)
            @test πest ≈ mcpiavx(a, b)
            @test πest ≈ mcpiavx_u4(a, b)
            @test πest ≈ mcpi_avx(a, b)
            @test πest ≈ mcpi_avx_u4(a, b)
        end

        if !(!Bool(LoopVectorization.VectorizationBase.has_feature(Val(:x86_64_avx2))) && T === Int32)
            @test dotloopinductvarpow(a) ≈ dotloopinductvarpowavx(a)
        end
        @test dot_from_n_to_100(a, b, 33) == @views mydotavx(a[33:100], b[33:100])

        a_re = rand(R, N); a_im = rand(R, N);
        b_re = rand(R, N); b_im = rand(R, N);
        ac = Complex.(a_re, a_im);
        bc = Complex.(b_re, b_im);

        @test mydot(ac, bc) ≈ complex_dot_soa(a_re, a_im, b_re, b_im)

        c_re1 = similar(a_re); c_im1 = similar(a_im);
        c_re2 = similar(a_re); c_im2 = similar(a_im);
        # b_re = rand(R, length(a_re) + 1); b_im = rand(R, length(a_im) + 1);
        complex_mul_with_index_offset!(c_re1, c_im1, a_re, a_im, b_re, b_im)
        complex_mul_with_index_offsetavx!(c_re2, c_im2, a_re, a_im, b_re, b_im)
        c_re1v, c_im1v, c_re2v, c_im2v = @views c_re1[1:end-1], c_im1[1:end-1], c_re2[1:end-1], c_im2[1:end-1];
        @test c_re1v ≈ c_re2v
        @test c_im1v ≈ c_im2v
        fill!(c_re2, -999999); fill!(c_im2, 99999999);
        complex_mul_with_index_offset_avx!(c_re2, c_im2, a_re, a_im, b_re, b_im)
        @test c_re1v ≈ c_re2v
        @test c_im1v ≈ c_im2v

        # Float32 is not accurate enough
        # Alternatively, loosen approx requirement?
        R == Float32 && continue
        A = rand(R, N, N, N);
        B = rand(R, N, N, N);
        @test mydot(A, B) ≈ mydotavx(A, B)
        # test CartesianIndices
        for i ∈ [3, :, 1:N-1], j ∈ [5, :, 1:N-2], k ∈ [:, 1:N-3]
            Av = view(A, i, j, k);
            Bv = view(B, i, j, k);
            # @show i, j, k
            @test mydot(Av, Bv) ≈ mydotavx(Av, Bv)
        end
    end
end


text:
using LinearAlgebra

mutable struct Bird <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    vel::NTuple{2,Float64}
    speed::Float64
    cohere_factor::Float64
    separation::Float64
    separate_factor::Float64
    match_factor::Float64
    visual_distance::Float64
end

"""
```julia
flocking(;
    n_birds = 100,
    speed = 1.0,
    cohere_factor = 0.25,
    separation = 4.0,
    separate_factor = 0.25,
    match_factor = 0.01,
    visual_distance = 5.0,
    dims = (100, 100),
)
```
Same as in [Flock model](@ref).
"""
function flocking(;
        n_birds = 100,
        speed = 1.0,
        cohere_factor = 0.25,
        separation = 4.0,
        separate_factor = 0.25,
        match_factor = 0.01,
        visual_distance = 5.0,
        dims = (100, 100),
    )

    space2d = ContinuousSpace(2; periodic = true, extend = dims)
    model = ABM(Bird, space2d, scheduler = random_activation)
    for _ in 1:n_birds
        vel = Tuple(rand(2) * 2 .- 1)
        add_agent!(
            model,
            vel,
            speed,
            cohere_factor,
            separation,
            separate_factor,
            match_factor,
            visual_distance,
        )
    end
    index!(model)
    return model, flocking_agent_step!, dummystep
end

function flocking_agent_step!(bird, model)
    ## Obtain the ids of neighbors within the bird's visual distance
    ids = nearby_ids(bird, model, bird.visual_distance)
    ## Compute velocity based on rules defined above
    bird.vel =
        (
            bird.vel .+ cohere(bird, model, ids) .+ separate(bird, model, ids) .+
            match(bird, model, ids)
        ) ./ 2
    bird.vel = bird.vel ./ norm(bird.vel)
    ## Move bird according to new velocity and speed
    move_agent!(bird, model, bird.speed)
end

distance(a1, a2) = sqrt(sum((a1.pos .- a2.pos) .^ 2))

get_heading(a1, a2) = a1.pos .- a2.pos

function cohere(bird, model, ids)
    N = max(length(ids), 1)
    birds = model.agents
    coherence = (0.0, 0.0)
    for id in ids
        coherence = coherence .+ get_heading(birds[id], bird)
    end
    return coherence ./ N .* bird.cohere_factor
end

function separate(bird, model, ids)
    seperation_vec = (0.0, 0.0)
    N = max(length(ids), 1)
    birds = model.agents
    for id in ids
        neighbor = birds[id]
        if distance(bird, neighbor) < bird.separation
            seperation_vec = seperation_vec .- get_heading(neighbor, bird)
        end
    end
    return seperation_vec ./ N .* bird.separate_factor
end

function match(bird, model, ids)
    match_vector = (0.0, 0.0)
    N = max(length(ids), 1)
    birds = model.agents
    for id in ids
        match_vector = match_vector .+ birds[id].vel
    end
    return match_vector ./ N .* bird.match_factor
end


text:
using Documenter, KroneckerTools

makedocs(sitename="KroneckerTools.jl")

deploydocs(
    repo = "github.com/DynareJulia/KroneckerTools.jl.git",
    devbranch = "main",
)


text:
export AbstractTrajectory, SART, SARTS, SARTSA, SLART, SLARTSL, SLARTSLA

"""
    AbstractTrajectory

A trajectory is used to record some useful information
during the interactions between agents and environments.
It behaves similar to a `NamedTuple` except that we extend it
with some optional methods.

Required Methods:

- `Base.getindex`
- `Base.keys`

Optional Methods:

- `Base.length`
- `Base.isempty`
- `Base.empty!`
- `Base.haskey`
- `Base.push!`
- `Base.pop!`
"""
abstract type AbstractTrajectory end

Base.haskey(t::AbstractTrajectory, s::Symbol) = s in keys(t)
Base.isempty(t::AbstractTrajectory) = all(k -> isempty(t[k]), keys(t))

function Base.empty!(t::AbstractTrajectory)
    for k in keys(t)
        empty!(t[k])
    end
end

function Base.push!(t::AbstractTrajectory; kwargs...)
    for (k, v) in kwargs
        push!(t[k], v)
    end
end

function Base.pop!(t::AbstractTrajectory)
    for k in keys(t)
        pop!(t[k])
    end
end

function Base.show(io::IO, t::AbstractTrajectory)
    println(io, "Trajectory of $(length(keys(t))) traces:")
    for k in keys(t)
        show(io, k)
        println(io, " $(summary(t[k]))")
    end
end

#####
# Common Keys
#####

const SART = (:state, :action, :reward, :terminal)
const SARTS = (:state, :action, :reward, :terminal, :next_state)
const SARTSA = (:state, :action, :reward, :terminal, :next_state, :next_action)
const SLART = (:state, :legal_actions_mask, :action, :reward, :terminal, :next_state)
const SLARTSL = (
    :state,
    :legal_actions_mask,
    :action,
    :reward,
    :terminal,
    :next_state,
    :next_legal_actions_mask,
)
const SLARTSLA = (
    :state,
    :legal_actions_mask,
    :action,
    :reward,
    :terminal,
    :next_state,
    :next_legal_actions_mask,
    :next_action,
)


text:
#Does a rough timing test of the bootstrapping calculation for the LCX Set

include("../helpers.jl")

function f1(boot_sample::Matrix, data::Matrix, numSamples::Int, a::Vector)
    Gamma::Float64 = 0.0
    for i = 1:numSamples
        randn!(a)
        a /= norm(a)
        const as = data * a
        const as_sample = boot_sample * a
        for b in as
            const lhs::Float64 = mean(max(as-b, 0)) - mean(max(as_sample-b, 0))   #bottleneck line
            if Gamma < lhs
                Gamma = lhs
            end
        end
    end
    Gamma
end

#Returns indx of last instance of val, assumes sorted
#0 if not a member
#Must have vector[start] >= val
function findlast_sort(val::Float64, vector::Vector; TOL::Float64=1e-10, start::Int64=1)
    i::Int64 = 1
    for i = int(start):length(vector)
        if abs(val - vector[i]) > TOL
            break
        end
    end
    i-1
end

function singlepass!(zetas::Vector, zetahats::Vector)
    sort!(zetas)
    sort!(zetahats)

    vstar::Float64 = mean(zetas) - zetas[1]
    vb::Float64    = mean(zetahats) - zetas[1]

    Gamma::Float64 = vstar - vb
    const N::Int64 = length(zetas)
    pbar::Float64    = 1.0
    hat_indx::Int64  = 1
    hat_indx_::Int64 = 0

    for k = 2:length(zetas)
        vstar += (zetas[k-1] - zetas[k]) * (N-k+1)/N
        hat_indx = findlast_sort(zetas[k-1], zetahats, start=hat_indx_ + 1)
        pbar -=  (hat_indx - hat_indx_)/N
        hat_indx_ = hat_indx
        vb  += (zetas[k-1] - zetas[k]) * pbar
        Gamma = max(Gamma, vstar - vb)
    end
    Gamma
end

#a::Vector is used as storage between bootstraps
function f2(boot_sample::Matrix, data::Matrix, numSamples::Int, a::Vector)
    Gamma::Float64 = 0.
    for i = 1:numSamples
        randn!(a)
        a /= norm(a)
        Gamma = max(Gamma, singlepass!(data*a, boot_sample*a))
    end
    Gamma
end

#Approximates the threshold by sampling a bunch of abs for each bootstrap rep
#VG Could be optimized
function calc_ab_thresh(data::Matrix, delta::Float64, numBoots::Int, numSamples::Int)
    const N = size(data, 1)
    const d = size(data, 2)
    a::Vector{Float64} = zeros(Float64, d)
    boot(data, f2, 1-delta, numBoots, data, numSamples, a)
end

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

data = randn(int(ARGS[1]), 2)
calc_ab_thresh(data, .1, 10, 10)  # one for the money!

Profile.init(10^6, .01)
tic()
@profile calc_ab_thresh(data, .1, 10, 1000)  #two for the show
println("Elapsed Time:\t", toc())
Profile.print()

#Baseline 47 s


text:
function Base.push!(A::AxisVector, item)
    StaticRanges.can_set_last(axes(A, 1)) || throw(MethodError(push!, (A, item)))
    push!(parent(A), item)
    grow_last!(axes(A, 1), 1)
    return A
end

function Base.push!(A::AxisVector, item::Pair)
    axis = axes(A, 1)
    StaticRanges.can_set_last(axis) || throw(MethodError(push!, (A, item)))
    push!(parent(A), last(item))
    push_key!(axis, first(item))
    return A
end

function Base.pushfirst!(A::AxisVector, item)
    can_change_size(A) || throw(MethodError(pushfirst!, (A, item)))
    pushfirst_axis!(axes(A, 1))
    pushfirst!(parent(A), item)
    return A
end

function pushfirst_axis!(axis::Axis, key)
    pushfirst!(keys(axis), key)
    grow_last!(parent(axis), 1)
    return nothing
end

function pushfirst_axis!(axis::AbstractAxis)
    grow_last!(parent(axis), 1)
    return nothing
end

function pushfirst_axis!(axis::Axis)
    grow_first!(keys(axis), 1)
    grow_last!(parent(axis), 1)
    return nothing
end

function Base.pushfirst!(A::AxisVector, item::Pair)
    can_change_size(A) || throw(MethodError(pushfirst!, (A, item)))
    axis = axes(A, 1)
    pushfirst_axis!(axis, first(item))
    pushfirst!(parent(A), last(item))
    return A
end

function Base.mapslices(f, a::AxisArray; dims, kwargs...)
    return reconstruct_reduction(a, Base.mapslices(f, parent(a); dims=dims, kwargs...), dims)
end

macro def_equals(f, X,Y)
    if X === :AxisArray
        if Y === :AxisArray
            esc(quote
                function Base.$f(x::$X, y::$Y)
                    if is_dense_wrapper(x) && is_dense_wrapper(y)
                        return Base.$f(parent(x), parent(y))
                    else
                        for (x_i,y_i) in zip(x,y)
                            Base.$f(x_i, y_i) || return false
                        end
                        return true
                    end
                end
            end)
        else
            esc(quote
                function Base.$f(x::$X, y::$Y)
                    if is_dense_wrapper(x)
                        return Base.$f(parent(x), y)
                    else
                        for (x_i,y_i) in zip(x,y)
                            Base.$f(x_i, y_i) || return false
                        end
                        return true
                    end
                end
            end)
        end
    else
        esc(quote
            function Base.$f(x::$X, y::$Y)
                if is_dense_wrapper(y)
                    return Base.$f(x, parent(y))
                else
                    for (x_i,y_i) in zip(x,y)
                        Base.$f(x_i, y_i) || return false
                    end
                    return true
                end
            end
        end)
    end
end

@def_equals(==, AxisArray, AxisArray)
@def_equals(==, AbstractArray, AxisArray)
@def_equals(==, AxisArray, AbstractArray)
@def_equals(==, AxisArray, AbstractAxis)
@def_equals(==, AbstractAxis, AxisArray)
@def_equals(==, AxisArray, GapRange)
@def_equals(==, GapRange, AxisArray)

@def_equals(isequal, AxisArray, AxisArray)
@def_equals(isequal, AbstractArray, AxisArray)
@def_equals(isequal, AxisArray, AbstractArray)
@def_equals(isequal, AxisArray, AbstractAxis)
@def_equals(isequal, AbstractAxis, AxisArray)

Base.isapprox(a::AxisArray, b::AxisArray; kw...) = isapprox(parent(a), parent(b); kw...)
Base.isapprox(a::AxisArray, b::AbstractArray; kw...) = isapprox(parent(a), b; kw...)
Base.isapprox(a::AbstractArray, b::AxisArray; kw...) = isapprox(a, parent(b); kw...)

Base.copy(A::AxisArray) = AxisArray(copy(parent(A)), map(copy, axes(A)))

for (tf, T, sf, S) in (
    (parent, :AxisArray, parent, :AxisArray),
    (parent, :AxisArray, identity, :AbstractArray),
    (identity, :AbstractArray, parent, :AxisArray))

    @eval function Base.cat(A::$T, B::$S, Cs::AbstractArray...; dims)
        p = cat($tf(A), $sf(B); dims=dims)
        return cat(AxisArray(p, cat_axes(A, B, p, dims)), Cs..., dims=dims)
    end
end

for (tf, T, sf, S) in (
    (parent, :AxisVecOrMat, parent, :AxisVecOrMat),
    (parent, :AxisArray, identity, :VecOrMat),
    (identity, :VecOrMat, parent, :AxisArray))
    @eval function Base.vcat(A::$T, B::$S, Cs::VecOrMat...)
        p = vcat($tf(A), $sf(B))
        return vcat(initialize_axis_array(p, vcat_axes(A, B, p)), Cs...)
    end

    @eval function Base.hcat(A::$T, B::$S, Cs::VecOrMat...)
        p = hcat($tf(A), $sf(B))
        return hcat(initialize_axis_array(p, hcat_axes(A, B, p)), Cs...)
    end
end

function Base.hcat(A::AxisArray{T,N}) where {T,N}
    if N === 1
        return AxisArray(hcat(parent(A)), (axes(A, 1), axes(A, 2)))
    else
        return A
    end
end

Base.vcat(A::AxisArray{T,N}) where {T,N} = A

Base.cat(A::AxisArray{T,N}; dims) where {T,N} = A


text:
include("testsystem.jl")
Random.seed!(11)
Y = Any[]
x = rand(StateObs(Gaussian(x0, P0), M.obs))
u = 1 => x[1]
Y = [1 => x[2]]
Y = trajectory(collecttrace(LinearObservation(LinearEvolution(Φ, Gaussian(b, Q)), H, R), Trace(Sample(u), Y, endtime(10))))
@show Y

Random.seed!(11)
x = rand(StateObs(Gaussian(x0, P0), M.obs))
X = trace(DynamicIterators.Sampled(M), 1 => x, endtime(10))

@test Trajectory(X.t, last.(X.x)) == Y


text:
# Tests to ensure datasets are available and correct

# Basic biosemi data file
data_path = joinpath(datadep"BioSemiTestFiles", "Newtest17-2048.bdf")
s = read_SSR(data_path)
@test samplingrate(s) == 2048
@test length(channelnames(s)) == 16
@test length(s.triggers["Index"]) == length(s.triggers["Code"])
@test length(s.triggers["Code"]) == length(s.triggers["Duration"])


text:
module DyestuffTest

println()
@info "Dyestuff example"

using GLMCopula, RDatasets, Test, Gurobi

# Dataframe with columns: Batch (Categorical), Yield (Int32)
dyestuff = dataset("lme4", "Dyestuff")
groups = unique(dyestuff[!, :Batch])
n, p, m = length(groups), 1, 1
gcs = Vector{GaussianCopulaVCObs{Float64}}(undef, n)
for (i, grp) in enumerate(groups)
    gidx = dyestuff[!, :Batch] .== grp
    ni = count(gidx)
    y = Float64.(dyestuff[gidx, :Yield])
    X = ones(ni, 1)
    V = [ones(ni, ni)]
    gcs[i] = GaussianCopulaVCObs(y, X, V)
end
gcm = GaussianCopulaVCModel(gcs);

# initialize β and τ from least square solution
@info "Initial point:"
init_β!(gcm)
@show gcm.β
# update σ2 and τ from β using the MM algorithm
fill!(gcm.Σ, 1)
# update_Σ!(gcm, 500, 1e-6, GurobiSolver(OutputFlag=0), true)
update_Σ!(gcm)
@show gcm.τ
@show gcm.Σ;
@test loglikelihood!(gcm, true, false) ≈ -164.00082379
@show gcm.∇β
@show gcm.∇τ
@show gcm.∇Σ

# fit model using NLP on profiled loglikelihood
@info "MLE:"
@time GLMCopula.fit!(gcm, IpoptSolver(print_level=5))
@show gcm.β
@show gcm.τ
@show gcm.Σ
@test loglikelihood!(gcm, true, false) ≈ -163.35545251
@show gcm.∇β
@show gcm.∇τ
@show gcm.∇Σ

end


text:
using PyPlot
using DelimitedFiles

R = readdlm("result.txt")
n = sqrt.([1, 4, 9, 16, 25])

close("all")
plot(n, R[1:5, 3], "or-", label = "1 Core")
plot(n[1:4], R[6:9, 3], "xg-", label = "2 Cores")
plot(n[1:3], R[10:12, 3], "^b-", label = "3 Cores")
plot(n, R[1:5, 5], "or--")
plot(n[1:4], R[6:9, 5], "xg--")
plot(n[1:3], R[10:12, 5], "^b--")
legend()
title("Weak Scaling for Forward Computation")
xlabel("\$N\$")
ylabel("Time (sec)")
savefig("forward_weak.png")

close("all")
plot(n, R[1:5, 4], "or-", label = "1 Core")
plot(n[1:4], R[6:9, 4], "xg-", label = "2 Cores")
plot(n[1:3], R[10:12, 4], "^b-", label = "3 Cores")
plot(n, R[1:5, 6], "or--")
plot(n[1:4], R[6:9, 6], "xg--")
plot(n[1:3], R[10:12, 6], "^b--")
legend()
title("Strong Scaling for Gradient Backpropagation")
xlabel("\$N\$")
ylabel("Time (sec)")
savefig("backward_weak.png")


text:
function ellipse_boundary_points(xc::SVector{2,T},
    theta::Real,
    a,
    b;
    N::Int = 200) where T
    # boday frame points
    thetaB = LinRange(0,2pi, N)
    X=a*cos.(thetaB)
    Y=b*sin.(thetaB)
    x=similar(X)
    y=similar(X)
    for i in 1:length(X)
        xy=body2lab(SVector{2,T}(X[i], Y[i]), xc, theta)
        x[i], y[i]=xy
    end
    return x,y
end


text:
@testset "GaloisField" begin

  for p in [31, 11, 101]
    F = GF(p)
    G, mG = unit_group(F)
    #Test generator
    g = mG(G[1])
    g1 = g
    for i = 2:p-2
      g1 *= g
      @test !isone(g1)
    end
    for i = 1:5
      x = rand(F)
      while iszero(x)
        x = rand(F)
      end
      elG = mG\x
      @test mG(elG) == x
    end

    G1, mG1 = unit_group(F, n_quo = 5)
    @test order(G1) == 5
    g = mG1(G1[1])
    g = g^(div(size(F)-1, 5))
    @test !isone(g)
    for i = 1:5
      x = rand(F)
      while iszero(x)
        x = rand(F)
      end
      c = mG1\x
      el = mG1(c)
      if el != x
        dl = mG\(el*inv(x))
        @test iszero(mod(dl[1], 5))
      end
    end
  end
end

@testset "GaloisFmpzField" begin

  for p in [31, 11, 101]
    F = GF(fmpz(p))
    G, mG = unit_group(F)
    #Test generator
    g = mG(G[1])
    g1 = g
    for i = 2:p-2
      g1 *= g
      @test !isone(g1)
    end
    for i = 1:5
      x = rand(F)
      while iszero(x)
        x = rand(F)
      end
      elG = mG\x
      @test mG(elG) == x
    end

    G1, mG1 = unit_group(F, n_quo = 5)
    @test order(G1) == 5

    g = mG1(G1[1])
    g = g^(div(size(F)-1, 5))
    @test !isone(g)
    for i = 1:5
      x = rand(F)
      while iszero(x)
        x = rand(F)
      end
      c = mG1\x
      el = mG1(c)
      if el != x
        dl = mG\(el*inv(x))
        @test iszero(mod(dl[1], 5))
      end
    end
  end
end

@testset "FqNmodFiniteField" begin

  for p in [31, 11, 101]
    F = FiniteField(p, 2)[1]
    G, mG = unit_group(F)
    #Test generator
    g = mG(G[1])
    g1 = g
    for i = 2:p^2-2
      g1 *= g
      @test !isone(g1)
    end
    for i = 1:5
      x = rand(F)
      while iszero(x)
        x = rand(F)
      end
      elG = mG\x
      @test mG(elG) == x
    end

    G1, mG1 = unit_group(F, n_quo = 5)
    @test order(G1) == 5
    g = mG1(G1[1])
    g = g^(div(size(F)-1, 5))
    @test !isone(g)

    for i = 1:5
      x = rand(F)
      while iszero(x)
        x = rand(F)
      end
      c = mG1\x
      el = mG1(c)
      if el != x
        dl = mG\(el*inv(x))
        @test iszero(mod(dl[1], 5))
      end
    end
  end
end

@testset "FqFiniteField" begin

  for p in [31, 11, 101]
    F = FiniteField(fmpz(p), 2, "a")[1]
    G, mG = unit_group(F)
    #Test generator
    g = mG(G[1])
    g1 = g
    for i = 2:p^2-2
      g1 *= g
      @test !isone(g1)
    end
    for i = 1:5
      x = rand(F)
      while iszero(x)
        x = rand(F)
      end
      elG = mG\x
      @test mG(elG) == x
    end

    G1, mG1 = unit_group(F, n_quo = 5)
    @test order(G1) == 5
    g = mG1(G1[1])
    g = g^(div(size(F)-1, 5))
    @test !isone(g)

    for i = 1:5
      x = rand(F)
      while iszero(x)
        x = rand(F)
      end
      c = mG1\x
      el = mG1(c)
      if el != x
        dl = mG\(el*inv(x))
        @test iszero(mod(dl[1], 5))
      end
    end
  end
end


text:
### A Pluto.jl notebook ###
# v0.14.6

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 839c36b1-1318-4892-84ab-389b40b2ea48
using PlutoUI, Plots

# ╔═╡ efaf0218-42bd-46dd-9f93-de666df37100
md"## Julia Lang

* A relatively new dynamic programming language (≈ 3 years since 1.0)
* General purpose and high level language
* Core paradigm is *Multiple Dispatch*
* Performance comparable to C/C++/Fortran
* For example: *for loops* don't have a performance penalty

### Reasons for Speed
* JIT - Just in Time compilation
* type stable functions
* because of the type system Julia can compile efficient code

### Helpfuls resources
* However, severalc differences to other languages
* Consult [performance tips](https://docs.julialang.org/en/v1/manual/performance-tips/)
* [Noteworthy differences to other language](https://docs.julialang.org/en/v1/manual/noteworthy-differences/)
* Ecosystem for image processing for microscopy data is small

### Some important packages
* [Bioformatsloader.jl](https://github.com/ahnlabb/BioformatsLoader.jl)
* [View5D.jl](https://github.com/RainerHeintzmann/View5D.jl)
* [DeconvOptim.jl](https://github.com/roflmaostc/DeconvOptim.jl)
* [Images.jl](https://github.com/JuliaImages/Images.jl/)
* [Napari.jl](https://github.com/mkitti/Napari.jl)
* [FourierTools.jl](https://github.com/bionanoimaging/FourierTools.jl)
* [TestImages.jl](https://github.com/JuliaImages/TestImages.jl)
"

# ╔═╡ ad907ec3-7f30-4508-9599-a734c18593af
md"#### sqrt example which fails in Julia because of the types"

# ╔═╡ f8d94633-d125-4da1-a446-bbf4e8fd43a9
sqrt(-1)

# ╔═╡ 8c83a90b-642a-4e9c-abe6-0e1ce0c819b5
sqrt(1) |> typeof

# ╔═╡ d78a922c-efa2-4935-881a-b5dccaa09caa
md"## Multiple Dispatch

* In general no type signatures are needed
* However we can specialize if we want
"

# ╔═╡ 8f247848-6258-4538-9bea-df4897d5716a
md"""
x = $(@bind x Slider(0:0.1:10; default=8, show_value=true))
z = $(@bind z Slider(1:0.01:10; default=8, show_value=true))
"""

# ╔═╡ 6b5f335d-15af-43f2-afc4-1b51acc693cc
# However, we can specialize if we want
function f(x::Int, y::Int)
	return 42
end

# ╔═╡ 8021d1d9-773d-46d4-9804-6fbe7978f4df
ϵ = 10

# ╔═╡ fc08b77e-565a-4edc-a2dd-0b7dc7acd155
md"## Arrays"

# ╔═╡ 26f41519-9f11-46e6-aa27-abd1142b00fb
arr = ones(Float32, (10, 10));

# ╔═╡ 708f34f7-1308-4946-9ace-998476606f9c
arr * arr # matrix multiplication

# ╔═╡ c24d72a7-e4f5-427e-a55a-6ca6b7ec9b9e
arr .* arr # point wise

# ╔═╡ 7183ba48-5785-41af-afc6-f80bc07916cf
md"## Small for loop array example

Consider an array where we want to sum up all positive numbers
"

# ╔═╡ aa72b577-18f6-4871-b232-d17c2bb5c7fd
Γ = randn((100_000_000, ));

# ╔═╡ 1e79336c-9f63-48ca-aa47-5c96bc7dfd9c
result1 = sum(Γ[ Γ .> 0.0])

# ╔═╡ dec2cd2e-2263-4139-b461-1c4c5c49976f
md"## Type Example"

# ╔═╡ 35e81c4d-8673-4b2a-aa0f-457aaefcddb7
abstract type Bit end

# ╔═╡ 6c57ae37-12d0-43c7-9043-3e5dbb0b85e3
struct One <: Bit end

# ╔═╡ e779870d-966a-4ed7-9ab4-54e534518642
struct Zero <: Bit end

# ╔═╡ 214b7930-763b-4034-bf50-05a0c369c02d
begin
	import Base.+
	+(::Bit, ::Bit) = Zero()
	+(::One, ::Zero) = One()
	+(::Zero, ::One) = One()
end

# ╔═╡ 3b8c572c-3d2e-4f68-96f3-b4552ac2d20b
# no type signatures needed
function f(x, y)
	return x + y
end

# ╔═╡ 6ea00fdb-7659-4bd5-83a1-09af4c6f31d7
f(x, z * 1im)

# ╔═╡ 853f3cd3-20c6-4479-8846-0c6b7a34e410
begin
	r = 0:0.005:10
	plot(r, sin.(x .* r .+ z))
end

# ╔═╡ 034676c5-8a30-488f-a685-d0d4c20d078b
f(ϵ, x + z)

# ╔═╡ b71710ee-09fe-4568-9aa7-fc5eaa4b0cb3
function sum_positive(Γ)
	res = zero(eltype(Γ))
	for entry in Γ
		if entry > 0
			res += entry
		end
	end
	return res
end

# ╔═╡ 76f78b1c-b796-4ec5-b857-12635be06a0d
result2 = sum_positive(Γ)

# ╔═╡ f9f7a2fe-c0e4-412d-bc88-22b2a8cf9f4c
result1 ≈ result2

# ╔═╡ b0c6f430-3ade-40ef-9468-9701391be788
One() + Zero()

# ╔═╡ e04b7ae2-3282-4093-a35e-1357d48cdff9
Zero() + Zero()

# ╔═╡ Cell order:
# ╠═839c36b1-1318-4892-84ab-389b40b2ea48
# ╟─efaf0218-42bd-46dd-9f93-de666df37100
# ╠═ad907ec3-7f30-4508-9599-a734c18593af
# ╠═f8d94633-d125-4da1-a446-bbf4e8fd43a9
# ╠═8c83a90b-642a-4e9c-abe6-0e1ce0c819b5
# ╠═d78a922c-efa2-4935-881a-b5dccaa09caa
# ╠═3b8c572c-3d2e-4f68-96f3-b4552ac2d20b
# ╟─8f247848-6258-4538-9bea-df4897d5716a
# ╠═6ea00fdb-7659-4bd5-83a1-09af4c6f31d7
# ╠═6b5f335d-15af-43f2-afc4-1b51acc693cc
# ╠═8021d1d9-773d-46d4-9804-6fbe7978f4df
# ╠═034676c5-8a30-488f-a685-d0d4c20d078b
# ╠═853f3cd3-20c6-4479-8846-0c6b7a34e410
# ╠═fc08b77e-565a-4edc-a2dd-0b7dc7acd155
# ╠═26f41519-9f11-46e6-aa27-abd1142b00fb
# ╠═708f34f7-1308-4946-9ace-998476606f9c
# ╠═c24d72a7-e4f5-427e-a55a-6ca6b7ec9b9e
# ╠═7183ba48-5785-41af-afc6-f80bc07916cf
# ╠═aa72b577-18f6-4871-b232-d17c2bb5c7fd
# ╠═1e79336c-9f63-48ca-aa47-5c96bc7dfd9c
# ╠═76f78b1c-b796-4ec5-b857-12635be06a0d
# ╠═b71710ee-09fe-4568-9aa7-fc5eaa4b0cb3
# ╠═f9f7a2fe-c0e4-412d-bc88-22b2a8cf9f4c
# ╠═dec2cd2e-2263-4139-b461-1c4c5c49976f
# ╠═35e81c4d-8673-4b2a-aa0f-457aaefcddb7
# ╠═6c57ae37-12d0-43c7-9043-3e5dbb0b85e3
# ╠═e779870d-966a-4ed7-9ab4-54e534518642
# ╠═214b7930-763b-4034-bf50-05a0c369c02d
# ╠═b0c6f430-3ade-40ef-9468-9701391be788
# ╠═e04b7ae2-3282-4093-a35e-1357d48cdff9


text:
import Base.Test: @test
using RNG.Xorshifts

stdout_ = STDOUT
pwd_ = pwd()
cd(joinpath(Pkg.dir("RNG"), "test/Xorshifts"))
mkpath("./actual")

for (rng_name, rng) in (
        (:Xorshift64, :(Xorshift64(123))),
        (:Xorshift64Star, :(Xorshift64Star(123))),
        (:Xorshift128, :(Xorshift128(123))),
        (:Xorshift128Star, :(Xorshift128Star(123))),
        (:Xorshift128Plus, :(Xorshift128Plus(123))),
        (:Xorshift1024, :(Xorshift1024(123))),
        (:Xorshift1024Star, :(Xorshift1024Star(123))),
        (:Xorshift1024Plus, :(Xorshift1024Plus(123))),
        (:Xoroshiro128, :(Xoroshiro128(123))),
        (:Xoroshiro128Star, :(Xoroshiro128Star(123))),
        (:Xoroshiro128Plus, :(Xoroshiro128Plus(123))),
    )

    outfile = open(string(
        "./actual/check-$(lowercase("$rng_name")).out"
    ), "w")
    redirect_stdout(outfile)

    @eval x = $rng_name()
    srand(x)
    srand(x, 1)
    @eval x = $rng

    for i in 1:100
        @printf "%.9f\n" rand(x)
    end

    close(outfile)
end
redirect_stdout(stdout_)

@test success(`diff -ru expected actual`)
cd(pwd_)


text:
"""
Creates Ohms constraints for DC branches

```
p[f_idx] + p[t_idx] == p * g[l] * (wdc[f_bus] - wdcr[f_bus,t_bus])
```
"""
function constraint_ohms_dc_branch(pm::_PM.AbstractBFQPModel, n::Int, f_bus, t_bus, f_idx, t_idx, r, p)
    l = f_idx[1];
    p_dc_fr = _PM.var(pm, n,  :p_dcgrid, f_idx)
    p_dc_to = _PM.var(pm, n,  :p_dcgrid, t_idx)
    ccm_dcgrid = _PM.var(pm, n,  :ccm_dcgrid, l)
    wdc_fr = _PM.var(pm, n,  :wdc, f_bus)
    wdc_to = _PM.var(pm, n,  :wdc, t_bus)

    JuMP.@constraint(pm.model, p_dc_fr + p_dc_to ==  r * p * ccm_dcgrid)
    JuMP.@constraint(pm.model, p_dc_fr^2 <= p^2 * wdc_fr * ccm_dcgrid)
    JuMP.@constraint(pm.model, wdc_to == wdc_fr - 2 * r * (p_dc_fr/p) + (r)^2 * ccm_dcgrid)
end

function constraint_ohms_dc_branch(pm::_PM.AbstractBFConicModel, n::Int, f_bus, t_bus, f_idx, t_idx, r, p)
    l = f_idx[1];
    p_dc_fr = _PM.var(pm, n,  :p_dcgrid, f_idx)
    p_dc_to = _PM.var(pm, n,  :p_dcgrid, t_idx)
    ccm_dcgrid = _PM.var(pm, n,  :ccm_dcgrid, l)
    wdc_fr = _PM.var(pm, n,  :wdc, f_bus)
    wdc_to = _PM.var(pm, n,  :wdc, t_bus)

    JuMP.@constraint(pm.model, p_dc_fr + p_dc_to ==  r * p * ccm_dcgrid)
    JuMP.@constraint(pm.model, [p*wdc_fr/sqrt(2), p*ccm_dcgrid/sqrt(2), p_dc_fr/sqrt(2), p_dc_fr/sqrt(2)] in JuMP.RotatedSecondOrderCone())
    JuMP.@constraint(pm.model, wdc_to == wdc_fr - 2 * r * (p_dc_fr/p) + (r)^2 * ccm_dcgrid)
end
"""
Model to approximate cross products of node voltages
```
wdcr[(i,j)] <= wdc[i]*wdc[j]
```
"""
function constraint_voltage_dc(pm::_PM.AbstractBFModel, n::Int = pm.cnw)
# do nothing
end

function variable_dcbranch_current(pm::_PM.AbstractBFModel; kwargs...)
    variable_dcbranch_current_sqr(pm; kwargs...)
end
function constraint_dc_branch_current(pm::_PM.AbstractBFModel, n::Int, f_bus, f_idx, ccm_max, p)
# do nothing
end

########## TNEP constraints #####################

"""
Creates Ohms constraints for DC branches

```
p[f_idx] + p[t_idx] == p * g[l] * (wdc[f_bus] - wdcr[f_bus,t_bus])
```
"""
function constraint_ohms_dc_branch_ne(pm::_PM.AbstractBFQPModel, n::Int, f_bus, t_bus, f_idx, t_idx, r, p)
    l = f_idx[1];
    p_dc_fr = _PM.var(pm, n, :p_dcgrid_ne, f_idx)
    p_dc_to = _PM.var(pm, n, :p_dcgrid_ne, t_idx)
    ccm_dcgrid = _PM.var(pm, n, :ccm_dcgrid_ne, l)
    z = _PM.var(pm, n, :branchdc_ne, l)

    wdc_to = []
    wdc_fr = []

    wdc_to, wdc_fr = contraint_ohms_dc_branch_busvoltage_structure_W(pm, n, f_bus, t_bus, wdc_to, wdc_fr)
    wdc_du_to = _PM.var(pm, n, :wdc_du_to, l)
    wdc_du_fr = _PM.var(pm, n, :wdc_du_fr, l)

    JuMP.@constraint(pm.model, p_dc_fr + p_dc_to ==  r * p * ccm_dcgrid)
    JuMP.@constraint(pm.model, p_dc_fr^2 <= p^2 * wdc_fr * ccm_dcgrid)

    # different type of on_off constraint
    JuMP.@constraint(pm.model, wdc_du_to  == wdc_du_fr - 2 * r * (p_dc_fr/p) + (r)^2 * ccm_dcgrid)
    JuMP.@constraint(pm.model, wdc_du_to <= wdc_to - JuMP.lower_bound(wdc_to)*(1-z))
    JuMP.@constraint(pm.model, wdc_du_to >= wdc_to - JuMP.upper_bound(wdc_fr)*(1-z))
    JuMP.@constraint(pm.model, wdc_du_fr <= wdc_fr - JuMP.lower_bound(wdc_to)*(1-z))
    JuMP.@constraint(pm.model, wdc_du_fr >= wdc_fr - JuMP.upper_bound(wdc_fr)*(1-z))
    JuMP.@constraint(pm.model, wdc_du_to <= z* JuMP.upper_bound(wdc_to))
    JuMP.@constraint(pm.model, wdc_du_to >= z* JuMP.lower_bound(wdc_to))
    JuMP.@constraint(pm.model, wdc_du_fr <= z* JuMP.upper_bound(wdc_fr))
    JuMP.@constraint(pm.model, wdc_du_fr >= z* JuMP.lower_bound(wdc_fr))

end

function constraint_ohms_dc_branch_ne(pm::_PM.AbstractBFConicModel, n::Int, f_bus, t_bus, f_idx, t_idx, r, p)
    l = f_idx[1];
    p_dc_fr = _PM.var(pm, n, :p_dcgrid_ne, f_idx)
    p_dc_to = _PM.var(pm, n, :p_dcgrid_ne, t_idx)
    ccm_dcgrid = _PM.var(pm, n, :ccm_dcgrid_ne, l)
    z = _PM.var(pm, n, :branchdc_ne, l)

    wdc_to = []
    wdc_fr = []

    wdc_to, wdc_fr = contraint_ohms_dc_branch_busvoltage_structure_W(pm, n, f_bus, t_bus, wdc_to, wdc_fr)
    wdc_du_to = _PM.var(pm, n, :wdc_du_to, l)
    wdc_du_fr = _PM.var(pm, n, :wdc_du_fr, l)

    JuMP.@constraint(pm.model, p_dc_fr + p_dc_to ==  r * p * ccm_dcgrid)
    JuMP.@constraint(pm.model, [p*wdc_fr/sqrt(2), p*ccm_dcgrid/sqrt(2), p_dc_fr/sqrt(2), p_dc_fr/sqrt(2)] in JuMP.RotatedSecondOrderCone())

    # different type of on_off constraint
    JuMP.@constraint(pm.model, wdc_du_to  == wdc_du_fr - 2 * r * (p_dc_fr/p) + (r)^2 * ccm_dcgrid)
    JuMP.@constraint(pm.model, wdc_du_to <= wdc_to - JuMP.lower_bound(wdc_to)*(1-z))
    JuMP.@constraint(pm.model, wdc_du_to >= wdc_to - JuMP.upper_bound(wdc_to)*(1-z))
    JuMP.@constraint(pm.model, wdc_du_fr <= wdc_fr - JuMP.lower_bound(wdc_fr)*(1-z))
    JuMP.@constraint(pm.model, wdc_du_fr >= wdc_fr - JuMP.upper_bound(wdc_fr)*(1-z))
    JuMP.@constraint(pm.model, wdc_du_to <= z* JuMP.upper_bound(wdc_to))
    JuMP.@constraint(pm.model, wdc_du_to >= z* JuMP.lower_bound(wdc_to))
    JuMP.@constraint(pm.model, wdc_du_fr <= z* JuMP.upper_bound(wdc_fr))
    JuMP.@constraint(pm.model, wdc_du_fr >= z* JuMP.lower_bound(wdc_fr))
end

"""
Model to approximate cross products of node voltages

```
wdcr[(i,j)] <= wdc[i]*wdc[j]
```
"""
function constraint_voltage_dc_ne(pm::_PM.AbstractBFModel, n::Int)
# do nothing
end
function variable_dcbranch_current_ne(pm::_PM.AbstractBFModel; kwargs...)
    variable_dcbranch_current_sqr_ne(pm; kwargs...)
end
function contraint_ohms_dc_branch_busvoltage_structure_W(pm::_PM.AbstractPowerModel, n::Int, f_bus, t_bus, wdc_to, wdc_fr)
    for i in _PM.ids(pm, n, :busdc_ne)
        if t_bus == i
            wdc_to = _PM.var(pm, n, :wdc_ne, t_bus)
        end
        if f_bus == i
            wdc_fr = _PM.var(pm, n, :wdc_ne, f_bus)
        end
    end
    for i in _PM.ids(pm, n, :busdc)
        if t_bus == i
            wdc_to = _PM.var(pm, n, :wdc, t_bus)
        end
        if f_bus == i
            wdc_fr = _PM.var(pm, n, :wdc, f_bus)
        end
    end
    return wdc_to, wdc_fr
end


text:
using ImageTransformations: box_extrapolation

@testset "deprecations" begin
    @info "deprecation warnings are expected"

    @testset "imrotate" begin
        # deprecate potionsal arguments in favor of keyword arguments
        img = repeat(range(0,stop=1,length=100), 1, 100)

        @test nearlysame(imrotate(img, π/4, Constant()),                imrotate(img, π/4, method=Constant()))
        @test nearlysame(imrotate(img, π/4, Constant(), 1),             imrotate(img, π/4, method=Constant(), fillvalue=1))
        @test nearlysame(imrotate(img, π/4, 1, Constant()),             imrotate(img, π/4, method=Constant(), fillvalue=1))
        @test nearlysame(imrotate(img, π/4, 1),                         imrotate(img, π/4, fillvalue=1))
        @test nearlysame(imrotate(img, π/4, Periodic()),                imrotate(img, π/4, fillvalue=Periodic()))

        @test nearlysame(imrotate(img, π/4, axes(img), Constant()),     imrotate(img, π/4, axes(img), method=Constant()))
        @test nearlysame(imrotate(img, π/4, axes(img), Constant(), 1),  imrotate(img, π/4, axes(img), method=Constant(), fillvalue=1))
        @test nearlysame(imrotate(img, π/4, axes(img), 1, Constant()),  imrotate(img, π/4, axes(img), method=Constant(), fillvalue=1))
        @test nearlysame(imrotate(img, π/4, axes(img), 1),              imrotate(img, π/4, axes(img), fillvalue=1))
        @test nearlysame(imrotate(img, π/4, axes(img), Periodic()),     imrotate(img, π/4, axes(img), fillvalue=Periodic()))
    end

    @testset "warp" begin
        # deprecate potionsal arguments in favor of keyword arguments
        img = repeat(range(0,stop=1,length=100), 1, 100)
        tfm = recenter(RotMatrix(pi/8), center(img))

        @test nearlysame(warp(img, tfm, Constant()),                warp(img, tfm, method=Constant()))
        @test nearlysame(warp(img, tfm, Constant(), 1),             warp(img, tfm, method=Constant(), fillvalue=1))
        @test nearlysame(warp(img, tfm, 1, Constant()),             warp(img, tfm, method=Constant(), fillvalue=1))
        @test nearlysame(warp(img, tfm, 1),                         warp(img, tfm, fillvalue=1))
        @test nearlysame(warp(img, tfm, Periodic()),                warp(img, tfm, fillvalue=Periodic()))

        @test nearlysame(warp(img, tfm, axes(img), Constant()),     warp(img, tfm, axes(img), method=Constant()))
        @test nearlysame(warp(img, tfm, axes(img), Constant(), 1),  warp(img, tfm, axes(img), method=Constant(), fillvalue=1))
        @test nearlysame(warp(img, tfm, axes(img), 1, Constant()),  warp(img, tfm, axes(img), method=Constant(), fillvalue=1))
        @test nearlysame(warp(img, tfm, axes(img), 1),              warp(img, tfm, axes(img), fillvalue=1))
        @test nearlysame(warp(img, tfm, axes(img), Periodic()),     warp(img, tfm, axes(img), fillvalue=Periodic()))
    end

    @testset "warpedview" begin
        # deprecate `warpedview` in favor of `WarpedView`
        img = repeat(range(0,stop=1,length=100), 1, 100)
        tfm = recenter(RotMatrix(pi/8), center(img))

        @test nearlysame(warpedview(img, tfm, Constant()),                WarpedView(img, tfm, method=Constant()))
        @test nearlysame(warpedview(img, tfm, Constant(), 1),             WarpedView(img, tfm, method=Constant(), fillvalue=1))
        @test nearlysame(warpedview(img, tfm, 1, Constant()),             WarpedView(img, tfm, method=Constant(), fillvalue=1))
        @test nearlysame(warpedview(img, tfm, 1),                         WarpedView(img, tfm, fillvalue=1))
        @test nearlysame(warpedview(img, tfm, Periodic()),                WarpedView(img, tfm, fillvalue=Periodic()))

        @test nearlysame(warpedview(img, tfm, axes(img), Constant()),     WarpedView(img, tfm, axes(img), method=Constant()))
        @test nearlysame(warpedview(img, tfm, axes(img), Constant(), 1),  WarpedView(img, tfm, axes(img), method=Constant(), fillvalue=1))
        @test nearlysame(warpedview(img, tfm, axes(img), 1, Constant()),  WarpedView(img, tfm, axes(img), method=Constant(), fillvalue=1))
        @test nearlysame(warpedview(img, tfm, axes(img), 1),              WarpedView(img, tfm, axes(img), fillvalue=1))
        @test nearlysame(warpedview(img, tfm, axes(img), Periodic()),     WarpedView(img, tfm, axes(img), fillvalue=Periodic()))
    end

    @testset "WarpedView" begin
        # deprecate potionsal arguments in favor of keyword arguments
        img = repeat(range(0,stop=1,length=100), 1, 100)
        tfm = recenter(RotMatrix(pi/8), center(img))

        @test nearlysame(WarpedView(img, tfm, Constant()),                WarpedView(img, tfm, method=Constant()))
        @test nearlysame(WarpedView(img, tfm, Constant(), 1),             WarpedView(img, tfm, method=Constant(), fillvalue=1))
        @test nearlysame(WarpedView(img, tfm, 1, Constant()),             WarpedView(img, tfm, method=Constant(), fillvalue=1))
        @test nearlysame(WarpedView(img, tfm, 1),                         WarpedView(img, tfm, fillvalue=1))
        @test nearlysame(WarpedView(img, tfm, Periodic()),                WarpedView(img, tfm, fillvalue=Periodic()))

        @test nearlysame(WarpedView(img, tfm, axes(img), Constant()),     WarpedView(img, tfm, axes(img), method=Constant()))
        @test nearlysame(WarpedView(img, tfm, axes(img), Constant(), 1),  WarpedView(img, tfm, axes(img), method=Constant(), fillvalue=1))
        @test nearlysame(WarpedView(img, tfm, axes(img), 1, Constant()),  WarpedView(img, tfm, axes(img), method=Constant(), fillvalue=1))
        @test nearlysame(WarpedView(img, tfm, axes(img), 1),              WarpedView(img, tfm, axes(img), fillvalue=1))
        @test nearlysame(WarpedView(img, tfm, axes(img), Periodic()),     WarpedView(img, tfm, axes(img), fillvalue=Periodic()))
    end

    @testset "invwarpedview" begin
        # deprecate potionsal arguments in favor of keyword arguments
        img = repeat(range(0,stop=1,length=100), 1, 100)
        tfm = recenter(RotMatrix(pi/8), center(img))

        @test nearlysame(invwarpedview(img, tfm, Constant()),                InvWarpedView(box_extrapolation(img; method=Constant()), tfm))
        @test nearlysame(invwarpedview(img, tfm, Constant(), 1),             InvWarpedView(box_extrapolation(img; method=Constant(), fillvalue=1), tfm))
        @test nearlysame(invwarpedview(img, tfm, 1, Constant()),             InvWarpedView(box_extrapolation(img; method=Constant(), fillvalue=1), tfm))
        @test nearlysame(invwarpedview(img, tfm, 1),                         InvWarpedView(box_extrapolation(img; fillvalue=1), tfm))
        @test nearlysame(invwarpedview(img, tfm, Periodic()),                InvWarpedView(box_extrapolation(img; fillvalue=Periodic()), tfm))

        @test nearlysame(InvWarpedView(img, tfm, Constant()),                InvWarpedView(img, tfm; method=Constant()))
        @test nearlysame(InvWarpedView(img, tfm, Constant(), 1),             InvWarpedView(img, tfm; method=Constant(), fillvalue=1))
        @test nearlysame(InvWarpedView(img, tfm, 1, Constant()),             InvWarpedView(img, tfm; method=Constant(), fillvalue=1))
        @test nearlysame(InvWarpedView(img, tfm, 1),                         InvWarpedView(img, tfm; fillvalue=1))
        @test nearlysame(InvWarpedView(img, tfm, Periodic()),                InvWarpedView(img, tfm; fillvalue=Periodic()))

        @test nearlysame(invwarpedview(img, tfm, axes(img), Constant()),                InvWarpedView(box_extrapolation(img; method=Constant()), tfm, axes(img)))
        @test nearlysame(invwarpedview(img, tfm, axes(img), Constant(), 1),             InvWarpedView(box_extrapolation(img; method=Constant(), fillvalue=1), tfm, axes(img)))
        @test nearlysame(invwarpedview(img, tfm, axes(img), 1, Constant()),             InvWarpedView(box_extrapolation(img; method=Constant(), fillvalue=1), tfm, axes(img)))
        @test nearlysame(invwarpedview(img, tfm, axes(img), 1),                         InvWarpedView(box_extrapolation(img; fillvalue=1), tfm, axes(img)))
        @test nearlysame(invwarpedview(img, tfm, axes(img), Periodic()),                InvWarpedView(box_extrapolation(img; fillvalue=Periodic()), tfm, axes(img)))

        @test nearlysame(InvWarpedView(img, tfm, axes(img), Constant()),                InvWarpedView(img, tfm, axes(img); method=Constant()))
        @test nearlysame(InvWarpedView(img, tfm, axes(img), Constant(), 1),             InvWarpedView(img, tfm, axes(img); method=Constant(), fillvalue=1))
        @test nearlysame(InvWarpedView(img, tfm, axes(img), 1, Constant()),             InvWarpedView(img, tfm, axes(img); method=Constant(), fillvalue=1))
        @test nearlysame(InvWarpedView(img, tfm, axes(img), 1),                         InvWarpedView(img, tfm, axes(img); fillvalue=1))
        @test nearlysame(InvWarpedView(img, tfm, axes(img), Periodic()),                InvWarpedView(img, tfm, axes(img); fillvalue=Periodic()))
    end
end


text:
"""
    $SIGNATURES

Perform upward inference using a state machine solution approach.

Notes:
- will call on values from children or parent cliques
- can be called multiple times
- Assumes all cliques in tree are being solved simultaneously and in similar manner.
- State machine rev.1 -- copied from first TreeBasedInitialization.jl.
- Doesn't do partial initialized state properly yet.
"""
function cliqInitSolveUpByStateMachine!(dfg::G,
                                        tree::AbstractBayesTree,
                                        cliq::TreeClique,
                                        timeout::Union{Nothing, <:Real}=nothing;
                                        N::Int=100,
                                        verbose::Bool=false,
                                        verbosefid=stdout,
                                        oldcliqdata::BayesTreeNodeData=BayesTreeNodeData(),
                                        drawtree::Bool=false,
                                        show::Bool=false,
                                        incremental::Bool=true,
                                        limititers::Int=-1,
                                        upsolve::Bool=true,
                                        downsolve::Bool=true,
                                        recordhistory::Bool=false,
                                        delay::Bool=false,
                                        injectDelayBefore::Union{Nothing,Pair{<:Function, <:Real}}=nothing,
                                        logger::SimpleLogger=SimpleLogger(Base.stdout)) where {G <: AbstractDFG, AL <: AbstractLogger}
  #
  children = getChildren(tree, cliq)#Graphs.out_neighbors(cliq, tree.bt)

  prnt = getParent(tree, cliq)

  destType = (G <: InMemoryDFGTypes) ? G : InMemDFGType

  csmc = CliqStateMachineContainer(dfg, initfg(destType, solverParams=getSolverParams(dfg)), tree, cliq, prnt, children, incremental, drawtree, downsolve, delay, getSolverParams(dfg), Dict{Symbol,String}(), oldcliqdata, logger)

  nxt = upsolve ? canCliqMargRecycle_StateMachine : (downsolve ? canCliqMargRecycle_StateMachine : error("must attempt either up or down solve"))

  csmiter_cb = getSolverParams(dfg).drawCSMIters ? ((st::StateMachine)->(cliq.attributes["xlabel"] = st.iter)) : ((st)->())

  statemachine = StateMachine{CliqStateMachineContainer}(next=nxt, name="cliq$(cliq.index)")

  # store statemachine and csmc in task
  if dfg.solverParams.dbg || recordhistory
    task_local_storage(:statemachine, statemachine)
    task_local_storage(:csmc, csmc)
  end

  while statemachine(csmc, timeout, verbose=verbose, verbosefid=verbosefid, verboseXtra=getCliqueStatus(csmc.cliq), iterlimit=limititers, recordhistory=recordhistory, housekeeping_cb=csmiter_cb, injectDelayBefore=injectDelayBefore); end
  statemachine.history
end

## ==============================================================================================
# Prepare CSM (based on FSM) entry points
## ==============================================================================================

function tryCliqStateMachineSolve!( dfg::AbstractDFG,
    treel::AbstractBayesTree,
    i::Int,
    timeout::Union{Nothing, <:Real}=nothing;
    verbose::Bool=false,
    verbosefid=stdout,
    N::Int=100,
    oldtree::AbstractBayesTree=emptyBayesTree(),
    drawtree::Bool=false,
    limititers::Int=-1,
    downsolve::Bool=false,
    incremental::Bool=false,
    injectDelayBefore::Union{Nothing,<:Pair{<:Function, <:Real}}=nothing,
    delaycliqs::Vector{Symbol}=Symbol[],
    recordcliqs::Vector{Symbol}=Symbol[])
#
clst = :na
cliq = getClique(treel, i)
syms = getCliqFrontalVarIds(cliq) # ids =
oldcliq = attemptTreeSimilarClique(oldtree, getCliqueData(cliq))
oldcliqdata = getCliqueData(oldcliq)
opts = getSolverParams(dfg)
# Base.rm(joinpath(opts.logpath,"logs/cliq$i"), recursive=true, force=true)
mkpath(joinpath(opts.logpath,"logs/cliq$i/"))
logger = SimpleLogger(open(joinpath(opts.logpath,"logs/cliq$i/log.txt"), "w+")) # NullLogger()
# global_logger(logger)
history = Vector{Tuple{DateTime, Int, Function, CliqStateMachineContainer}}()
recordthiscliq = length(intersect(recordcliqs,syms)) > 0
delaythiscliq = length(intersect(delaycliqs,syms)) > 0
try
history = cliqInitSolveUpByStateMachine!(dfg, treel, cliq, timeout, N=N,
              verbose=verbose, verbosefid=verbosefid, drawtree=drawtree,
              oldcliqdata=oldcliqdata,
              injectDelayBefore=injectDelayBefore,
              limititers=limititers, downsolve=downsolve, recordhistory=recordthiscliq, incremental=incremental, delay=delaythiscliq, logger=logger )
#
if getSolverParams(dfg).dbg || length(history) >= limititers && limititers != -1
@info "writing logs/cliq$i/csm.txt"
# @save "/tmp/cliqHistories/cliq$i.jld2" history
fid = open(joinpath(opts.logpath,"logs/cliq$i/csm.txt"), "w")
printCliqHistorySummary(fid, history)
close(fid)
end
flush(logger.stream)
close(logger.stream)
# clst = getCliqueStatus(cliq)
# clst = cliqInitSolveUp!(dfg, treel, cliq, drawtree=drawtree, limititers=limititers )
catch err
## TODO -- use this format instead
# io = IOBuffer()
# showerror(io, ex, catch_backtrace())
# err = String(take!(io))
# msg = "Error while packing '$(f.label)' as '$fnctype', please check the unpacking/packing converters for this factor - \r\n$err"
# error(msg)

## OLD format
bt = catch_backtrace()
println()
showerror(stderr, err, bt)
# @warn "writing /tmp/caesar/logs/cliq$i/*.txt"
fid = open(joinpath(opts.logpath,"logs/cliq$i/stacktrace.txt"), "w")
showerror(fid, err, bt)
close(fid)
fid = open(joinpath(opts.logpath,"logs/cliq$(i)_stacktrace.txt"), "w")
showerror(fid, err, bt)
close(fid)
# @save "/tmp/cliqHistories/$(cliq.label).jld2" history
fid = open(joinpath(opts.logpath,"logs/cliq$i/csm.txt"), "w")
printCliqHistorySummary(fid, history)
close(fid)
fid = open(joinpath(opts.logpath,"logs/cliq$(i)_csm.txt"), "w")
printCliqHistorySummary(fid, history)
close(fid)
flush(logger.stream)
close(logger.stream)
# error(err)
rethrow()
end
# if !(clst in [:upsolved; :downsolved; :marginalized])
#   error("Clique $(cliq.index), initInferTreeUp! -- cliqInitSolveUp! did not arrive at the desired solution statu: $clst")
# end
return history
end

"""
    $SIGNATURES

Perform tree based initialization of all variables not yet initialized in factor graph as non-blocking method.

Notes:
- To simplify debugging, this method does not include the usual `@ sync` around all the state machine async processes.
- Extract the error stack with a `fetch` on the failed process return by this function.

Related

initInferTreeUp!
"""
function asyncTreeInferUp!( dfg::AbstractDFG,
                            treel::AbstractBayesTree,
                            timeout::Union{Nothing, <:Real}=nothing;
                            oldtree::AbstractBayesTree=emptyBayesTree(),
                            verbose::Bool=false,
                            verbosefid=stdout,
                            drawtree::Bool=false,
                            N::Int=100,
                            limititers::Int=-1,
                            downsolve::Bool=false,
                            incremental::Bool=false,
                            limititercliqs::Vector{Pair{Symbol, Int}}=Pair{Symbol, Int}[],
                            injectDelayBefore::Union{Nothing,Vector{<:Pair{Int,<:Pair{<:Function,<:Real}}}}=nothing,
                            skipcliqids::Vector{Symbol}=Symbol[],
                            delaycliqs::Vector{Symbol}=Symbol[],
                            recordcliqs::Vector{Symbol}=Symbol[] )
  #
  resetTreeCliquesForUpSolve!(treel)
  if drawtree
    pdfpath = joinLogPath(dfg,"bt.pdf")
    drawTree(treel, show=false, filepath=pdfpath)
  end

  # queue all the tasks
  alltasks = Vector{Task}(undef, length(getCliques(treel)))
  # cliqHistories = Dict{Int,Vector{Tuple{DateTime, Int, Function, CliqStateMachineContainer}}}()
  if !isTreeSolved(treel, skipinitialized=true)
    # @sync begin
      # duplicate int i into async (important for concurrency)
      for i in 1:length(getCliques(treel))
        scsym = getCliqFrontalVarIds(getClique(treel, i))
        if length(intersect(scsym, skipcliqids)) == 0
          limthiscsm = filter(x -> (x[1] in scsym), limititercliqs)
          limiter = 0<length(limthiscsm) ? limthiscsm[1][2] : limititers
          injDelay = if injectDelayBefore === nothing
            nothing
          else
            idb = filter((x)->x[1]==i,injectDelayBefore)
            length(idb) == 1 ? idb[1][2] : nothing
          end
          alltasks[i] = @async tryCliqStateMachineSolve!(dfg, treel, i, timeout, oldtree=oldtree, verbose=verbose, verbosefid=verbosefid, drawtree=drawtree, limititers=limiter, downsolve=downsolve, delaycliqs=delaycliqs, recordcliqs=recordcliqs, injectDelayBefore=injDelay, incremental=incremental, N=N)
        end # if
      end # for
    # end # sync
  end # if

  return alltasks #, cliqHistories
end

"""
    $SIGNATURES

Perform tree based initialization of all variables not yet initialized in factor graph.

Related

asyncTreeInferUp!
"""
function initInferTreeUp!(dfg::AbstractDFG,
                          treel::AbstractBayesTree,
                          timeout::Union{Nothing, <:Real}=nothing;
                          oldtree::AbstractBayesTree=emptyBayesTree(),
                          verbose::Bool=false,
                          verbosefid=stdout,
                          drawtree::Bool=false,
                          N::Int=100,
                          limititers::Int=-1,
                          downsolve::Bool=false,
                          incremental::Bool=false,
                          limititercliqs::Vector{Pair{Symbol, Int}}=Pair{Symbol, Int}[],
                          injectDelayBefore::Union{Nothing,Vector{<:Pair{Int,<:Pair{<:Function,<:Real}}}}=nothing,
                          skipcliqids::Vector{Symbol}=Symbol[],
                          recordcliqs::Vector{Symbol}=Symbol[],
                          delaycliqs::Vector{Symbol}=Symbol[],
                          alltasks::Vector{Task}=Task[],
                          runtaskmonitor::Bool=true)
  #
  # revert :downsolved status to :initialized in preparation for new upsolve
  resetTreeCliquesForUpSolve!(treel)
  if drawtree
    pdfpath = joinLogPath(dfg,"bt.pdf")
    drawTree(treel, show=false, filepath=pdfpath)
  end

  # queue all the tasks
  resize!(alltasks,length(getCliques(treel)))
  cliqHistories = Dict{Int,Vector{Tuple{DateTime, Int, Function, CliqStateMachineContainer}}}()
  if !isTreeSolved(treel, skipinitialized=true)
    @sync begin
      runtaskmonitor ? (global monitortask = monitorCSMs(treel, alltasks; forceIntExc = true)) : nothing
      # duplicate int i into async (important for concurrency)
      for i in 1:length(getCliques(treel))
        scsym = getCliqFrontalVarIds(getClique(treel, i))
        if length(intersect(scsym, skipcliqids)) == 0
          limthiscsm = filter(x -> (x[1] in scsym), limititercliqs)
          limiter = 0<length(limthiscsm) ? limthiscsm[1][2] : limititers
          injDelay = if injectDelayBefore === nothing
            nothing
          else
            idb = filter((x)->x[1]==i,injectDelayBefore)
            length(idb) == 1 ? idb[1][2] : nothing
          end
          alltasks[i] = @async tryCliqStateMachineSolve!(dfg, treel, i, timeout, oldtree=oldtree, verbose=verbose, verbosefid=verbosefid, drawtree=drawtree, limititers=limiter, downsolve=downsolve, incremental=incremental, delaycliqs=delaycliqs, injectDelayBefore=injDelay, recordcliqs=recordcliqs,  N=N)
        end # if
      end # for
    end # sync
  end # if

  # if record cliques is in use, else skip computational delay
  0 == length(recordcliqs) ? nothing : fetchCliqHistoryAll!(alltasks, cliqHistories)

  return alltasks, cliqHistories
end


text:
@testset "Read MATPOWER" begin
    path = "./test/data/case6ww.m"
    network = PowerFlowNetwork(path, "MATPOWERM")
    @test network.baseMVA == 100
    bus = zeros(6, 17)
    bus[:, 1:13] = [1  3  0   0   0  0  1  1.05  0  230  1  1.05  1.05;
                    2  2  0   0   0  0  1  1.05  0  230  1  1.05  1.05;
                    3  2  0   0   0  0  1  1.07  0  230  1  1.07  1.07;
                    4  1  70  70  0  0  1  1     0  230  1  1.05  0.95;
                    5  1  70  70  0  0  1  1     0  230  1  1.05  0.95;
                    6  1  70  70  0  0  1  1     0  230  1  1.05  0.95;]
    gen = zeros(3, 25)
    gen[:, 1:21] = [1  0   0  100  -100  1.05  100  1  200  50    0  0  0  0  0  0  0  0  0  0  0;
                    2  50  0  100  -100  1.05  100  1  150  37.5  0  0  0  0  0  0  0  0  0  0  0;
                    3  60  0  100  -100  1.07  100  1  180  45    0  0  0  0  0  0  0  0  0  0  0;]
    branch = zeros(11, 21)
    branch[:, 1:13] = [1  2  0.1   0.2   0.04  40  40  40  0  0  1  -360  360;
                       1  4  0.05  0.2   0.04  60  60  60  0  0  1  -360  360;
                       1  5  0.08  0.3   0.06  40  40  40  0  0  1  -360  360;
                       2  3  0.05  0.25  0.06  40  40  40  0  0  1  -360  360;
                       2  4  0.05  0.1   0.02  60  60  60  0  0  1  -360  360;
                       2  5  0.1   0.3   0.04  30  30  30  0  0  1  -360  360;
                       2  6  0.07  0.2   0.05  90  90  90  0  0  1  -360  360;
                       3  5  0.12  0.26  0.05  70  70  70  0  0  1  -360  360;
                       3  6  0.02  0.1   0.02  80  80  80  0  0  1  -360  360;
                       4  5  0.2   0.4   0.08  20  20  20  0  0  1  -360  360;
                       5  6  0.1   0.3   0.06  40  40  40  0  0  1  -360  360;]
    colnames = Dict(
        "bus" => ["ID", "TYPE", "PD", "QD", "GS", "BS", "BUS_AREA", "VM", "VA", "BASE_KV", "ZONE", "VMAX", "VMIN", "LAM_P",
                  "LAM_Q", "MU_VMAX", "MU_VMIN"],
        "gen" => ["ID", "PG", "QG", "QMAX", "QMIN", "VG", "MBASE", "GEN_STATUS", "PMAX", "PMIN", "PC1", "PC2", "QC1MIN",
                  "QC1MAX", "QC2MIN", "QC2MAX", "RAMP_AGC", "RAMP_10", "RAMP_30", "RAMP_Q", "APF", "MU_PMAX", "MU_PMIN",
                  "MU_QMAX", "MU_QMIN"],
        "branch" => ["SRC", "DST", "BR_R", "BR_X", "BR_B", "RATE_A", "RATE_B", "RATE_C", "TAP", "SHIFT", "BR_STATUS",
                     "ANGMIN", "ANGMAX", "PF", "QF", "PT", "QT", "MU_SF", "MU_ST", "MU_ANGMIN", "MU_ANGMAX"],
        "gencost" => ["MODEL", "STARTUP", "SHUTDOWN", "NCOST", "COST"]
    )
    @test network.bus == DataFrame(bus, colnames["bus"])
    @test network.gen == DataFrame(gen, colnames["gen"])
    @test network.branch == DataFrame(branch, colnames["branch"])
    @test nbus(network) == nbus(path, "MATPOWERM")
    @test nbranch(network) == nbranch(path, "MATPOWERM")
    @test nbranch(network; distinct_pair=true) == nbranch(path, "MATPOWERM"; distinct_pair=true)
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
@testset "Gini coefficient" begin
@testset "random data" begin
    n = 20
    data = rand(n)
    curve = CalculatedABC.ABCcurve(data)

    answer = 0.3008352645541903
    @test isapprox(CalculatedABC.gini_coeff(data), answer)
    @test isapprox(CalculatedABC.gini_coeff(curve), answer)
end

@testset "SwissInhabitants data" begin
    data = readdlm("swissinhabitants")[:,1]
    curve = CalculatedABC.ABCcurve(data)

    answer = 0.6294073840902157
    @test isapprox(CalculatedABC.gini_coeff(data), answer)
    @test isapprox(CalculatedABC.gini_coeff(curve), answer)
end
end


text:
using JuMP

function bench_add(n)
    v = JuMP.IndexedVector(Float64, n)
    for i in 1:2:n
        JuMP.addelt!(v, i, 1.0)
    end
    for i in 2:2:n
        JuMP.addelt!(v, i, 1.0)
    end
    for i in 1:3:n
        JuMP.addelt!(v, i, 1.0)
    end
end

bench_add(10)
gc_enable(false)
@time bench_add(1000000)
@time bench_add(1000000)
@time bench_add(1000000)
gc_enable(true)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text: include("quasinewton.jl")


text:
binomial_rv4 = (n, p) -> length(filter(x -> x < p, [rand() for i in range(1, n)]))

head = 1/2 < rand() ? true : false

function pay()
  s = 0
  for i in range(1, 10)
    1/2 < rand() ? s += 1 : s = 0
  end
  s >= 3 ? true : false
end

#println(head)
println(pay())

function ts(t, alpha)
  if t > 0
    a = ts(t-1, alpha)
    return push!(a, alpha*a[end] + randn())
  else
    return [0.0]
  end
end

function ts2(T, alpha)
  series = zeros(T+1)
  for i in 1:T
    series[i+1] = alpha * series[i] + randn()
  end
  return series
end

println(ts(200, 0.9))
println(ts2(200, 0.9))


text:
@testset "The Stiefel-Manifold V(3,4)" begin
  import Random: seed!
  seed!(42);
  M = Stiefel(3,4)
  x = StPoint([1. 0. 0.; 0. 1. 0.; 0. 0. 1.; 0. 0. 0.])
  y = StPoint([0. 0. 0.; 1. 0. 0.; 0. 1. 0.; 0. 0. 1.])
  z = StPoint([-0.851712 0.0231522 0.502726; 0.509415 -0.0752401 0.741671; -0.0955954 -0.918939 -0.212126; -0.0771071 0.386462 -0.39012])
  η = StTVector([0.0 0.359296 -0.387003; -0.359296 0.0 0.0453968; 0.387003 -0.0453968 0.0; 0.0037522 -0.176039 0.638079])
  ω = StTVector([0.0 -0.162363 -0.639113; 0.162363 0.0 0.153497; 0.639113 -0.153497 0.0; 0.133987 0.155025 -0.203072])
  ξ = StTVector([0.0 -0.306164 0.31943; 0.306164 0.0 -0.467336; -0.31943 0.467336 0.0; -0.0815496 0.193606 -0.357094])
  ν = randomTVector(M,x)
  x2 = addNoise(M,x,0.5)
  x3 = addNoise(M,x,500.0)
  w = randomMPoint(M)
  s = pi

  # Test Stiefel
  @test_throws ErrorException Stiefel(4,3)
  # Test Minus operator
  @test getValue( (-η) ) == getValue(-η)
  # Test Cauchy-Schwarz-Inequation
  @test dot(M,x,ω,η)^2 <= norm(M,x,ω)^2*norm(M,x,η)^2
  @test dot(M,x,ω,η) <= norm(M,x,ω)*norm(M,x,η)
  # Test dot
  @test dot(M,x,ω,η) ≈ dot(M,x,η,ω) atol = 10.0^(-16)
  @test dot(M,x,ω,η+ξ) ≈ dot(M,x,ω,η) + dot(M,x,ω,ξ) atol = 10.0^(-15)
  @test dot(M,x,ω,s*η) ≈ s*dot(M,x,ω,η) atol = 10.0^(-14)
  # Test norm
  @test norm(M,x,ω)^2 ≈ dot(M,x,ω,ω) atol = 10.0^(-14)
  @test norm(M,x,ω) ≈ sqrt(dot(M,x,ω,ω)) atol = 10.0^(-16)
  # Test randomMPoint and randomTVector
  @test norm(transpose(getValue(w))*getValue(w) - one(transpose(getValue(w))*getValue(w))) ≈ 0 atol = 10.0^(-14)
  @test norm(transpose(getValue(x))*getValue(ν) + transpose(transpose(getValue(x))*getValue(ν))) ≈ 0 atol = 10.0^(-15)
  # Test addNoise
  @test norm(transpose(getValue(x2))*getValue(x2) - one(transpose(getValue(x2))*getValue(x2))) ≈ 0 atol=10.0^(-15)
  @test norm(transpose(getValue(x3))*getValue(x3) - one(transpose(getValue(x3))*getValue(x3))) ≈ 0 atol=10.0^(-12)
  #Test retraction
  @test norm(getValue(inverseRetractionQR(M,x,retractionQR(M,x,ω))) - getValue(ω)) ≈ 0 atol = 10.0^(-14)
  @test norm(getValue(inverseRetraction(M,x,retraction(M,x,ω))) - getValue(ω)) ≈ 0 atol = 10.0^(-14)
  @test norm(getValue(inverseRetractionPolar(M,x,retractionPolar(M,x,ω))) - getValue(ω)) ≈ 0 atol = 10.0^(-14)
  @test norm(transpose(getValue(retractionQR(M,x,ω))) * getValue(retractionQR(M,x,ω)) - one(transpose(getValue(x))*getValue(x))) ≈ 0 atol = 10.0^(-15)
  @test norm(transpose(getValue(retraction(M,x,ω))) * getValue(retraction(M,x,ω)) - one(transpose(getValue(x))*getValue(x))) ≈ 0 atol = 10.0^(-15)
  @test norm(transpose(getValue(retractionPolar(M,x,ω))) * getValue(retractionPolar(M,x,ω)) - one(transpose(getValue(x))*getValue(x))) ≈ 0 atol = 10.0^(-15)
  #Test manifoldDimension
  @test manifoldDimension(M) == manifoldDimension(x)
  @test manifoldDimension(M) == 6
  @test manifoldDimension(x) == 6
  # Test distance – not yet implemented
  @test_throws DomainError distance(M,x,y)
  # Test parallelTransport
  @test norm(getValue(parallelTransport(M,x,z,η)) - getValue(projection(M,z,getValue(η)))) ≈ 0 atol = 10.0^(-16)
  # Test zeroTVector
  @test norm(M,x,zeroTVector(M,x)) ≈ 0 atol = 10.0^(-16)
  # Test validateMPoint and validateTVector
  ynot1 = StPoint([1. 0. 0.; 0. 1. 0.; 0. 0. 1.; 0. 0. 0.; 0. 0. 0.])
  ynot2 = StPoint([1. 0. 0. 0.; 0. 1. 0. 0.; 0. 0. 1. 0.; 0. 0. 0. 1.])
  ynot3 = StPoint([1. 0. 0.; 0. 1. 0.; 0. 0. 0; 0. 0. 0.])
  ynot4 = StPoint([1. 2. 0.; 0. 1. 0.; 0. 4. 0; 5. 0. 6.])
  ξnot1 = StTVector([1. 0. 0.; 0. 1. 0.; 0. 0. 1.; 0. 0. 0.; 0. 0. 0.])
  ξnot2 = StTVector([1. 0. 0. 0.; 0. 1. 0. 0.; 0. 0. 1. 0.; 0. 0. 0. 1.])
  ξnot3 = StTVector([ 1. 2. 3.; 4. 5. 6.; 7. 8. 9.; 0. 0. 0.])
  @test_throws ErrorException validateMPoint(M,ynot1)
  @test_throws ErrorException validateMPoint(M,ynot2)
  @test_throws ErrorException validateMPoint(M,ynot3)
  @test_throws ErrorException validateMPoint(M,ynot4)
  @test_throws ErrorException validateTVector(M,x,ξnot1)
  @test_throws ErrorException validateTVector(M,x,ξnot2)
  @test_throws ErrorException validateTVector(M,x,ξnot3)

  N = Stiefel{Complex{Float64}}(3,4)
  xcompl = StPoint{Complex{Float64}}([-0.761286+0.0462087im -0.18395+0.0566532im 0.525627+0.265715im; -0.104276+0.262847im -0.55139+0.157097im -0.582766+0.202596im; -0.324463+0.199309im -0.041674-0.427848im -0.202927+0.0456294im; -0.327554-0.29335im -0.078241+0.66583im -0.232484-0.418496im])
  ycompl = StPoint{Complex{Float64}}(Complex{Float64}[-0.503249+0.0451954im -0.106808-0.271355im -0.422033-0.674388im; 0.263371-0.361562im -0.212757-0.435859im -0.208488+0.233779im; 0.598422-0.153485im 0.289058+0.159355im -0.254219-0.321635im; -0.240499-0.32419im 0.252482+0.712056im -0.314129+0.0470639im])
  zcompl = StPoint{Complex{Float64}}(Complex{Float64}[-0.385214-0.00723968im -0.0221375-0.419541im -0.518744-0.430094im; -0.286931+0.256567im -0.174173+0.01588im 0.00148874-0.527176im; 0.516436+0.382863im 0.592404-0.321075im -0.062389-0.136172im; 0.0362203-0.537401im 0.322842+0.484405im -0.49128-0.0649267im])
  ηcompl = StTVector{Complex{Float64}}(Complex{Float64}[-0.144443+0.050361im 0.0883068-0.262202im -0.0739903+0.145955im; -0.117226-0.00340324im 0.0298396+0.182582im -0.0763624-0.314719im; 0.423676-0.0952191im 0.0747197-0.168339im -0.448123+0.0629099im; -0.227158+0.141742im 0.437945-0.0236934im 0.0141424+0.170016im])
  ωcompl = StTVector{Complex{Float64}}(Complex{Float64}[0.0340014+0.025003im -0.175668-0.373111im 0.0973745-0.27843im; 0.0580946-0.252544im -0.13527+0.203748im -0.121749+0.270144im; 0.0677755-0.266317im -0.198544+0.306502im 0.0923344+0.0506648im; -0.151337-0.418159im 0.208301+0.0321241im 0.243117+0.0715284im])
  ξcompl = StTVector{Complex{Float64}}(Complex{Float64}[0.0416013-0.190526im 0.114072-0.0248113im 0.0703455+0.367732im; 0.124414+0.38413im -0.225464-0.120304im 0.39404-0.00423936im; 0.0805494-0.322981im 0.2325+0.234871im -0.151441-0.20267im; -0.0148454-0.129968im -0.314709+0.00379301im -0.0887356-0.128297im])
  νcompl = randomTVector(N,xcompl)
  wcompl = randomMPoint(N)
  x2compl = addNoise(N,xcompl,0.5)
  x3compl = addNoise(N,xcompl,500.0)

  # Test Stiefel
  @test_throws ErrorException Stiefel{Complex{Float64}}(4,3)
  # Test Minus operator
  @test getValue( (-ηcompl) ) == getValue(-ηcompl)
  # Test Cauchy-Schwarz-Inequation
  @test abs(dot(N,xcompl,ωcompl,ηcompl)^2) <= abs(norm(N,xcompl,ωcompl)^2*norm(N,xcompl,ηcompl)^2)
  @test dot(N,xcompl,ωcompl,ηcompl) <= norm(N,xcompl,ωcompl)*norm(N,xcompl,ηcompl)
  # Test dot
  @test dot(N,xcompl,ωcompl,ηcompl) ≈ dot(N,xcompl,ηcompl,ωcompl) atol = 10.0^(-16)
  @test dot(N,xcompl,ωcompl,ηcompl+ξcompl) ≈ dot(N,xcompl,ωcompl,ηcompl) + dot(N,xcompl,ωcompl,ξcompl) atol = 10.0^(-15)
  @test dot(N,xcompl,ωcompl,s*ηcompl) ≈ s*dot(N,xcompl,ωcompl,ηcompl) atol = 10.0^(-14)
  # Test norm
  @test norm(N,xcompl,ωcompl)^2 ≈ dot(N,xcompl,ωcompl,ωcompl) atol = 10.0^(-14)
  @test norm(N,xcompl,ωcompl) ≈ sqrt(dot(N,xcompl,ωcompl,ωcompl)) atol = 10.0^(-14)
  # Test randomMPoint and randomTVector
  @test norm(getValue(wcompl)'*getValue(wcompl) - one(getValue(wcompl)'*getValue(wcompl))) ≈ 0 atol = 10.0^(-14)
  @test norm(getValue(xcompl)'*getValue(νcompl) + (getValue(xcompl)'*getValue(νcompl))') ≈ 0 atol = 10.0^(-5)
  # Test addNoise
  @test norm(getValue(x2compl)'*getValue(x2compl) - one(getValue(x2compl)'*getValue(x2compl))) ≈ 0 atol=10.0^(-5)
  @test norm(getValue(x3compl)'*getValue(x3compl) - one(getValue(x3compl)'*getValue(x3compl))) ≈ 0 atol=10.0^(-5)
  # Test parallelTransport
  @test norm(getValue(parallelTransport(N,xcompl,zcompl,ηcompl)) - getValue(projection(N,zcompl,getValue(ηcompl)))) ≈ 0 atol = 10.0^(-16)

  #Test manifoldDimension
  @test manifoldDimension(N) == manifoldDimension(wcompl)
  @test manifoldDimension(N) == 15
  @test manifoldDimension(xcompl) == 15
end


text:
abstract type AbstractMouseState end

"""
    MouseState{T<:AbstractMouseState}

Describes a mouse state change.
Fields:
- `typ`: Symbol describing the mouse state
- `t`: Time of the event
- `pos`: Mouse position
- `tprev`: Time of previous event
- `prev`: Previous mouse position
"""
struct MouseState{T<:AbstractMouseState}
    typ::T
    t::Float64
    pos::Point2f0
    tprev::Float64
    prev::Point2f0
end

mousestates = (:MouseOut, :MouseEnter, :MouseOver,
    :MouseLeftDown, :MouseRightDown, :MouseMiddleDown,
    :MouseLeftUp, :MouseRightUp, :MouseMiddleUp,
    :MouseLeftDragStart, :MouseRightDragStart, :MouseMiddleDragStart,
    :MouseLeftDrag, :MouseRightDrag, :MouseMiddleDrag,
    :MouseLeftDragStop, :MouseRightDragStop, :MouseMiddleDragStop,
    :MouseLeftClick, :MouseRightClick, :MouseMiddleClick,
    :MouseLeftDoubleclick, :MouseRightDoubleclick, :MouseMiddleDoubleclick,
    :MouseDownOutside
    )

for statetype in mousestates
    onfunctionname = Symbol("on" * lowercase(String(statetype)))
    @eval begin
        struct $statetype <: AbstractMouseState end

        """
        Executes the function f whenever the `Node{MouseState}` statenode transitions
        to `$($statetype)`.
        """
        function $onfunctionname(f, statenode::Node{MouseState})
            on(statenode) do state
                if state.typ isa $statetype
                    f(state)
                end
            end
        end
        export $onfunctionname
    end
end

function Base.show(io::IO, ms::MouseState{T}) where T
    print(io, "$T(t: $(ms.t), pos: $(ms.pos[1]), $(ms.pos[2]), tprev: $(ms.tprev), prev: $(ms.prev[1]), $(ms.prev[2]))")
end

"""
    addmousestate!(scene, elements...)

Returns an `Observable{MouseState}` which is triggered by all mouse
interactions with the `scene` and optionally restricted to all given
plot objects in `elements`.

To react to mouse events, use the onmouse... handlers.

Example:

```
mousestate = addmousestate!(scene, scatterplot)

onmouseleftclick(mousestate) do state
    # do something with the mousestate
end
```
"""
function addmousestate!(scene, elements...)

    Mouse = AbstractPlotting.Mouse
    dblclick_max_interval = 0.2

    mousestate = Node{MouseState}(MouseState(MouseOut(), 0.0, Point2f0(0, 0), 0.0, Point2f0(0, 0)))

    is_mouse_over_relevant_area() = isempty(elements) ? AbstractPlotting.is_mouseinside(scene) : mouseover(scene, elements...)

    # initialize state variables
    last_mousestate = Ref{Mouse.DragEnum}(events(scene).mousedrag[])
    prev = Ref(mouseposition(AbstractPlotting.rootparent(scene)))
    mouse_downed_inside = Ref(false)
    mouse_downed_button = Ref{Optional{Mouse.Button}}(nothing)
    drag_ongoing = Ref(false)
    mouse_was_inside = Ref(false)
    tprev = Ref(0.0)
    t_last_click = Ref(0.0)
    b_last_click = Ref{Optional{Mouse.Button}}(nothing)
    last_click_was_double = Ref(false)

    # react to mouse position changes
    on(events(scene).mouseposition) do mp

        t = time()
        pos = mouseposition(AbstractPlotting.rootparent(scene))
        mouse_inside = is_mouse_over_relevant_area()

        # movement while mouse is pressed
        if last_mousestate[] == Mouse.pressed
            # must have been a registered drag (otherwise could have come from outside)
            if drag_ongoing[]
                event = @match mouse_downed_button[] begin
                    Mouse.left => MouseLeftDrag()
                    Mouse.right => MouseRightDrag()
                    Mouse.middle => MouseMiddleDrag()
                    x => error("No recognized mouse button $x")
                end
                mousestate[] = MouseState(event, t, pos, tprev[], prev[])
            end
        # mouse moved while just having been pressed down
        elseif last_mousestate[] == Mouse.down
            # mouse must have been downed inside
            # that means a drag started
            if mouse_downed_inside[]
                drag_ongoing[] = true
                event = @match mouse_downed_button[] begin
                    Mouse.left => MouseLeftDragStart()
                    Mouse.right => MouseRightDragStart()
                    Mouse.middle => MouseMiddleDragStart()
                    x => error("No recognized mouse button $x")
                end
                mousestate[] = MouseState(event, t, pos, tprev[], prev[])

                event = @match mouse_downed_button[] begin
                    Mouse.left => MouseLeftDrag()
                    Mouse.right => MouseRightDrag()
                    Mouse.middle => MouseMiddleDrag()
                    x => error("No recognized mouse button $x")
                end
                mousestate[] = MouseState(event, t, pos, tprev[], prev[])
            end
        else
            if mouse_inside
                if mouse_was_inside[]
                    mousestate[] = MouseState(MouseOver(), t, pos, tprev[], prev[])
                else
                    mousestate[] = MouseState(MouseEnter(), t, pos, tprev[], prev[])
                end
            else
                if mouse_was_inside[]
                    mousestate[] = MouseState(MouseOut(), t, pos, tprev[], prev[])
                end
            end
        end

        mouse_was_inside[] = mouse_inside
        prev[] = pos
        tprev[] = t
    end

    # react to mouse button changes
    on(events(scene).mousedrag) do mousedrag

        t = time()
        pos = prev[]

        pressed_buttons = events(scene).mousebuttons[]

        # we only need to handle mousedown and mouseup
        # pressed and not pressed are redundant events with mouse position changes
        if mousedrag == Mouse.down
            if length(pressed_buttons) == 1
                button = only(pressed_buttons)
                mouse_downed_button[] = button

                if mouse_was_inside[]
                    event = @match mouse_downed_button[] begin
                        Mouse.left => MouseLeftDown()
                        Mouse.right => MouseRightDown()
                        Mouse.middle => MouseMiddleDown()
                        x => error("No recognized mouse button $x")
                    end
                    mousestate[] = MouseState(event, t, pos, tprev[], prev[])
                    mouse_downed_inside[] = true
                else
                    mouse_downed_inside[] = false
                    mousestate[] = MouseState(MouseDownOutside(), t, pos, tprev[], prev[])
                end
            end
        elseif mousedrag == Mouse.up
            # only register up events and clicks if the upped button matches
            # the recorded downed one
            # and it can't be nothing (if the first up event comes from outside)

            downed_button_missing_from_pressed = !(mouse_downed_button[] in pressed_buttons)
            some_mouse_button_had_been_downed = !isnothing(mouse_downed_button[])

            if downed_button_missing_from_pressed && some_mouse_button_had_been_downed

                if drag_ongoing[]
                    event = @match mouse_downed_button[] begin
                        Mouse.left => MouseLeftDragStop()
                        Mouse.right => MouseRightDragStop()
                        Mouse.middle => MouseMiddleDragStop()
                        x => error("No recognized mouse button $x")
                    end
                    mousestate[] = MouseState(event, t, pos, tprev[], prev[])
                    drag_ongoing[] = false

                    if mouse_was_inside[]
                        # up after drag done over element
                        event = @match mouse_downed_button[] begin
                            Mouse.left => MouseLeftUp()
                            Mouse.right => MouseRightUp()
                            Mouse.middle => MouseMiddleUp()
                            x => error("No recognized mouse button $x")
                        end

                        mousestate[] = MouseState(event, t, pos, tprev[], prev[])
                    else
                        # mouse could be not over elements after drag is over
                        mousestate[] = MouseState(MouseOut(), t, pos, tprev[], prev[])
                    end
                else
                    if mouse_was_inside[]
                        dt_last_click = t - t_last_click[]
                        t_last_click[] = t

                        # guard against mouse coming in from outside, then mouse upping
                        if mouse_downed_inside[]
                            if dt_last_click < dblclick_max_interval && !last_click_was_double[] &&
                                    mouse_downed_button[] == b_last_click[]

                                event = @match mouse_downed_button[] begin
                                    Mouse.left => MouseLeftDoubleclick()
                                    Mouse.right => MouseRightDoubleclick()
                                    Mouse.middle => MouseMiddleDoubleclick()
                                    x => error("No recognized mouse button $x")
                                end
                                mousestate[] = MouseState(event, t, pos, tprev[], prev[])
                                last_click_was_double[] = true
                            else
                                event = @match mouse_downed_button[] begin
                                    Mouse.left => MouseLeftClick()
                                    Mouse.right => MouseRightClick()
                                    Mouse.middle => MouseMiddleClick()
                                    x => error("No recognized mouse button $x")
                                end
                                mousestate[] = MouseState(event, t, pos, tprev[], prev[])
                                last_click_was_double[] = false
                            end
                            # save what type the last downed button was
                            b_last_click[] = mouse_downed_button[]
                        end
                        mouse_downed_inside[] = false

                        # up after click
                        event = @match mouse_downed_button[] begin
                            Mouse.left => MouseLeftUp()
                            Mouse.right => MouseRightUp()
                            Mouse.middle => MouseMiddleUp()
                            x => error("No recognized mouse button $x")
                        end

                        mousestate[] = MouseState(event, t, pos, tprev[], prev[])
                    end
                end

            end
        end

        last_mousestate[] = mousedrag
        tprev[] = t
    end

    mousestate
end


text:
##############################   Coulomb  ###################################
# mutable struct Coulomb <: EmpiricalPotential
#     q_1::Float64
#     q_2::Float64
#     ϵ0::Float64
# end

# function Coulomb()
#     #ToDO
#     return Coulomb(1.0, 1.0, 1.0)
# end

# function get_trainable_params(c::Coulomb)
#     p = Parameter{}(())
#     return p
# end

# function get_nontrainable_params(c::Coulomb)
#     return (q_1 = c.q_1, q_2 = c.q_2, ϵ0 = c.ϵ0)
# end

# ############################# Energies ##########################################

# function potential_energy(r::Vector{<:Real}, p::Coulomb)
#     return p.q_1 * p.q_2 / (4.0 * π * p.ϵ0 * norm(r))
# end

# ############################### Forces ##########################################

# function force(r::Vector{<:Real}, p::Coulomb)
#     d = norm(r)
#     return p.q_1 * p.q_2 / (4.0 * π * p.ϵ0 * d^2) .* r ./ d
# end

# ##############################   Gradients  ###################################

# function grad_potential_energy(r::Vector{<:Real}, p::Coulomb)
#     println("The Coulomb potential has no trainable parameters")
#     return (dpdnull = 0.0)
# end
# function grad_force(r:: Vector{<:Real}, p::Coulomb)
#     println("The Coulomb potential has no trainable parameters")
#     return (dfdnull =  0 .* r ./ d )
# end

# function grad_virial(r::Vector{<:Real}, p::Coulomb)
#     println("The Coulomb potential has no trainable parameters")
#     return (dvdnull =  0 .* r ./ d )
# end


text:
#=
Gaussian model:
\[\vec{d}(\vec{p}) = -\mathrm{i}\left(\frac{1}{\pi\alpha}\right)^{3/4}
\frac{\vec{p}}{\alpha}e^{-p^2/2\alpha}\]
=#
d_gaussian(p,α) = -im*(1/(π*α))^(3/4) * p/α*exp(p^2/(2α))

#=
Hydrogenlike atom:
\[\vec{d}(\vec{p}) = -\mathrm{i}\left(\frac{2^{7/2}\alpha^{5/4}}{\pi}\right)
\frac{\vec{p}}{(p^2+\alpha)^3}\]
=#
d_hyd(p,α) = -im*(2^(7/2)*α^(5/4)/π)*p/((norm(p)^2+α)^3)
hyd_α(Ip) = 2Ip

hydrogen_like(Ip) = Ip, p -> d_hyd(p, hyd_α(Ip))
hydrogen = hydrogen_like(0.5)

function yukawa_fourier(q, a, b, λ)
    # Eq. (47) of
    #
    # - D B Milo\vsevi\'c, Paulus, G. G., Bauer, D., & Becker,
    #   W. (2006). Above-Threshold Ionization By Few-Cycle Pulses. Journal
    #   of Physics B: Atomic, Molecular and Optical Physics, 39(14),
    #   203–262. http://dx.doi.org/10.1088/0953-4075/39/14/r01
    C = norm(q)^2 + λ^2
    -(2b*λ + a*C)/(2π^2*C^2)
end

coulomb_fourier(q, Z=1) = yukawa_fourier(q, Z, 0, 0)

export d_gaussian, d_hyd, hyd_α, hydrogen_like, hydrogen,
    yukawa_fourier, coulomb_fourier


text:
"""
# module Inputs

# Examples

```julia
julia>
```
"""
module Inputs

using AbInitioSoftwareBase.Inputs: Input, InputEntry, Namelist, Card, Setter, groupname
using Compat: only, isnothing

export optionof,
    optionpool,
    groupname,
    required_namelists,
    optional_namelists,
    required_cards,
    optional_cards,
    allnamelists,
    allcards

"""
    dropdefault(nml::Namelist)

Return an `AbstractDict` of non-default values of a `Namelist`.
"""
function dropdefault(nml::Namelist)
    default = typeof(nml)()  # Create a `Namelist` with all default values
    # Compare `default` with `nml`, discard the same values
    result = Iterators.filter(item -> item.second != getfield(default, item.first), nml)
    # for (drivingarg, passivearg) in _coupledargs(typeof(nml))
    # rule
    # end
    if isempty(result)
        @info "Every entry in the namelist is the default value!"
    end
    return Dict{Symbol,Any}(result)
end

"""
    optionof(x::Card)

Return a `String` representing the option for `Card` `x`.

!!! warning
    Do not use `x.option` to access a `Card`'s `option`.
"""
optionof(card::Card) = hasfield(typeof(card), :option) ? getfield(card, :option) : nothing

"""
    optionpool(T::Type{<:Card})

Return the allowed options for `Card` `T`.

# Examples
```julia
julia> optionpool(AtomicPositionsCard)
("alat", "bohr", "angstrom", "crystal", "crystal_sg")

julia> optionpool(CellParametersCard)
("alat", "bohr", "angstrom")

julia> optionpool(SpecialKPointsCard)
("tpiba", "crystal", "tpiba_b", "crystal_b", "tpiba_c", "crystal_c")
```
"""
function optionpool end

"Represent input files of executables (such as `pw.x` and `cp.x`)."
abstract type QuantumESPRESSOInput <: Input end

function allnamelists end

function allcards end

function required_namelists end

function optional_namelists end

function required_cards end

function optional_cards end

struct VerbositySetter <: Setter
    v::String
    function VerbositySetter(v)
        @assert v in ("high", "low")
        return new(v)
    end
end

include("crystallography.jl")
include("PWscf/PWscf.jl")
# include("CP/CP.jl")
include("PHonon/PHonon.jl")

end


text:
using Base.Test
using Lint

lintpkg("VirtualTestDrive")

using VirtualTestDrive
using LibExpat

include("test_io.jl")
include("test_virescsvs.jl")

max_stars_repo_name: sisl/VirtualTestDrive

max_issues_repo_name: sisl/VirtualTestDrive

max_forks_repo_name: sisl/VirtualTestDrive