text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# macro wrappers for various reflection functions

import Base: typesof, insert!

separate_kwargs(args...; kwargs...) = (args, kwargs.data)

"""
Transform a dot expression into one where each argument has been replaced by a
variable "xj" (with j an integer from 1 to the returned i).
The list `args` contains the original arguments that have been replaced.
"""
function recursive_dotcalls!(ex, args, i=1)
    if !(ex isa Expr) || ((ex.head !== :. || !(ex.args[2] isa Expr)) &&
                          (ex.head !== :call || string(ex.args[1])[1] != '.'))
        newarg = Symbol('x', i)
        if Meta.isexpr(ex, :...)
            push!(args, only(ex.args))
            return Expr(:..., newarg), i+1
        else
            push!(args, ex)
            return newarg, i+1
        end
    end
    (start, branches) = ex.head === :. ? (1, ex.args[2].args) : (2, ex.args)
    for j in start:length(branches)
        branch, i = recursive_dotcalls!(branches[j], args, i)
        branches[j] = branch
    end
    return ex, i
end

function gen_call_with_extracted_types(__module__, fcn, ex0, kws=Expr[])
    if isa(ex0, Expr)
        if ex0.head === :do && Meta.isexpr(get(ex0.args, 1, nothing), :call)
            if length(ex0.args) != 2
                return Expr(:call, :error, "ill-formed do call")
            end
            i = findlast(a->(Meta.isexpr(a, :kw) || Meta.isexpr(a, :parameters)), ex0.args[1].args)
            args = copy(ex0.args[1].args)
            insert!(args, (isnothing(i) ? 2 : i+1), ex0.args[2])
            ex0 = Expr(:call, args...)
        end
        if ex0.head === :. || (ex0.head === :call && ex0.args[1] !== :.. && string(ex0.args[1])[1] == '.')
            codemacro = startswith(string(fcn), "code_")
            if codemacro && ex0.args[2] isa Expr
                # Manually wrap a dot call in a function
                args = Any[]
                ex, i = recursive_dotcalls!(copy(ex0), args)
                xargs = [Symbol('x', j) for j in 1:i-1]
                dotfuncname = gensym("dotfunction")
                dotfuncdef = Expr(:local, Expr(:(=), Expr(:call, dotfuncname, xargs...), ex))
                return quote
                    $(esc(dotfuncdef))
                    local args = typesof($(map(esc, args)...))
                    $(fcn)($(esc(dotfuncname)), args; $(kws...))
                end
            elseif !codemacro
                fully_qualified_symbol = true # of the form A.B.C.D
                ex1 = ex0
                while ex1 isa Expr && ex1.head === :.
                    fully_qualified_symbol = (length(ex1.args) == 2 &&
                                              ex1.args[2] isa QuoteNode &&
                                              ex1.args[2].value isa Symbol)
                    fully_qualified_symbol || break
                    ex1 = ex1.args[1]
                end
                fully_qualified_symbol &= ex1 isa Symbol
                if fully_qualified_symbol
                    return quote
                        local arg1 = $(esc(ex0.args[1]))
                        if isa(arg1, Module)
                            $(if string(fcn) == "which"
                                  :(which(arg1, $(ex0.args[2])))
                              else
                                  :(error("expression is not a function call"))
                              end)
                        else
                            local args = typesof($(map(esc, ex0.args)...))
                            $(fcn)(Base.getproperty, args)
                        end
                    end
                else
                    return Expr(:call, :error, "dot expressions are not lowered to "
                                * "a single function call, so @$fcn cannot analyze "
                                * "them. You may want to use Meta.@lower to identify "
                                * "which function call to target.")
                end
            end
        end
        if any(a->(Meta.isexpr(a, :kw) || Meta.isexpr(a, :parameters)), ex0.args)
            return quote
                local arg1 = $(esc(ex0.args[1]))
                local args, kwargs = $separate_kwargs($(map(esc, ex0.args[2:end])...))
                $(fcn)(Core.kwfunc(arg1),
                       Tuple{typeof(kwargs), Core.Typeof(arg1), map(Core.Typeof, args)...};
                       $(kws...))
            end
        elseif ex0.head === :call
            return Expr(:call, fcn, esc(ex0.args[1]),
                        Expr(:call, typesof, map(esc, ex0.args[2:end])...),
                        kws...)
        elseif ex0.head === :(=) && length(ex0.args) == 2
            lhs, rhs = ex0.args
            if isa(lhs, Expr)
                if lhs.head === :(.)
                    return Expr(:call, fcn, Base.setproperty!,
                                Expr(:call, typesof, map(esc, lhs.args)..., esc(rhs)), kws...)
                elseif lhs.head === :ref
                    return Expr(:call, fcn, Base.setindex!,
                                Expr(:call, typesof, esc(lhs.args[1]), esc(rhs), map(esc, lhs.args[2:end])...), kws...)
                end
            end
        elseif ex0.head === :vcat || ex0.head === :typed_vcat
            if ex0.head === :vcat
                f, hf = Base.vcat, Base.hvcat
                args = ex0.args
            else
                f, hf = Base.typed_vcat, Base.typed_hvcat
                args = ex0.args[2:end]
            end
            if any(a->isa(a,Expr) && a.head === :row, args)
                rows = Any[ (isa(x,Expr) && x.head === :row ? x.args : Any[x]) for x in args ]
                lens = map(length, rows)
                return Expr(:call, fcn, hf,
                            Expr(:call, typesof,
                                 (ex0.head === :vcat ? [] : Any[esc(ex0.args[1])])...,
                                 Expr(:tuple, lens...),
                                 map(esc, vcat(rows...))...), kws...)
            else
                return Expr(:call, fcn, f,
                            Expr(:call, typesof, map(esc, ex0.args)...), kws...)
            end
        else
            for (head, f) in (:ref => Base.getindex, :hcat => Base.hcat, :(.) => Base.getproperty, :vect => Base.vect, Symbol("'") => Base.adjoint, :typed_hcat => Base.typed_hcat, :string => string)
                if ex0.head === head
                    return Expr(:call, fcn, f,
                                Expr(:call, typesof, map(esc, ex0.args)...), kws...)
                end
            end
        end
    end
    if isa(ex0, Expr) && ex0.head === :macrocall # Make @edit @time 1+2 edit the macro by using the types of the *expressions*
        return Expr(:call, fcn, esc(ex0.args[1]), Tuple{#=__source__=#LineNumberNode, #=__module__=#Module, Any[ Core.Typeof(a) for a in ex0.args[3:end] ]...}, kws...)
    end

    ex = Meta.lower(__module__, ex0)
    if !isa(ex, Expr)
        return Expr(:call, :error, "expression is not a function call or symbol")
    end

    exret = Expr(:none)
    if ex.head === :call
        if any(e->(isa(e, Expr) && e.head === :(...)), ex0.args) &&
            (ex.args[1] === GlobalRef(Core,:_apply_iterate) ||
             ex.args[1] === GlobalRef(Base,:_apply_iterate))
            # check for splatting
            exret = Expr(:call, ex.args[2], fcn,
                        Expr(:tuple, esc(ex.args[3]),
                            Expr(:call, typesof, map(esc, ex.args[4:end])...)))
        else
            exret = Expr(:call, fcn, esc(ex.args[1]),
                         Expr(:call, typesof, map(esc, ex.args[2:end])...), kws...)
        end
    end
    if ex.head === :thunk || exret.head === :none
        exret = Expr(:call, :error, "expression is not a function call, "
                                  * "or is too complex for @$fcn to analyze; "
                                  * "break it down to simpler parts if possible. "
                                  * "In some cases, you may want to use Meta.@lower.")
    end
    return exret
end

"""
Same behaviour as gen_call_with_extracted_types except that keyword arguments
of the form "foo=bar" are passed on to the called function as well.
The keyword arguments must be given before the mandatory argument.
"""
function gen_call_with_extracted_types_and_kwargs(__module__, fcn, ex0)
    kws = Expr[]
    arg = ex0[end] # Mandatory argument
    for i in 1:length(ex0)-1
        x = ex0[i]
        if x isa Expr && x.head === :(=) # Keyword given of the form "foo=bar"
            if length(x.args) != 2
                return Expr(:call, :error, "Invalid keyword argument: $x")
            end
            push!(kws, Expr(:kw, x.args[1], x.args[2]))
        else
            return Expr(:call, :error, "@$fcn expects only one non-keyword argument")
        end
    end
    return gen_call_with_extracted_types(__module__, fcn, arg, kws)
end

for fname in [:which, :less, :edit, :functionloc]
    @eval begin
        macro ($fname)(ex0)
            gen_call_with_extracted_types(__module__, $(Expr(:quote, fname)), ex0)
        end
    end
end

macro which(ex0::Symbol)
    ex0 = QuoteNode(ex0)
    return :(which($__module__, $ex0))
end

for fname in [:code_warntype, :code_llvm, :code_native]
    @eval begin
        macro ($fname)(ex0...)
            gen_call_with_extracted_types_and_kwargs(__module__, $(Expr(:quote, fname)), ex0)
        end
    end
end

macro code_typed(ex0...)
    thecall = gen_call_with_extracted_types_and_kwargs(__module__, :code_typed, ex0)
    quote
        local results = $thecall
        length(results) == 1 ? results[1] : results
    end
end

macro code_lowered(ex0...)
    thecall = gen_call_with_extracted_types_and_kwargs(__module__, :code_lowered, ex0)
    quote
        local results = $thecall
        length(results) == 1 ? results[1] : results
    end
end

"""
    @functionloc

Applied to a function or macro call, it evaluates the arguments to the specified call, and
returns a tuple `(filename,line)` giving the location for the method that would be called for those arguments.
It calls out to the `functionloc` function.
"""
:@functionloc

"""
    @which

Applied to a function or macro call, it evaluates the arguments to the specified call, and
returns the `Method` object for the method that would be called for those arguments. Applied
to a variable, it returns the module in which the variable was bound. It calls out to the
`which` function.
"""
:@which

"""
    @less

Evaluates the arguments to the function or macro call, determines their types, and calls the `less`
function on the resulting expression.
"""
:@less

"""
    @edit

Evaluates the arguments to the function or macro call, determines their types, and calls the `edit`
function on the resulting expression.
"""
:@edit

"""
    @code_typed

Evaluates the arguments to the function or macro call, determines their types, and calls
[`code_typed`](@ref) on the resulting expression. Use the optional argument `optimize` with

    @code_typed optimize=true foo(x)

to control whether additional optimizations, such as inlining, are also applied.
"""
:@code_typed

"""
    @code_lowered

Evaluates the arguments to the function or macro call, determines their types, and calls
[`code_lowered`](@ref) on the resulting expression.
"""
:@code_lowered

"""
    @code_warntype

Evaluates the arguments to the function or macro call, determines their types, and calls
[`code_warntype`](@ref) on the resulting expression.
"""
:@code_warntype

"""
    @code_llvm

Evaluates the arguments to the function or macro call, determines their types, and calls
[`code_llvm`](@ref) on the resulting expression.
Set the optional keyword arguments `raw`, `dump_module`, `debuginfo`, `optimize`
by putting them and their value before the function call, like this:

    @code_llvm raw=true dump_module=true debuginfo=:default f(x)
    @code_llvm optimize=false f(x)

`optimize` controls whether additional optimizations, such as inlining, are also applied.
`raw` makes all metadata and dbg.* calls visible.
`debuginfo` may be one of `:source` (default) or `:none`,  to specify the verbosity of code comments.
`dump_module` prints the entire module that encapsulates the function.
"""
:@code_llvm

"""
    @code_native

Evaluates the arguments to the function or macro call, determines their types, and calls
[`code_native`](@ref) on the resulting expression.

Set the optional keyword argument `debuginfo` by putting it before the function call, like this:

    @code_native debuginfo=:default f(x)

`debuginfo` may be one of `:source` (default) or `:none`, to specify the verbosity of code comments.
"""
:@code_native

max_stars_repo_name: udesou/julia

max_issues_repo_name: udesou/julia

max_forks_repo_name: udesou/julia


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

__precompile__(true)

module IntervalArithmetic

import CRlibm

using FastRounding
using SetRounding

using Markdown

import Base:
    +, -, *, /, //, fma,
    <, >, ==, !=, ⊆, ^, <=,
    zero, one, eps, abs, abs2, real, min, max,
    sqrt, exp, log, sin, cos, tan, cot, inv, cbrt, csc, hypot, sec,
    exp2, exp10, expm1, log2, log10, log1p,
    asin, acos, atan, acot,
    sinh, cosh, tanh, coth, csch, sech, asinh, acosh, atanh, acoth,
    sinpi, cospi,
    in, ∈, intersect, ∩, isempty, ⊆, ⊇,
    convert, promote_rule, float, real, big,
    eps, precision, setrounding, setprecision,
    floor, ceil, trunc, sign, round, copysign,
    show, hash, valtype

import .Broadcast: broadcasted

export
    AbstractInterval, Interval, ValidInterval, EmptyInterval,
    interval, emptyinterval, ∅, nai, entireinterval, widened, hull,
    inf, sup, diam, radius, mid, mag, mig, midpoint_radius,
    isempty, isnai, isentire, isthin, isatomic, iscommon,
    isinterior, ⪽, precedes, strictprecedes, ≼, ≺, ⊂, ⊃, ⊇,
    cancelminus, cancelplus, nthroot, mul_rev_to_pair,
    RoundTiesToEven, RoundTiesToAway

# export forbidden for type piracy: ±

if VERSION >= v"1.5.0-DEV.124"
    import Base: isdisjoint
else
    export isdisjoint
end

export show

## Includes
include("intervals.jl")
include(joinpath("internals", "special.jl"))
include(joinpath("internals", "rounding_macros.jl"))
include(joinpath("internals", "rounding.jl"))
include(joinpath("internals", "conversion.jl"))
include(joinpath("operations", "boolean.jl"))
include(joinpath("operations", "numeric.jl"))
include(joinpath("operations", "set_operations.jl"))
include(joinpath("operations", "arithmetic.jl"))
include(joinpath("operations", "round.jl"))
include(joinpath("operations", "powers.jl"))
include(joinpath("operations", "exp_log.jl"))
include(joinpath("operations", "trigonometric.jl"))
include(joinpath("operations", "hyperbolic.jl"))
include(joinpath("complementary", "complementary.jl"))

# These definitions has been put there because generated functions must be
# defined after all methods they use.
@generated function interval(::Type{T}, x::I) where {T <: Real, I <: Irrational}
    res = unsafe_atomic(T, x())  # Precompute the interval
    return :(return $res)  # Set body of the function to return the precomputed result
end

end # module IntervalArithmetic

ext: jl


text:
################################################################################
#                              TEST: Quaternions
################################################################################

for i = 1:samples
    # Rotations
    # =========

    # X Axis
    # ------

    # Sample an angle.
    ang = -pi + 2pi*rand()

    # Create a quaternion that rotates about X axis.
    q = angle_to_quat(ang, 0.0, 0.0, :XYZ)

    # Create a vector that does not have X component.
    v = [0;randn();randn()]

    # Rotate the reference using the quaternion.
    v_r = vect(inv_rotation(q)*v*q)

    # Get the sine of the angle between the representations.
    sin_ang = (cross(v_r, v)/norm(v)^2)[1]

    # Test the angle between the two representations.
    @test abs(sin(ang)-sin_ang) < 1e-10

    # Test the same algorithm using different functions.
    v_r = imag((1/q)*(v*q))

    # Get the sine of the angle between the representations.
    sin_ang = (cross(v_r, v)/norm(v)^2)[1]

    # Test the angle between the two representations.
    @test abs(sin(ang)-sin_ang) < 1e-10

    # Y Axis
    # ------

    # Sample an angle.
    ang = -pi + 2pi*rand()

    # Create a quaternion that rotates about Y axis.
    q = angle_to_quat(ang, 0.0, 0.0, :YXZ)

    # Create a vector that does not have Y component.
    v = [randn();0;randn()]

    # Rotate the reference using the quaternion.
    v_r = vect(inv_rotation(q)*v*q)

    # Get the sine of the angle between the representations.
    sin_ang = (cross(v_r, v)/norm(v)^2)[2]

    # Test the angle between the two representations.
    @test abs(sin(ang)-sin_ang) < 1e-10

    # Test the same algorithm using different functions.
    v_r = imag(q\v*q)

    # Get the sine of the angle between the representations.
    sin_ang = (cross(v_r, v)/norm(v)^2)[2]

    # Test the angle between the two representations.
    @test abs(sin(ang)-sin_ang) < 1e-10

    # Z Axis
    # ------

    # Sample an angle.
    ang = -pi + 2pi*rand()

    # Create a quaternion that rotates about Z axis.
    q = angle_to_quat(ang, 0.0, 0.0, :ZXY)

    # Create a vector that does not have Z component.
    v = [randn();randn();0]

    # Rotate the reference using the quaternion.
    v_r = vect(inv_rotation(q)*v*q)

    # Get the sine of the angle between the representations.
    sin_ang = (cross(v_r, v)/norm(v)^2)[3]

    # Test the angle between the two representations.
    @test abs(sin(ang)-sin_ang) < 1e-10

    # Test the same algorithm using different functions.
    v_r = imag(I/q*v*q)

    # Get the sine of the angle between the representations.
    sin_ang = (cross(v_r, v)/norm(v)^2)[3]

    # Test the angle between the two representations.
    @test abs(sin(ang)-sin_ang) < 1e-10

    # Other functions
    # ===============

    # eye
    # ---

    # Create a random quaternion with Float32 type.
    q = Quaternion{Float32}(randn(), randn(), randn(), randn())
    q = q/norm(q)

    # Creata an identity quaternion using many methods.
    q_eye_1 = Quaternion{Float64}(I)
    q_eye_2 = Quaternion(1.f0I)
    q_eye_3 = Quaternion(I,q)

    # Check if the types are correct.
    @test eltype(q_eye_1.q0) != eltype(real(q))
    @test eltype(q_eye_2.q0) == eltype(real(q))
    @test eltype(q_eye_3.q0) == eltype(real(q))

    # Check if the value of the quaternion is correct.
    @test norm(vect(q_eye_3)) == 0.0
    @test real(q_eye_3) == 1.0

    # zeros
    # -----

    # Create a random quaternion with Float32 type.
    q = Quaternion{Float32}(randn(), randn(), randn(), randn())
    q = q/norm(q)

    # Creata a zero quaternion using many methods.
    q_zeros_1 = zeros(Quaternion)
    q_zeros_2 = zeros(Quaternion{Float32})
    q_zeros_3 = zeros(q)

    # Check if the types are correct.
    @test eltype(q_zeros_1.q0) != eltype(real(q))
    @test eltype(q_zeros_2.q0) == eltype(real(q))
    @test eltype(q_zeros_3.q0) == eltype(real(q))

    # Check if the value of the quaternion is correct.
    @test norm(q_zeros_3) == 0.0

    # Create quaternion from vectors
    # ------------------------------

    # Vector
    let
        v = rand(3)
        q = Quaternion(v)

        @test q.q0 == 0
        @test q.q1 == v[1]
        @test q.q2 == v[2]
        @test q.q3 == v[3]

        v = rand(4)
        q = Quaternion(v)

        @test q.q0 == v[1]
        @test q.q1 == v[2]
        @test q.q2 == v[3]
        @test q.q3 == v[4]
    end

    # SVector
    let
        v = @SVector rand(3)
        q = Quaternion(v)

        @test q.q0 == 0
        @test q.q1 == v[1]
        @test q.q2 == v[2]
        @test q.q3 == v[3]

        v = @SVector rand(4)
        q = Quaternion(v)

        @test q.q0 == v[1]
        @test q.q1 == v[2]
        @test q.q2 == v[3]
        @test q.q3 == v[4]
    end

    # MVector
    let
        v = @MVector rand(3)
        q = Quaternion(v)

        @test q.q0 == 0
        @test q.q1 == v[1]
        @test q.q2 == v[2]
        @test q.q3 == v[3]

        v = @MVector rand(4)
        q = Quaternion(v)

        @test q.q0 == v[1]
        @test q.q1 == v[2]
        @test q.q2 == v[3]
        @test q.q3 == v[4]
    end
end

# Operations that were not tested yet.
# ====================================

let
    q = Quaternion(1,1,1,1)
    @test typeof(q) == Quaternion{Int}

    q = Quaternion(1,1,1,1.f0)
    @test typeof(q) == Quaternion{Float32}

    q = Quaternion(1.,0,0.f0,0)
    @test typeof(q) == Quaternion{Float64}

    q = Quaternion(SVector{3}(1.,2.,3.))
    @test q.q0 == 0
    @test q.q1 == 1
    @test q.q2 == 2
    @test q.q3 == 3

    q = Quaternion(SVector{4}(1.,2.,3.,4.))
    @test q.q0 == 1
    @test q.q1 == 2
    @test q.q2 == 3
    @test q.q3 == 4

    q = 3*q
    @test q.q0 == 3
    @test q.q1 == 6
    @test q.q2 == 9
    @test q.q3 == 12

    qc = conj(q)
    @test qc.q0 == +3
    @test qc.q1 == -6
    @test qc.q2 == -9
    @test qc.q3 == -12

    q = Quaternion(1.0,1.0,1.0,1.0)
    @test inv_rotation(q) != inv(q)
    @test (q/q)[:] ≈ [1;0;0;0]
    @test (q\q)[:] ≈ [1;0;0;0]

    # Uniform scaling.
    @test Quaternion(1,1,1,1) +    I === Quaternion{Int}(2,1,1,1)
    @test Quaternion(1,1,1,1) + 1.0I === Quaternion{Float64}(2,1,1,1)

    @test Quaternion(1,1,1,1) -    I === Quaternion{Int}(0,1,1,1)
    @test Quaternion(1,1,1,1) - 1.0I === Quaternion{Float64}(0,1,1,1)

    @test Quaternion(1,1,1,1)*I      === Quaternion{Int}(1,1,1,1)
    @test Quaternion(1,1,1,1)*1.0I   === Quaternion{Float64}(1,1,1,1)
    @test I*Quaternion(1,1,1,1)      === Quaternion{Int}(1,1,1,1)
    @test 1.0I*Quaternion(1,1,1,1)   === Quaternion{Float64}(1,1,1,1)

    @test Quaternion(1,1,1,1)/I      === Quaternion{Float64}(1,1,1,1)
    @test Quaternion(1,1,1,1)/1.0I   === Quaternion{Float64}(1,1,1,1)
    @test I/Quaternion(1,1,1,1)      === inv(Quaternion(1,1,1,1))
    @test 1.0I/Quaternion(1,1,1,1)   === inv(Quaternion(1,1,1,1))

    @test Quaternion(1,1,1,1)\I      === inv(Quaternion(1,1,1,1))
    @test Quaternion(1,1,1,1)\1.0I   === inv(Quaternion(1,1,1,1))
    @test I\Quaternion(1,1,1,1)      === Quaternion{Float64}(1,1,1,1)
    @test 1.0I\Quaternion(1,1,1,1)   === Quaternion{Float64}(1,1,1,1)

end

# Test printing
# =============

# Color test is not working on AppVeyor. Hence, if the OS is windows, then we
# just skip.
if !Sys.iswindows()
    # With colors.
    expected = """
Quaternion{Float64}:
  + 0.11111111 + 0.22222222.\e[1mi\e[0m + 0.33333333.\e[1mj\e[0m + 0.44444444.\e[1mk\e[0m"""

    result = sprint((io,quat)->show(io, MIME"text/plain"(), quat),
                    Quaternion(0.11111111,0.22222222,0.33333333,0.44444444);
                    context = :color => true)
    @test result == expected
end

# Without colors.
expected = """
Quaternion{Float64}:
  + 0.11111111 + 0.22222222.i + 0.33333333.j + 0.44444444.k"""

result = sprint((io,quat)->show(io, MIME"text/plain"(), quat),
                Quaternion(0.11111111,0.22222222,0.33333333,0.44444444);
                context = :color => false)
@test result == expected

# Inline.
expected = """
Quaternion{Float64}: + 0.11111111 + 0.22222222.i + 0.33333333.j + 0.44444444.k"""

result = sprint(print, Quaternion(0.11111111,0.22222222,0.33333333,0.44444444))
@test result == expected

# Test exceptions
# ===============

@test_throws ArgumentError Quaternion([1])
@test_throws ArgumentError Quaternion([1;2])
@test_throws ArgumentError Quaternion([1;2;3;4;5])
@test_throws ArgumentError Quaternion([1;2;3;4;5;6])
@test_throws ArgumentError dquat(Quaternion(I), [1])
@test_throws ArgumentError dquat(Quaternion(I), [1;2])
@test_throws ArgumentError dquat(Quaternion(I), [1;2;3;4])
@test_throws ArgumentError dquat(Quaternion(I), [1;2;3;4;5])


text:
import Base: *

"""
    CircularConvolution{M, N, T}

A preplanned, circular convolution operator on an M × N matrix of data of type T

# Fields
- `Ĝ`: DFT coefficients of the convolution kernel
- `F`: preplanned rFFT operator
- `F⁻¹`: preplanned irFFT operator
- `paddedSpace`: scratch space to zero-pad the input matrix
- `Â`: scratch space to store the DFT coefficients of the zero-padded input matrix

# Constructors:

- `CircularConvolution(G::Matrix{T})`

# Example:
```jldoctest
julia> G = repeat(1.0:3,1,4)
3×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0
 3.0  3.0  3.0  3.0

julia> C = CircularConvolution(G)
Circular convolution on a 3 × 4 matrix of data type Float64

julia> C*reshape(1:12, 3, 4)
3×4 Array{Int64,2}:
 164  164  164  164
 130  130  130  130
 148  148  148  148
```
"""
struct CircularConvolution{M, N, T, K, KI}
    Ĝ::Matrix{ComplexF64}
    F::K
    F⁻¹::KI

    paddedSpace::Matrix{T}
    Â::Matrix{ComplexF64}
end

function Base.show(io::IO, c::CircularConvolution{M, N, T}) where {M, N, T}
    print(io, "Circular convolution on a $M × $N matrix of data type $T")
end

function CircularConvolution(G::AbstractMatrix{T},fftw_flags = FFTW.ESTIMATE; dtype = Float64) where {T}
    FFTW.set_num_threads(2)

    M, N = size(G)
    #paddedSpace = Matrix{Float64}(undef, 2M-1, 2N-1)
    paddedSpace = Matrix{dtype}(undef, 2M, 2N)

    if dtype == ComplexF64
      F = FFTW.plan_fft(paddedSpace, flags = fftw_flags)
    else
      F = FFTW.plan_rfft(paddedSpace, flags = fftw_flags)
    end

    mirror!(paddedSpace, G)
    Ĝ = F * paddedSpace

    Â = similar(Ĝ)
    #F⁻¹ = FFTW.plan_irfft(Â, 2M - 1, flags = fftw_flags)
    if dtype == ComplexF64
      F⁻¹ = FFTW.plan_ifft(Â, flags = fftw_flags)
    else
      F⁻¹ = FFTW.plan_irfft(Â, 2M, flags = fftw_flags)
    end

    CircularConvolution{M, N, dtype, typeof(F), typeof(F⁻¹)}(Ĝ, F, F⁻¹, paddedSpace, Â)
end

function mul!(out, C::CircularConvolution{M, N, T}, B) where {M, N, T}
    FFTW.set_num_threads(2)
    @assert size(out) == size(B) == (M, N)

    inds = CartesianIndices((M,N))
    fill!(C.paddedSpace, 0)
    copyto!(C.paddedSpace, inds, B, inds)
    mul!(C.Â, C.F, C.paddedSpace)

    C.Â .*= C.Ĝ

    mul!(C.paddedSpace, C.F⁻¹, C.Â)

    #copyto!(out, inds, C.paddedSpace, CartesianIndices((M:2M-1,N:2N-1)))
    copyto!(out, inds, C.paddedSpace, CartesianIndices((M+1:2M,N+1:2N)))

end

C::CircularConvolution * B = mul!(similar(B), C, B)

function mirror!(A, a::AbstractArray{T,2}) where {T}
    Nr, Nc = size(a)
    #A[1:Nr-1, 1:Nc-1] .= a[Nr:-1:2, Nc:-1:2]
    #A[1:Nr-1, Nc:end] .= a[Nr:-1:2, 1:Nc]
    #A[Nr:end, 1:Nc-1] .= a[1:Nr, Nc:-1:2]
    #A[Nr:end, Nc:end] .= a
    A .= 0
    A[2:Nr, 2:Nc] .= a[Nr:-1:2, Nc:-1:2]
    A[2:Nr, Nc+1:end] .= a[Nr:-1:2, 1:Nc]
    A[Nr+1:end, 2:Nc] .= a[1:Nr, Nc:-1:2]
    A[Nr+1:end, Nc+1:end] .= a[1:Nr, 1:Nc]
    A
end

function mirror(a::AbstractArray{T,2}) where {T}
    Nr, Nc = size(a)
    mirror!(zeros(T, 2Nr-1, 2Nc-1), a)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## async event notifications

"""
    AsyncCondition()

Create a async condition that wakes up tasks waiting for it
(by calling [`wait`](@ref) on the object)
when notified from C by a call to `uv_async_send`.
Waiting tasks are woken with an error when the object is closed (by [`close`](@ref)).
Use [`isopen`](@ref) to check whether it is still active.
"""
mutable struct AsyncCondition
    handle::Ptr{Cvoid}
    cond::ThreadSynchronizer
    isopen::Bool
    set::Bool

    function AsyncCondition()
        this = new(Libc.malloc(_sizeof_uv_async), ThreadSynchronizer(), true, false)
        iolock_begin()
        associate_julia_struct(this.handle, this)
        err = ccall(:uv_async_init, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
            eventloop(), this, @cfunction(uv_asynccb, Cvoid, (Ptr{Cvoid},)))
        if err != 0
            #TODO: this codepath is currently not tested
            Libc.free(this.handle)
            this.handle = C_NULL
            throw(_UVError("uv_async_init", err))
        end
        finalizer(uvfinalize, this)
        iolock_end()
        return this
    end
end

"""
    AsyncCondition(callback::Function)

Create a async condition that calls the given `callback` function. The `callback` is passed one argument,
the async condition object itself.
"""
function AsyncCondition(cb::Function)
    async = AsyncCondition()
    @async while _trywait(async)
            cb(async)
            isopen(async) || return
        end
    return async
end

## timer-based notifications

"""
    Timer(delay; interval = 0)

Create a timer that wakes up tasks waiting for it (by calling [`wait`](@ref) on the timer object).

Waiting tasks are woken after an initial delay of `delay` seconds, and then repeating with the given
`interval` in seconds. If `interval` is equal to `0`, the timer is only triggered once. When
the timer is closed (by [`close`](@ref)) waiting tasks are woken with an error. Use [`isopen`](@ref)
to check whether a timer is still active.
"""
mutable struct Timer
    handle::Ptr{Cvoid}
    cond::ThreadSynchronizer
    isopen::Bool
    set::Bool

    function Timer(timeout::Real; interval::Real = 0.0)
        timeout ≥ 0 || throw(ArgumentError("timer cannot have negative timeout of $timeout seconds"))
        interval ≥ 0 || throw(ArgumentError("timer cannot have negative repeat interval of $interval seconds"))
        timeout = UInt64(round(timeout * 1000)) + 1
        interval = UInt64(round(interval * 1000))
        loop = eventloop()

        this = new(Libc.malloc(_sizeof_uv_timer), ThreadSynchronizer(), true, false)
        associate_julia_struct(this.handle, this)
        iolock_begin()
        err = ccall(:uv_timer_init, Cint, (Ptr{Cvoid}, Ptr{Cvoid}), loop, this)
        @assert err == 0
        finalizer(uvfinalize, this)
        ccall(:uv_update_time, Cvoid, (Ptr{Cvoid},), loop)
        err = ccall(:uv_timer_start, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, UInt64, UInt64),
            this, @cfunction(uv_timercb, Cvoid, (Ptr{Cvoid},)),
            timeout, interval)
        @assert err == 0
        iolock_end()
        return this
    end
end

unsafe_convert(::Type{Ptr{Cvoid}}, t::Timer) = t.handle
unsafe_convert(::Type{Ptr{Cvoid}}, async::AsyncCondition) = async.handle

function _trywait(t::Union{Timer, AsyncCondition})
    set = t.set
    if !set
        t.handle == C_NULL && return false
        iolock_begin()
        set = t.set
        if !set
            preserve_handle(t)
            lock(t.cond)
            try
                set = t.set
                if !set
                    if t.handle != C_NULL
                        iolock_end()
                        set = wait(t.cond)
                        unlock(t.cond)
                        iolock_begin()
                        lock(t.cond)
                    end
                end
            finally
                unlock(t.cond)
                unpreserve_handle(t)
            end
        end
        iolock_end()
    end
    t.set = false
    return set
end

function wait(t::Union{Timer, AsyncCondition})
    _trywait(t) || throw(EOFError())
    nothing
end

isopen(t::Union{Timer, AsyncCondition}) = t.isopen

function close(t::Union{Timer, AsyncCondition})
    iolock_begin()
    if t.handle != C_NULL && isopen(t)
        t.isopen = false
        ccall(:jl_close_uv, Cvoid, (Ptr{Cvoid},), t)
    end
    iolock_end()
    nothing
end

function uvfinalize(t::Union{Timer, AsyncCondition})
    iolock_begin()
    lock(t.cond)
    try
        if t.handle != C_NULL
            disassociate_julia_struct(t.handle) # not going to call the usual close hooks
            if t.isopen
                t.isopen = false
                ccall(:jl_close_uv, Cvoid, (Ptr{Cvoid},), t)
            end
            t.handle = C_NULL
            notify(t.cond, false)
        end
    finally
        unlock(t.cond)
    end
    iolock_end()
    nothing
end

function _uv_hook_close(t::Union{Timer, AsyncCondition})
    lock(t.cond)
    try
        t.isopen = false
        t.handle = C_NULL
        notify(t.cond, t.set)
    finally
        unlock(t.cond)
    end
    nothing
end

function uv_asynccb(handle::Ptr{Cvoid})
    async = @handle_as handle AsyncCondition
    lock(async.cond)
    try
        async.set = true
        notify(async.cond, true)
    finally
        unlock(async.cond)
    end
    nothing
end

function uv_timercb(handle::Ptr{Cvoid})
    t = @handle_as handle Timer
    lock(t.cond)
    try
        t.set = true
        if ccall(:uv_timer_get_repeat, UInt64, (Ptr{Cvoid},), t) == 0
            # timer is stopped now
            close(t)
        end
        notify(t.cond, true)
    finally
        unlock(t.cond)
    end
    nothing
end

"""
    sleep(seconds)

Block the current task for a specified number of seconds. The minimum sleep time is 1
millisecond or input of `0.001`.
"""
function sleep(sec::Real)
    sec ≥ 0 || throw(ArgumentError("cannot sleep for $sec seconds"))
    wait(Timer(sec))
    nothing
end

# timer with repeated callback
"""
    Timer(callback::Function, delay; interval = 0)

Create a timer that wakes up tasks waiting for it (by calling [`wait`](@ref) on the timer object) and
calls the function `callback`.

Waiting tasks are woken and the function `callback` is called after an initial delay of `delay` seconds,
and then repeating with the given `interval` in seconds. If `interval` is equal to `0`, the timer
is only triggered once. The function `callback` is called with a single argument, the timer itself.
When the timer is closed (by [`close`](@ref)) waiting tasks are woken with an error. Use [`isopen`](@ref)
to check whether a timer is still active.

# Examples

Here the first number is printed after a delay of two seconds, then the following numbers are printed quickly.

```julia-repl
julia> begin
           i = 0
           cb(timer) = (global i += 1; println(i))
           t = Timer(cb, 2, interval=0.2)
           wait(t)
           sleep(0.5)
           close(t)
       end
1
2
3
```
"""
function Timer(cb::Function, timeout::Real; interval::Real=0.0)
    timer = Timer(timeout, interval=interval)
    @async while _trywait(timer)
            cb(timer)
            isopen(timer) || return
        end
    return timer
end

"""
    timedwait(callback::Function, timeout::Real; pollint::Real=0.1)

Waits until `callback` returns `true` or `timeout` seconds have passed, whichever is earlier.
`callback` is polled every `pollint` seconds. The minimum value for `timeout` and `pollint`
is `0.001`, that is, 1 millisecond.

Returns :ok or :timed_out
"""
function timedwait(testcb::Function, timeout::Real; pollint::Real=0.1)
    pollint >= 1e-3 || throw(ArgumentError("pollint must be ≥ 1 millisecond"))
    start = time_ns()
    ns_timeout = 1e9 * timeout
    done = Channel(1)
    function timercb(aw)
        try
            if testcb()
                put!(done, (:ok, nothing))
            elseif (time_ns() - start) > ns_timeout
                put!(done, (:timed_out, nothing))
            end
        catch e
            put!(done, (:error, CapturedException(e, catch_backtrace())))
        finally
            isready(done) && close(aw)
        end
        nothing
    end

    try
        testcb() && return :ok
    catch e
        throw(CapturedException(e, catch_backtrace()))
    end

    t = Timer(timercb, pollint, interval = pollint)
    ret, e = fetch(done)
    close(t)

    ret === :error && throw(e)

    return ret
end

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
if VERSION >= v"1.1.0-DEV.215"
function dce!(ci, mi)
    argtypes = Core.Compiler.matching_cache_argtypes(mi, nothing)[1]
    ir = Compiler.inflate_ir(ci, sptypes_from_meth_instance(mi),
                             argtypes)
    compact = Core.Compiler.IncrementalCompact(ir, true)
    # Just run through the iterator without any processing
    Core.Compiler.foreach(x -> nothing, compact)
    ir = Core.Compiler.finish(compact)

    Core.Compiler.replace_code_newstyle!(ci, ir, length(argtypes)-1)
end
else
function dce!(ci, mi)
end
end

function preprocess!(::Consumer, ref::Reflection, optimize)
    if optimize
        # if the optimizer hasn't run, the IR hasn't been converted
        # to SSA form yet and dce is not legal
        dce!(ref.CI, ref.mi)
    end
end


text:
using Flux, Test
using DiffEqFlux: destructure, restructure

model = Chain(Dense(10, 5, relu), Dense(5, 2))

p = destructure(model)

m2 = restructure(model, p)

@test Flux.isleaf(m2[1].W)

x = rand(10)

@test model(x) == m2(x)


text:
using LearnBase
using DataFrames
import Base: cat, vcat, hcat, _cat, lastindex, getindex

abstract type AbstractNode end
abstract type AbstractProductNode <: AbstractNode end
abstract type AbstractBagNode <: AbstractNode end

# FIXME: this alias would better be Union{AbstractVector{T}, Tuple{Vararg{T}}}
# and method signatures should do AbstractVecOrTuple{<:T} when they want covariance,
# but that solution currently fails (see #27188 and #27224)
AbstractVecOrTuple{T} = Union{AbstractVector{<:T}, Tuple{Vararg{T}}}

"""
    data(x::AbstractNode)

    return data hold by the datanode
"""
data(x::AbstractNode) = x.data
data(x) = x

"""
    catobs(as...)

    concatenates `as...` into a single datanode while preserving their structure
"""
catobs(as...) = reduce(catobs, collect(as))

# reduction of common datatypes the way we like it
reduce(::typeof(catobs), as::Vector{<: AbstractMatrix}) = reduce(hcat, as)
@adjoint function reduce(::typeof(catobs), as::Vector{<: AbstractMatrix})
  sz = cumsum(size.(as, 2))
  return reduce(hcat, as), Δ -> (nothing, map(n -> Zygote.pull_block_horz(sz[n], Δ, as[n]), eachindex(as)))
end
reduce(::typeof(catobs), as::Vector{<: AbstractVector}) = reduce(vcat, as)
reduce(::typeof(catobs), as::Vector{<: DataFrame}) = reduce(vcat, as)
reduce(::typeof(catobs), as::Vector{<: Missing}) = missing
reduce(::typeof(catobs), as::Vector{<: Nothing}) = nothing
reduce(::typeof(catobs), as::Vector{<: Union{Missing, Nothing}}) = nothing
function reduce(::typeof(catobs), as::Vector{T}) where {T <: Union{Missing, AbstractNode}}
    reduce(catobs, [a for a in as if !ismissing(a)])
end

function reduce(::typeof(catobs), as::Vector{<: Any})
    isempty(as) && return(as)
    T = mapreduce(typeof, typejoin, as)
    T == Any && @error "cannot reduce Any"
    reduce(catobs, Vector{T}(as))
end

Base.cat(as::AbstractNode...; dims = :) = reduce(catobs, collect(as))

_cattrees(as::Vector{T}) where T <: Union{Tuple, Vector}  = tuple([reduce(catobs, [a[i] for a in as]) for i in 1:length(as[1])]...)
function _cattrees(as::Vector{T}) where T <: NamedTuple
    ks = keys(as[1])
    vs = [k => reduce(catobs, [a[k] for a in as]) for k in ks]
    (;vs...)
end

mapdata(f, x) = f(x)

# functions to make datanodes compatible with getindex and with MLDataPattern
Base.getindex(x::T, i::BitArray{1}) where T <: AbstractNode = x[findall(i)]
Base.getindex(x::T, i::Vector{Bool}) where T <: AbstractNode = x[findall(i)]
Base.getindex(x::AbstractNode, i::Int) = x[i:i]
Base.lastindex(ds::AbstractNode) = nobs(ds)
MLDataPattern.getobs(x::AbstractNode, i) = x[i]
MLDataPattern.getobs(x::AbstractNode, i, ::LearnBase.ObsDim.Undefined) = x[i]
MLDataPattern.getobs(x::AbstractNode, i, ::LearnBase.ObsDim.Last) = x[i]

#subset of common datatypes the way we like them
subset(x::AbstractMatrix, i) = x[:, i]
subset(x::AbstractVector, i) = x[i]
subset(x::AbstractNode, i) = x[i]
subset(x::DataFrame, i) = x[i, :]
subset(::Missing, i) = missing
subset(::Nothing, i) = nothing
subset(xs::Tuple, i) = tuple(map(x -> x[i], xs)...)
subset(xs::NamedTuple, i) = (; [k => xs[k][i] for k in keys(xs)]...)

include("arraynode.jl")

# definitions needed for all types of bag nodes
_len(a::UnitRange) = max(a.stop - a.start + 1, 0)
_len(a::Vector) = length(a)
LearnBase.nobs(a::AbstractBagNode) = length(a.bags)
LearnBase.nobs(a::AbstractBagNode, ::Type{ObsDim.Last}) = nobs(a)
Base.ndims(x::AbstractBagNode) = Colon()

include("bagnode.jl")
include("weighted_bagnode.jl")
include("ngrams.jl")
include("productnode.jl")
include("missingnode.jl")
include("lazynode.jl")


text:
module BenchPartitionBy

using Statistics
using BenchmarkTools
using Transducers

function manual_partition_by(xs)
    acc = zero(eltype(xs))
    buffer = similar(xs, 0)
    pval = xs[1] > 0
    push!(buffer, xs[1])
    for x in @view xs[2:end]
        val = x > 0
        if val == pval
            push!(buffer, x)
        else
            if mean(abs, buffer) < 1.0
                acc += prod(buffer)
            end
            empty!(buffer)
            push!(buffer, x)
            pval = val
        end
    end
    if mean(abs, buffer) < 1.0
        acc += prod(buffer)
    end
    return acc
end

suite = BenchmarkGroup()

xf = PartitionBy(x -> x > 0) |>
    Filter(xs -> mean(abs, xs) < 1.0) |>
    Map(prod)

n = 10^5

let xs = randn(100)
    @assert manual_partition_by(xs) == mapfoldl(xf, +, xs, init=0.0)
end

suite["xf"] = @benchmarkable(
    mapfoldl($xf, +,  xs, init=0.0),
    setup=(xs = randn($n)))
suite["man"] = @benchmarkable(
    manual_partition_by(xs),
    setup=(xs = randn($n)))

end  # module
BenchPartitionBy.suite


text:
struct PrepDataFrame
    taxa::Vector{<:AbstractString}
    df::DataFrame

    function PrepDataFrame(df::DataFrame)
        if names(df)[1] != "taxon"
            error("first column of dataframe must be \"$(TAXON_NAME)\"")
        end

        for i = 2:size(df, 2)
            col = df[!, i]
            if !(eltype(col) <: Float64)
                error("all columns (except the first) must have elements of type Float64")
            end
        end

        return new(df[!, 1], df[!, 2:end])
    end
end

import Base.getindex

function getindex(pdf::PrepDataFrame, args...)
    return getindex(pdf.df, args...)
end

import Base.setindex!

function setindex!(pdf::PrepDataFrame, args...)
    return setindex!(pdf.df, args...)
end


text:
import LinearAlgebra
import FileIO
import GeometryBasics
import Observables
import Makie
import GLMakie

using Porta

startframe = 1
FPS = 60
resolution = (3840, 2160)
segments = 48
speed = 1
factor = 0.0025
scale = 1.0
radius = 0.005
number = 720
frames = 3600
basemapradius = 0.325
basemapsegments = 2segments
basepointradius = 0.01
basepointsegments = 10
basepointtransparency = true
exportmode = ["gif", "frames", "video"]
exportmode = exportmode[2]
modelname = "spinor2"

# The scene object that contains other visual objects
scene = Makie.Scene(backgroundcolor = :black,
                               show_axis = false,
                               resolution = resolution)

fmap(b::S²) = b

"""
    getpointonpath(t)

Get a point on a path on S² with the given parameter `t`. t ∈ [0, 1].
"""
getpointonpath(t::Real) = begin
    t = t * 12π
    d = ℯ^cos(t) - 2cos(4t) - sin(t / 12)^5
    x, y = sin(t) * d, cos(t) * d
    ComplexLine(x + im * y)
end

"""
    getbutterflycurve(points)

Get butterfly curve by Temple H. Fay (1989) with the given number of `points`.
"""
function getbutterflycurve(points::Int)
    array = Array{ComplexLine,1}(undef, points)
    # 0 ≤ t ≤ 12π
    for i in 1:points
        t = (i - 1) / points * 12π
        d = ℯ^cos(t) - 2cos(4t) - sin(t / 12)^5
        x, y = sin(t) * d, cos(t) * d
        array[i] = ComplexLine(x + im * y)
    end
    array
end

quaternionicrgb(p::Quaternion, q::Quaternion, s::Real, v::Real) = begin
    hue = acos(dot(p, q)) / π
    hsvtorgb([hue * 360; s; v])
end

"""
    getpoints(z [, segments])

Get a circle of point in the base space with the given center `z` in the Riemann
sphere. The optional argument `segments` determines how many points should the
circle have.
"""
function getpoints(z::S²; segments::Int = 30)
    g = Geographic(z)
    r, ϕ, θ = vec(g)
    lspace = collect(range(float(-pi), stop = float(pi), length = segments))
    array = Array{Geographic,1}(undef, segments)
    for i in 1:segments
        ϕ′, θ′ = g.ϕ + factor * cos(lspace[i]), g.θ + factor * sin(lspace[i])
        array[i] = Geographic(r, ϕ′, θ′)
    end
    array
end

A(s, p::S²) = begin
    z = s(p)
    z₀, z₁ = vec(ComplexPlane(z))
    r = vec(z) # radial
    f = [-r[2]; r[1]; -r[4]; r[3]] # fiber tangent
    X₀, X₁ = vec(ComplexPlane(Quaternion(f)))
    1/2 * (conj(z₀) * X₀ - z₀ * conj(X₀) + conj(z₁) * X₁ - z₁ * conj(X₁))
end

"""
    getflower(;N=4, A=.5, B=-pi/7, P=pi/2, Q=0, number=300)

Calculates the x, y and z points of a flower in the base space.
with the given number of petals N, the fattness of the petals A,
the height of the petals B, the latitude of the flower P,
the rotation of the flower Q, and the total number of points in the grid.
"""
function getflower(;N=4, A=.5, B=-pi/7, P=pi/2, Q=0, number=300)
    N = 6
    A = .5
    B = -pi/7
    P = π / 2 - π / 4
    Q = 0
    t = range(0, stop = 2pi, length = number)
    az = 2pi .* t + A .* cos.(N .* 2pi .* t) .+ Q
    po = B .* sin.(N .* 2pi .* t) .+ P
    x = cos.(az).*sin.(po)
    y = sin.(az).*sin.(po)
    z = cos.(po)
    points = Array{Geographic}(undef, number)
    for i in 1:number
        points[i] = Geographic(Cartesian(x[i], y[i], z[i]))
    end
    points
end

#basemapconfig = Biquaternion(ℝ³(0, -1.25, 0))
basemapconfig = Biquaternion(Quaternion(1, 0, 0, 0), ℝ³(-0.9 + basemapradius / 2, -√2 - basemapradius / 2, 0))
basemapcolor = Makie.RGBAf(0.8, 0.8, 0.8, 0.2)
#basemapcolor = load("gallery/mars_1k_color.jpg")
basemaptransparency = true
basemap = Sphere(basemapconfig,
                 scene,
                 radius = basemapradius,
                 segments = basemapsegments,
                 color = basemapcolor,
                 transparency = basemaptransparency)

curvepoints = getflower(number = number)
solidcolors = Array{Makie.RGBAf,1}(undef, number)
ghostcolors = similar(solidcolors)
basepoints = Array{Sphere,1}(undef, number)
for i in 1:number
    config = Biquaternion(ℝ³(Cartesian(curvepoints[i])) * basemapradius)
    q₁ = λ⁻¹map(curvepoints[i])
    q₂ = Quaternion(1, 0, 0, 0)
    rgb = quaternionicrgb(q₁, q₂, 1, 1)
    # hue = (i - 1) / number * 360
    # rgb = hsvtorgb([hue, 1, 1])
    solidcolors[i] = Makie.RGBAf(rgb..., 0.7)
    ghostcolors[i] = Makie.RGBAf(rgb..., 0.1)
    sphere = Sphere(basemapconfig * config,
                    scene,
                    radius = basepointradius,
                    segments = basepointsegments,
                    color = solidcolors[i],
                    transparency = basepointtransparency)
    basepoints[i] = sphere
end

#bundleconfig = Biquaternion(ℝ³(0, 1.25, 0))
bundleconfig = Biquaternion(ℝ³(0, 0, 0))
s3rotation = Quaternion(1, 0, 0, 0)
solidfibers = Array{Fiber,1}(undef, number)
ghostfibers = Array{Fiber,1}(undef, number)
solidtop = U1(0)
solidbottom = U1(2π)
ghosttop = U1(0)
ghostbottom = U1(2π)
for i in 1:number
    point = curvepoints[i]
    transparency = false
    solidfiber = Fiber(scene,
                       point,
                       λ⁻¹map,
                       fmap,
                       radius = radius,
                       top = solidtop,
                       bottom = solidbottom,
                       s3rotation = s3rotation,
                       config = bundleconfig,
                       segments = segments,
                       color = solidcolors[i],
                       transparency = transparency,
                       scale = scale)
    solidfibers[i] = solidfiber
    transparency = true
    ghostfiber = Fiber(scene,
                       point,
                       λ⁻¹map,
                       fmap,
                       radius = radius,
                       top = ghosttop,
                       bottom = ghostbottom,
                       s3rotation = s3rotation,
                       config = bundleconfig,
                       segments = segments,
                       color = ghostcolors[i],
                       transparency = transparency,
                       scale = scale)
    ghostfibers[i] = ghostfiber
end

basespacepoints = Array{ComplexLine,1}(undef, number)
phases = Array{U1,2}(undef, number, 2)
s3rotations = Array{Quaternion,1}(undef, number)

len = 0.1
width = 0.01
transparency = false
tail, head = ℝ³(rand(3)), ℝ³(rand(3))
arrows = []
for i in 1:number
    push!(arrows, Arrow(tail, head, scene, width = width, color = solidcolors[i]))
end

group2algebra(A::SU2, r::Float64) = begin
    τ₃ = -0.5 .* [im 0; 0 -im]
    (-2 * r) .* (A.a * τ₃ * LinearAlgebra.inv(A.a))
end

exponentiate(M::Array{<:Complex,2}, t::Float64) = begin
    eigenvalues = LinearAlgebra.eigvals(M)
    eigenvectors = LinearAlgebra.eigvecs(M)
    S = eigenvectors
    S⁻¹ = LinearAlgebra.inv(S)
    e = LinearAlgebra.Diagonal(eigenvalues) .* t
    S = convert.(Complex{BigFloat}, S)
    S⁻¹ = convert.(Complex{BigFloat}, S⁻¹)
    e = convert.(Complex{BigFloat}, e)
    G = S * ℯ.^e * S⁻¹
    z₁, z₂ = G[1,1], G[1,2]
    v = LinearAlgebra.normalize([real(z₁); imag(z₁); real(z₂); imag(z₂)])
    z₁, z₂ = Complex{Float64}(v[1] + im * v[2]), Complex{Float64}(v[3] + im * v[4])
    G = [z₁ -conj(z₂); z₂ conj(z₁)]
    SU2(G)
end

B = SU2(Quaternion(π / 8, ℝ³(Cartesian(Geographic(1, π / 3, π / 4)))))
# ComplexPlane(0.9238795325112867 + 0.13529902503654928im, 0.2343447855778369 + 0.27059805007309845im)

L₁ = [Complex(0) Complex(1); Complex(1) Complex(0)]
L₂ = [Complex(0) Complex(-im); Complex(im) Complex(0)]
L₃ = [Complex(1) Complex(0); Complex(0) Complex(-1)]
rotation = 2π
frames = 360
ϵ = rotation / frames
g = L₃ # group2algebra(B, ϵ)
q₂ = g

for i in 1:frames
    h = group2algebra(exponentiate(g, i * ϵ - π), 1.0)
    q = rotate(h, ComplexPlane(B))
    q₂ = q
    println(q)
    if i % (frames / 4) == 0.0
        println("")
    end
end

# ComplexPlane(-0.010686830141843868 + 0.7131691405919512im, -0.612279177876824 + 0.34115945964803246im)
# ComplexPlane(-0.6123681646060238 + 0.356191878306672im, -0.0022871329322613106 - 0.7057813725427761im)
# ComplexPlane(-0.35355339059327384 - 0.6123724356957946im, -0.7071067811865475 + 4.336808689942018e-19im)
# ComplexPlane(0.6123680128847738 - 0.3562382883042756im, 0.002327401029800577 + 0.7057579485532132im)
# ComplexPlane(0.35355339059327384 + 0.6123724356957946im, 0.7071067811865475 - 8.673617379884035e-19im)
# ComplexPlane(-0.6123680128847738 + 0.3562382883042756im, -0.0023274010298005776 - 0.7057579485532132im)

println(q₀)
println(Quaternion(q))

"""
    animate(i)

Update the state of observables with the given frame number `i`.
"""
function animate(i::Int)
    step = (i - 1) / frames
    τ = step * rotation
    global q = rotate(q, e)

    for (index, item) in enumerate(curvepoints)
        z = item
        h = rotate(λ⁻¹map(z), q)
        z = πmap(h)
        q₂ = Quaternion(1, 0, 0, 0)
        rgb = quaternionicrgb(h, q₂, 0.8 + 0.2 * rand(), 0.8 + 0.2 * rand())
        solidcolors[index] = Makie.RGBAf(rgb..., 0.5 + 0.3 * rand())
        ghostcolors[index] = Makie.RGBAf(rgb..., 0.1 + 0.1 * rand())

        basespacepoints[index] = z
        phases[index, 1] = U1(0)
        phases[index, 2] = U1(τ - 2π)
    end
    for (index, item) in enumerate(basepoints)
        z = basespacepoints[index]
        p = ℝ³(Cartesian(z)) * basemapradius
        config = Biquaternion(p)
        update(item, basemapconfig * config)
        update(item, solidcolors[index])
    end
    for (index, item) in enumerate(solidfibers)
        top = phases[index, 1]
        bottom = phases[index, 2]
        update(item, curvepoints[index], q, radius, top, bottom)
        update(item, solidcolors[index])
    end
    for (index, item) in enumerate(ghostfibers)
        top = U1(0)
        bottom = U1(2π)
        update(item, curvepoints[index], q, radius, top, bottom)
        update(item, ghostcolors[index])
    end

    # update tangent vectors
    for (index, item) in enumerate(arrows)
        w = basespacepoints[index]
        tail = ℝ³(Cartesian(w)) * basemapradius
        initial = ℝ³(0, 0, 1) * basemapradius
        q = getrotation(normalize(tail), normalize(initial))
        α = phases[index, 2]
        head = initial + ℝ³(1, 0, 0) * (1 / 10)
        head = rotate(head, q)
        q = Quaternion(angle(α), normalize(tail))
        head = rotate(head, q)
        tail, head = applyconfig([tail; head], basemapconfig)
        update(item, tail, head - tail)
        update(item, ghostcolors[index])
    end
end

n = ℝ³(1, 0, 0)
v = ℝ³(0, 0, 2) * √2
# update eye position
# scene.camera.eyeposition.val
upvector = GeometryBasics.Vec3f0(vec(n)...)
eyeposition = GeometryBasics.Vec3f0(vec(v)...)
lookat = GeometryBasics.Vec3f0(0, 0, 0)
Makie.update_cam!(scene, eyeposition, lookat, upvector)
scene.center = false # prevent scene from recentering on display
if exportmode ∈ ["gif", "video"]
    outputextension = exportmode == "gif" ? "gif" : "mkv"
    Makie.record(scene, "gallery/$modelname.$outputextension",
                            framerate = FPS) do io
        for i in startframe:frames
            animate(i) # animate the scene
            Makie.recordframe!(io) # record a new frame
            step = (i - 1) / frames
            println("Completed step $(100step).\n")
        end
    end
elseif exportmode == "frames"
    directory = joinpath("gallery", modelname)
    !isdir(directory) && mkdir(directory)
    directory = joinpath("gallery", modelname)
    !isdir(directory) && mkdir(directory)
    for i in startframe:frames
        start = time()
        animate(i) # animate the scene
        elapsed = time() - start
        sleep(0.1)
        println("Generating frame $i took $elapsed (s).")

        start = time()

        paddingnumber = length(digits(frames))
        imageid = lpad(i, paddingnumber, "0")
        imagename = "$(modelname)_$(resolution[2])p_$(FPS)fps_$imageid.jpeg"
        filepath = joinpath(directory, imagename)
        FileIO.save(filepath,
                    scene;
                    resolution = resolution,
                    pt_per_unit = 400.0,
                    px_per_unit = 400.0)
        elapsed = time() - start
        sleep(0.1)
        println("Saving file $filepath took $elapsed (s).")
        step = (i - 1) / frames
        println("Completed step $(100step).\n")

        stitch() = begin
            part = i ÷ (frames / 3)
            exportdir = joinpath(directory, "export")
            !isdir(exportdir) && mkdir(exportdir)
            WxH = "$(resolution[1])x$(resolution[2])"
            #WxH = "1920x1080"
            commonpart = "$(modelname)_$(resolution[2])p_$(FPS)fps"
            inputname = "$(commonpart)_%0$(paddingnumber)d.jpeg"
            inputpath = joinpath(directory, inputname)
            outputname = "$(commonpart)_$(part).mp4"
            outputpath = joinpath(exportdir, outputname)
            command1 = `ffmpeg -y -f image2 -framerate $FPS -i $inputpath -s $WxH -pix_fmt yuvj420p $outputpath`
            run(command1)
        end

        if isapprox(i % (frames / 3), 0)
            stitch()
            sleep(1)
        end
    end
end

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
"""
    struct CoreModelCoupled <: MetabolicModel

The linear model with additional coupling constraints in the form
```
    cₗ ≤ C x ≤ cᵤ
```
"""
mutable struct CoreModelCoupled <: MetabolicModel
    lm::CoreModel
    C::SparseMat
    cl::SparseVec
    cu::SparseVec

    function CoreModelCoupled(
        lm::MetabolicModel,
        C::M,
        cl::V,
        cu::V,
    ) where {V<:VecType,M<:MatType}

        length(cu) == length(cl) ||
            throw(DimensionMismatch("`cl` and `cu` need to have the same size"))
        size(C) == (length(cu), n_reactions(lm)) ||
            throw(DimensionMismatch("wrong dimensions of `C`"))

        new(convert(CoreModel, lm), sparse(C), sparse(cl), sparse(cu))
    end
end

"""
    reactions(a::CoreModelCoupled)

Extract reactions from [`CoreModelCoupled`](@ref) (uses the internal
[`CoreModel`](@ref)).
"""
reactions(a::CoreModelCoupled) = reactions(a.lm)

"""
    metabolites(a::CoreModelCoupled)

Extract metabolites from [`CoreModelCoupled`](@ref) (uses the internal
[`CoreModel`](@ref)).
"""
metabolites(a::CoreModelCoupled) = metabolites(a.lm)

"""
    stoichiometry(a::CoreModelCoupled)

Extract stoichiometry from [`CoreModelCoupled`](@ref) (uses the internal
[`CoreModel`](@ref)).
"""
stoichiometry(a::CoreModelCoupled) = stoichiometry(a.lm)

"""
    bounds(a::CoreModelCoupled)

Extract bounds from [`CoreModelCoupled`](@ref) (uses the internal
[`CoreModel`](@ref)).
"""
bounds(a::CoreModelCoupled) = bounds(a.lm)

"""
    balance(a::CoreModelCoupled)

Extract balance from [`CoreModelCoupled`](@ref) (uses the internal
[`CoreModel`](@ref)).
"""
balance(a::CoreModelCoupled) = balance(a.lm)

"""
    objective(a::CoreModelCoupled)

Extract objective from [`CoreModelCoupled`](@ref) (uses the internal
[`CoreModel`](@ref)).
"""
objective(a::CoreModelCoupled) = objective(a.lm)

"""
    coupling(a::CoreModelCoupled)::SparseMat

Coupling constraint matrix for a `CoreModelCoupled`.
"""
function coupling(a::CoreModelCoupled)::SparseMat
    a.C
end

"""
    n_coupling_constraints(a::CoreModelCoupled)::Int

The number of coupling constraints in a `CoreModelCoupled`.
"""
function n_coupling_constraints(a::CoreModelCoupled)::Int
    return size(a.C, 1)
end

"""
    coupling_bounds(a::CoreModelCoupled)::Tuple{SparseVec,SparseVec}

Coupling bounds for a `CoreModelCoupled`.
"""
function coupling_bounds(a::CoreModelCoupled)::Tuple{SparseVec,SparseVec}
    (a.cl, a.cu)
end

"""
    reaction_stoichiometry(model::CoreModelCoupled, rxn_id::String)::Dict{String, Float64}

Return the reaction equation of reaction with id `rxn_id` in model. The reaction
equation maps metabolite ids to their stoichiometric coefficients.
"""
function reaction_stoichiometry(m::CoreModelCoupled, rxn_id::String)::Dict{String,Float64}
    reaction_stoichiometry(m.lm, rxn_id)
end

"""
    reaction_stoichiometry(model::CoreModelCoupled, rxn_ind::Int)::Dict{String, Float64}

Return the reaction equation of reaction with id `rxn_ind` in model. The reaction
equation maps metabolite ids to their stoichiometric coefficients.
"""
function reaction_stoichiometry(m::CoreModelCoupled, rxn_ind::Int)::Dict{String,Float64}
    reaction_stoichiometry(m.lm, rxn_ind)
end

"""
    Base.convert(::Type{CoreModelCoupled}, mm::MetabolicModel)

Make a `CoreModelCoupled` out of any compatible model type.
"""
function Base.convert(::Type{CoreModelCoupled}, mm::MetabolicModel)
    if typeof(mm) == CoreModelCoupled
        return mm
    end

    (cl, cu) = coupling_bounds(mm)
    CoreModelCoupled(convert(CoreModel, mm), coupling(mm), cl, cu)
end


text:
module Precompile2

using Precompile1
test() = Precompile1.Foo((b=1, a=2)).nt

end # module


text:
export InitVaguePipelineStage

# TODO, wont' work for nodes with different variate_form on different edges

"""
    InitVaguePipelineStage{D} <: AbstractPipelineStage

Uses the `start_with()` operator from `Rocket.jl` library and `vague` function to initialise first updates for a `stream` in the form of vague messages.

# Arguments
- `dimensionality`: (optional) dimensionality of initial vague messages

See also: [`AbstractPipelineStage`](@ref), [`apply_pipeline_stage`](@ref), [`EmptyPipelineStage`](@ref), [`CompositePipelineStage`](@ref)
"""
struct InitVaguePipelineStage{D} <: AbstractPipelineStage
    dimensionality :: D
end

InitVaguePipelineStage() = InitVaguePipelineStage(nothing)

getdimensionality(stage::InitVaguePipelineStage) = stage.dimensionality

variate_form(stage::InitVaguePipelineStage) = variate_form(stage, getdimensionality(stage))

variate_form(::InitVaguePipelineStage, ::Nothing)         = Univariate
variate_form(::InitVaguePipelineStage, ::Int)             = Multivariate
variate_form(::InitVaguePipelineStage, ::Tuple{Int, Int}) = Matrixvariate

apply_pipeline_stage(stage::InitVaguePipelineStage, factornode, tag, stream) = __apply_vague_pipeline_stage(variate_form(stage), stage, factornode, tag, stream)

function __apply_vague_pipeline_stage(::Type{ Univariate }, stage::InitVaguePipelineStage, factornode, tag, stream)
    return stream |> start_with(Message(vague(conjugate_type(functionalform(factornode), tag)), false, true))
end

function __apply_vague_pipeline_stage(::Type{ Multivariate }, stage::InitVaguePipelineStage, factornode, tag, stream)
    return stream |> start_with(Message(vague(conjugate_type(functionalform(factornode), tag), getdimensionality(stage)), false, true))
end

function __apply_vague_pipeline_stage(::Type{ Matrixvariate }, stage::InitVaguePipelineStage, factornode, tag, stream)
    return stream |> start_with(Message(vague(conjugate_type(functionalform(factornode), tag), getdimensionality(stage)...), false, true))
end


text:
@testset "Copy sequence" begin
    # Test copy from sequence to sequence
    function test_copy(A, str)
        seq = LongSequence{A}(str)
        str2 = String(copy(seq))
        @test (str == str2 == String(seq))
    end

    for len in [0, 1, 15, 33]
        test_copy(DNAAlphabet{4}, random_dna(len))
        test_copy(RNAAlphabet{2}, random_rna(len, [0.25, 0.25, 0.25, 0.25]))
        test_copy(AminoAcidAlphabet, random_aa(len))
        test_copy(CharAlphabet, random_aa(len))
    end

    # Test copying shared seq
    seq = LongSequence{DNAAlphabet{2}}(random_dna(50, [0.25, 0.25, 0.25, 0.25]))
    seq2 = seq[5:44]
    @test String(seq)[5:44] == String(seq2)
end # testset

@testset "Copy! sequence" begin
    function test_copy!(A, srctxt)
        src = LongSequence{A}(srctxt)
        dst = LongSequence{A}(0)
        for len in [max(0, length(src) - 3), length(src), length(src) + 4]
            resize!(dst, len)
            copy!(dst, src)
            @test String(dst) == String(src)
        end
    end

    test_copy!(DNAAlphabet{4}, random_dna(14))
    test_copy!(RNAAlphabet{2}, random_rna(55, [0.25, 0.25, 0.25, 0.25]))
    test_copy!(AminoAcidAlphabet, random_aa(18))
    test_copy!(CharAlphabet, random_aa(9))

    #  Also works across nucleotide types!
    for N in (2,4)
        src = LongSequence{DNAAlphabet{N}}(random_dna(33, [0.25, 0.25, 0.25, 0.25]))
        dst = LongSequence{RNAAlphabet{N}}(random_rna(31, [0.25, 0.25, 0.25, 0.25]))
        copy!(dst, src)
        @test String(typeof(dst)(src)) == String(dst)
        resize!(dst, 16)
        copy!(src, dst)
        @test String(typeof(dst)(src)) == String(dst)
    end

    # Doesn't work for wrong types
    @test_throws MethodError copy!(LongDNASeq("TAG"), LongAminoAcidSeq("WGM"))
    @test_throws MethodError copy!(LongSequence{DNAAlphabet{2}}("TAG"), LongRNASeq("UGM"))
end

@testset "Copyto! sequence" begin
    function test_copyto!1(A1, dst, A2, src)
        dst_ = LongSequence{A1}(dst)
        src_ = LongSequence{A2}(src)
        copyto!(dst_, src_)
        @test String(src_) == String(dst_[1:length(src_)])
    end

    test_copyto!1(DNAAlphabet{4}, random_dna(19), DNAAlphabet{4}, random_dna(17))
    test_copyto!1(RNAAlphabet{2}, random_rna(31, [0.25, 0.25, 0.25, 0.25]),
                  RNAAlphabet{2}, random_rna(11, [0.25, 0.25, 0.25, 0.25]))
    test_copyto!1(AminoAcidAlphabet, random_aa(61), AminoAcidAlphabet, random_aa(61))

    function test_copyto2!(A, F)
        for len in [10, 17, 51]
            start = rand(1:3)
            N = len - rand(1:4) - start
            src = LongSequence{A}(F(len))
            dst = LongSequence{A}(F(len))
            copyto!(dst, start, src, start + 1, N)
            @test String(dst[start:start+N-1]) == String(src[start+1:start+N])
        end
    end

    test_copyto2!(DNAAlphabet{2}, len -> random_dna(len, [0.25, 0.25, 0.25, 0.25]))
    test_copyto2!(RNAAlphabet{4}, random_rna)
    test_copyto2!(AminoAcidAlphabet, random_aa)
    test_copyto2!(CharAlphabet, random_aa)

    # Test bug when copying to self
    src = LongDNASeq("A"^16 * "C"^16 * "A"^16)
    copyto!(src, 17, src, 1, 32)
    @test String(src) == "A"^32 * "C"^16
end

@testset "Copy! data" begin
    function test_copy!(seq, src)
        @test String(src) == String(copy!(seq, src))
    end
    # Needed because conversion to String truncates vector.
    function test_copy!(seq, src::Vector)
        @test String(copy(src)) == String(copy!(copy(seq), src))
    end

    probs = [0.25, 0.25, 0.25, 0.25, 0.00]
    dna2 = LongSequence{DNAAlphabet{2}}(6)
    dna4 = LongSequence{DNAAlphabet{4}}(6)
    rna2 = LongSequence{RNAAlphabet{2}}(6)
    rna4 = LongSequence{RNAAlphabet{4}}(6)
    aa = LongSequence{AminoAcidAlphabet}(6)
    charseq = LongSequence{CharAlphabet}(6)
    for dtype in [Vector{UInt8}, Vector{Char}, String, Test.GenericString]
        for len in [0, 1, 10, 16, 32, 100, 5]
            test_copy!(dna2, dtype(random_dna(len, probs)))
            test_copy!(dna4, dtype(random_dna(len)))
            test_copy!(rna2, dtype(random_rna(len, probs)))
            test_copy!(rna4, dtype(random_rna(len)))
            test_copy!(aa, dtype(random_aa(len)))
            test_copy!(charseq, dtype(random_aa(len)))
        end
    end
    test_copy!(charseq, "ϐʌ⨝W")
end

@testset "Copyto! data" begin
    function test_twoarg_copyto!(seq, src)
        copyto!(seq, src)
        @test String(src[1:length(src)]) == String(src)
    end
    # Needed because conversion to String truncates vector.
    function test_twoarg_copyto!(seq, src::Vector)
        copyto!(seq, src)
        @test String(src[1:length(src)]) == String(copy(src))
    end

    probs = [0.25, 0.25, 0.25, 0.25, 0.00]
    dna2 = LongSequence{DNAAlphabet{2}}(50)
    dna4 = LongSequence{DNAAlphabet{4}}(50)
    rna2 = LongSequence{RNAAlphabet{2}}(50)
    rna4 = LongSequence{RNAAlphabet{4}}(50)
    aa = LongSequence{AminoAcidAlphabet}(50)
    charseq = LongSequence{CharAlphabet}(50)
    for dtype in [Vector{UInt8}, Vector{Char}, String, Test.GenericString]
        for len in [0, 1, 10, 16, 32, 5]
            test_twoarg_copyto!(dna2, dtype(random_dna(len, probs)))
            test_twoarg_copyto!(dna4, dtype(random_dna(len)))
            test_twoarg_copyto!(rna2, dtype(random_rna(len, probs)))
            test_twoarg_copyto!(rna4, dtype(random_rna(len)))
            test_twoarg_copyto!(aa, dtype(random_aa(len)))
            test_twoarg_copyto!(charseq, dtype(random_aa(len)))
        end
    end
    copyto!(charseq, "ϐʌ⨝W")
    @test collect(charseq[1:4]) == collect("ϐʌ⨝W")

    # Five-arg copyto!
    function test_fivearg_copyto!(seq, src)
        for soff in (1, 3)
            for doff in (1, 5)
                for N in (0, 5, 18, 30)
                    copyto!(seq, doff, src, soff, N)
                    @test String(seq[doff:doff+N-1]) == String(src[soff:soff+N-1])
                end
            end
        end
    end

    for dtype in [Vector{UInt8}, Vector{Char}, String, Test.GenericString]
        test_fivearg_copyto!(dna2, dtype(random_dna(60, probs)))
        test_fivearg_copyto!(dna4, dtype(random_dna(60)))
        test_fivearg_copyto!(rna2, dtype(random_rna(60, probs)))
        test_fivearg_copyto!(rna4, dtype(random_rna(60)))
        test_fivearg_copyto!(aa, dtype(random_aa(60)))
        test_fivearg_copyto!(charseq, dtype(random_aa(60)))
    end

end

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

@testset "Concatenation" begin
    function test_concatenation(A, chunks)
        parts = UnitRange{Int}[]
        for i in 1:lastindex(chunks)
            start = rand(1:length(chunks[i]))
            stop = rand(start:length(chunks[i]))
            push!(parts, start:stop)
        end
        str = string([chunk[parts[i]] for (i, chunk) in enumerate(chunks)]...)
        seq = *([LongSequence{A}(chunk)[parts[i]] for (i, chunk) in enumerate(chunks)]...)
        @test convert(String, seq) == uppercase(str)
    end

    for _ in 1:100
        n = rand(1:10)
        chunks = [random_dna(rand(1:100)) for _ in 1:n]
        test_concatenation(DNAAlphabet{4}, chunks)

        chunks = [random_rna(rand(1:100)) for _ in 1:n]
        test_concatenation(RNAAlphabet{4}, chunks)

        chunks = [random_aa(rand(1:100)) for _ in 1:n]
        test_concatenation(AminoAcidAlphabet, chunks)

        probs = [0.25, 0.25, 0.25, 0.25, 0.00]
        chunks = [random_dna(rand(1:100), probs) for _ in 1:n]
        test_concatenation(DNAAlphabet{2}, chunks)

        chunks = [random_rna(rand(1:100), probs) for _ in 1:n]
        test_concatenation(RNAAlphabet{2}, chunks)
    end
end

@testset "Repetition" begin
    function test_repetition(A, chunk)
        start = rand(1:length(chunk))
        stop = rand(start:length(chunk))
        n = rand(1:10)
        str = chunk[start:stop] ^ n
        seq = LongSequence{A}(chunk)[start:stop] ^ n
        @test convert(String, seq) == uppercase(str)
    end

    for _ in 1:10
        chunk = random_dna(rand(1:100))
        test_repetition(DNAAlphabet{4}, chunk)

        chunk = random_rna(rand(1:100))
        test_repetition(RNAAlphabet{4}, chunk)

        chunk = random_aa(rand(1:100))
        test_repetition(AminoAcidAlphabet, chunk)

        probs = [0.25, 0.25, 0.25, 0.25, 0.00]
        chunk = random_dna(rand(1:100), probs)
        test_repetition(DNAAlphabet{2}, chunk)

        chunk = random_rna(rand(1:100), probs)
        test_repetition(RNAAlphabet{2}, chunk)
    end
end

@testset "Length" begin
    for len in [0, 1, 2, 3, 10, 16, 32, 1000, 10000]
        seq = LongDNASeq(random_dna(len))
        @test length(seq) === lastindex(seq) === len

        seq = LongRNASeq(random_rna(len))
        @test length(seq) === lastindex(seq) === len

        seq = LongAminoAcidSeq(random_aa(len))
        @test length(seq) === lastindex(seq) === len
    end

    @test length(char"いろはabc") === 6
end

@testset "Access" begin
    dna_seq = dna"ACTG"

    @test dna_seq[1] === DNA_A
    @test dna_seq[2] === DNA_C
    @test dna_seq[3] === DNA_T
    @test dna_seq[4] === DNA_G

    # Access indexes out of bounds
    @test_throws BoundsError dna_seq[-1]
    @test_throws BoundsError dna_seq[0]
    @test_throws BoundsError dna_seq[5]

    @test dna"ACTGNACTGN"[1:5] == dna"ACTGN"
    @test dna"ACTGNACTGN"[5:1] == dna""

    rna_seq = rna"ACUG"
    @test rna_seq[1] === RNA_A
    @test rna_seq[2] === RNA_C
    @test rna_seq[3] === RNA_U
    @test rna_seq[4] === RNA_G

    # Access indexes out of bounds
    @test_throws BoundsError rna_seq[-1]
    @test_throws BoundsError rna_seq[0]
    @test_throws BoundsError rna_seq[5]

    @test rna"ACUGNACUGN"[1:5] == rna"ACUGN"
    @test rna"ACUGNACUGN"[5:1] == rna""

    @test aa"KSAAV"[3] == AA_A
    @test char"いろはにほ"[3] == 'は'
end

@testset "Equality" begin
    @testset "DNA" begin
        a = dna"ACTGN"
        b = dna"ACTGN"
        @test a == b
        @test dna"ACTGN" == dna"ACTGN"
        @test dna"ACTGN" != dna"ACTGA"
        @test dna"ACTGN" != dna"ACTG"
        @test dna"ACTG"  != dna"ACTGN"

        a = dna"ACGTNACGTN"
        b = dna"""
        ACGTN
        ACGTN
        """
        @test a == b
    end

    @testset "RNA" begin
        a = rna"ACUGN"
        b = rna"ACUGN"
        @test a == b
        @test rna"ACUGN" == rna"ACUGN"
        @test rna"ACUGN" != rna"ACUGA"
        @test rna"ACUGN" != rna"ACUG"
        @test rna"ACUG"  != rna"ACUGN"

        a = rna"ACUGNACUGN"
        b = rna"""
        ACUGN
        ACUGN
        """
        @test a == b
    end

    @testset "AminoAcid" begin
        a = aa"ARNDCQEGHILKMFPSTWYVX"
        b = aa"ARNDCQEGHILKMFPSTWYVX"
        @test a == b
        @test a != aa"ARNDCQEGHILKMFPSTWYXV"
        @test a != aa"ARNDCQEGHLKMFPSTWYVX"

        b = aa"""
        ARNDCQEGHI
        LKMFPSTWYV
        X
        """
        @test a == b
    end
end


text:
"""
    Experimental

This module contains experimental stuff that does not respect SemVer
"""
module Experimental

using Graphs: AbstractGraph, weights, nv, vertices, outneighbors, DijkstraState, DefaultDistance
using DataStructures: PriorityQueue, dequeue!
using SimpleValueGraphs: ValGraph, edgevals_type, outedgevals
using SimpleWeightedGraphs: SimpleWeightedGraph, weighttype

struct InternalEdgeVals{key} end

internaledgevals_or_weights(g::AbstractGraph) = weights(g)

internaledgevals_or_weights(g::SimpleWeightedGraph) = InternalEdgeVals{nothing}()

internaledgevals_or_weights(g::ValGraph, key=1) = InternalEdgeVals{key}()

outneighbors_and_edgevals(g::AbstractGraph, weights, u) =
    ((v, @inbounds weights[u, v]) for v in outneighbors(g, u))

outneighbors_and_edgevals(g::ValGraph, ::InternalEdgeVals{key}, u) where key =
    zip(outneighbors(g, u), outedgevals(g, u, key))

function outneighbors_and_edgevals(g::SimpleWeightedGraph, ::InternalEdgeVals, u)

    w = g.weights
    zip(outneighbors(g, u), view(w.nzval, w.colptr[1]:w.colptr[u+1]-1))
end

function _edgevals_type(g::ValGraph, ::InternalEdgeVals{key}) where key

    return edgevals_type(g, key)
end

function _edgevals_type(g::SimpleWeightedGraph, ::InternalEdgeVals)

    return weighttype(g)
end

function dijkstra_shortest_paths(g::AbstractGraph, srcs::Vector{<:Integer}, distmx=internaledgevals_or_weights(g);
    allpaths=false,
    trackvertices=false)

    srcs ⊆ vertices(g) || error("srcs must be vertices of g")

    V = eltype(g)
    W = if distmx isa InternalEdgeVals
        _edgevals_type(g, distmx)
    else
        size(distmx, 1) >= nv(g) || error("distmx is too small")
        size(distmx, 2) >= nv(g) || error("distmx is too small")

        eltype(distmx)
    end

    nvg = Int(nv(g))
    dists = fill(typemax(W), nvg)
    parents = zeros(V, nvg)
    visited = falses(nvg)
    pathcounts = zeros(nvg)
    closest_vertices = sizehint!(Vector{V}(), nvg)
    preds = fill(Vector{V}(), nvg)

    q = PriorityQueue{V, W}()

    for src in srcs
        @inbounds dists[src] = zero(W)
        @inbounds visited[src] = true
        @inbounds pathcounts[src] = one(Float64)
        q[src] = zero(W)
    end

    while !isempty(q)
        u = dequeue!(q)

        if trackvertices
            push!(closest_vertices, u)
        end

        @inbounds d = dists[u]
        @inbounds for (v, w) in outneighbors_and_edgevals(g, distmx, u)
            alt = d + w

            if !visited[v]
                visited[v] = true
                dists[v] = alt
                parents[v] = u

                pathcounts[v] = pathcounts[u]
                if allpaths
                    preds[v] = u
                end
                q[v] = alt
            elseif alt < dists[v]
                dists[v] = alt
                parents[v] = u
                pathcounts[v] = pathcounts[u]
                if allpaths
                    resize!(preds[v], 1)
                    preds[v][1] = u
                end
                q[v] = alt
            elseif alt == dists[v]
                pathcounts[v] += pathcounts[u]
                if allpaths
                    push!(preds[v], u)
                end
            end
        end
    end

    for src in srcs
        @inbounds pathcounts[src] = one(Float64)
        @inbounds parents[src] = 0
        @inbounds empty!(preds[src])
    end

    return DijkstraState{W, V}(parents, dists, preds, pathcounts, closest_vertices)
end

dijkstra_shortest_paths(g::AbstractGraph, src::Integer, distmx=internaledgevals_or_weights(g); kwargs...) =
    dijkstra_shortest_paths(g, [src], distmx; kwargs...)

end # end module


text:
module Benchmark
# ----------------------------------------------------
# --- Benchmark module
# ----------------------------------------------------
# This module uses UHDBindings and can be used to benchmark output rate of the SDR
# Benchmark.main(5e6,args) returns a tuple (ideally (5e6,5e6)) where
# Input
# - First parmaeter is the desired radio rate
# - Second parameter is the radio parameter (USRP device)
# Output
# First parameter is the actual radio rate
# Second parameter is the calculated rate

# To benchmark a range one can do Benchmark.main.(1e6:1e6:16e6)

# --- External modules
using UHDBindings
# --- Functions
"""
Calculate rate based on Julia timing
"""
function getRate(tInit,tFinal,nbSamples)
	return nbSamples / (tFinal-tInit);
end

"""
Main call to monitor Rx rate
"""
function main(samplingRate,args,duration=2)
	# ----------------------------------------------------
	# --- Physical layer and RF parameters
	# ----------------------------------------------------
	# --- Create the radio object in function
	carrierFreq		= 770e6;
	gain			= 50.0;
    radio			= openUHD(carrierFreq,samplingRate,gain,args=args)
	# --- Print the configuration
	print(radio);
	# --- Init parameters
	# Get the radio size for buffer pre-allocation
	nbSamples 		= radio.rx.packetSize;
    @show nbSamples
	# We will get complex samples from recv! method
    sig		  = zeros(Complex{Cfloat},nbSamples);
	# --- Targeting 2 seconds acquisition
	# Init counter increment
	nS		  = 0;
	# Max counter definition
	nbBuffer  = samplingRate*duration;
	# --- Timestamp init
	p 			= recv!(sig,radio);
	nS			+= p;
	timeInit  	= time();
	while true
		# --- Direct call to avoid allocation
		p = recv!(sig,radio);
		# # --- Ensure packet is OK
		# err 	= getError(radioRx);
		# --- Update counter
		nS		+= p;
		# --- Interruption
		if nS > nbBuffer
			break
		end
	end
	# --- Last timeStamp and rate
	timeFinal = time();
	# --- Getting effective rate
	radioRate	  = radio.rx.samplingRate;
    effectiveRate = getRate(timeInit,timeFinal,nS);
	# --- Free all and return
	close(radio);
	return (radioRate,effectiveRate);
    end
end

max_stars_repo_name: RGerzaguet/UHDBindings

max_issues_repo_name: RGerzaguet/UHDBindings


text:
using Test
using OpenQuantumSystems
using LinearAlgebra
using SparseArrays

# https://github.com/JuliaIO/Suppressor.jl/blob/master/src/Suppressor.jl
macro suppress(block)
    quote
        if ccall(:jl_generating_output, Cint, ()) == 0
            original_stdout = stdout
            out_rd, out_wr = redirect_stdout()
            out_reader = @async read(out_rd, String)

            original_stderr = stderr
            err_rd, err_wr = redirect_stderr()
            err_reader = @async read(err_rd, String)

            # approach adapted from https://github.com/JuliaLang/IJulia.jl/pull/667/files
            logstate = Base.CoreLogging._global_logstate
            logger = logstate.logger
            if logger.stream == original_stderr
                new_logstate =
                    Base.CoreLogging.LogState(typeof(logger)(err_wr, logger.min_level))
                Core.eval(Base.CoreLogging, Expr(:(=), :(_global_logstate), new_logstate))
            end
        end

        try
            $(esc(block))
        finally
            if ccall(:jl_generating_output, Cint, ()) == 0
                redirect_stdout(original_stdout)
                close(out_wr)

                redirect_stderr(original_stderr)
                close(err_wr)

                if logger.stream == stderr
                    Core.eval(Base.CoreLogging, Expr(:(=), :(_global_logstate), logstate))
                end
            end
        end
    end
end

@testset "aggregate" begin

    # D(op1::Array, op2::Array) = abs(norm(op1 - op2))

    mode1 = Mode(0.2, 1.0)
    mode2 = Mode(0.3, 2.0)
    Energy = [0.0, 200.0]
    mol1 = Molecule([mode1], 3, [2.0, 200.0])
    mol2 = Molecule([mode2], 3, [3.0, 300.0])
    aggCore = AggregateCore([mol1, mol2])
    aggCore.coupling[2, 3] = 50
    aggCore.coupling[3, 2] = 50
    aggTools = AggregateTools(aggCore)
    aggOperators = AggregateOperators(aggCore, aggTools; groundEnergy=true)

    agg = setupAggregate(aggCore; groundEnergy=true)

    @test agg.core == aggCore
    @test agg.tools == aggTools
    @test agg.operators == aggOperators

    agg = Aggregate(aggCore, nothing, nothing)
    agg_ = Aggregate(aggCore)
    @test agg_ == agg

    agg = setupAggregate(aggCore)
    aggTools_ = AggregateTools(agg)
    @test aggTools_ == aggTools

    aggOperators_ = AggregateOperators(agg)
    @test aggOperators_ == aggOperators

    agg = setupAggregate(aggCore)
    agg_ = Aggregate(aggCore, nothing, nothing)
    setupAggregate!(agg_)
    @test agg_ == agg

end

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems

max_forks_repo_name: detrin/OpenQuantumSystems


text:
## File for normal stacks. Use for test and experiments.

type NormalStack{T,D}
  # IO
  input::IOStream # Pointer on a file to avoid memory consumption
  output::Nullable{IOStream} # Pointer to an optional (Nullable) output file
  # Functions defining the behavior of the stack in case of push/pop
  push_condition::Function
  push_action::Function
  pop_condition::Function
  pop_action::Function
  # (Explicit) data and context
  data::Vector{Data{D}} # Data type is in base.jl
  context::Nullable{T}
  index::Int
end

function NormalStack(input::IOStream, context_type::DataType,
  data_type::DataType, push_action::Function, push_condition::Function,
  pop_action::Function, pop_condition::Function;
  index = 0, context = Nullable{context_type}(), output = Nullable{IOStream}())

  data = Vector{Data{data_type}}()
  NormalStack(input, output, push_condition, push_action, pop_condition,
  pop_action, data, context, index)
end

## Normal Stack from input file
function NormalStack(name::AbstractString, pop_action::Function,
  push_action::Function, pop_condition::Function, push_condition::Function,
  context_type::DataType, data_type::DataType)

  input = open(name, "r")
  output_name = name * "_out"
  output = open(output_name, "w")
  settings = get_settings(input)

  NormalStack(input, context_type, data_type, push_action, push_condition, pop_action,
  pop_condition; output= Nullable{IOStream}(output))
end


text:
using Documenter, GPlot

makedocs(
    modules = [GPlot],
    format = Documenter.HTML(
        # Use clean URLs, unless built as a "local" build
        prettyurls = !("local" in ARGS),
        assets = ["assets/custom.css"],
        ),
    sitename = "GPlot.jl",
    authors  = "Thibaut Lienart",
    pages    = [
        "Home" => "index.md",
        "Manual" => [
            "Installation"         => "man/installation.md",
            "Quick start"          => "man/quickstart.md",
            "Line & Scatter plots" => "man/line-scatter.md",
            "Hist & Bar plots"     => "man/hist-bar.md",
            "Legend"               => "man/legend.md",
            "Axes and Axis"        => "man/figure-axes-axis.md",
            "Fill between lines"   => "man/fill-between.md",
            "Styling"              => "man/styling.md",
            "Boxplot"              => "man/boxplot.md",
            "Heatmap"              => "man/heatmap.md",
            "Annotations"          => "man/annotations.md",
            "Global settings"      => "man/global-settings.md",
            ],
        "Appendix" => [
            "Fonts" => "appendix/fonts.md",
            ],
        "Library" => [
            "Public"    => "lib/public.md",
            "Internals" => "lib/internals.md",
            ],
        ], # end pages
)

deploydocs(
    repo = "github.com/tlienart/GPlot.jl.git"
)


text:
# # Utility functions used in `demo_trsp.ipynb`

using MeshArrays, FortranFiles

# #### Function that reads files that were generated by `trsp_prep`

function trsp_read(mygrid::gcmgrid,dirIn::String)
    #read data from files generated by trsp_prep
    TrspX=mygrid.read(dirIn*"TrspX.bin",MeshArray(mygrid,Float32))
    TrspY=mygrid.read(dirIn*"TrspY.bin",MeshArray(mygrid,Float32))
    TauX=mygrid.read(dirIn*"TauX.bin",MeshArray(mygrid,Float32))
    TauY=mygrid.read(dirIn*"TauY.bin",MeshArray(mygrid,Float32))
    SSH=mygrid.read(dirIn*"SSH.bin",MeshArray(mygrid,Float32))

    return TrspX, TrspY, TauX, TauY, SSH
end

# #### Function that generates small binary files (2D) from large netcdf ones (4D)

function trsp_prep(mygrid::gcmgrid,GridVariables::Dict,dirOut::String="")

    #wind stress
    fileName="nctiles_climatology/oceTAUX/oceTAUX"
    oceTAUX=Main.read_nctiles(fileName,"oceTAUX",mygrid)
    fileName="nctiles_climatology/oceTAUY/oceTAUY"
    oceTAUY=Main.read_nctiles(fileName,"oceTAUY",mygrid)
    oceTAUX=mask(oceTAUX,0.0)
    oceTAUY=mask(oceTAUY,0.0)

    #sea surface height anomaly
    fileName="nctiles_climatology/ETAN/ETAN"
    ETAN=Main.read_nctiles(fileName,"ETAN",mygrid)
    fileName="nctiles_climatology/sIceLoad/sIceLoad"
    sIceLoad=Main.read_nctiles(fileName,"sIceLoad",mygrid)
    rhoconst=1029.0
    myssh=(ETAN+sIceLoad./rhoconst)
    myssh=mask(myssh,0.0)

    #seawater transports
    fileName="nctiles_climatology/UVELMASS/UVELMASS"
    U=Main.read_nctiles(fileName,"UVELMASS",mygrid)
    fileName="nctiles_climatology/VVELMASS/VVELMASS"
    V=Main.read_nctiles(fileName,"VVELMASS",mygrid)
    U=mask(U,0.0)
    V=mask(V,0.0)

    #time averaging and vertical integration
    TrspX=similar(GridVariables["DXC"])
    TrspY=similar(GridVariables["DYC"])
    TauX=similar(GridVariables["DXC"])
    TauY=similar(GridVariables["DYC"])
    SSH=similar(GridVariables["XC"])

    for i=1:mygrid.nFaces
        tmpX=mean(U.f[i],dims=4)
        tmpY=mean(V.f[i],dims=4)
        for k=1:length(GridVariables["RC"])
            tmpX[:,:,k]=tmpX[:,:,k].*GridVariables["DYG"].f[i]
            tmpX[:,:,k]=tmpX[:,:,k].*GridVariables["DRF"][k]
            tmpY[:,:,k]=tmpY[:,:,k].*GridVariables["DXG"].f[i]
            tmpY[:,:,k]=tmpY[:,:,k].*GridVariables["DRF"][k]
        end
        TrspX.f[i]=dropdims(sum(tmpX,dims=3),dims=(3,4))
        TrspY.f[i]=dropdims(sum(tmpY,dims=3),dims=(3,4))
        TauX.f[i]=dropdims(mean(oceTAUX.f[i],dims=3),dims=3)
        TauY.f[i]=dropdims(mean(oceTAUY.f[i],dims=3),dims=3)
        SSH.f[i]=dropdims(mean(myssh.f[i],dims=3),dims=3)
    end

    if !isempty(dirOut)
        write_bin(TrspX,dirOut*"TrspX.bin")
        write_bin(TrspY,dirOut*"TrspY.bin")
        write_bin(TauX,dirOut*"TauX.bin")
        write_bin(TauY,dirOut*"TauY.bin")
        write_bin(SSH,dirOut*"SSH.bin")
    end

    return TrspX, TrspY, TauX, TauY, SSH
end

# #### Function that writes a `MeshArray` object to a binary file using `FortranFile`

function write_bin(inFLD::MeshArray,filOut::String)
    recl=prod(inFLD.grid.ioSize)*4
    tmp=Float32.(convert2gcmfaces(inFLD))
    println("saving to file: "*filOut)
    f =  FortranFile(filOut,"w",access="direct",recl=recl,convert="big-endian")
    write(f,rec=1,tmp)
    close(f)
end

max_stars_repo_name: lmilechin/MeshArrayNotebooks

max_issues_repo_name: lmilechin/MeshArrayNotebooks

max_forks_repo_name: lmilechin/MeshArrayNotebooks


text:
"""
This function generates an Array of floats where each entry represents the RHS of the nodal balance equations. The corresponding values are the net-load values for each node and each time-step
"""
function active_timeseries_netinjection(sys::PowerSystems.PowerSystem)

    tsnetinjection =  zeros(Float64, length(sys.buses), sys.time_periods)

    # Note: This function has separate loops for load and resources such that it is possible to break the loops when there is no
    # fix loads or fix resources in the bus and in that way avoid looping through all the buses*time_steps`

    for source in sys.generators

         typeof(source) <: Array{<:ThermalGen} ? continue : (isa(source, Nothing) ? continue : true)

         for b in sys.buses

             for t in 1:sys.time_periods

                 fixed_source = [fs.tech.installedcapacity * values(fs.scalingfactor)[t] for fs in source if (fs.bus == b && isa(fs,fix_resource) && fs.available)]

                 isempty(fixed_source) ? break : fixed_source = tsnetinjection[b.number,t] -= sum(fixed_source)

             end

         end

     end

     for b in sys.buses

             for t in 1:sys.time_periods

             staticload = [sl.maxactivepower * values(sl.scalingfactor)[t] for sl in sys.loads if sl.bus == b]

             isempty(staticload) ? break : tsnetinjection[b.number,t] += sum(staticload)

         end
     end

    return  tsnetinjection
end

function reactive_timeseries_netinjection(sys::PowerSystems.PowerSystem)

    tsnetinjection =  zeros(Float64, length(sys.buses), sys.time_periods)

    # Note: This function has separate loops for load and resources such that it is possible to break the loops when there is no
    # fix loads or fix resources in the bus and in that way avoid looping through all the buses*time_steps`

    for source in sys.generators

         typeof(source) <: Array{<:ThermalGen} ? continue : (isa(source, Nothing) ? continue : true)

         for b in sys.buses

             for t in 1:sys.time_periods

                 fixed_source = [fs.tech.installedcapacity*sin(acos(fs.tech.powerfactor))* values(fs.scalingfactor)[t] for fs in source if (fs.bus == b && isa(fs,fix_resource))]

                 isempty(fixed_source) ? break : fixed_source = tsnetinjection[b.number,t] -= sum(fixed_source)

             end

         end

     end

     for b in sys.buses

             for t in 1:sys.time_periods

             staticload = [sl.maxreactivepower * values(sl.scalingfactor)[t] for sl in sys.loads if sl.bus == b]

             isempty(staticload) ? break : tsnetinjection[b.number,t] += sum(staticload)

         end
     end

    return  tsnetinjection
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Scenario-decomposition spconstructs #
# =================================== #
function EWS(stochasticprogram::StochasticProgram, structure::ScenarioDecompositionStructure)
    # Restore structure if optimization has been run before
    restore_structure!(optimizer(stochasticprogram))
    # Dispatch on scenarioproblems
    return EWS_horizontal(scenarioproblems(structure))
end
function EWS_horizontal(scenarioproblems::ScenarioProblems)
    return outcome_mean(subproblems(scenarioproblems), probability.(scenarios(scenarioproblems)))
end
function EWS_horizontal(scenarioproblems::DistributedScenarioProblems)
    partial_ews = Vector{Float64}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_ews[i] = remotecall_fetch(
                w,
                scenarioproblems[w-1]) do sp
                    scenarioproblems = fetch(sp)
                    num_scenarios(scenarioproblems) == 0 && return 0.0
                    return outcome_mean(subproblems(scenarioproblems),
                                        probability.(scenarios(scenarioproblems)))
                end
        end
    end
    return sum(partial_ews)
end

function statistical_EWS(::StochasticProgram, structure::ScenarioDecompositionStructure)
    return statistical_EWS_horizontal(scenarioproblems(structure))
end
function statistical_EWS_horizontal(scenarioproblems::ScenarioProblems)
    # Welford algorithm on WS subproblems
    return welford(subproblems(scenarioproblems),
                   probability.(scenarios(scenarioproblems)))
end
function statistical_EWS_horizontal(scenarioproblems::DistributedScenarioProblems)
    partial_welfords = Vector{Tuple{Float64,Float64,Float64,Int}}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_welfords[i] = remotecall_fetch(
                w,
                scenarioproblems[w-1]) do sp
                    scenarioproblems = fetch(sp)
                    num_scenarios(scenarioproblems) == 0 && return 0.0, 0.0, 0.0, 0
                    return welford(subproblems(scenarioproblems),
                                   probability.(scenarios(scenarioproblems)))
                end
        end
    end
    𝔼WS, σ², _ = reduce(aggregate_welford, partial_welfords)
    return 𝔼WS, sqrt(σ²)
end


text:
using ProgressMeter
using Dates

function learn_basic()
    n = 50
    p = Progress(n, dt=0.1, desc="Computing...", color=:green)
    for i in 1:n
        sleep(0.1)
        next!(p)
    end
    finish!(p)
end

function learn_showvalues()
    n = 50
    p = Progress(n, dt=0.1, desc="Computing...", color=:green)
    for i in 1:n
        sleep(0.1)
        next!(p, showvalues = [("Time",Dates.now())])
    end
    finish!(p)
end

function learn_update()
    n = 50
    p = Progress(n, dt=0.1, desc="Computing...", color=:green)
    for i in 1:n
        sleep(0.1)
        ProgressMeter.update!(p, i, showvalues = [("Time",Dates.now())])
    end
    finish!(p)
end

#learn_basic()
#learn_showvalues()
learn_update()

nothing