text:
module Iterators
using Base

import Base: start, next, done, count, take, eltype, length

export
    count,
    take,
    takestrict,
    drop,
    cycle,
    repeated,
    repeatedly,
    chain,
    product,
    distinct,
    partition,
    groupby,
    imap,
    subsets,
    iterate,
    @itr

# Infinite counting

immutable Count{S<:Number}
    start::S
    step::S
end

eltype{S}(it::Count{S}) = S

count(start::Number, step::Number) = Count(promote(start, step)...)
count(start::Number)               = Count(start, one(start))
count()                            = Count(0, 1)

start(it::Count) = it.start
next(it::Count, state) = (state, state + it.step)
done(it::Count, state) = false

# Iterate through the first n elements

immutable Take{I}
    xs::I
    n::Int
end

eltype(it::Take) = eltype(it.xs)

take(xs, n::Int) = Take(xs, n)

start(it::Take) = (it.n, start(it.xs))

function next(it::Take, state)
    n, xs_state = state
    v, xs_state = next(it.xs, xs_state)
    return v, (n - 1, xs_state)
end

function done(it::Take, state)
    n, xs_state = state
    return n <= 0 || done(it.xs, xs_state)
end

# Iterate through the first n elements, throwing an exception if
# fewer than n items ar encountered.

immutable TakeStrict{I}
    xs::I
    n::Int
end

eltype(it::TakeStrict) = eltype(it.xs)

takestrict(xs, n::Int) = TakeStrict(xs, n)

start(it::TakeStrict) = (it.n, start(it.xs))

function next(it::TakeStrict, state)
    n, xs_state = state
    v, xs_state = next(it.xs, xs_state)
    return v, (n - 1, xs_state)
end

function done(it::TakeStrict, state)
    n, xs_state = state
    if n <= 0
        return true
    elseif done(it.xs, xs_state)
        error("In takestrict(xs, n), xs had fewer than n items to take.")
    else
        return false
    end
end

function length(it::TakeStrict)
    return it.n
end

# Iterator through all but the first n elements

immutable Drop{I}
    xs::I
    n::Int
end

eltype(it::Drop) = eltype(it.xs)

drop(xs, n::Int) = Drop(xs, n)

function start(it::Drop)
    xs_state = start(it.xs)
    for i in 1:it.n
        if done(it.xs, xs_state)
            break
        end

        _, xs_state = next(it.xs, xs_state)
    end
    xs_state
end

next(it::Drop, state) = next(it.xs, state)
done(it::Drop, state) = done(it.xs, state)

# Cycle an iterator forever

immutable Cycle{I}
    xs::I
end

eltype(it::Cycle) = eltype(it.xs)

cycle(xs) = Cycle(xs)

function start(it::Cycle)
    s = start(it.xs)
    return s, done(it.xs, s)
end

function next(it::Cycle, state)
    s, d = state
    if done(it.xs, s)
        s = start(it.xs)
    end
    v, s = next(it.xs, s)
    return v, (s, false)
end

done(it::Cycle, state) = state[2]

# Repeat an object n (or infinitely many) times.

immutable Repeat{O}
    x::O
    n::Int
end

eltype{O}(it::Repeat{O}) = O
length(it::Repeat) = it.n

repeated(x, n) = Repeat(x, n)

start(it::Repeat) = it.n
next(it::Repeat, state) = (it.x, state - 1)
done(it::Repeat, state) = state <= 0

immutable RepeatForever{O}
    x::O
end

eltype{O}(r::RepeatForever{O}) = O

repeated(x) = RepeatForever(x)

start(it::RepeatForever) = nothing
next(it::RepeatForever, state) = (it.x, nothing)
done(it::RepeatForever, state) = false

# Repeat a function application n (or infinitely many) times.

immutable RepeatCall
    f::Function
    n::Int
end

length(it::RepeatCall) = it.n
repeatedly(f, n) = RepeatCall(f, n)

start(it::RepeatCall) = it.n
next(it::RepeatCall, state) = (it.f(), state - 1)
done(it::RepeatCall, state) = state <= 0

immutable RepeatCallForever
    f::Function
end

repeatedly(f) = RepeatCallForever(f)

start(it::RepeatCallForever) = nothing
next(it::RepeatCallForever, state) = (it.f(), nothing)
done(it::RepeatCallForever, state) = false

# Concatenate the output of n iterators

immutable Chain
    xss::Vector{Any}
    function Chain(xss...)
        new(Any[xss...])
    end
end

function eltype(it::Chain)
    try
        typejoin([eltype(xs) for xs in it.xss]...)
    catch
        Any
    end
end

chain(xss...) = Chain(xss...)

function start(it::Chain)
    i = 1
    xs_state = nothing
    while i <= length(it.xss)
        xs_state = start(it.xss[i])
        if !done(it.xss[i], xs_state)
            break
        end
        i += 1
    end
    return i, xs_state
end

function next(it::Chain, state)
    i, xs_state = state
    v, xs_state = next(it.xss[i], xs_state)
    while done(it.xss[i], xs_state)
        i += 1
        if i > length(it.xss)
            break
        end
        xs_state = start(it.xss[i])
    end
    return v, (i, xs_state)
end

done(it::Chain, state) = state[1] > length(it.xss)

# Cartesian product as a sequence of tuples

immutable Product
    xss::Vector{Any}
    function Product(xss...)
        new(Any[xss...])
    end
end

eltype(p::Product) = tuple(map(eltype, p.xss)...)
length(p::Product) = prod(map(length, p.xss))

product(xss...) = Product(xss...)

function start(it::Product)
    n = length(it.xss)
    js = Any[start(xs) for xs in it.xss]
    if n == 0
        return js, nothing
    end
    for i = 1:n
        if done(it.xss[i], js[i])
            return js, nothing
        end
    end
    vs = Array(Any, n)
    for i = 1:n
        vs[i], js[i] = next(it.xss[i], js[i])
    end
    return js, vs
end

function next(it::Product, state)
    js = copy(state[1])
    vs = copy(state[2])
    ans = tuple(vs...)

    n = length(it.xss)
    for i in 1:n
        if !done(it.xss[i], js[i])
            vs[i], js[i] = next(it.xss[i], js[i])
            break
        elseif i == n
            vs = nothing
            break
        end

        js[i] = start(it.xss[i])
        vs[i], js[i] = next(it.xss[i], js[i])
    end
    return ans, (js, vs)
end

done(it::Product, state) = state[2] === nothing

# Filter out reccuring elements.

immutable Distinct{I}
    xs::I

    # Map elements to the index at which it was first seen, so given an iterator
    # state (index) we can test if an element has previously been observed.
    seen::Dict{Any, Int}

    Distinct(xs) = new(xs, Dict{Any, Int}())
end

eltype(it::Distinct) = eltype(it.xs)

distinct{I}(xs::I) = Distinct{I}(xs)

function start(it::Distinct)
    start(it.xs), 1
end

function next(it::Distinct, state)
    s, i = state
    x, s = next(it.xs, s)
    it.seen[x] = i
    i += 1

    while !done(it.xs, s)
        y, t = next(it.xs, s)
        if !haskey(it.seen, y) || it.seen[y] >= i
            break
        end
        s = t
        i += 1
    end

    x, (s, i)
end

done(it::Distinct, state) = done(it.xs, state[1])

# Group output from at iterator into tuples.
# E.g.,
#   partition(count(1), 2) = (1,2), (3,4), (5,6) ...
#   partition(count(1), 2, 1) = (1,2), (2,3), (3,4) ...
#   partition(count(1), 2, 3) = (1,2), (4,5), (7,8) ...

immutable Partition{I}
    xs::I
    n::Int
    step::Int
end

eltype(it::Partition) = tuple(fill(eltype(it.xs),it.n)...)

function partition(xs, n::Int)
    Partition(xs, n, n)
end

function partition(xs, n::Int, step::Int)
    if step < 1
        error("Partition step must be at least 1.")
    end

    Partition(xs, n, step)
end

function start(it::Partition)
    p = Array(eltype(it.xs), it.n)
    s = start(it.xs)
    for i in 1:(it.n - 1)
        if done(it.xs, s)
            break
        end
        (p[i], s) = next(it.xs, s)
    end
    (s, p)
end

function next(it::Partition, state)
    (s, p0) = state
    (x, s) = next(it.xs, s)
    ans = p0; ans[end] = x

    p = similar(p0)
    overlap = max(0, it.n - it.step)
    for i in 1:overlap
        p[i] = ans[it.step + i]
    end

    # when step > n, skip over some elements
    for i in 1:max(0, it.step - it.n)
        if done(it.xs, s)
            break
        end
        (x, s) = next(it.xs, s)
    end

    for i in (overlap + 1):(it.n - 1)
        if done(it.xs, s)
            break
        end

        (x, s) = next(it.xs, s)
        p[i] = x
    end

    (tuple(ans...), (s, p))
end

done(it::Partition, state) = done(it.xs, state[1])

# Group output from an iterator based on a key function.
# Consecutive entries from the iterator with the same
# key value will be returned in a single array.
# Inspired by itertools.groupby in python.
# E.g.,
#   x = ["face", "foo", "bar", "book", "baz", "zzz"]
#   groupby(x, z -> z[1]) =
#       ["face", "foo"]
#       ["bar", "book", "baz"]
#       ["zzz"]
immutable GroupBy{I}
    xs::I
    keyfunc::Function
end

eltype{I}(it::GroupBy{I}) = I
eltype{I<:Ranges}(it::GroupBy{I}) = Array{eltype(it.xs),}

function groupby(xs, keyfunc)
    GroupBy(xs, keyfunc)
end

function start(it::GroupBy)
    s = start(it.xs)
    prev_value = nothing
    prev_key = nothing
    return (s, (prev_key, prev_value))
end

function next(it::GroupBy, state)
    (s, (prev_key, prev_value)) = state
    values = Array(eltype(it.xs), 0)
    # We had a left over value from the last time the key changed.
    if prev_value != nothing || prev_key != nothing
        push!(values, prev_value)
    end
    prev_value = nothing
    while !done(it.xs, s)
        (x, s) = next(it.xs, s)
        key = it.keyfunc(x)
        # Did the key change?
        if prev_key != nothing && key != prev_key
            prev_key = key
            prev_value = x
            break
        else
            push!(values, x)
        end
        prev_key = key
    end
    # We either reached the end of the input or the key changed,
    # either way emit what we have so far.
    return (values, (s, (prev_key, prev_value)))
end

function done(it::GroupBy, state)
    return state[2][2] == nothing && done(it.xs, state[1])
end

# Like map, except returns the output as an iterator.  The iterator
# is done when any of the input iterators have been exhausted.
# E.g.,
#   imap(+, count(), [1, 2, 3]) = 1, 3, 5 ...
immutable IMap
    mapfunc::Base.Callable
    xs::Vector{Any}
end

function imap(mapfunc, it1, its...)
    IMap(mapfunc, Any[it1, its...])
end

function start(it::IMap)
    map(start, it.xs)
end

function next(it::IMap, state)
    next_result = map(next, it.xs, state)
    return (
        it.mapfunc(map(x -> x[1], next_result)...),
        map(x -> x[2], next_result)
    )
end

function done(it::IMap, state)
    any(map(done, it.xs, state))
end

# Iterate over all subsets of a collection

immutable Subsets
    xs
end

eltype(it::Subsets) = Array{eltype(it.xs),1}
length(it::Subsets) = 1 << length(it.xs)

function subsets(xs)
    Subsets(xs)
end

function start(it::Subsets)
    # one extra bit to indicated that we are at the end
    BitVector(length(it.xs) + 1)
end

function next(it::Subsets, state)
    ss = Array(eltype(it.xs), 0)
    for i = 1:length(it.xs)
        if state[i]
            push!(ss, it.xs[i])
        end
    end

    state = copy(state)
    state[1] = !state[1]
    for i in 2:length(state)
        if !state[i-1]
            state[i] = !state[i]
        else
            break
        end
    end

    (ss, state)
end

function done(it::Subsets, state)
    state[end]
end

# Unfolding (anamorphism)
# Outputs the stream: seed, f(seed), f(f(seed)), ...

immutable Iterate{T}
    f::Function
    seed::T
end

iterate(f, seed) = Iterate(f, seed)
start(it::Iterate) = it.seed
next(it::Iterate, state) = (state, it.f(state))
done(it::Iterate, state) = (state==None)

using Base.Meta

## @itr macro for auto-inlining in for loops
#
# it dispatches on macros defined below

macro itr(ex)
    isexpr(ex, :for) || error("@itr macro expects a for loop")
    isexpr(ex.args[1], :(=)) || error("malformed or unsupported for loop in @itr macro")
    isexpr(ex.args[1].args[2], :call) || error("@itr macro expects an iterator call, e.g. @itr for (x,y) = zip(a,b)")
    iterator = ex.args[1].args[2].args[1]
    ex.args[1].args[2] = Expr(:tuple, ex.args[1].args[2].args[2:end]...)
    if iterator == :zip
        rex = :(@zip($(esc(ex))))
    elseif iterator == :enumerate
        rex = :(@enumerate($(esc(ex))))
    elseif iterator == :take
        rex = :(@take($(esc(ex))))
    elseif iterator == :takestrict
        rex = :(@takestrict($(esc(ex))))
    elseif iterator == :drop
        rex = :(@drop($(esc(ex))))
    elseif iterator == :chain
        rex = :(@chain($(esc(ex))))
    else
        error("unknown or unsupported iterator $iterator in @itr macro")
    end
    return rex
end

macro zip(ex)
    @assert ex.head == :for
    @assert ex.args[1].head == :(=)
    isexpr(ex.args[1].args[1], :tuple) || error("@zip macro needs explicit tuple arguments")
    isexpr(ex.args[1].args[2], :tuple) || error("@zip macro needs explicit tuple arguments")
    n = length(ex.args[1].args[1].args)
    length(ex.args[1].args[2].args) == n || error("unequal tuple sizes in @zip macro")
    body = esc(ex.args[2])
    vars = map(esc, ex.args[1].args[1].args)
    iters = map(esc, ex.args[1].args[2].args)
    states = [gensym("s") for i=1:n]
    as = Any[Expr(:call, :(Base.start), iters[i]) for i=1:n]
    startex = Expr(:(=), Expr(:tuple, states...), Expr(:tuple, as...))
    ad = Any[Expr(:call, :(Base.done), iters[i], states[i]) for i = 1:n]
    notdoneex = Expr(:call, :(!), Expr(:||, ad...))
    nextex = Expr(:(=), Expr(:tuple, Any[Expr(:tuple, vars[i], states[i]) for i=1:n]...),
                        Expr(:tuple, Any[Expr(:call, :(Base.next), iters[i], states[i]) for i=1:n]...))
    Expr(:let, Expr(:block,
        startex,
        Expr(:while, notdoneex, Expr(:block,
            nextex,
            body))),
        states...)
end

macro enumerate(ex)
    @assert ex.head == :for
    @assert ex.args[1].head == :(=)
    isexpr(ex.args[1].args[1], :tuple) || error("@enumerate macro needs an explicit tuple argument")
    length(ex.args[1].args[1].args) == 2 || error("lentgh of tuple must be 2 in @enumerate macro")
    body = esc(ex.args[2])
    vars = map(esc, ex.args[1].args[1].args)
    if isexpr(ex.args[1].args[2], :tuple) && length(ex.args[1].args[2].args) == 1
        ex.args[1].args[2] = ex.args[1].args[2].args[1]
    end
    iter = esc(ex.args[1].args[2])
    ind = gensym("i")
    startex = Expr(:(=), ind, 0)
    forex = Expr(:(=), vars[2], iter)
    index = Expr(:(=), vars[1], ind)
    increx = Expr(:(=), ind, Expr(:call, :(+), ind, 1))
    Expr(:let, Expr(:block,
        startex,
        Expr(:for, forex, Expr(:block,
            increx,
            index,
            body))),
        ind)
end

# both @take and @takestrict use @_take

macro _take(ex, strict)
    mname = strict ? "takestrict" : "take"
    @assert ex.head == :for
    @assert ex.args[1].head == :(=)
    isexpr(ex.args[1].args[2], :tuple) || error("@$(mname) macro needs an explicit tuple argument")
    length(ex.args[1].args[2].args) == 2 || error("length of tuple must be 2 in @$(mname) macro")
    body = esc(ex.args[2])
    var = esc(ex.args[1].args[1])
    iter = esc(ex.args[1].args[2].args[1])
    n = esc(ex.args[1].args[2].args[2])
    ind = gensym("i")
    state = gensym("s")
    startex = Expr(:block,
        Expr(:(=), ind, 0),
        Expr(:(=), state, Expr(:call, :(Base.start), iter)))
    notdoneex = Expr(:call, :(!), Expr(:||,
        Expr(:call, :(>=), ind, n),
        Expr(:call, :(Base.done), iter, state)))
    nextex = Expr(:block,
        Expr(:(=), Expr(:tuple, var, state),
                   Expr(:call, :(Base.next), iter, state)),
        Expr(:(=), ind, Expr(:call, :(+), ind, 1)))
    if strict
        checkex = Expr(:if, Expr(:call, :(<), ind, n),
            Expr(:call, :error, "in takestrict(xs, n), xs had fewer than n items to take."))
    else
        checkex = :nothing
    end

    Expr(:let, Expr(:block,
        startex,
        Expr(:while, notdoneex, Expr(:block,
            nextex,
            body)),
        checkex),
        ind, state)
end

macro take(ex)
    :(@_take($(esc(ex)), false))
end

macro takestrict(ex)
    :(@_take($(esc(ex)), true))
end

macro drop(ex)
    @assert ex.head == :for
    @assert ex.args[1].head == :(=)
    isexpr(ex.args[1].args[2], :tuple) || error("@drop macro needs an explicit tuple argument")
    length(ex.args[1].args[2].args) == 2 || error("length of tuple must be 2 in @drop macro")
    body = esc(ex.args[2])
    var = esc(ex.args[1].args[1])
    iter = esc(ex.args[1].args[2].args[1])
    n = esc(ex.args[1].args[2].args[2])
    ind = gensym("i")
    state = gensym("s")
    startex = Expr(:block,
        Expr(:(=), ind, 0),
        Expr(:(=), state, Expr(:call, :(Base.start), iter)))
    notdoneex1 = Expr(:call, :(!), Expr(:||,
        Expr(:call, :(>=), ind, n),
        Expr(:call, :(Base.done), iter, state)))
    nextex1 = Expr(:block,
        Expr(:(=), Expr(:tuple, var, state),
                   Expr(:call, :(Base.next), iter, state)),
        Expr(:(=), ind, Expr(:call, :(+), ind, 1)))
    notdoneex2 = Expr(:call, :(!), Expr(:call, :(Base.done), iter, state))
    nextex2 = Expr(:(=), Expr(:tuple, var, state), Expr(:call, :(Base.next), iter, state))

    Expr(:let, Expr(:block,
        startex,
        Expr(:while, notdoneex1, nextex1),
        Expr(:while, notdoneex2, Expr(:block,
            nextex2,
            body))),
        ind, state)
end

macro chain(ex)
    @assert ex.head == :for
    @assert ex.args[1].head == :(=)
    isexpr(ex.args[1].args[2], :tuple) || error("@chain macro needs explicit tuple arguments")
    n = length(ex.args[1].args[2].args)
    body = esc(ex.args[2])
    var = esc(ex.args[1].args[1])
    iters = map(esc, ex.args[1].args[2].args)
    states = [gensym("s") for i=1:n]

    cycleex = [Expr(:block,
        Expr(:(=), states[i], Expr(:call, :(Base.start), iters[i])),
        Expr(:while, Expr(:call, :(!), Expr(:call, :(Base.done), iters[i], states[i])), Expr(:block,
            Expr(:(=), Expr(:tuple, var, states[i]), Expr(:call, :(Base.next), iters[i], states[i])),
            body))) for i = 1:n]

    Expr(:let, Expr(:block, cycleex...), states...)
end

end # module Iterators


text:
#Came from GLWindow/core.jl

"""
Create a new rectangle with x,y == 0,0 while taking the widths from the original
Rectangle
"""
zeroposition(r::SimpleRectangle{T}) where {T} = SimpleRectangle(zero(T), zero(T), r.w, r.h)

ext: jl


text:
function make_model(f, ϵsub, ϵ, cellL, thickness, order, lb, ub, filename)
    # TODO: call out to Python
    pts = chebpoints(order, lb, ub)
end

function get_model(order, lb, ub, filename)
    f = open(filename)
    function val(line)
        dat = split(line)
        parse(Float64, dat[1]) + parse(Float64, dat[2]) * im
    end
    vals = [val(line) for line in eachline(f)]
    chebinterp(vals, lb, ub)
end

# rrule for Chebyshev polynomial functor. TODO: support chebjacobian (or explicitly don't support it)
# TODO: support x real
function rrule(c::ChebPoly, x::AbstractVector)
    project_x = ProjectTo(x)
    y, Δy = chebgradient(c, x)
    pullback(∂y) = NoTangent(), project_x(∂y * Δy')
    y, pullback
end


text:
"""
    LennardJones <: AbstractPairPotential

Holds constants for Lennard-Jones potential
    E = C12/R^12 - C6/R^6
"""
mutable struct LennardJones <: AbstractPairPotential
    C6::Float64
    C12::Float64
    LennardJones() = new(0.0,0.0)
    LennardJones(c6,c12) = new(c6,c12)
end

function Base.show(io::IO, potential::LennardJones; energy_unit="cm^-1")
    if potential.C6 > 0 && potential.C12 > 0
        ϵ = energy_to(potential.C6^2/(4*potential.C12), energy_unit)
        σ = (potential.C12/potential.C6)^(1/6)
        print(io, "$(typeof(potential)),  ϵ=$(ϵ) $(energy_unit), σ=$(σ) Å")
    else
        print(io, "$(typeof(potential)),  C6=$(potential.C6), C12=$(potential.C12)")
    end
end

function potential_variables(potential::LennardJones,
                             c::AbstractCluster, indices)
    r=distances(c, indices[1], indices[2])
    return [r^-6 r^-12]
end

function get_potential!(potential::LennardJones, x6, x12)
    potential.C6 = -x6
    potential.C12 = x12
    return potential
end

function (p::LennardJones)(r::Number)::Float64
    r6 = r^-6
    muladd(r6, p.C12, -p.C6) * r6
end


text:
############################################################################################
include("abstractarray.jl")
include("tuple.jl")
include("namedtuple.jl")

############################################################################################
# Export


text:
path = dirname(@__FILE__)
using DSGE

m = AnSchorfheide()

@testset "Testing getting MeansBands Input and Output Files" begin
    @test get_meansbands_input_file(m, :mode, :none, :histobs) ==
        joinpath(rawpath(m, "forecast"), "histobs_cond=none_para=mode_vint=" * data_vintage(m) * ".jld2")
    @test get_meansbands_input_file("a", ["b"], :mode, :none, :histobs, fileformat = :jld2) ==
        joinpath("a", "histobs_b_cond=none_para=mode.jld2")

    @test get_meansbands_output_file(m, :mode, :none, :histobs) ==
        joinpath(workpath(m, "forecast"), "mbhistobs_cond=none_para=mode_vint=" * data_vintage(m) * ".jld2")
    @test get_meansbands_output_file("a", ["b"], :mode, :none, :histobs, fileformat = "jld2") ==
        joinpath("a", "mbhistobs_b_cond=none_para=mode.jld2")

    run(`rm -r a`)

    @test Matrix(read_mb("$path/../reference/mbhistobs.jld2").means) == Matrix(load("$path/../reference/mbhistobs.jld2", "mb").means)
    @test Matrix(read_mb("$path/../reference/mbhistobs.jld2", "$path/../reference/mbhistobs.jld2").means) ==
                 Matrix(load("$path/../reference/mbhistobs.jld2", "mb").means)
    @test_throws AssertionError read_bdd_and_unbdd_mb("", "")
    @test Matrix(read_bdd_and_unbdd_mb("$path/../reference/mbbddforecastobs.jld2", "$path/../reference/mbforecastobs.jld2").means) == Matrix(load("$path/../reference/bdd_unbdd.jld2", "bdd_unbdd"))
    @test Matrix(read_mb("$path/../reference/mbbddforecastobs.jld2", "$path/../reference/mbbddforecastobs.jld2").means) ==
                 Matrix(load("$path/../reference/mbbddforecastobs.jld2", "mb").means)
end

@testset "Testing add_requisite_output_vars_meansbands" begin
    @test DSGE.add_requisite_output_vars_meansbands([:histobs, :shockdecpseudo]) ==
        [:histobs, :shockdecpseudo, :dettrendpseudo, :trendpseudo, :histforecastpseudo]
    @test DSGE.add_requisite_output_vars_meansbands([:histobs, :shockdecobs]) ==
        [:histobs, :shockdecobs, :dettrendobs, :trendobs, :histforecastobs]
end


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

# Set up environment for Julia OSX binary distribution
let
    ROOT = abspath(JULIA_HOME,"..")
    ENV["PATH"]="$JULIA_HOME:$(ENV["PATH"])"
    ENV["FONTCONFIG_PATH"] = joinpath(ROOT, "etc", "fonts")
    ENV["TK_LIBRARY"] = "/System/Library/Frameworks/Tk.framework/Versions/8.5/Resources/Scripts"
end


text:
using Unpack
using BenchmarkTools

function unpack_bc(w::Array{<:Tuple})
    Tuple((v->v[i]).(w) for i=1:length(w[1]))
end

g(a, b) = a+1im*b, a*b, convert(Int, round(a-b))
packed = g.(rand(100,100,1), rand(1,1,100))
println("Broadcast")
@btime unpack_bc($packed)
println("Cartesian")
@btime unpack($packed)


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -58.37772561872491 <= q <= 51.748312049639644)

                     add_NL_constraint(m, :(softplus(-0.3976822916974827 + -0.11888104611683081*softplus(-0.09559544092481564 + 0.23201220435158865*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + -0.9126138494455787*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + 0.24324641666113456*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + 0.1690500900766101*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + 0.20833731206579253*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))) + -0.8277675607937618*softplus(-0.41653739657549194 + 0.7547021537813423*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + 0.5787262183714494*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + -0.11492508246438904*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + 0.30062983013792843*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + -0.1276250114011983*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))) + 0.37430250255255215*softplus(0.9557400948867283 + -0.3128610884258398*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + 0.29518845150992634*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + -0.18188148649835556*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + 0.8330182900675651*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + -0.8804319372999063*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))) + 0.9159784725268092*softplus(-0.32868009819829513 + 0.7041500634234112*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + -0.44147758504740464*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + -0.5937281428689025*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + -0.3198554528784734*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + -0.047330375221362964*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))) + -0.21239133839693736*softplus(-0.6560670298284208 + 0.11579976545709414*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + 0.7227941593613254*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + 0.5041216754022764*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + -0.9501104307252297*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + 0.7033943599032555*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5])))) + softplus(-0.962366690533889 + 0.21288820713294498*softplus(-0.09559544092481564 + 0.23201220435158865*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + -0.9126138494455787*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + 0.24324641666113456*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + 0.1690500900766101*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + 0.20833731206579253*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))) + 0.04621026197884115*softplus(-0.41653739657549194 + 0.7547021537813423*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + 0.5787262183714494*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + -0.11492508246438904*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + 0.30062983013792843*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + -0.1276250114011983*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))) + 0.3579305330052933*softplus(0.9557400948867283 + -0.3128610884258398*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + 0.29518845150992634*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + -0.18188148649835556*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + 0.8330182900675651*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + -0.8804319372999063*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))) + -0.40247732796428703*softplus(-0.32868009819829513 + 0.7041500634234112*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + -0.44147758504740464*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + -0.5937281428689025*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + -0.3198554528784734*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + -0.047330375221362964*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))) + -0.3937461182921762*softplus(-0.6560670298284208 + 0.11579976545709414*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + 0.7227941593613254*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + 0.5041216754022764*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + -0.9501104307252297*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + 0.7033943599032555*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5])))) + softplus(-0.9189915189234159 + -0.5125130010222394*softplus(-0.09559544092481564 + 0.23201220435158865*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + -0.9126138494455787*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + 0.24324641666113456*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + 0.1690500900766101*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + 0.20833731206579253*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))) + -0.3502499398907637*softplus(-0.41653739657549194 + 0.7547021537813423*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + 0.5787262183714494*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + -0.11492508246438904*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + 0.30062983013792843*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + -0.1276250114011983*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))) + -0.8274375080497247*softplus(0.9557400948867283 + -0.3128610884258398*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + 0.29518845150992634*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + -0.18188148649835556*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + 0.8330182900675651*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + -0.8804319372999063*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))) + -0.7757421130805153*softplus(-0.32868009819829513 + 0.7041500634234112*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + -0.44147758504740464*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + -0.5937281428689025*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + -0.3198554528784734*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + -0.047330375221362964*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))) + 0.6602598595684599*softplus(-0.6560670298284208 + 0.11579976545709414*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + 0.7227941593613254*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + 0.5041216754022764*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + -0.9501104307252297*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + 0.7033943599032555*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5])))) + softplus(-0.760456357412739 + 0.8974231726544617*softplus(-0.09559544092481564 + 0.23201220435158865*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + -0.9126138494455787*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + 0.24324641666113456*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + 0.1690500900766101*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + 0.20833731206579253*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))) + -0.7589867374120325*softplus(-0.41653739657549194 + 0.7547021537813423*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + 0.5787262183714494*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + -0.11492508246438904*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + 0.30062983013792843*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + -0.1276250114011983*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))) + -0.3014899851721653*softplus(0.9557400948867283 + -0.3128610884258398*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + 0.29518845150992634*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + -0.18188148649835556*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + 0.8330182900675651*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + -0.8804319372999063*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))) + 0.03263132145111802*softplus(-0.32868009819829513 + 0.7041500634234112*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + -0.44147758504740464*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + -0.5937281428689025*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + -0.3198554528784734*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + -0.047330375221362964*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))) + 0.8460059521417822*softplus(-0.6560670298284208 + 0.11579976545709414*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + 0.7227941593613254*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + 0.5041216754022764*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + -0.9501104307252297*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + 0.7033943599032555*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5])))) + softplus(-0.1364886291423204 + -0.7689732133011344*softplus(-0.09559544092481564 + 0.23201220435158865*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + -0.9126138494455787*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + 0.24324641666113456*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + 0.1690500900766101*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + 0.20833731206579253*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))) + 0.9275671307049969*softplus(-0.41653739657549194 + 0.7547021537813423*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + 0.5787262183714494*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + -0.11492508246438904*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + 0.30062983013792843*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + -0.1276250114011983*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))) + -0.311285022430249*softplus(0.9557400948867283 + -0.3128610884258398*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + 0.29518845150992634*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + -0.18188148649835556*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + 0.8330182900675651*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + -0.8804319372999063*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))) + -0.2373305013184268*softplus(-0.32868009819829513 + 0.7041500634234112*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + -0.44147758504740464*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + -0.5937281428689025*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + -0.3198554528784734*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + -0.047330375221362964*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))) + 0.7676258102518863*softplus(-0.6560670298284208 + 0.11579976545709414*softplus(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])) + 0.7227941593613254*softplus(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])) + 0.5041216754022764*softplus(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])) + -0.9501104307252297*softplus(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])) + 0.7033943599032555*softplus(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using ProbNumDiffEq
using ProbNumDiffEq: stack
using DifferentialEquations
using Statistics
using Plots

C1, C2 = "#107D79", "#FF9933"

# ODE definition as in DifferentialEquations.jl
function fitz(u, p, t)
    a, b, c = p
    return [
        c * (u[1] - u[1]^3 / 3 + u[2])
        -(1 / c) * (u[1] - a - b * u[2])
    ]
end
u0 = [-1.0; 1.0]
tspan = (0.0, 20.0)
p = (0.2, 0.2, 3.0)
prob = ODEProblem(fitz, u0, tspan, p)

appxsol = solve(prob, abstol=1e-12, reltol=1e-12)

############################################################################################
# Solve with integrator interface and make animation
integ = init(prob, EK0(order=1), adaptive=false, dt=7e-2)

anim = @animate for i in 1:(prob.tspan[2]/integ.dt)
    step!(integ)

    plot(
        integ.sol;
        color=[C1 C2],
        ribbon_width=3,
        xlims=prob.tspan,
        ylims=(-2.5, 2.5),
        ylabel="u(t)",
        xlabel="t",
        label=["u₁(t)" "u₂(t)"],
    )

    times = integ.t:0.1:prob.tspan[2]
    post = integ.sol(times).u
    plot!(
        times,
        stack(mean(post)),
        ribbon=3stack(std(post)),
        linestyle=:dot,
        color=[C1 C2],
        fillalpha=0.2,
        label="",
    )
    scatter!(integ.sol.t, stack(integ.sol.u), color=[C1 C2], label="")
    plot!(appxsol, linestyle=:dash, color=:black, label="")
    vline!([integ.t], label="", color=:black)
end

gif(anim, "./fitzhughnagumo_solve.gif")


text:
@recipe function f(ints::Vector{T}) where {T<:AbstractInterval}
    Dict(0=>ints)
end

# Plot intervals as a persistance diagram or barcode
@recipe function f(dgm::Dict{Int, Vector{T}};
                   maxoutdim=1,
                   skipzero = false) where {T<:AbstractInterval}
    # set default plot type
    seriestype --> :diagram

    # get interval dimensions
    vals = vcat(([i.b, i.d] for i in vcat(values(dgm)...))...)
    dims = [d for d in sort!(filter!(d->d<=maxoutdim, collect(keys(dgm))))]
    # compute range of intervals
    minval, maxval = extrema(filter!(!isinf, vals))

    if plotattributes[:seriestype] == :barcode # Barcode
        #xticks := minval:maxval
        xlims  := (minval, maxval*1.01)
        yticks := nothing
        ylims := (0.0, 1.05)
        legend := :none
        layout := (length(dims),1)
        for (pi, d) in enumerate(dims)
            dmaxval = max(maxval, foldl(max, filter(!isinf, map(i->i.d, dgm[d])); init=-Inf))
            dminval = min(minval, foldl(min, filter(!isinf, map(i->i.b, dgm[d])); init=Inf))
            points = map(i -> i.b == i.d, dgm[d])
            step = 1 / (skipzero ? sum(.!points) : length(dgm[d]))
            y = step
            idxs = sortperm(dgm[d])
            for (i, ispoint) in zip(dgm[d][idxs], points)
                xs, ys = [i.b, isinf(i.d) ? maxval : i.d], [y, y]
                ispoint && skipzero && continue
                @series begin
                    title := "Degree $d"
                    titlefontsize --> 12
                    linewidth --> 2
                    subplot := pi
                    primary := true
                    arrow := isinf(i.d)
                    seriescolor --> :blue
                    seriestype := ispoint ? :scatter : :path
                    markersize --> 2
                    xs, ys
                end
                y += step
            end
        end
    else # Persistance Diagram
        padding = (maxval - minval)*0.01 # 3%
        xs = map(l->map(i->i.b,l), dgm[d] for d in dims)
        ys = map(l->map(i->isinf(i.d) ? maxval : i.d,l), dgm[d] for d in dims)
        mkr = map(l->map(i->isinf(i.d) ? :rect : :circle,l), dgm[d] for d in dims)
        legend --> :bottomright
        legendtitle --> "Homology"
        # label --> ["Degree $(i-1)" for i in 1:length(dgm)]
        # xticks := minval:maxval
        xlims := (minval-padding, maxval+padding)
        xguide := "birth"
        # yticks := 0:maxval
        yguide := "death"
        ylims := (minval, maxval*1.03)

        @series begin
            primary := false
            seriestype := :path
            linecolor := :black
            [0, maxval], [0, maxval]
        end

        for i in 1:length(xs)
            idxs = skipzero ? xs[i] .!== ys[i] : BitArray(true for i in 1:length(xs[i]))
            @series begin
                primary := true
                seriestype := :scatter
                markershape := mkr[i][idxs]
                markersize --> 2
                label --> "Degree $(i-1)"
                # markerstrokecolor := col
                xs[i][idxs], ys[i][idxs]
            end
        end
    end
end


text:
abstract type Person
end

abstract type Musician <: Person
end

mutable struct Rockstar <: Musician
	name::String
	instrument::String
	bandName::String
	headbandColor::String
end

struct ClassicMusician <: Musician
	name::String
	instrument::String
end

mutable struct Physicist <: Person
	name::String
	sleepHours::Float64
	favouriteLanguage::String
end
#%%
aure = Physicist("Aurelio", 6, "Julia")
aure.sleepHours = 8
aure_musician = ClassicMusician("Aurelio", "Violin")
aure_musician.instrument="Cello"

ricky = Rockstar("Riccardo", "Voice", "Black Lotus", "red", 2)

ricky.headbandColor

#%%
function introduceMe(person::Person)
    println("Hello, my name is $(person.name).")
end

function introduceMe(person::Musician)
    println("Hello, my name is $(person.name) and I play $(person.instrument).")
end

function introduceMe(person::Rockstar)
	if person.instrument == "Voice"
		println("Hello, my name is $(person.name) and I sing.")
	else
		println("Hello, my name is $(person.name) and I play $(person.instrument).")
	end

	println("My band name is $(person.bandName) and my favourite headband colour is $(person.headbandColor)!")
end
#%%
introduceMe(aure)

introduceMe(aure_musician)

introduceMe(ricky)

#%%
mutable struct MyData
	x::Float64
	x2::Float64
	y::Float64
	z::Float64
	function MyData(x::Float64, y::Float64)
		x2=x^2
		z = sin(x2+y)
		new(x, x2, y, z)
	end
end

mutable struct MyData2{T<:Real}
	x::T
	x2::T
	y::T
	z::Float64
	function MyData2{T}(x::T, y::T) where {T<:Real}
		x2=x^2
		z = sin(x2+y)
		new(x, x2, y, z)
	end
end

#%%
MyData2{Float64}(2.0,3.0)
MyData2{Int}(2,3)

#%%
module TestModuleTypes

export Circle, computePerimeter, computeArea, printCircleEquation

mutable struct Circle{T<:Real}
    radius::T
    perimeter::Float64
    area::Float64

    function Circle{T}(radius::T) where {T<:Real}
# we initialize perimeter and area to -1.0, which is not a possible value
        new(radius, -1.0, -1.0)
    end
end

@doc raw"""
	computePerimeter(circle::Circle)

Compute the perimeter of `circle` and store the value.
"""
function computePerimeter(circle::Circle)
    circle.perimeter = 2 * π * circle.radius
    return circle.perimeter
end

@doc raw"""
	computeArea(circle::Circle)

Compute the area of `circle` and store the value.
"""
function computeArea(circle::Circle)
    circle.area = π * circle.radius^2
    return circle.area
end

@doc raw"""
	printCircleEquation(xc::Real, yc::Real, circle::Circle )

Print the equation of a cricle with center at (xc, yc) and radius given by circle.
"""
function printCircleEquation(xc::Real, yc::Real, circle::Circle)
    println("(x - $xc)^2 + (y - $yc)^2 = $(circle.radius^2)")
    return
end
end # end module
#%%
using .TestModuleTypes

circle1 = Circle{Float64}(5.0)

computePerimeter(circle1)
circle1.perimeter

computeArea(circle1)
circle1.area

printCircleEquation(2, 3, circle1)

max_stars_repo_name: JeffreySarnoff/techytok-examples

max_issues_repo_name: JeffreySarnoff/techytok-examples

max_forks_repo_name: JeffreySarnoff/techytok-examples


text:
function myFunction (a{T};
                     x::Float64=sin(1))::Bool
    return false
end


text:
"""
    @reduce() do (acc₁ [= init₁]; x₁), ..., (accₙ [= initₙ]; xₙ)
        ...
    end
    @reduce(acc₁ op₁= x₁, ..., accₙ opₙ= xₙ)
    @reduce(acc₁ = op₁(init₁, x₁), ..., accₙ = opₙ(initₙ, xₙ))

Declare how accumulators are updated in the sequential basecase and
how the resulting accumulators from two basecases are combined.

The arguments `accᵢ` and `xᵢ` must be symbols except for `xᵢ` of the
last two forms in which an expression can be used at `xᵢ`.

In the first form,

```julia
function ((acc₁, acc₂, ..., accₙ), (x₁, x₂, ..., xₙ))
    ...  # body of the `do` block
    return (acc₁, acc₂, ..., accₙ)
end
```

should be an associative function.

In the last two forms, every `opᵢ` should be an associative function.

If `initᵢ` is specified, the tuple `(init₁, init₂, ..., initₙ)` should
be the identify of the related associative function.  `accᵢ = initᵢ`
is evaluated for each basecase (each `Task`) in the beginning.

Consider a loop with the following form

```julia
@floop for ...
    # code computing (x₁, x₂, ..., xₙ)
    @reduce() do (acc₁ = init₁; x₁), ..., (accₙ = initₙ; xₙ)
        # code updating (acc₁, acc₂, ..., accₙ) using (x₁, x₂, ..., xₙ)
    end
end
```

This is converted to

```julia
acc₁ = init₁
...
accₙ = initₙ
for ...
    # code computing (x₁, x₂, ..., xₙ)
    # code updating (acc₁, acc₂, ..., accₙ) using (x₁, x₂, ..., xₙ)
end
```

for computing `(acc₁, acc₂, ..., accₙ)` of each basecase.  The
accumulators `accᵢ` of two basecases are combined using "code updating
`(acc₁, acc₂, ..., accₙ)` using `(x₁, x₂, ..., xₙ)`" where `(x₁, x₂,
..., xₙ)` are replaced with `(acc₁, acc₂, ..., accₙ)` of the next
basecase.  Note that "code computing `(x₁, x₂, ..., xₙ)`" is not used
for combining the basecases.

# Examples
```julia
@reduce() do (vmax=-Inf; v), (imax=0; i)
    if isless(vmax, v)
        vmax = v
        imax = i
    end
end

@reduce(s += y, p *= y)

@reduce(xs = append!!(EmptyVector(), x), ys = append!!(EmptyVector(), y))
```
"""
macro reduce(args...)
    # We insert a custom non-`Expr` object `ReduceOpSpec` in the AST
    # so that the argument to `@reduce` can be later analyzed by
    # `floop_parallel` below.  We do this because `floop_parallel` is
    # run after `macroexpand`.  We insert `throw` of `ReduceOpSpec` so
    # that unprocessed `ReduceOpSpec` becomes an error at runtime
    # rather than silently failing.
    :(throw($(ReduceOpSpec(collect(Any, args)))))
end
# TODO: detect free variables in `do` blocks

struct ReduceOpSpec
    args::Vector{Any}
    visible::Vector{Symbol}
end

ReduceOpSpec(args::Vector{Any}) = ReduceOpSpec(args, Symbol[])

"""
    @init begin
        pv₁ = init₁
        ...
        pvₙ = initₙ
    end

Initialize private variables `pvᵢ` with initializer expression `initᵢ` for
each task. This can be used for mutating objects in a data race-free manner.
"""
macro init(ex)
    :(throw($(initspec(ex))))
end

struct InitSpec
    expr::Expr
    lhs::Vector{Union{Symbol,Expr}}
    rhs::Vector{Any}
    visible::Vector{Symbol}
end

InitSpec(expr::Expr, lhs::Vector{Union{Symbol,Expr}}, rhs::Vector{Any}) =
    InitSpec(expr, lhs, rhs, Symbol[])

initspec(@nospecialize x) = invalid_at_init(x)
initspec(ex::Expr) = InitSpec(ex, collect_assignments(ex)...)

invalid_at_init(@nospecialize ex) =
    error("`@init` requires an assignment or a sequence of assignments; got:\n", ex)

# TODO: merge with `assigned_vars`?
collect_assignments(ex) = collect_assignments!(Union{Symbol,Expr}[], [], ex)
collect_assignments!(_, _, @nospecialize(x)) = invalid_at_init(x)
collect_assignments!(lhs, rhs, ::LineNumberNode) = (lhs, rhs)
function collect_assignments!(lhs, rhs, ex::Expr)
    @match ex begin
        Expr(:(=), l::Union{Symbol,Expr}, r) => begin
            push!(lhs, l)
            push!(rhs, r)
        end
        Expr(:block, args...) => begin
            for a in args
                collect_assignments!(lhs, rhs, a)
            end
        end
        # TODO: should we support other things, like side effects, here?
        _ => invalid_at_init(ex)
    end
    return (lhs, rhs)
end

""" Recursively copy `expr::Expr` but not `ReduceOpSpec` or `InitSpec`. """
copyexpr(expr::Expr) = Expr(expr.head, (copyexpr(a) for a in expr.args)...)
copyexpr(@nospecialize x) = x

function analyze_loop_local_variables!(spec::Union{ReduceOpSpec,InitSpec}, scopes)
    @assert isempty(spec.visible)
    append!(spec.visible, (var.name for sc in scopes for var in sc.bounds))
    unique!(spec.visible)
    nothing
end
analyze_loop_local_variables!(expr::Expr, scopes) =
    @match expr begin
        Expr(:scoped, sc, inner) => begin
            push!(scopes, sc)
            analyze_loop_local_variables!(inner, scopes)
            pop!(scopes)
            nothing
        end
        _ => begin
            for a in expr.args
                analyze_loop_local_variables!(a, scopes)
            end
        end
    end
analyze_loop_local_variables!(@nospecialize(_), _) = nothing

""" Fill `visible` field of `ReduceOpSpec` and `InitSpec` in-place. """
function fill_loop_local_variables!(expr)
    ex = copyexpr(expr)  # deep copy expr but not ReduceOpSpec or InitSpec
    analyze_loop_local_variables!(solve_from_local!(simplify_ex(ex)), Any[])
    return expr
end
#
# Currently, we overestimate variables visible to `@reduce` and `@init`. For
# example, we treat that `a` is visible to `init` clause of `@reduce` (i.e.,
# `f(a)`) in the following example:
#
#     @reduce y = op(f(a), x)
#     a = 1
#
# This is for (relatively easily) supporting the uses of `let` like this:
#
#     julia> let
#                let
#                    a = 1
#                end
#                a
#            end
#     ERROR: UndefVarError: a not defined
#     ...
#
#     julia> let
#                let
#                    a = 1
#                end
#                @show a
#                a = 0
#                a
#            end
#     a = 1
#     1
#
# The latter is the example where the visibility of `a` at `@show a` is
# "retroactively" changed by the assignment after `@show a`; i.e., we cannot
# use simple single forward-pass algorithm to compute the set of accessible
# variables at given points.

function unpack_kwargs(;
    otherwise = donothing,
    on_expr = otherwise,
    on_init = otherwise,
    kwargs...
)
    @assert isempty(kwargs)
    return (otherwise, on_expr, on_init)
end

function on_reduce_op_spec(on_spec, ex; kwargs...)
    (otherwise, on_expr, on_init) = unpack_kwargs(; kwargs...)
    @match ex begin
        Expr(:call, throw′, spec::ReduceOpSpec) => on_spec(spec)
        Expr(:call, throw′, spec::InitSpec) => on_init(spec)
        Expr(head, args...) => begin
            new_args = map(args) do x
                on_reduce_op_spec(on_spec, x; kwargs...)
            end
            on_expr(head, new_args...)
        end
        _ => otherwise(ex)
    end
end

on_reduce_op_spec_reconstructing(
    on_spec,
    ex;
    otherwise = identity,
    on_init = otherwise,
) = on_reduce_op_spec(
    on_spec,
    ex;
    on_expr = Expr,
    otherwise = otherwise,
    on_init = on_init,
)

is_parallel(ex) = on_reduce_op_spec(
    _ -> true,
    ex;
    on_init = _ -> true,
    on_expr = (_, args...) -> any(args),
    otherwise = _ -> false,
)

function floop_parallel(ctx::MacroContext, ex::Expr, simd, executor = nothing)
    if !isexpr(ex, :for, 2)
        error("expected a `for` loop; got:\n", ex)
    end
    iterspec, body = ex.args
    parallel_loop_ex = @match iterspec begin
        Expr(:block, loop_axes...) => begin
            rf_arg, coll = transform_multi_loop(loop_axes)
            as_parallel_loop(ctx, rf_arg, coll, body, simd, executor)
        end
        Expr(:(=), rf_arg, coll) => begin
            as_parallel_loop(ctx, rf_arg, coll, body, simd, executor)
        end
    end
    return parallel_loop_ex
end

function analyze_rf_args(ex::Expr)
    @assert isexpr(ex, :tuple)
    accs = []
    inits = []
    inputs = []
    for arg in ex.args
        @match arg begin
            Expr(:block, acc_init, x) || Expr(:block, acc_init, ::LineNumberNode, x) =>
                begin
                    push!(inputs, x)
                    @match acc_init begin
                        Expr(:(=), a, i) => begin
                            push!(accs, a)
                            push!(inits, i)
                        end
                        a => push!(accs, a)
                    end
                end
            Expr(:tuple, a, x) => begin
                throw(ArgumentError("got `($a, $x)` use `($a; $x)` instead"))
            end
        end
    end
    if !isempty(inits) && length(inits) != length(inputs)
        # TODO: Remove this restriction.  If not all `init`s are
        # specified, use it only for the empty case.
        throw(ArgumentError(string(
            "[NOT IMPLEMENTED]",
            " Currently, initial value should be specified for all accumulators",
            " when it is specified for at least one accumulator.",
        )))
    end
    if isempty(inits)
        inits = nothing
    end
    return accs, inits, inputs
end

function verify_unique_symbols(all_vars, kind)
    if length(Set(all_vars)) != length(all_vars)
        counts = Dict{Any,Int}()
        for var in all_vars
            counts[var] = get(counts, var, 0) + 1
        end
        dups = sort!([(n, var) for (var, n) in counts if n > 1])
        msg = sprint() do io
            print(io, "Same $kind variable used multiple times.")
            for (n, var) in dups
                println(io)
                print(io, "* `", var, "` used ", n, " times")
            end
        end
        throw(ArgumentError(msg))
    end
end

# To allow something like `@reduce(c += 1)` and `@reduce(c = 0 + 1)`,
# assign the right (second) argument to a temporary variable:
function extract_pre_updates(raw_inputs)
    inputs = []
    pre_updates = []
    for x in raw_inputs
        if x isa Symbol
            push!(inputs, x)
        else
            @gensym tmp
            push!(pre_updates, :($tmp = $x))
            push!(inputs, tmp)
        end
    end
    return (inputs, pre_updates)
end

function uniquify_inputs(inputs)
    uniquified = empty(inputs)
    pre_updates = Expr[]
    seen = Set{eltype(inputs)}()
    for x in inputs
        if x in seen
            y = gensym(x)
            push!(pre_updates, :($y = $x))
        else
            push!(seen, x)
            y = x
        end
        push!(uniquified, y)
    end
    return uniquified, pre_updates
end

function inject_spec_expr_for_analysis(ex::Expr)
    function on_init(spec::InitSpec)
        quote
            $(spec.expr)
            throw($spec)
        end
    end
    function on_reduce(spec::ReduceOpSpec)
        # If we want to make
        #     @reduce(a = op(a0, x))
        #     @reduce(b = op(f(a), x))
        # work, we need to do something here. But, for the moment, let us
        # pretend that the output of `@reduce` is not visible outside the
        # `@reduce` block.
        :(throw($spec))
    end
    return on_reduce_op_spec_reconstructing(on_reduce, ex; on_init = on_init)
end

function as_parallel_loop(ctx::MacroContext, rf_arg, coll, body0::Expr, simd, executor)
    body0 = deepcopy(body0)  # To be mutated by `fill_loop_local_variables!`.
    dummy_loop_body = quote  # Dummy expression that simulates the loop body for JuliaVariables.
        $rf_arg = nothing
        $(inject_spec_expr_for_analysis(body0))
    end
    fill_loop_local_variables!(dummy_loop_body)

    accs_symbols = Symbol[]
    inputs_symbols = Symbol[]
    init_exprs = []
    combine_bodies = []
    is_init = Bool[]  # `is_init[i]` is true iff `i`th accumulator is due to `@init`
    all_rf_inits = []
    all_rf_accs = []
    all_rf_inputs = []

    function check_invariance()
        num_state_groups = length(accs_symbols)
        @assert length(inputs_symbols) == num_state_groups
        @assert length(init_exprs) == num_state_groups
        @assert length(is_init) == num_state_groups
        @assert length(all_rf_inits) == num_state_groups
        @assert length(all_rf_accs) == num_state_groups
        @assert length(all_rf_inputs) == num_state_groups
        @assert length(combine_bodies) == num_state_groups

        nums_grouped_states = map(length, all_rf_accs)
        function check_num_states(xs)
            for (x, n) in zip(xs, nums_grouped_states)
                x === nothing && continue
                @assert length(x) == n
            end
        end
        check_num_states(all_rf_inits)
        check_num_states(all_rf_inputs)
    end

    function on_init(spec::InitSpec)
        @gensym grouped_private_states
        push!(accs_symbols, grouped_private_states)
        push!(inputs_symbols, :_)

        accs = spec.lhs
        push!(is_init, true)
        push!(all_rf_inits, nothing)
        push!(all_rf_accs, accs)
        push!(all_rf_inputs, nothing)
        verify_unique_symbols(accs, "private")

        # The corresponding combine function is "keep left" (i.e., do nothing):
        push!(combine_bodies, nothing)

        @gensym init_allocator
        scratch = quote
            $(spec.expr)  # the expression from `@init $initializer`; sets `accs`
            $Base.@inline $init_allocator() = tuple($(spec.rhs...))
            $ScratchSpace($init_allocator, ($(accs...),))
        end

        if isempty(intersect(spec.visible, unbound_rhs(spec.expr)))
            # Hoisting out `@init`, since it is not accessing variables used
            # inside the loop body.
            push!(init_exprs, scratch)
            return quote
                # Just in case it is demoted to `EmptyScratchSpace`:
                $grouped_private_states = $allocate($grouped_private_states)
                ($(accs...),) = $grouped_private_states.value
            end
        else
            push!(init_exprs, _FLoopInit())
            return quote
                if $grouped_private_states isa $_FLoopInit
                    # `$scratch` sets `$accs`.
                    # Assigning to `grouped_private_states` for reusing it next time.
                    $grouped_private_states = $scratch
                else
                    # Just in case it is demoted to `EmptyScratchSpace`:
                    $grouped_private_states = $allocate($grouped_private_states)
                    # After the initialization, just carry it over to the next iteration:
                    ($(accs...),) = $grouped_private_states.value
                end
            end
        end
    end

    body1 = on_reduce_op_spec_reconstructing(body0; on_init = on_init) do spec
        opspecs = spec.args
        @gensym grouped_accs grouped_inputs
        push!(accs_symbols, grouped_accs)
        push!(inputs_symbols, grouped_inputs)
        if length(opspecs) == 1 && is_function(opspecs[1])
            # handle: @reduce() do ...
            rf_ex, = opspecs
            # rf_ex = :(((acc1; input1), ..., (accN; inputN)) -> rf_body)
            accs, inits, inputs = analyze_rf_args(rf_ex.args[1])
            rf_body = rf_ex.args[2]
            pre_updates = []
            updaters = [rf_body]
        else
            if all(is_rebinding_update, opspecs)
                # handle: @reduce(acc₁ op₁= x₁, ..., accₙ opₙ= xₙ)
                ops = [Symbol(String(x.head)[1:end-1]) for x in opspecs]
                accs = [x.args[1] for x in opspecs]
                inits = nothing
                (inputs, pre_updates) = extract_pre_updates([x.args[2] for x in opspecs])
            elseif all(x -> isexpr(x, :(=), 2) && isexpr(x.args[2], :call, 3), opspecs)
                # handle: @reduce(acc₁ = op₁(init₁, x₁), ..., accₙ = opₙ(initₙ, xₙ))
                ops = [x.args[2].args[1] for x in opspecs]
                accs = [x.args[1] for x in opspecs]
                inits = [x.args[2].args[2] for x in opspecs]
                (inputs, pre_updates) =
                    extract_pre_updates([x.args[2].args[3] for x in opspecs])
            else
                error(join(vcat(["unsupported:"], opspecs), "\n"))
            end
            inputs, pre_updates2 = uniquify_inputs(inputs)
            append!(pre_updates, pre_updates2)
            updaters = [:($a = $op($a, $x)) for (op, a, x) in zip(ops, accs, inputs)]
        end
        push!(is_init, false)
        push!(all_rf_inits, inits)
        push!(all_rf_accs, accs)
        push!(all_rf_inputs, inputs)
        verify_unique_symbols(accs, "accumulator")
        verify_unique_symbols(inputs, "input")
        # TODO: input symbols just have to be unique within a
        # `@reduce` block.  This restriction (unique across all
        # `@reduce`) can be removed.
        use_oninit = false
        if inits === nothing
            initializers = [:($a = $x) for (a, x) in zip(accs, inputs)]
            updaters0 = []
            push!(init_exprs, _FLoopInit())
        elseif any(!isempty(intersect(spec.visible, unbound_rhs(ex))) for ex in inits)
            initializers = [:($a = $x) for (a, x) in zip(accs, inits)]
            updaters0 = updaters
            push!(init_exprs, _FLoopInit())
        else
            # Hoisting out `init` clauses, since it is not accessing variables
            # used inside the loop body.
            use_oninit = true
            initializers = updaters0 = nothing
            push!(init_exprs, :(tuple($(inits...))))
        end
        function rf_body_with_init(pre_updates = [])
            if use_oninit
                return quote
                    $(pre_updates...)
                    ($(accs...),) = $grouped_accs
                    $(updaters...)
                    $grouped_accs = ($(accs...),)
                end
            end
            quote
                $(pre_updates...)
                if $grouped_accs isa $_FLoopInit
                    $(initializers...)
                    $(updaters0...)
                else
                    ($(accs...),) = $grouped_accs
                    $(updaters...)
                end
                $grouped_accs = ($(accs...),)
            end
        end
        combine_body0 = quote
            ($(inputs...),) = $grouped_inputs
            $(rf_body_with_init())
        end
        combine_body = if use_oninit
            combine_body0
        else
            quote
                if $grouped_inputs isa $_FLoopInit
                else
                    $combine_body0
                end
            end
        end
        push!(combine_bodies, combine_body)
        return rf_body_with_init(pre_updates)
    end
    check_invariance()

    body2, info = transform_loop_body(body1, accs_symbols)

    @gensym oninit_function reducing_function combine_function result
    if ctx.module_ === Main
        # Ref: https://github.com/JuliaLang/julia/issues/39895
        oninit_function = Symbol(:__, oninit_function)
        reducing_function = Symbol(:__, reducing_function)
        combine_function = Symbol(:__, combine_function)
    end

    unpackers = map(
        enumerate(zip(is_init, all_rf_accs, all_rf_inits)),
    ) do (i, (nounpack, accs, inits))
        @gensym grouped_accs
        if nounpack
            # This accumulator is from `@init`.
            nothing
        elseif inits === nothing
            quote
                $grouped_accs = $result[$i]
                # Assign to accumulator only if it is updated at least once:
                if $grouped_accs isa $_FLoopInit
                else
                    ($(accs...),) = $grouped_accs
                end
            end
        else
            quote
                $grouped_accs = $result[$i]
                ($(accs...),) = if $grouped_accs isa $_FLoopInit
                    $unreachable_floop()
                    # ($(inits...),)
                else
                    $grouped_accs
                end
            end
        end
    end
    unpack_state = Expr(:block, unpackers...)
    gotos = gotos_for(info.external_labels, unpack_state, result)
    mkdecl(x) = x |> NotA(Nothing) |> Cat() |> Map(a -> :(local $a)) |> collect
    accs_declarations = mkdecl(all_rf_accs)
    inputs_declarations = mkdecl(all_rf_inputs)

    return quote
        $Base.@inline function $oninit_function()
            $(accs_declarations...)
            return tuple($(init_exprs...))
        end
        $Base.@inline function $reducing_function(($(accs_symbols...),), $rf_arg)
            $(accs_declarations...)
            $body2
            return ($(accs_symbols...),)
        end
        $combine_function(_, b::$(Union{Goto,Return})) = b
        function $combine_function(($(accs_symbols...),), ($(inputs_symbols...),))
            $(accs_declarations...)
            $(inputs_declarations...)
            $(combine_bodies...)
            return ($(accs_symbols...),)
        end
        $_verify_no_boxes($reducing_function)
        $result = $_fold(
            $wheninit(
                $oninit_function,
                $whencombine($combine_function, $reducing_function),
            ),
            $coll,
            $executor,
            $(Val(simd)),
        )
        $result isa $Return && return $result.value
        $(gotos...)
        $unpack_state
        nothing
    end
end

struct _FLoopInit end

@noinline unreachable_floop() = error("unrechable reached (FLoops.jl bug)")

@inline _fold(rf::RF, coll, ::Nothing, simd) where {RF} =
    _fold(rf, coll, PreferParallel(), simd)
@inline _fold(rf::RF, coll, exc::Executor, simd) where {RF} = unreduced(
    transduce(IdentityTransducer(), rf, DefaultInit, coll, maybe_set_simd(exc, simd)),
)

function Base.showerror(io::IO, opspecs::ReduceOpSpec)
    print(io, "`@reduce(")
    join(io, opspecs.args, ", ")
    print(io, ")` used outside `@floop`")
end

function Base.showerror(io::IO, spec::InitSpec)
    ex = spec.expr
    print(io, "`@init", ex, "` used outside `@floop`")
end


text:
mutable struct State{T,dim} <: AbstractState
  t::Float64
  step::Int
  dt::Float64
  sol::Array{T,dim}
end

State(T::DataType, sz::Tuple, dt) = State(0.0, 0, dt, zeros(T, sz))

mutable struct DualState{Tc,Tr,dimc,dimr} <: AbstractState
  t::Float64
  step::Int
  dt::Float64
  solc::Array{Tc,dimc}
  solr::Array{Tr,dimr}
end

DualState(Tc::DataType, Tr::DataType, sizec, sizer, dt) = DualState{Tc,Tr,length(sizec),length(sizer)}(
  0.0, 0, dt, zeros(Tc, sizec), zeros(Tr, sizer))

"""
This type defines the linear implicit and explicit components of an equation.
The linear implicit part of an equation is defined by an array of coefficients
which multiply the solution. The explicit part of an equation is calculated
by a function that may define linear and nonlinear parts.
"""
struct Equation{T,dim} <: AbstractEquation
  LC::Array{T,dim} # Coeffs of the eqn's implicit linear part
  calcN!::Function # Function that calcs linear & nonlinear parts
end

struct DualEquation{T,dimc,dimr} <: AbstractEquation
  LCc::Array{T,dimc}
  LCr::Array{T,dimr}
  calcN!::Function
end

"""
    Problem(g, v, p, eq, ts)

Initialize a FourierFlows problem on grid g, with variables v, parameters p,
equation eq, and timestepper ts.
"""
mutable struct Problem <: AbstractProblem
  grid::AbstractGrid
  vars::AbstractVars
  params::AbstractParams
  eqn::AbstractEquation
  ts::AbstractTimeStepper
  state::AbstractState
  t::Float64
  step::Int
end

Problem(g, v, p, eq, ts, st) = Problem(g, v, p, eq, ts, st, st.t, st.step)

function Problem(g, v, p, eq::Equation, ts; cxsol=true)
  Tsol = cxsol ? cxeltype(eq.LC) : eltype(eq.LC)
  Problem(g, v, p, eq, ts, State(Tsol, size(eq.LC), ts.dt))
end

function Problem(g, v, p, eq::DualEquation, ts; cxsolc=true, cxsolr=true)
  Tc = cxsolc ? cxeltype(eq.LCc) : eltype(eq.LCc)
  Tr = cxsolr ? cxeltype(eq.LCr) : eltype(eq.LCr)
  Problem(g, v, p, eq, ts, DualState(Tc, Tr, size(eq.LCc), size(eq.LCr), ts.c.dt))
end

unpack(prob::Problem) = prob.grid, prob.vars, prob.params, prob.state

#=
For Julia v1.0 release:
  The `t` and `step` properties can be removed from Problem, instead
  overloading the `getfield` function to intercept attempts to access
  `t` and `step`, which will be redirected to prob.state.t and
  prob.state.step (and perhaps sol as well). This'll make lots of things
  just a little bit nicer.

import Base: getfield

function getfield(prob::AbstractProblem, name)
  if name ∈ [:t, :step]
    return getfield(prob.state, name)
  else
    return getfield(prob, name)
  end
end
=#

#=
For Julia v1.0 release:
  The `t` and `step` properties can be removed from Problem, instead
  overloading the `getfield` function to intercept attempts to access
  `t` and `step`, which will be redirected to prob.state.t and
  prob.state.step (and perhaps sol as well). This'll make lots of things
  just a little bit nicer.

import Base: getfield

function getfield(prob::AbstractProblem, name)
  if name ∈ [:t, :step]
    return getfield(prob.state, name)
  else
    return getfield(prob, name)
  end
end
=#


text:
module Hmm

  type BaumWelch
    A::Array{Float64,2}
    B::Array{Float64,2}
    ρ::Array{Float64,2}
    state_num::Int64
    symbol_num::Int64
    symbol::Array{Int64,2}
    pA::Array{Float64,2}
    pB::Array{Float64,2}
    pρ::Array{Float64,2}
    A_numer::Array{Float64,2}
    A_denom::Array{Float64,2}
    B_numer::Array{Float64,2}
    B_denom::Array{Float64,2}
    newρ::Array{Float64,2}
    α::Array{Float64,2}
    β::Array{Float64,2}
    c::Array{Float64,1}
  end

  function hmm_initialization(A, B, ρ)

    state_num = length(A[1,:])
    symbol_num = length(B[2,:])
    symbol = [0 1]
    # pA = A
    # pB = B
    # pρ = ρ
    pA = zeros(state_num, state_num)
    pB = zeros(state_num, symbol_num)
    pρ = zeros(1,state_num)
    A_numer = zeros(state_num, state_num)
    A_denom = zeros(state_num, state_num)
    B_numer = zeros(state_num, symbol_num)
    B_denom = zeros(state_num, symbol_num)
    newρ = zeros(1,state_num)
    α = zeros(0, state_num)
    β = zeros(0, state_num)
    c = zeros(0)

    hmm = BaumWelch(A, B, ρ, state_num, symbol_num, symbol, pA, pB, pρ,
                    A_numer, A_denom, B_numer, B_denom, newρ, α, β, c)

    return hmm
  end

  function init_variables(hmm)
    hmm.A_numer = zeros(hmm.state_num, hmm.state_num)
    hmm.A_denom = zeros(hmm.state_num, hmm.state_num)
    hmm.B_numer = zeros(hmm.state_num, hmm.symbol_num)
    hmm.B_denom = zeros(hmm.state_num, hmm.symbol_num)
    hmm.newρ = zeros(1,hmm.state_num)
  end

  function forward(hmm,obs)
    # 観測系列の長さ
    n = length(obs)
    # scaled forwardアルゴリズム
    # 変数の初期化
    α = zeros(n, hmm.state_num)
    c = zeros(n)

    # 初期化
    α[1, :] = hmm.ρ[:] .* hmm.B[:, Int(obs[1])+1]
    c[1] = 1.0 / sum(α[1, :])
    α[1, :] = c[1] * α[1, :]

    # 再帰的計算
    for t in 2:n
      α[t, :] = (α[t-1, :]' * hmm.A)' .* hmm.B[:, Int(obs[t])+1]
      c[t] = 1.0 / sum(α[t, :])
      α[t, :] = c[t] * α[t, :]
    end

    hmm.α = α
    hmm.c = c
  end

  function backward(hmm, obs)
    # 観測系列の長さ
    n = length(obs)
    # scaled backwardアルゴリズム
    # 変数の初期化
    β = zeros(n, hmm.state_num)

    # 初期化
    β[n, :] = hmm.c[n]

    # 再帰的計算
    for t in n:-1:2
      β[t-1, :] = hmm.A * (hmm.B[:, Int(obs[t]+1)] .* β[t, :])
      β[t-1, :] = hmm.c[t-1] * β[t-1, :]
    end

    hmm.β = β
  end

  function maximization_step(hmm, obs)
    # 観測系列の長さ
    n = length(obs)

    # update A
    for i in 1:hmm.state_num
      for j in 1:hmm.state_num
        A_numer = A_denom = 0.0
        for t in 1:n-1
          A_numer += hmm.α[t,i] * hmm.A[i,j] * hmm.B[j, Int(obs[t+1])+1] * hmm.β[t+1,j]
          A_denom += hmm.α[t,i] * hmm.β[t,i] / hmm.c[t]
        end
        hmm.A_numer[i, j] += A_numer
        hmm.A_denom[i, j] += A_denom
      end
    end

    # update B
    for j in 1:hmm.state_num
      for k in 1:hmm.symbol_num
        B_numer = B_denom = 0.0
        for t in 1:n
          B_numer += (obs[t] == hmm.symbol[k]) * hmm.α[t, j] * hmm.β[t, j] / hmm.c[t]
          B_denom +=  hmm.α[t, j] * hmm.β[t, j] / hmm.c[t]
        end
        hmm.B_numer[j, k] += B_numer
        hmm.B_denom[j, k] += B_denom
      end
    end

    # update ρ
    hmm.newρ += hmm.α[1, :]' .* hmm.β[1, :]' / hmm.c[1]

  end

  function check_convergence(hmm, eps)
    diff = 0.0
    diff += sum((hmm.A - hmm.pA).^2)
    diff += sum((hmm.B - hmm.pB).^2)
    diff += sum((hmm.ρ - hmm.pρ).^2)

    return sqrt(diff) < eps
  end

  function train(hmm, obs, eps = 1e-9, max_iter = 10000)

    # init
    seq_num = length(obs)
    loglik = zeros(seq_num)
    iter = 0

    for count in 1:max_iter
      if count % 10 == 0
        println("iter: [", count, "]", )
      end

      init_variables(hmm)

      # calc α, β, c each sequences
      for s in 1:seq_num
        s_obs = obs[s]
        # E-Step
        forward(hmm, s_obs)
        backward(hmm, s_obs)
        # M-Step
        maximization_step(hmm, s_obs)
        # calc log-likelihood
        loglik[s] = -sum(log.(hmm.c))
      end

      # update parameter
      hmm.A = hmm.A_numer ./ hmm.A_denom
      hmm.B = hmm.B_numer ./ hmm.B_denom
      hmm.ρ = hmm.newρ ./ seq_num

      # convergence check
      if check_convergence(hmm, eps)
        println("convergence !!")
        break
      end

      # save previous parameter
      hmm.pA = hmm.A
      hmm.pB = hmm.B
      hmm.pρ = hmm.ρ

      iter = count
    end

    print("iteration: ", iter)
    println(" log-likelihood: ", mean(loglik))

  end
end


text:
export topoSort

function topoSort(g)
  graph=copy(g)
  # Initialization
  indegree = sum(graph,1); # Incoming link counter
  indegree0 = indegree; # Initial incoming link counter

  # Topological sort
  G = [1:size(graph,1);]   # Set of nodes to sort
  L = [];     # Empty list for the sorted elem.
  S = [];     # Empty Set for nodes with no incoming edges
  level = ones(size(G));  # Level of the node in L
  dag = graph;    # Directed Acyclic Graph

  while !isempty(G) # G is non-empty
      # Search for free nodes
      for n = G
          if indegree[n]==0 # Incoming link counter of n = 0
              S = [S;n];
          end
      end

      if !isempty(S) # S is non-empty, perform ordinary top-sort
        n = S[1];  # Remove a node n from S and G
        S = S[S.!=n];
        G = G[G.!=n];
        L = [L;n]; # Insert n into L
        for m = find(graph[n,:]) # Foreach node m with an edge e from n to m
            graph[n,m] = 0; # Remove e
            level[m] = max(level[m], level[n]+1);
            indegree[m] = indegree[m] -1; # Decrement incoming link counter from m
        end
      end
    end
    return L
end


text:
"""
    MatShell{T}(obj, m, n)

Create a `m×n` PETSc shell matrix object wrapping `obj`.

If `obj` is a `Function`, then the multiply action `obj(y,x)`; otherwise it calls `mul!(y, obj, x)`.
This can be changed by `PETSc._mul!`.

"""
mutable struct MatShell{T,A} <: AbstractMat{T}
  ptr::CMat
  comm::MPI.Comm
  obj::A
end

struct MatOp{T,Op}
end

function _mul!(y,mat::MatShell{T,F},x) where {T, F<:Function}
  mat.obj(y, x)
end

function _mul!(y,mat::MatShell{T},x) where {T}
  LinearAlgebra.mul!(y, mat.obj, x)
end

MatShell{T}(obj, m, n) where {T} = MatShell{T}(obj, MPI.COMM_SELF, m, n, m, n)

@for_libpetsc begin
  function MatShell{$PetscScalar}(obj::A, comm::MPI.Comm, m, n, M, N) where {A}
    mat = MatShell{$PetscScalar,A}(C_NULL, comm, obj)
    # we use the MatShell object itsel
    ctx = pointer_from_objref(mat)
    @chk ccall((:MatCreateShell, $libpetsc), PetscErrorCode,
      (MPI.MPI_Comm,$PetscInt,$PetscInt,$PetscInt,$PetscInt,Ptr{Cvoid},Ptr{CMat}),
      comm, m, n, M, N, ctx, mat)

    mulptr = @cfunction(MatOp{$PetscScalar, MATOP_MULT}(), $PetscInt, (CMat, CVec, CVec))
    @chk ccall((:MatShellSetOperation, $libpetsc), PetscErrorCode, (CMat, MatOperation, Ptr{Cvoid}), mat, MATOP_MULT, mulptr)
    return mat
  end

  function (::MatOp{$PetscScalar, MATOP_MULT})(M::CMat,cx::CVec,cy::CVec)::$PetscInt
    #try
      r_ctx = Ref{Ptr{Cvoid}}()
      @chk ccall((:MatShellGetContext, $libpetsc), PetscErrorCode, (CMat, Ptr{Ptr{Cvoid}}), M, r_ctx)
      ptr = r_ctx[]
      mat = unsafe_pointer_to_objref(ptr)

      x = unsafe_localarray($PetscScalar, cx)
      y = unsafe_localarray($PetscScalar, cy)

      _mul!(y,mat,x)

      Base.finalize(y)
      Base.finalize(x)
      return $PetscInt(0)
    # catch e
    #   if e isa PetscErrorCode
    #     return e.code
    #   else
    #     return
    #   end
    # end
  end

end


text:
using MoistAir, Unitful,ThermoState
using Test
const λ = VariableSpec()

include("test_hyland83.jl")
include("test_hyland83a.jl")
include("test_user.jl")


text:
"""
Original author: Joe Ninan
Converted to Julia and optimized by Christian Gilbertson & Eric Ford
Returns a cubit interpolator for windowed sinc Filter curve.
no_of_points: number of intepolation points to use in cubic inteprolator
"""
function create_filter_curve_orig(no_of_points::Int; filter_size::Int=23, kaiserB::Int=13)
	x = range(-Int(floor(filter_size / 2)), stop=Int(floor(filter_size / 2)), length=no_of_points)
	Sinc = sinc.(x)
	Window = kaiser(length(x), kaiserB / π)
	FilterResponse = Window .* Sinc
	x_step = (x[end-1] - x[2]) / (no_of_points - 1)
	# append 0 to both ends far at the next node for preventing cubic spline
	# from extrapolating spurious values
	x = range(x[1] - x_step, stop=x[end] + x_step, length=no_of_points+2)
	#return CubicSplineInterpolation(x, multiple_append!([0.], FilterResponse, [0.]); extrapolation_bc=Line())
	return CubicSplineInterpolation(x, vcat(0., FilterResponse, 0.); extrapolation_bc=Line())
end

function create_filter_curve(no_of_points::Int; filter_size::Int=23, kaiserB::Int=13)
	ret_t = Float64
	x = range(-floor(Int,filter_size // 2), stop=floor(Int,filter_size // 2), length=no_of_points)
	# append 0 to both ends far at the next node for preventing cubic spline
	# from extrapolating spurious values
	Sinc = Array{ret_t,1}(undef,no_of_points+2)
	Window = Array{ret_t,1}(undef,no_of_points+2)
	Sinc[1] = zero(ret_t)
	Sinc[2:end-1] .= sinc.(x)
	Sinc[end] = zero(ret_t)
	Window[1] = zero(ret_t)
	Window[2:end-1] .= kaiser(length(x), kaiserB / π)
	Window[end] = zero(ret_t)
	FilterResponse = Window .* Sinc
	x_step = (x[end-1] - x[2]) / (no_of_points - 1)
	x = range(x[1] - x_step, stop=x[end] + x_step, length=no_of_points+2)
	return CubicSplineInterpolation(x, FilterResponse; extrapolation_bc=Line())
end

function next_min_arg(minarg::Integer, minvalue::Real, len::Integer)
	local sign = minvalue < 0
	local next_min_arg = minarg +sign -~sign
	if next_min_arg < 1       next_min_arg = len-next_min_arg end
	if next_min_arg > len     next_min_arg = next_min_arg-len end
	return next_min_arg
end

"""
Original author: Joe Ninan
Converted to Julia and optimized by Christian Gilbertson
Further adapted/optimized by Eric Ford
Additional optimizations possible by preallocating arrasy for minargs, Nminargs, minvalues, FilterValues and OldYCoords.
"""
function spectra_interpolate(
	newX::V1,
	oldX::V2,
	oldY::V3;
	PeriodicBoundary::Bool=false,
	filter_size::Int=23,
	kaiserB::Int=13,
	Filter::V4 = create_filter_curve(filter_size*21; filter_size=filter_size, kaiserB=kaiserB),
	 assume_sorted::Bool=false, verbose::Bool = false) where {
	 	T1<:Real, V1<:AbstractVector{T1}, T2<:Real, V2<:AbstractVector{T2}, T3<:Real, V3<:AbstractVector{T3}, T4<:Real, V4<:AbstractVector{T4}  }
	if !assume_sorted
		@assert issorted(newX)
		@assert issorted(oldX)
	end

	#Filter = create_filter_curve(filter_size*21; filter_size=filter_size, kaiserB=kaiserB)
	#pixarray = collect(-floor(Int,filter_size//2):floor(Int,filter_size//2))
	pixarray = -floor(Int,filter_size//2):floor(Int,filter_size//2)

	# oXsize = length(oldX)
    # First generate a 2D array of difference in pixel values
    # OldXminusNewX = np.array(oldX)[:,np.newaxis] - np.array(newX)
    # Find the minimum position to find nearest pixel for each each newX
	minargs = RvSpectMLBase.searchsortednearest(oldX, newX)
    # Pickout the those minumum values from 2D array
	minvalues = oldX[minargs] - newX
	Nminargs = next_min_arg.(minargs,minvalues,length(oldX))
	# EBF: used above to replace next few due to 1-based arrays, simplification and reduced temporaries
	#=
	sign = minvalues .< 0  # True means new X is infront of nearest old X
	# coordinate of the next adjacent bracketing point
	Nminargs = minargs .+sign -.~sign
	#Nminargs = Nminargs .% oXsize  # Periodic boundary  # WARN: Not right for 0-based arrays
	# In terms of pixel coordinates the shift values will be
	=#
	# EBF: consolidated to help reduce temporaries.  I'm not sure if if mattered.
	#=
	shiftvalues = minvalues./abs.(oldX[minargs]-oldX[Nminargs])
	# Coordinates to calculate the Filter values
	FilterCoords = shiftvalues .+ pixarray'
	FilterValues = Filter.(FilterCoords)
	=#
	FilterValues = Filter.(minvalues./abs.(oldX[minargs]-oldX[Nminargs]) .+ pixarray')
	# Coordinates to pick the values to be multiplied with Filter and summed
	OldYCoords = minargs .+ pixarray'
	if PeriodicBoundary
		# OldYCoords = OldYCoords .% oXsize  # Periodic boundary  # EBF: changed for 1-based arrays, TODO: WARN: need to test line below
		for i in 1:length(OldYCoords)
			if OldYCoords[i] < 1               		 OldYCoords[i] = length(oldX)-OldYCoords[i]   end
			if OldYCoords[i] > length(OldYCoords)+1  OldYCoords[i] = OldYCoords[i]-length(oldX)   end
		end

	else  # Extrapolate the last value till end..
		OldYCoords[OldYCoords .>= length(oldX) + 1] .= length(oldX)
		OldYCoords[OldYCoords .< 1] .= 1
	end
	# EBF consolidated to reduce temporaries.
	#OldYSlices = oldY[OldYCoords]
	#return vec(sum(OldYSlices .* FilterValues, dims=2))
	return vec(sum(oldY[OldYCoords] .* FilterValues, dims=2))
end


text:
using ADCME

PWD = pwd()
change_directory("$(@__DIR__)/CustomOps/build")
require_file("build.ninja") do
    ADCME.cmake()
end
require_library("ADSeismic") do
    ADCME.make()
end
cd(PWD)


text:
abstract type Pet end

struct Dog <: Pet
    name::AbstractString
end

struct Cat <: Pet
    name::AbstractString
end

"""
    encounter(a, b)

Simulate an encounter between `a` and `b`.
"""
encounter(a, b) = "$(name(a)) meets $(name(b)) and $(meets(a, b))."

"""
    name(p::Pet)

Return the name of the pet `p`.
"""
name(p::Pet) = p.name

# fallbacks
meets(a::Pet, b::Pet) = "is cautious"
meets(a::Pet, b) = "runs away"

"""
    meets(a, b)

Return a string describing what happens when `a` and `b` meet.
"""
meets(a, b) = "nothing happens"

# specific types
meets(a::Dog, b::Dog) = "sniffs"
meets(a::Dog, b::Cat) = "chases"
meets(a::Cat, b::Dog) = "hisses"
meets(a::Cat, b::Cat) = "slinks"