text:
ff_promote_rule(::Type{T1}, ::Type{T2}, ::Val{W}) where {T1,T2,W} = promote_type(T1,T2)
Base.@pure function _ff_promote_rule(::Type{T1}, ::Type{T2}, ::Val{W}) where {T1, T2, W}
    T_canon = promote_type(T1,T2)
    pick_vector_width(T_canon) < W ? T1 : T_canon
end
ff_promote_rule(::Type{T1}, ::Type{T2}, ::Val{W}) where {T1 <: Integer, T2 <: Integer,W} = _ff_promote_rule(T1,T2,Val{W}())
ff_promote_rule(::Type{T1}, ::Type{T2}, ::Val{W}) where {T1 <: FloatingTypes, T2<:FloatingTypes,W} = _ff_promote_rule(T1,T2,Val{W}())
@generated function ff_promote_rule(::Type{T1}, ::Type{T2}, ::Val{W}) where {T1 <: IntegerTypes, T2 <: FloatingTypes,W}
    T_canon = promote_type(T1,T2)
    (sizeof(T_canon) * W ≤ REGISTER_SIZE) && return T_canon
    @assert sizeof(T1) * W ≤ REGISTER_SIZE
    @assert sizeof(T1) == 4
    Float32
    # N, r = 4W ÷ REGISTER_SIZE
    # @assert iszero(r)
end

Base.promote_rule(::Type{V}, ::Type{T2}) where {W,T1,T2<:NativeTypes,V<:AbstractSIMDVector{W,T1}} = Vec{W,ff_promote_rule(T1,T2,Val{W}())}
@generated function Base.promote_rule(::Type{V1}, ::Type{V2}) where {W,T1,T2,V1<:AbstractSIMDVector{W,T1},V2<:AbstractSIMDVector{W,T2}}
    T = promote_type(T1,T2)
    if pick_vector_width(T) ≥ W
        return :(Vec{$W,$T})
    end
    if T === Float64 || T === Float32
        N, r1 = (sizeof(T) * W) ÷ REGISTER_SIZE
        Wnew, r2 = divrem(W, N)
        @assert iszero(r)

        return :(VecUnroll{$(N-1),$Wnew,$T,Vec{$Wnew,$T}})
        # Should we demote `Float64` -> `Float32`?
        # return :(Vec{$W,$T})
        # don't demote to smaller than `Float32`
        # return :(Vec{$W,Float32})
    end
    # They're both of integers
    V1MM = V1 <: MM
    V2MM = V2 <: MM
    if V1MM ⊻ V2MM
        V1MM ? :(Vec{$W,$T2}) : :(Vec{$W,$T1})
    else
        I = pick_integer(W, sizeof(Int), 4 - 3V1MM) #
        T <: Unsigned ? :(Vec{$W,$(unsigned(I))}) : :(Vec{$W,$I})
    end
end

@generated function Base.promote_rule(
    ::Type{<:VecUnroll{Nm1,Wsplit,T,V1}}, ::Type{V2}
) where {Nm1,Wsplit,T,T2,V1,W,V2<:AbstractSIMDVector{W,T2}}
    N = Nm1 + 1
    @assert N * Wsplit == W
    V3 = if V2 <: Mask
        Mask{Wsplit,mask_type(Wsplit)}
    else
        Vec{Wsplit,T2}
    end
    _assemble_vec_unroll(Val{Nm1}(), promote_type(V1,V3))
end

_assemble_vec_unroll(::Val{N}, ::Type{V}) where {N,W,T,V<:AbstractSIMDVector{W,T}} = VecUnroll{N,W,T,V}
Base.promote_rule(::Type{VecUnroll{N,W,T1,V}}, ::Type{T2}) where {N,W,T1,V,T2<:NativeTypes} = _assemble_vec_unroll(Val{N}(), promote_type(V,T2))
Base.promote_rule(::Type{<:VecUnroll{N,W,T,V1}}, ::Type{V2}) where {N,W,T,V1,V2<:AbstractSIMDVector{W}} = _assemble_vec_unroll(Val{N}(), promote_type(V1,V2))
Base.promote_rule(::Type{VecUnroll{N,W,T1,V1}}, ::Type{VecUnroll{N,W,T2,V2}}) where {N,W,T1,T2,V1,V2} = _assemble_vec_unroll(Val{N}(), promote_type(V1,V2))

Base.promote_rule(::Type{Bit}, ::Type{T}) where {T <: Number} = T

issigned(x) = issigned(typeof(x))
issigned(::Type{<:Signed}) = Val{true}()
issigned(::Type{<:Unsigned}) = Val{false}()
issigned(::Type{<:AbstractSIMD{<:Any,T}}) where {T} = issigned(T)
issigned(::Type{T}) where {T} = nothing
"""
 Promote, favoring <:Signed or <:Unsigned of first arg.
"""
@inline promote_div(x::Union{Integer,AbstractSIMD{<:Any,<:Integer}}, y::Union{Integer,AbstractSIMD{<:Any,<:Integer}}) = promote_div(x, y, issigned(x))
@inline promote_div(x, y) = promote(x, y)
@inline promote_div(x, y, ::Nothing) = promote(x, y) # for Integers that are neither Signed or Unsigned, e.g. Bool
@inline function promote_div(x::T1, y::T2, ::Val{true}) where {T1,T2}
    T = promote_type(T1, T2)
    signed(x % T), signed(y % T)
end
@inline function promote_div(x::T1, y::T2, ::Val{false}) where {T1,T2}
    T = promote_type(T1, T2)
    unsigned(x % T), unsigned(y % T)
end
itosize(i::Union{I,AbstractSIMD{<:Any,I}}, ::Type{J}) where {I,J} = signorunsign(i % J, issigned(I))
signorunsign(i, ::Val{true}) = signed(i)
signorunsign(i, ::Val{false}) = unsigned(i)

# Base.promote_rule(::Type{VecTile{M,N,W,T1}}, ::Type{T2}) where {M,N,W,T1,T2<:NativeTypes} = VecTile{M,N,W,promote_rule(T1,T2)}
# Base.promote_rule(::Type{VecTile{M,N,W,T1}}, ::Type{Vec{W,T2}}) where {M,N,W,T1,T2} = VecTile{M,N,W,promote_rule(T1,T2)}
# Base.promote_rule(::Type{VecTile{M,N,W,T1}}, ::Type{VecUnroll{M,W,T2}}) where {M,N,W,T1,T2} = VecTile{M,N,W,promote_rule(T1,T2)}
# Base.promote_rule(::Type{VecTile{M,N,W,T1}}, ::Type{VecTile{M,N,W,T2}}) where {M,N,W,T1,T2} = VecTile{M,N,W,promote_rule(T1,T2)}


text:
# Generate initial conditions
# Alex Tam, 12/10/2020

"Initial condition for actin filaments"
function actin_ic(s::State, parN, parA, Lxx, Lyy)
    af = Vector{Actin_Filament}(undef, parN.nA)
    for i = 1:length(af)
        af[i] = Actin_Filament(s, i, parA, Lxx, Lyy); # Generate actin filament
    end
    return af, s
end

"Initial condition for myosin motors"
function myosin_ic(s::State, mm, parN, xl, Lxx, Lxy, Lyx, Lyy)
    xl_indices = sample(1:length(xl), parN.nM, replace = false); # Sample intersections without replacement
    for i = 1:parN.nM
        push!(mm, Myosin_Motor(s, xl[xl_indices[i]], i, Lxx, Lxy, Lyx, Lyy)); # Generate myosin motor
    end
    deleteat!(xl, sort(xl_indices)); # Remove cross-links at intersections with a motor
    return mm, xl, s
end

ext: jl


text:
@testset "Stack Operators" begin
    @testset "encode/decode num" begin
        tests = [([0x75, 0x39, 0xd7, 0xf9], -2044148085),
        ([0x1b, 0x44, 0x8d, 0x17, 0x61, 0x36, 0x45, 0x51, 0x44, 0x38, 0x15, 0x20, 0xeb, 0x3f, 0xd1, 0xbc], -80840166210881125725200074314630382619),
        ([0xd3, 0x2c, 0xd5, 0xb9, 0x8b, 0xfa, 0x37, 0xd5, 0x24, 0xd0, 0xd2, 0xe8, 0x81, 0x8c, 0xcb, 0x15, 0xe5, 0xfe, 0xdd, 0x95, 0x55, 0xe4, 0x3d, 0x01, 0x82, 0x67, 0xc2, 0x23, 0xc5, 0xd6, 0x0c, 0x90, 0xe9, 0x72, 0xff, 0x9a, 0xd1, 0xcf, 0x16, 0xfa, 0x83], -8494629442587523471751141882831378743107401446683360409874466988135963066279425110931949604908243),
        ([0x7e, 0x5e, 0x14, 0x7b, 0x73, 0x6b, 0x57, 0x64, 0x67, 0x6d, 0x0b, 0x70, 0xa1, 0x6e, 0x71], 588973400481137931526979715425656446),
        ([0x7b], 123)]
        for i in 1:length(tests)
            @test BitcoinPrimitives.decode_num(tests[i][1]) == tests[i][2]
            @test BitcoinPrimitives.encode_num(tests[i][2]) == tests[i][1]
        end
    end
    @testset "Push value onto stack" begin
        test_stack = [[0x75, 0x39, 0xd7, 0xf9],
        [0x86, 0x8f, 0xf9, 0x73, 0x0a, 0x3c, 0x0c, 0xb9],
        [0xe2, 0xd6],
        [0x7b],
        [0x92, 0x8b, 0xae, 0xed, 0x1e, 0xfb, 0x04, 0x58],
        [0xdb, 0xaf, 0x44],
        [0x84, 0xd7, 0x35, 0xb7],
        [0xaf, 0x80, 0x91]]
        empty_stack = UInt8[]
        @testset "op_0" begin
            stack = test_stack
            want = stack
            push!(want, [0x00])
            @test BitcoinPrimitives.op_0(stack) == true
            @test stack == want
        end
        @testset "op_1negate" begin
            stack = test_stack
            want = stack
            push!(want, [0x81])
            @test BitcoinPrimitives.op_1negate(stack) == true
            @test stack == want
        end
        @testset "op_2" begin
            stack = test_stack
            want = stack
            push!(want, [0x02])
            @test BitcoinPrimitives.op_2(stack) == true
            @test stack == want
        end
        @testset "op_3" begin
            stack = test_stack
            want = stack
            push!(want, [0x03])
            @test BitcoinPrimitives.op_3(stack) == true
            @test stack == want
        end
        @testset "op_4" begin
            stack = test_stack
            want = stack
            push!(want, [0x01])
            @test BitcoinPrimitives.op_4(stack) == true
            @test stack == want
        end
        @testset "op_5" begin
            stack = test_stack
            want = stack
            push!(want, [0x05])
            @test BitcoinPrimitives.op_5(stack) == true
            @test stack == want
        end
        @testset "op_6" begin
            stack = test_stack
            want = stack
            push!(want, [0x06])
            @test BitcoinPrimitives.op_6(stack) == true
            @test stack == want
        end
        @testset "op_7" begin
            stack = test_stack
            want = stack
            push!(want, [0x07])
            @test BitcoinPrimitives.op_7(stack) == true
            @test stack == want
        end
        @testset "op_8" begin
            stack = test_stack
            want = stack
            push!(want, [0x08])
            @test BitcoinPrimitives.op_8(stack) == true
            @test stack == want
        end
        @testset "op_9" begin
            stack = test_stack
            want = stack
            push!(want, [0x09])
            @test BitcoinPrimitives.op_9(stack) == true
            @test stack == want
        end
        @testset "op_10" begin
            stack = test_stack
            want = stack
            push!(want, [0x0a])
            @test BitcoinPrimitives.op_10(stack) == true
            @test stack == want
        end
        @testset "op_11" begin
            stack = test_stack
            want = stack
            push!(want, [0x0b])
            @test BitcoinPrimitives.op_11(stack) == true
            @test stack == want
        end
        @testset "op_12" begin
            stack = test_stack
            want = stack
            push!(want, [0x0c])
            @test BitcoinPrimitives.op_12(stack) == true
            @test stack == want
        end
        @testset "op_13" begin
            stack = test_stack
            want = stack
            push!(want, [0x0d])
            @test BitcoinPrimitives.op_13(stack) == true
            @test stack == want
        end
        @testset "op_14" begin
            stack = test_stack
            want = stack
            push!(want, [0x0e])
            @test BitcoinPrimitives.op_14(stack) == true
            @test stack == want
        end
        @testset "op_15" begin
            stack = test_stack
            want = stack
            push!(want, [0x0f])
            @test BitcoinPrimitives.op_15(stack) == true
            @test stack == want
        end
        @testset "op_16" begin
            stack = test_stack
            want = stack
            push!(want, [0x10])
            @test BitcoinPrimitives.op_16(stack) == true
            @test stack == want
        end
    end
    @testset "Conditional control flow" begin
        @testset "op_nop" begin
            @test BitcoinPrimitives.op_nop([UInt8[]]) == true
            @test BitcoinPrimitives.op_nop([[0x64],
                [0x86, 0x8f, 0xf9, 0x73, 0x0a, 0x3c, 0x0c, 0xb9],
                [0xe2, 0xd6]]) == true
        end
        # TODO review those tests
        @testset "op_if" begin
            want = [[0x64],
                [0x86, 0x8f, 0xf9, 0x73, 0x0a, 0x3c, 0x0c, 0xb9],
                [0xe2, 0xd6]]
            stack = Vector{UInt8}[]
            items = [0x63,0x0c,0xe2]
            @test_throws AssertionError BitcoinPrimitives.op_if(stack,items)
            @test stack == []
            @test items == [0x63,0x0c,0xe2]

            stack = want
            @test BitcoinPrimitives.op_if(stack,items) == false
            @test stack == want
            @test items == UInt8[]

            items == [0x67,0x0d,0x69]
            @test BitcoinPrimitives.op_if(stack,items) == false
            @test stack == want
            @test items == UInt8[]

            items == [0x63,0x68,0x6b]
            @test BitcoinPrimitives.op_if(stack,items) == false
            @test stack == want
            @test items == UInt8[]

            items == [0x0d,0x0f,0x7a]
            @test BitcoinPrimitives.op_if(stack,items) == false
            @test stack == want
            @test items == UInt8[]

            push!(stack,UInt8[])
            items = [0x68,0xc4,0x07]
            @test BitcoinPrimitives.op_if(stack,items) == true
            @test stack == want
            @test items == [0xc4,0x07]

            items = [0x68,0xc4,0x07]
            @test BitcoinPrimitives.op_if(stack,items) == true
            pop!(want)
            @test stack == want
            @test items == [0xc4,0x07]
        end
        # TODO create those tests
        @testset "op_notif" begin
        end
        @testset "op_verify" begin
            stack = Vector{UInt8}[]
            @test_throws AssertionError BitcoinPrimitives.op_verify(stack)
            @test stack ==  Vector{UInt8}[]
            stack = [[0x64], [0x86, 0x8f, 0xf9], [0xe2, 0xd6]]
            @test BitcoinPrimitives.op_verify(stack) == true
            @test stack ==  [[0x64], [0x86, 0x8f, 0xf9]]
            stack = [[0x64], [0x86, 0x8f, 0xf9], [0xe2, 0xd6], UInt8[]]
            @test BitcoinPrimitives.op_verify(stack) == false
            @test stack ==  [[0x64], [0x86, 0x8f, 0xf9], [0xe2, 0xd6]]
        end
        @testset "op_return" begin
            @test BitcoinPrimitives.op_return([UInt8[]]) == false
            @test BitcoinPrimitives.op_return([[0x64],
                [0x86, 0x8f, 0xf9, 0x73, 0x0a, 0x3c, 0x0c, 0xb9],
                [0xe2, 0xd6]]) == false
        end
    end
    @testset "Timelock operations" begin
        empty = Vector{UInt8}[]
        stack = Vector{UInt8}[[0x03], [0x04]]
        @testset "op_checklocktimeverify" begin
            @test BitcoinPrimitives.op_checklocktimeverify(stack, 1, 0xffffffff) == false
            @test BitcoinPrimitives.op_checklocktimeverify(empty, 100, 1) == false
            push!(stack, [0x81])
            @test BitcoinPrimitives.op_checklocktimeverify(stack, 100, 1) == false
            pop!(stack)
            @test BitcoinPrimitives.op_checklocktimeverify(stack, 500000001, 1) == false
            @test BitcoinPrimitives.op_checklocktimeverify(stack, 1, 1) == false
            @test BitcoinPrimitives.op_checklocktimeverify(stack, 100, 1) == true
        end
        # TODO
        @testset "op_checksequenceverify" begin
        end
    end
    @testset "Stack operations" begin
        @testset "op_toaltstack" begin
            stack = Vector{UInt8}[]
            altstack = stack
            @test_throws AssertionError BitcoinPrimitives.op_toaltstack(stack, altstack)
            stack = [[0x04], [0x06, 0x0f, 0x09], [0x02, 0x06]]
            altstack = [[0xf4], [0xf6, 0xff, 0xf9], [0xf2, 0xf6]]
            @test BitcoinPrimitives.op_toaltstack(stack, altstack) == true
            @test stack == [[0x04], [0x06, 0x0f, 0x09]] && altstack == [[0xf4], [0xf6, 0xff, 0xf9], [0xf2, 0xf6], [0x02, 0x06]]
        end
        @testset "op_fromaltstack" begin
            stack = Vector{UInt8}[]
            altstack = stack
            @test_throws AssertionError BitcoinPrimitives.op_fromaltstack(stack, altstack)
            stack = [[0x04], [0x06, 0x0f, 0x09], [0x02, 0x06]]
            altstack = [[0xf4], [0xf6, 0xff, 0xf9], [0xf2, 0xf6]]
            @test BitcoinPrimitives.op_fromaltstack(stack, altstack) == true
            @test stack == [[0x04], [0x06, 0x0f, 0x09], [0x02, 0x06],[0xf2, 0xf6]] && altstack == [[0xf4], [0xf6, 0xff, 0xf9]]
        end
        @testset "op_2drop" begin
            stack = [[0x04]]
            @test_throws AssertionError BitcoinPrimitives.op_2drop(stack)
            stack = [[0x04], [0x06, 0x0f, 0x09], [0x02, 0x06]]
            @test BitcoinPrimitives.op_2drop(stack) == true
            @test stack == [[0x04]]
        end
        @testset "op_2dup" begin
            stack = [[0x04]]
            @test_throws AssertionError BitcoinPrimitives.op_2dup(stack)
            stack = [[0x04], [0x06, 0x0f, 0x09], [0x02, 0x06]]
            @test BitcoinPrimitives.op_2dup(stack) == true
            @test stack == [[0x04], [0x06, 0x0f, 0x09], [0x02, 0x06], [0x06, 0x0f, 0x09], [0x02, 0x06]]
        end
        @testset "op_3dup" begin
            stack = [[0x04], [0x06, 0x0f, 0x09]]
            @test_throws AssertionError BitcoinPrimitives.op_3dup(stack)
            stack = [[0x04], [0x06, 0x0f, 0x09], [0x02, 0x06]]
            @test BitcoinPrimitives.op_3dup(stack) == true
            @test stack == [[0x04], [0x06, 0x0f, 0x09], [0x02, 0x06], [0x04], [0x06, 0x0f, 0x09], [0x02, 0x06]]
        end
        @testset "op_2over" begin
            stack = [[0x01], [0x02], [0x03]]
            @test_throws AssertionError BitcoinPrimitives.op_2over(stack)
            stack = [[0x01], [0x02], [0x03], [0x04], [0x05]]
            @test BitcoinPrimitives.op_2over(stack) == true
            @test stack == [[0x01], [0x02], [0x03], [0x04], [0x05], [0x02], [0x03]]
        end
        @testset "op_2rot" begin
            stack = [[0x01], [0x02], [0x03], [0x04], [0x05]]
            @test BitcoinPrimitives.op_2rot(stack) == false
            stack = [[0x01], [0x02], [0x03], [0x04], [0x05], [0x06], [0x07]]
            @test BitcoinPrimitives.op_2rot(stack) == true
            @test stack == [[0x01], [0x02], [0x03], [0x04], [0x05], [0x06], [0x07], [0x02], [0x03]]
        end
        @testset "op_2swap" begin
            stack = [[0x01], [0x02], [0x03]]
            @test_throws AssertionError BitcoinPrimitives.op_2swap(stack)
            stack = [[0x01], [0x02], [0x03], [0x04], [0x05]]
            @test BitcoinPrimitives.op_2swap(stack) == true
            @test stack == [[0x01], [0x04], [0x05], [0x02], [0x03]]
        end
        @testset "op_ifdup" begin
            stack = Vector{UInt8}[]
            @test_throws AssertionError BitcoinPrimitives.op_ifdup(stack)
            stack = [[0x01], [0x00]]
            @test BitcoinPrimitives.op_ifdup(stack) == true
            @test stack == [[0x01], [0x00]]
            stack = [[0x01], [0x02]]
            @test BitcoinPrimitives.op_ifdup(stack) == true
            @test stack == [[0x01], [0x02], [0x02]]
        end
        @testset "op_depth" begin
            stack = Vector{UInt8}[]
            @test BitcoinPrimitives.op_depth(stack) == true
            @test stack == [[0x00]]
            stack = [[0x01], [0x02]]
            @test BitcoinPrimitives.op_depth(stack) == true
            @test stack == [[0x01], [0x02], [0x02]]
        end
        @testset "op_dup" begin
            stack = Vector{UInt8}[]
            @test_throws AssertionError BitcoinPrimitives.op_dup(stack)
            stack = [[0x01], [0x02], [0x03]]
            @test BitcoinPrimitives.op_dup(stack) == true
            @test stack == [[0x01], [0x02], [0x03], [0x03]]
        end
        @testset "op_nip" begin
            stack = Vector{UInt8}[[0x01]]
            @test_throws AssertionError BitcoinPrimitives.op_nip(stack)
            stack = [[0x01], [0x02], [0x03]]
            @test BitcoinPrimitives.op_nip(stack) == true
            @test stack == [[0x01], [0x03]]
        end
        @testset "op_over" begin
            stack = [[0x01]]
            @test_throws AssertionError BitcoinPrimitives.op_over(stack)
            stack = [[0x01], [0x02], [0x03], [0x04]]
            @test BitcoinPrimitives.op_over(stack) == true
            @test stack == [[0x01], [0x02], [0x03], [0x04], [0x03]]
        end
        @testset "op_pick" begin
            stack = Vector{UInt8}[[0x01]]
            @test BitcoinPrimitives.op_pick(stack) == false
            stack = [[0x01], [0x0b], [0x03], [0x05]]
            @test BitcoinPrimitives.op_pick(stack) == false
            stack = [[0x01], [0x0b], [0x03], [0x02]]
            @test BitcoinPrimitives.op_pick(stack) == true
            @test stack == [[0x01], [0x0b], [0x03], [0x0b]]
        end
        @testset "op_roll" begin
            stack = Vector{UInt8}[]
            @test_throws AssertionError BitcoinPrimitives.op_roll(stack)
            stack = [[0x01], [0x02], [0x03], [0x04]]
            @test BitcoinPrimitives.op_roll(stack) == false
            @test stack == [[0x01], [0x02], [0x03]]
            stack = [[0x01], [0x0b], [0x03], [0x00]]
            @test BitcoinPrimitives.op_roll(stack) == true
            @test stack == [[0x01], [0x0b], [0x03]]
            stack = [[0x01], [0x0b], [0x03], [0x02]]
            @test BitcoinPrimitives.op_roll(stack) == true
            @test stack == [[0x01], [0x03], [0x0b]]
        end
        @testset "op_rot" begin
            stack = [[0x01], [0x02]]
            @test_throws AssertionError BitcoinPrimitives.op_rot(stack)
            stack = [[0x01], [0x02], [0x03], [0x04]]
            @test BitcoinPrimitives.op_rot(stack) == true
            @test stack == [[0x01], [0x03], [0x04], [0x02]]
        end
        @testset "op_swap" begin
            stack = [[0x01]]
            @test_throws AssertionError BitcoinPrimitives.op_swap(stack)
            stack = [[0x01], [0x02], [0x03], [0x04]]
            @test BitcoinPrimitives.op_swap(stack) == true
            @test stack == [[0x01], [0x02], [0x04], [0x03]]
        end
        @testset "op_tuck" begin
            stack = [[0x01]]
            @test_throws AssertionError BitcoinPrimitives.op_tuck(stack)
            stack = [[0x01], [0x02], [0x03], [0x04]]
            @test BitcoinPrimitives.op_tuck(stack) == true
            @test stack == [[0x01], [0x02], [0x04], [0x03], [0x04]]
        end
    end
    @testset "String splice operations" begin
        @testset "op_size" begin
            stack = Vector{UInt8}[]
            @test_throws AssertionError BitcoinPrimitives.op_size(stack)
            stack = [[0x01], [0x02], [0x01, 0x08]]
            @test BitcoinPrimitives.op_size(stack) == true
            @test stack == [[0x01], [0x02], [0x01, 0x08], [0x02]]
        end
    end
    @testset "Binary arithmetic and conditionals" begin
        @testset "op_equal" begin
            stack = Vector{UInt8}[]
            @test_throws AssertionError BitcoinPrimitives.op_equal(stack)
            stack = [[0x01], [0x02], [0x03], [0x04]]
            @test BitcoinPrimitives.op_equal(stack) == true
            @test stack == [[0x01], [0x02], [0x00]]
            stack = [[0x01], [0x02], [0x03], [0x03]]
            @test BitcoinPrimitives.op_equal(stack) == true
            @test stack == [[0x01], [0x02], [0x01]]
        end
        @testset "op_equalverify" begin
            stack = Vector{UInt8}[]
            @test_throws AssertionError BitcoinPrimitives.op_equalverify(stack)
            stack = [[0x01], [0x02], [0x03], [0x04]]
            @test BitcoinPrimitives.op_equalverify(stack) == false
            stack = [[0x01], [0x02], [0x03], [0x03]]
            @test BitcoinPrimitives.op_equalverify(stack) == true
        end
    end
    @testset "Numeric operators" begin
        empty = Vector{UInt8}[]
        stack = Vector{UInt8}[[0x03]]
        @testset "op_1add" begin
            @test_throws AssertionError BitcoinPrimitives.op_1add(empty)
            @test BitcoinPrimitives.op_1add(stack) == true
            @test stack == [[0x04]]
        end
        @testset "op_1sub" begin
            @test_throws AssertionError BitcoinPrimitives.op_1sub(empty)
            @test BitcoinPrimitives.op_1sub(stack) == true
            @test stack == [[0x03]]
        end
        @testset "op_negate" begin
            @test_throws AssertionError BitcoinPrimitives.op_negate(empty)
            @test BitcoinPrimitives.op_negate(stack) == true
            @test stack == [[0x83]]
        end
        @testset "op_abs" begin
            @test_throws AssertionError BitcoinPrimitives.op_abs(empty)
            @test BitcoinPrimitives.op_abs(stack) == true
            @test stack == [[0x03]]
            @test BitcoinPrimitives.op_abs(stack) == true
            @test stack == [[0x03]]
        end
        @testset "op_not" begin
            @test_throws AssertionError BitcoinPrimitives.op_not(empty)
            @test BitcoinPrimitives.op_not(stack) == true
            @test stack == [[0x00]]
            @test BitcoinPrimitives.op_not(stack) == true
            @test stack == [[0x01]]
        end
        @testset "op_0notequal" begin
            @test_throws AssertionError BitcoinPrimitives.op_0notequal(empty)
            stack = [[0x1a]]
            @test BitcoinPrimitives.op_0notequal(stack) == true
            @test stack == [[0x01]]
            stack = [[0x00]]
            @test BitcoinPrimitives.op_0notequal(stack) == true
            @test stack == [[0x00]]
        end
        @testset "op_add" begin
            @test_throws AssertionError BitcoinPrimitives.op_add(empty)
            stack = [[0x09], [0x02], [0x03]]
            @test BitcoinPrimitives.op_add(stack) == true
            @test stack == [[0x09], [0x05]]
        end
        @testset "op_sub" begin
            @test_throws AssertionError BitcoinPrimitives.op_sub(empty)
            @test BitcoinPrimitives.op_sub(stack) == true
            @test stack == [[0x04]]
        end
        @testset "op_booland" begin
            @test_throws AssertionError BitcoinPrimitives.op_booland(empty)
            stack = [[0x00], [0x00]]
            @test BitcoinPrimitives.op_booland(stack) == true
            @test stack == [[0x00]]
            stack = [[0x01], [0x00]]
            @test BitcoinPrimitives.op_booland(stack) == true
            @test stack == [[0x00]]
            stack = [[0x00], [0x01]]
            @test BitcoinPrimitives.op_booland(stack) == true
            @test stack == [[0x00]]
            stack = [[0x01], [0x01]]
            @test BitcoinPrimitives.op_booland(stack) == true
            @test stack == [UInt8[0x01]]
        end
        @testset "op_boolor" begin
            @test_throws AssertionError BitcoinPrimitives.op_boolor(empty)
            stack = [[0x00], [0x00]]
            @test BitcoinPrimitives.op_boolor(stack) == true
            @test stack == [[0x00]]
            stack = [[0x00], [0x01]]
            @test BitcoinPrimitives.op_boolor(stack) == true
            @test stack == [UInt8[0x01]]
            stack = [[0x01], [0x00]]
            @test BitcoinPrimitives.op_boolor(stack) == true
            @test stack == [UInt8[0x01]]
            stack = [[0x01], [0x01]]
            @test BitcoinPrimitives.op_boolor(stack) == true
            @test stack == [UInt8[0x01]]
        end
        @testset "op_numequal" begin
            @test_throws AssertionError BitcoinPrimitives.op_numequal(empty)
            stack = [[0x01], [0x02], [0x00], [0x00]]
            @test BitcoinPrimitives.op_numequal(stack) == true
            @test stack == [[0x01], [0x02], [0x01]]
            @test BitcoinPrimitives.op_numequal(stack) == true
            @test stack == [[0x01], [0x00]]
        end
        @testset "op_numequalverify" begin
            @test_throws AssertionError BitcoinPrimitives.op_numequalverify(empty)
            stack = [[0x01], [0x02], [0x00], [0x00]]
            @test BitcoinPrimitives.op_numequalverify(stack) == true
            @test stack == [[0x01], [0x02]]
            @test BitcoinPrimitives.op_numequalverify(stack) == false
            @test stack == empty
        end
        @testset "op_numnotequal" begin
            @test_throws AssertionError BitcoinPrimitives.op_numnotequal(empty)
            stack = [[0x01], [0x02], [0x03], [0x03]]
            @test BitcoinPrimitives.op_numnotequal(stack) == true
            @test stack == [[0x01], [0x02], [0x00]]
            stack = [[0x01], [0x02], [0x03], [0x04]]
            @test BitcoinPrimitives.op_numnotequal(stack) == true
            @test stack == [[0x01], [0x02], [0x01]]
        end
        @testset "op_lessthan" begin
            @test_throws AssertionError BitcoinPrimitives.op_lessthan(empty)
            stack = [[0x01], [0x01], [0x03], [0x04]]
            @test BitcoinPrimitives.op_lessthan(stack) == true
            @test stack == [[0x01], [0x01], [0x01]]
            @test BitcoinPrimitives.op_lessthan(stack) == true
            @test stack == [[0x01], [0x00]]
        end
        @testset "op_greaterthan" begin
            @test_throws AssertionError BitcoinPrimitives.op_greaterthan(empty)
            stack = [[0x04], [0x01], [0x02], [0x01]]
            @test BitcoinPrimitives.op_greaterthan(stack) == true
            @test stack == [[0x04], [0x01], [0x01]]
            @test BitcoinPrimitives.op_greaterthan(stack) == true
            @test stack == [[0x04], [0x00]]
        end
        @testset "op_lessthanorequal" begin
            @test_throws AssertionError BitcoinPrimitives.op_lessthanorequal(empty)
            stack = [[0x01], [0x02], [0x04], [0x04]]
            @test BitcoinPrimitives.op_lessthanorequal(stack) == true
            @test stack == [[0x01], [0x02], [0x01]]
            @test BitcoinPrimitives.op_lessthanorequal(stack) == true
            @test stack == [[0x01], [0x00]]
        end
        @testset "op_greaterthanorequal" begin
            @test_throws AssertionError BitcoinPrimitives.op_greaterthanorequal(empty)
            stack = [[0x04], [0x00], [0x02], [0x02]]
            @test BitcoinPrimitives.op_greaterthanorequal(stack) == true
            @test stack == [[0x04], [0x00], [0x01]]
            @test BitcoinPrimitives.op_greaterthanorequal(stack) == true
            @test stack == [[0x04], [0x00]]
        end
        @testset "op_min" begin
            @test_throws AssertionError BitcoinPrimitives.op_min(empty)
            stack = [[0x01], [0x08], [0x03], [0x04]]
            @test BitcoinPrimitives.op_min(stack) == true
            @test stack == [[0x01], [0x08], [0x03]]
            @test BitcoinPrimitives.op_min(stack) == true
            @test stack == [[0x01], [0x03]]
        end
        @testset "op_max" begin
            @test_throws AssertionError BitcoinPrimitives.op_max(empty)
            stack = [[0x01], [0x08], [0x03], [0x04]]
            @test BitcoinPrimitives.op_max(stack) == true
            @test stack == [[0x01], [0x08], [0x04]]
            @test BitcoinPrimitives.op_max(stack) == true
            @test stack == [[0x01], [0x08]]
        end
        @testset "op_within" begin
            @test_throws AssertionError BitcoinPrimitives.op_within(empty)
            stack = [[0xa1], [0x02], [0x02], [0x04]]
            @test BitcoinPrimitives.op_within(stack) == true
            @test stack == [[0xa1], [0x01]]
            stack = [[0xa1], [0x08], [0x02], [0x04]]
            @test BitcoinPrimitives.op_within(stack) == true
            @test stack == [[0xa1], [0x00]]
        end
    end
    @testset "Cryptographic and hashing operations" begin
        empty = Vector{UInt8}[]
        stack = Vector{UInt8}[[0x03], [0x04]]
        @testset "op_ripemd160" begin
            @test_throws AssertionError BitcoinPrimitives.op_ripemd160(empty)
            @test BitcoinPrimitives.op_ripemd160(stack) == true
            @test stack == Vector{UInt8}[[0x03], [0x44, 0x9b, 0x34, 0xb6, 0xa3, 0x41, 0x19, 0x43, 0xe3, 0x3a, 0x25, 0x87, 0xeb, 0xf2, 0x81, 0xca, 0xff, 0x16, 0x74, 0x98]]
        end
        stack = Vector{UInt8}[[0x03], [0x04]]
        @testset "op_sha1" begin
            @test_throws AssertionError BitcoinPrimitives.op_sha1(empty)
            @test BitcoinPrimitives.op_sha1(stack) == true
            @test stack == Vector{UInt8}[[0x03], [0xa4, 0x2c, 0x6c, 0xf1, 0xde, 0x3a, 0xbf, 0xde, 0xa9, 0xb9, 0x5f, 0x34, 0x68, 0x7c, 0xbb, 0xe9, 0x2b, 0x9a, 0x73, 0x83]]
        end
        stack = Vector{UInt8}[[0x03], [0x04]]
        @testset "op_sha256" begin
            @test_throws AssertionError BitcoinPrimitives.op_sha256(empty)
            @test BitcoinPrimitives.op_sha256(stack) == true
            @test stack == Vector{UInt8}[[0x03], [0xe5, 0x2d, 0x9c, 0x50, 0x8c, 0x50, 0x23, 0x47, 0x34, 0x4d, 0x8c, 0x07, 0xad, 0x91, 0xcb, 0xd6, 0x06, 0x8a, 0xfc, 0x75, 0xff, 0x62, 0x92, 0xf0, 0x62, 0xa0, 0x9c, 0xa3, 0x81, 0xc8, 0x9e, 0x71]]
        end
        stack = Vector{UInt8}[[0x03], [0x04]]
        @testset "op_hash160" begin
            @test_throws AssertionError BitcoinPrimitives.op_hash160(empty)
            @test BitcoinPrimitives.op_hash160(stack) == true
            @test stack == Vector{UInt8}[[0x03], [0x6d, 0x4c, 0x0a, 0xa9, 0x72, 0xc3, 0x14, 0x84, 0x0a, 0xc0, 0x7b, 0xe9, 0x6c, 0x5d, 0xde, 0x9c, 0x71, 0x4c, 0x9c, 0xa4]]
        end
        stack = Vector{UInt8}[[0x03], [0x04]]
        @testset "op_hash256" begin
            @test_throws AssertionError BitcoinPrimitives.op_hash256(empty)
            @test BitcoinPrimitives.op_hash256(stack) == true
            @test stack == Vector{UInt8}[[0x03], [0x21, 0x4e, 0x63, 0xbf, 0x41, 0x49, 0x0e, 0x67, 0xd3, 0x44, 0x76, 0x77, 0x8f, 0x67, 0x07, 0xaa, 0x6c, 0x8d, 0x2c, 0x8d, 0xcc, 0xdf, 0x78, 0xae, 0x11, 0xe4, 0x0e, 0xe9, 0xf9, 0x1e, 0x89, 0xa7]]
        end
        # TODO
        @testset "op_codeseparator" begin
        end
        # TODO
        @testset "op_checksig" begin
            z     = big"0x7c076ff316692a3d7eb3c3bb0f8b1488cf72e1afcd929e29307032997a838a3d"
            sec   = hex2bytes("04887387e452b8eacc4acfde10d9aaf7f6d9a0f975aabb10d006e4da568744d06c61de6d95231cd89026e286df3b6ae4a894a3378e393e93a0f45b666329a0ae34")
            sig   = hex2bytes("3045022000eff69ef2b1bd93a66ed5219add4fb51e11a840f404876325a1e8ffe0529a2c022100c7207fee197d27c618aea621406f6bf5ef6fca38681d82b2f06fddbdce6feab601")
            stack = [sig, sec]
            @test BitcoinPrimitives.op_checksig(stack, z) == true
            @test BitcoinPrimitives.decode_num(stack[1]) == 1
        end
        # TODO
        @testset "op_checksigverify" begin
        end
        @testset "op_checkmultisig" begin
            z     = big"0xe71bfa115715d6fd33796948126f40a8cdd39f187e4afb03896795189fe1423c"
            sig1  = hex2bytes("3045022100dc92655fe37036f47756db8102e0d7d5e28b3beb83a8fef4f5dc0559bddfb94e02205a36d4e4e6c7fcd16658c50783e00c341609977aed3ad00937bf4ee942a8993701")
            sig2  = hex2bytes("3045022100da6bee3c93766232079a01639d07fa869598749729ae323eab8eef53577d611b02207bef15429dcadce2121ea07f233115c6f09034c0be68db99980b9a6c5e75402201")
            sec1  = hex2bytes("022626e955ea6ea6d98850c994f9107b036b1334f18ca8830bfff1295d21cfdb70")
            sec2  = hex2bytes("03b287eaf122eea69030a0e9feed096bed8045c8b98bec453e1ffac7fbdbd4bb71")
            stack = [UInt8[], sig1, sig2, [0x02], sec1, sec2, [0x02]]
            @test BitcoinPrimitives.op_checkmultisig(stack, z)
            @test BitcoinPrimitives.decode_num(stack[1]) == 1
        end
        # TODO
        @testset "op_checkmultisigverify" begin
        end
    end
    @testset "Nonoperators" begin
    end
    @testset "Reserved OP codes for internal use by the parser operators" begin
    end
end


text:
#=
Shows

@author : Spencer Lyon <spencer.lyon@nyu.edu>
          Victoria Gregory <victoria.gregory@nyu.edu>

@date: 2014-08-05

References
----------

http://quant-econ.net/career.html
=#

srand(41)  # reproducible results
wp = CareerWorkerProblem()
v_init = fill(100.0, wp.N, wp.N)
func(x) = update_bellman(wp, x)
v = compute_fixed_point(func, v_init, max_iter=500, verbose=false)

surface(wp.theta, wp.epsilon, v', alpha=0.5,
        linewidth=0.25, zlims=(150, 200),
        xlabel="theta", ylabel="epsilon", cbar=false)


text:
using .IntervalArithmetic

# Method used for Taylor1{Interval{T}}^n
function ^(a::Taylor1{T}, n::Integer) where {T<:Interval}
    n == 0 && return one(a)
    n == 1 && return copy(a)
    n == 2 && return square(a)
    n < 0 && return a^float(n)
    return power_by_squaring(a, n)
end

function evaluate(a::Taylor1, dx::Interval)
    order = a.order
    uno = one(dx)
    dx2 = dx^2
    if iseven(order)
        kend = order-2
        @inbounds sum_even = a[end]*uno
        @inbounds sum_odd = a[end-1]*zero(dx)
    else
        kend = order-3
        @inbounds sum_odd = a[end]*uno
        @inbounds sum_even = a[end-1]*uno
    end
    @inbounds for k in kend:-2:0
        sum_odd = sum_odd*dx2 + a[k+1]
        sum_even = sum_even*dx2 + a[k]
    end
    return sum_even + sum_odd*dx
end

function evaluate(a::TaylorN, dx::IntervalBox{N,T}) where {N,T}

    @assert N == get_numvars()
    a_length = length(a)
    suma = zero(constant_term(a)) + Interval{T}(0, 0)
    @inbounds for homPol in length(a):-1:1
        suma += evaluate(a.coeffs[homPol], dx)
    end

    return suma
end

function evaluate(a::HomogeneousPolynomial, dx::IntervalBox{N,T}) where {T, N}

    @assert N == get_numvars()
    dx == IntervalBox(-1..1, Val(N)) && return _evaluate(a, dx, Val(true))
    dx == IntervalBox( 0..1, Val(N)) && return _evaluate(a, dx, Val(false))

    return evaluate(a, dx...)
end

function _evaluate(a::HomogeneousPolynomial, dx::IntervalBox{N,T}, ::Val{true} ) where {T, N}
    a.order == 0 && return a[1] + Interval{T}(0, 0)

    ct = coeff_table[a.order+1]
    @inbounds suma = a[1]*Interval{T}(0,0)

    Ieven = Interval{T}(0,1)
    for (i,a_coeff) in enumerate(a.coeffs)
        iszero(a_coeff) && continue
        if isodd(sum(ct[i]))
            suma += sum(a_coeff) * dx[1]
            continue
        end
        @inbounds tmp = iseven(ct[i][1]) ? Ieven : dx[1]
        for n in 2:N
            @inbounds vv = iseven(ct[i][n]) ? Ieven : dx[1]
            tmp *= vv
        end
        suma += a_coeff * tmp
    end
    return suma
end

function _evaluate(a::HomogeneousPolynomial, dx::IntervalBox{N,T}, ::Val{false} ) where {T, N}
    a.order == 0 && return a[1] + Interval{T}(0, 0)

    @inbounds suma = zero(a[1])*dx[1]
    @inbounds for homPol in a.coeffs
        suma += homPol*dx[1]
    end
    return suma
end

"""
    normalize_taylor(a::Taylor1, I::Interval, symI::Bool=true)

Normalizes `a::Taylor1` such that the interval `I` is mapped
by an affine transformation to the interval `-1..1` (`symI=true`)
or to `0..1` (`symI=false`).
"""
normalize_taylor(a::Taylor1, I::Interval{T}, symI::Bool=true) where {T} =
    _normalize(a, I, Val(symI))

"""
    normalize_taylor(a::TaylorN, I::IntervalBox, symI::Bool=true)

Normalize `a::TaylorN` such that the intervals in `I::IntervalBox`
are mapped by an affine transformation to the intervals `-1..1`
(`symI=true`) or to `0..1` (`symI=false`).
"""
normalize_taylor(a::TaylorN, I::IntervalBox{N,T}, symI::Bool=true) where {N,T} =
    _normalize(a, I, Val(symI))

#  I -> -1..1
function _normalize(a::Taylor1, I::Interval{T}, ::Val{true}) where {T}
    order = get_order(a)
    t = Taylor1(T, order)
    tnew = mid(I) + t*radius(I)
    return a(tnew)
end

#  I -> 0..1
function _normalize(a::Taylor1, I::Interval{T}, ::Val{false}) where {T}
    order = get_order(a)
    t = Taylor1(T, order)
    tnew = inf(I) + t*diam(I)
    return a(tnew)
end

#  I -> IntervalBox(-1..1, Val(N))
function _normalize(a::TaylorN, I::IntervalBox{N,T}, ::Val{true}) where {N,T}
    order = get_order(a)
    x = Vector{typeof(a)}(undef, N)
    for ind in eachindex(x)
        x[ind] = mid(I[ind]) + TaylorN(ind, order=order)*radius(I[ind])
    end
    return a(x)
end

#  I -> IntervalBox(0..1, Val(N))
function _normalize(a::TaylorN, I::IntervalBox{N,T}, ::Val{false}) where {N,T}
    order = get_order(a)
    x = Vector{typeof(a)}(undef, N)
    for ind in eachindex(x)
        x[ind] = inf(I[ind]) + TaylorN(ind, order=order)*diam(I[ind])
    end
    return a(x)
end


text:
# Machine learning API for gradient boosting.
module ML

importall GradientBoost.LossFunctions
importall GradientBoost.GB
importall GradientBoost.GBDecisionTree
importall GradientBoost.GBBaseLearner

export GBLearner,
       fit!,
       predict!,
       LossFunction,
       LeastSquares,
       LeastAbsoluteDeviation,
       BinomialDeviance,
       GBDT,
       GBBL,
       learner_fit,
       learner_predict

# Gradient boosting learner as defined by ML API.
type GBLearner
  algorithm::GBAlgorithm
  output::Symbol
  model

  function GBLearner(algorithm, output=:regression)
    new(algorithm, output, nothing)
  end
end

function fit!(gbl::GBLearner, instances, labels)
  error("Instance type: $(typeof(instances))
    and label type: $(typeof(labels)) together is currently not supported.")
end
function predict!(gbl::GBLearner, instances)
  error("Instance type: $(typeof(instances)) is currently not supported.")
end

function fit!(gbl::GBLearner,
  instances::Matrix{Float64}, labels::Vector{Float64})

  # No special processing required.
  gbl.model = fit(gbl.algorithm, instances, labels)
end

function predict!(gbl::GBLearner,
  instances::Matrix{Float64})

  # Predict with GB algorithm
  predictions = predict(gbl.model, instances)

  # Postprocess according to output and loss function
  predictions = postprocess_pred(
    gbl.output, gbl.algorithm.loss_function, predictions
  )

  predictions
end

# Postprocesses predictions according to
# output and loss function.
function postprocess_pred(
  output::Symbol, lf::LossFunction, predictions::Vector{Float64})

  if output == :class && typeof(lf) <: BinomialDeviance
    return round(logistic(predictions))
  elseif output == :class_prob && typeof(lf) <: BinomialDeviance
    return logistic(predictions)
  elseif output == :regression && !(typeof(lf) <: BinomialDeviance)
    return predictions
  else
    error("Cannot handle $(output) and $(typeof(lf)) together.")
  end
end

# Logistic function.
function logistic(x)
  1 ./ (1 .+ exp(-x))
end

end # module


text:
module Capture

#using Debug

import CompilerTools
import ..API
import ..operators
import ..binary_operators
import ..rename_if_needed
import ParallelAccelerator

import CompilerTools.DebugMsg
DebugMsg.init()

const binary_operator_set = Set(binary_operators)

#data_source_num = 0

ref_assign_map = Dict{Symbol, Symbol}(
    :(+=) => :pa_api_add,
    :(-=) => :pa_api_sub,
    :(*=) => :pa_api_mul,
    :(/=) => :pa_api_div,
    :(.+=) => :pa_api_elem_add,
    :(.-=) => :pa_api_elem_sub,
    :(.*=) => :pa_api_elem_mul,
    :(./=) => :pa_api_elem_div
)

type process_node_state
    array
    dim
end

function process_node(node::Symbol, state, top_level_number, is_top_level, read)
    if node == :end
        if state == nothing
            @dprintln(3, "Found :end symbol in AST but had no array name to generate replacement length expression.")
        else
            @assert isa(state, process_node_state) "process_node state was not nothing nor a process_node_state."
            return Expr(:call, GlobalRef(Base, :arraysize), state.array, state.dim)
        end
    end

    CompilerTools.AstWalker.ASTWALK_RECURSE
end

"""
At macro level, we translate function calls and operators that matches operator names
in our API module to direct call to those in the API module.
"""
function process_node(node::Expr, state, top_level_number, is_top_level, read)
    head = node.head
    @dprintln(3, "api-capture process_node ", node, " head = ", head)
    if head == :comparison
        # Surprise! Expressions like x > y are not :call in AST at macro level
        opr = node.args[2]
        if isa(opr, Symbol) && in(opr, operators)
            node.args[2] = GlobalRef(API, opr)
        end
    elseif head == :call
        # f(...)
        opr = node.args[1]
        process_operator(node, opr)
    elseif head == :(=)
        process_assignment(node, node.args[1], node.args[2])
    elseif haskey(ref_assign_map, head) && isa(node.args[1], Expr) && node.args[1].head == :ref
        # x[...] += ...
        lhs = node.args[1].args[1]
        idx = node.args[1].args[2:end]
        @dprintln(3, "idx before end substitution in process_node = ", idx)
        for i = 1:length(idx)
            idx[i] = CompilerTools.AstWalker.AstWalk(idx[i], process_node, process_node_state(lhs, i))
        end
        @dprintln(3, "idx after end substitution in process_node = ", idx)
        rhs = node.args[2]
        tmpvar = gensym()
        node.head = :block
        node.args = Any[
        Expr(:(=), tmpvar, Expr(:call, GlobalRef(API, ref_assign_map[head]),
        Expr(:call, GlobalRef(API, :getindex), lhs, idx...),
        rhs)),
        Expr(:call, GlobalRef(API, :setindex!), lhs, tmpvar, idx...),
        tmpvar]
    elseif haskey(ref_assign_map, head)
        # x += ...
        lhs = node.args[1]
        rhs = node.args[2]
        tmpvar = gensym()
        node.head = ((string(head)[1] == '.') && VERSION >= v"0.5.0-dev+5381") ? :(.=) : :(=)
        node.args = Any[ lhs, Expr(:call, GlobalRef(API, ref_assign_map[head]), lhs, rhs) ]
    elseif node.head == :ref
        node.head = :call
        @dprintln(3, "ref array = ", node.args[1], " args = ", node.args[2:end])
        node.args[1] = CompilerTools.AstWalker.AstWalk(node.args[1], process_node, nothing)
        for i = 2:length(node.args)
            node.args[i] = CompilerTools.AstWalker.AstWalk(node.args[i], process_node, process_node_state(node.args[1], i-1))
        end
        @dprintln(3, "after args = ", node.args[2:end])
        node.args = Any[GlobalRef(API, :getindex), node.args...]
        return node
    end

    CompilerTools.AstWalker.ASTWALK_RECURSE
end

function process_node(node::Any, state, top_level_number, is_top_level, read)
    CompilerTools.AstWalker.ASTWALK_RECURSE
end

function process_operator(node::Expr, opr::Symbol)
    rename_opr = rename_if_needed(opr)
    api_opr = GlobalRef(API, rename_opr)
    if in(opr, operators)
        node.args[1] = api_opr
    end
    if in(opr, binary_operator_set) && length(node.args) > 3
        # we'll turn multiple arguments into pairwise
        expr = foldl((a, b) -> Expr(:call, api_opr, a, b), node.args[2:end])
        node.args = expr.args
    end
end

function process_operator(node::Expr, opr::Any)
end

function process_assignment(node, lhs::Expr, rhs::Any)
    if lhs.head == :ref
        # x[...] = ...
        lhs = node.args[1].args[1]
        idx = node.args[1].args[2:end]
        @dprintln(3, "idx before end substitution in process_assignment = ", idx)
        for i = 1:length(idx)
            idx[i] = CompilerTools.AstWalker.AstWalk(idx[i], process_node, process_node_state(lhs, i))
        end
        @dprintln(3, "idx after end substitution in process_assignment = ", idx)
        rhs = node.args[2]
        tmpvar = gensym()
        node.head = :block
        node.args = Any[
        Expr(:(=), tmpvar, rhs),
        Expr(:call, GlobalRef(API, :setindex!), lhs, tmpvar, idx...),
        tmpvar]
    end
end

function process_assignment(node, lhs::Any, rhs::Any)
end

#function get_unique_data_source_num()
#    global data_source_num
#    data_source_num += 1
#    return data_source_num
#end

function translate_par(args)
    na = length(args)
    @assert (na > 0) "Expect a for loop as argument to @par"
    redVars = Array{Symbol}(na-1)
    redOps = Array{Symbol}(na-1)
    loop = args[end]
    if !isa(loop,Expr) || !(loop.head === :for)
        error("malformed @par loop")
    end
    if !isa(loop.args[1], Expr)
        error("maltformed for loop")
    end
    for i = 1:na-1
        @assert (isa(args[i], Expr) && args[i].head == :call) "Expect @par reduction in the form of var(op), but got " * string(args[i])
        v = args[i].args[1]
        op = args[i].args[2]
        #println("got reduction variable ", v, " with function ", op)
        @assert (isa(v, Symbol)) "Expect reduction variable to be symbols, but got " * string(v)
        @assert (isa(op, Symbol)) "Expect reduction operator to be symbols, but got " * string(op)
        redVars[i] = v
        redOps[i] = op
    end
    if loop.args[1].head == :block
        loopheads = loop.args[1].args
    else
        loopheads = Any[loop.args[1]]
    end
    ndim = length(loopheads)
    body = loop.args[2]
    params = Array{Symbol}(ndim)
    indices = Array{Symbol}(ndim)
    ranges = Array{Any}(ndim)
    headers = Array{Any}(ndim)
    for i = 1:ndim
        r = loopheads[i]
        assert(r.head == :(=))
        indices[i] = r.args[1]
        params[i] = gensym(string(indices[i]))
        ranges[i] = r.args[2]
        #headers[i] = Expr(:(=), indices[i], Expr(:call, GlobalRef(Base, :getindex), ranges[i], params[i]))
        #first(v) + (i-1)*step(v)
        first = GlobalRef(Base, :first)
        step = GlobalRef(Base, :step)
        headers[i] = :($(indices[i]) = $first($(ranges[i])) + ($(params[i]) - 1) * $step($(ranges[i])))
        # if range is 1:N, generate simple assignment instead of formula
        # this helps array index analysis (motivated by Kernel score in HPAT)
        if isa(ranges[i], Expr) && isStart1UnitRange(ranges[i])
            headers[i] = :($(indices[i]) = $(params[i]))
        end
    end
    args = Expr(:tuple, params...)
    dims = Expr(:tuple, [ Expr(:call, GlobalRef(Base, :length), r) for r in ranges ]...)
    tmpret = gensym("tmp")
    tmpinits = [ :($idx = 1) for idx in params ]
    #println("indices = ", indices)
    #println("params = ", params)
    #println("ranges = ", ranges)
    #println("headers = ", headers)
    #println("tmpinits = ", tmpinits)
    #println("body = ", body)
    #if na==1
    #    thecall = :(pfor($(make_pfor_body(var, body, :therange)), length(therange)))
    #else
    #    thecall = :(preduce($(esc(reducer)),
    #                        $(make_preduce_body(reducer, var, body, :therange)), length(therange)))
    #end
    #localize_vars(quote therange = $(esc(r)); $thecall; end)
    ast =Expr(:call, GlobalRef(API, :cartesianmapreduce),
                    :($args -> let $(headers...); $body; return end), :($dims))
    for i=1:na-1
        redVar = redVars[i]
        redArg = gensym(string(redVar))
        opr = redOps[i]
        if in(opr, operators)
            if haskey(API.rename_forward, opr)
                opr = API.rename_forward[opr]
            end
            opr = GlobalRef(API, opr)
        end
        redBody = Expr(:(=), redVars[i], Expr(:call, opr, redVar, redArg))
        # redVals = Expr(:call, GlobalRef(Base, :copy), redVar)
        push!(ast.args, :($redArg -> begin $redBody; return $redVar; end, $redVar))
    end
    #println("ast = ", ast)
    ast
    #args -> esc(loop)
end

function isStart1UnitRange(node::Expr)
    if node.head==:(:) && length(node.args)==2 && node.args[1]==1
        return true
    end
    return false
end

function process_par_macro(node::Expr, state, top_level_number, is_top_level, read)
    if node.head === :macrocall
        #println("got par macro, args = ", node.args)
        if node.args[1] == Symbol("@par")
            ast = translate_par(node.args[2:end])
            node.head = ast.head
            node.args = ast.args
            node.typ  = ast.typ
        end
    end
    CompilerTools.AstWalker.ASTWALK_RECURSE
end

function process_par_macro(node::ANY, state, top_level_number, is_top_level, read)
    CompilerTools.AstWalker.ASTWALK_RECURSE
end

end # module Capture


text:
"""
A LAS "variable length record" - the generic way to store extra user or
organization defined binary metadata in LAS files.
"""
struct LasVariableLengthRecord
    reserved::UInt16
    user_id::AbstractString
    record_id::UInt16
    description::AbstractString
    data  # anything with read+write+sizeof methods, like GeoKeys or Vector{UInt8}
end

# Read a variable length metadata record from a stream.
#
# If `extended` is true, the VLR is one of the extended VLR types specified in
# the LAS 1.4 spec which can be larger and come after the point data.
function Base.read(io::IO, ::Type{LasVariableLengthRecord}, extended::Bool=false)
    # `reserved` is meant to be 0 according to the LAS spec 1.4, but earlier
    # versions set it to 0xAABB.  Whatever, I guess we just store&ignore for now.
    # See https://groups.google.com/forum/#!topic/lasroom/SVtNBA2y9iI
    reserved = read(io, UInt16)
    user_id = readstring(io, 16)
    record_id = read(io, UInt16)
    record_data_length::Int = extended ? read(io, UInt64) : read(io, UInt16)
    description = readstring(io, 32)
    data = read_vlr_data(io, record_id, record_data_length)
    LasVariableLengthRecord(
        reserved,
        user_id,
        record_id,
        description,
        data
    )
end

function Base.write(io::IO, vlr::LasVariableLengthRecord, extended::Bool=false)
    write(io, vlr.reserved)
    writestring(io, vlr.user_id, 16)
    write(io, vlr.record_id)
    record_data_length = extended ? UInt64(sizeof(vlr.data)) : UInt16(sizeof(vlr.data))
    write(io, record_data_length)
    writestring(io, vlr.description, 32)
    write(io, vlr.data)
    nothing
end

# size of a VLR in bytes
# assumes it is not extended VLR
Base.sizeof(vlr::LasVariableLengthRecord) = 54 + sizeof(vlr.data)


text:
module DReal

using AbstractDomains
using Compat
using MathProgBase

@windows_only error("Windows not supported")

# @show joinpath(dirname(@__FILE__),"..","deps","deps.jl")
# if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
#     include("../deps/deps.jl")
# else
#     error("DReal not properly installed. Please run Pkg.build(\"DReal\")")
# end

VERSION < v"0.4-" && using Docile

deps_dir = joinpath(joinpath(Pkg.dir("DReal"),"deps"))
prefix = joinpath(deps_dir,"usr")
src_dir = joinpath(prefix,"src")
bin_dir = joinpath(prefix,"bin")
lib_dir = joinpath(prefix,"lib")

try
  @osx? begin
    @compat Libdl.dlopen(joinpath(lib_dir, "libdreal.dylib"), Libdl.RTLD_LAZY|Libdl.RTLD_DEEPBIND|Libdl.RTLD_GLOBAL)
  end : begin
    @compat Libdl.dlopen(joinpath(lib_dir, "libprim.so"), Libdl.RTLD_LAZY|Libdl.RTLD_DEEPBIND|Libdl.RTLD_GLOBAL)
    @compat Libdl.dlopen(joinpath(lib_dir, "libCoinUtils.so"), Libdl.RTLD_LAZY|Libdl.RTLD_DEEPBIND|Libdl.RTLD_GLOBAL)
    @compat Libdl.dlopen(joinpath(lib_dir, "libClp.so"), Libdl.RTLD_LAZY|Libdl.RTLD_DEEPBIND|Libdl.RTLD_GLOBAL)
    @compat Libdl.dlopen(joinpath(lib_dir, "libcapd.so"), Libdl.RTLD_LAZY|Libdl.RTLD_DEEPBIND|Libdl.RTLD_GLOBAL)
    @compat Libdl.dlopen(joinpath(lib_dir, "libibex.so"), Libdl.RTLD_LAZY|Libdl.RTLD_DEEPBIND|Libdl.RTLD_GLOBAL)
    @compat Libdl.dlopen(joinpath(lib_dir, "libdreal.so"), Libdl.RTLD_LAZY|Libdl.RTLD_DEEPBIND|Libdl.RTLD_GLOBAL)
  end
catch e
  println("Could not load required shared libraries")
  rethrow(e)
end

import Base: print, show
import Base:  abs, exp, log,
              ^, sin, cos,
              tan, asin, acos,
              atan, sinh, cosh,
              tanh, atan2, sqrt

import Base: ifelse
import Base: convert, push!
import Base: ==, >=, <=, >, <, /, -, +, *, &, |, !

export
  Context,
  Logic,
  model,
  init_dreal!,
  set_precision!,
  init_dreal!,
  Var,
  ForallVar,
  Ex,
  add!,
  is_satisfiable,
  global_context,
  set_verbosity_level!,
  push_ctx!,
  pop_ctx!,
  reset_ctx!,
  delete_ctx!,
  reset_global_ctx!,
  →,
  implies,
  minimize,
  default_global_context

## Global defaults
const DEFAULT_PRECISION = 0.001

include("wrap_capi.jl")
include("logic.jl")
include("context.jl")
include("environment.jl")
include("expression.jl")
include("construct.jl")
include("optimize.jl")
include("ode.jl")
include("SolverInterface.jl")
include("util.jl")

export
  qf_uf,         # Uninterpreted Functions
  qf_nra,        # Non-Linear Real Arithmetic
  qf_nra_ode,    # Non-Linear Real Arithmetic
  qf_bv,         # BitVectors
  qf_rdl,        # Real difference logics
  qf_idl,        # Integer difference logics
  qf_lra,        # Real linear arithmetic
  qf_lia,        # Integer linear arithmetic
  qf_ufidl,      # UF + IDL
  qf_uflra,      # UF + LRA
  qf_bool,       # Only booleans
  qf_ct          # Cost

init_dreal!()
create_global_ctx!()

end


text:
using Go
using Base.Test

function random_position(seed, board_size, n_moves)
    srand(seed)
    p = Position{board_size}()
    for i in 1:n_moves
        for t in (black, white)
            ok = false
            while !ok
                try
                    p = move!(Position(p), t, rand(1:board_size), rand(1:board_size))
                    ok = true
                catch e
                    if !isa(e, IllegalMove)
                        rethrow()
                    end
                end
            end
        end
    end
    p
end

function random_expression(n)
    @assert n > 7
    e = rand(UInt8, n)
    e[1:4] = map(UInt8, collect("goxp"))
    e[5] = 0x00
    e[6:7] = reinterpret(UInt8, UInt16[n])
    e
end

function random_population(n, m)
    [random_expression(n) for i in 1:m]
end

function compare(path, objs...)
    print(path)
    s = join([sprint(print, x) for x in objs], "\n") * "\n"
    if !exists(path)
        print(" (w)")
        open(path, "w") do io
            print(io, s)
        end
    end
    @test s == open(readall, path, "r")
    println(" ok")
end

include("position.jl")
include("expression.jl")
include("evolution.jl")


text:
using Mocking

cancel_job_1_patch = @patch function cancel_job(a...; b...)
    return Dict{String, Any}()
end

create_compute_environment_1_patch = @patch function create_compute_environment(a...; b...)
    return Dict{String, Any}("computeEnvironmentName" => "C4OnDemand", "computeEnvironmentArn" => "arn:aws:batch:us-east-1:012345678910:compute-environment/C4OnDemand")
end

create_compute_environment_2_patch = @patch function create_compute_environment(a...; b...)
    return Dict{String, Any}("computeEnvironmentName" => "M4Spot", "computeEnvironmentArn" => "arn:aws:batch:us-east-1:012345678910:compute-environment/M4Spot")
end

create_job_queue_1_patch = @patch function create_job_queue(a...; b...)
    return Dict{String, Any}("jobQueueName" => "LowPriority", "jobQueueArn" => "arn:aws:batch:us-east-1:012345678910:job-queue/LowPriority")
end

create_job_queue_2_patch = @patch function create_job_queue(a...; b...)
    return Dict{String, Any}("jobQueueName" => "HighPriority", "jobQueueArn" => "arn:aws:batch:us-east-1:012345678910:job-queue/HighPriority")
end

delete_compute_environment_1_patch = @patch function delete_compute_environment(a...; b...)
    return Dict{String, Any}()
end

delete_job_queue_1_patch = @patch function delete_job_queue(a...; b...)
    return Dict{String, Any}()
end

deregister_job_definition_1_patch = @patch function deregister_job_definition(a...; b...)
    return Dict{String, Any}()
end

describe_compute_environments_1_patch = @patch function describe_compute_environments(a...; b...)
    return Dict{String, Any}("computeEnvironments" => Any[Dict{String, Any}("state" => "ENABLED", "statusReason" => "ComputeEnvironment Healthy", "status" => "VALID", "computeEnvironmentName" => "P2OnDemand", "serviceRole" => "arn:aws:iam::012345678910:role/AWSBatchServiceRole", "computeEnvironmentArn" => "arn:aws:batch:us-east-1:012345678910:compute-environment/P2OnDemand", "type" => "MANAGED", "ecsClusterArn" => "arn:aws:ecs:us-east-1:012345678910:cluster/P2OnDemand_Batch_2c06f29d-d1fe-3a49-879d-42394c86effc", "computeResources" => Dict{String, Any}("instanceTypes" => Any["p2"], "minvCpus" => 0, "subnets" => Any["subnet-220c0e0a", "subnet-1a95556d", "subnet-978f6dce"], "desiredvCpus" => 48, "ec2KeyPair" => "id_rsa", "maxvCpus" => 128, "securityGroupIds" => Any["sg-cf5093b2"], "tags" => Dict{String, Any}("Name" => "Batch Instance - P2OnDemand"), "type" => "EC2", "instanceRole" => "ecsInstanceRole"))])
end

describe_job_definitions_1_patch = @patch function describe_job_definitions(a...; b...)
    return Dict{String, Any}("jobDefinitions" => Any[Dict{String, Any}("revision" => 1, "status" => "ACTIVE", "containerProperties" => Dict{String, Any}("ulimits" => Any[], "volumes" => Any[], "image" => "busybox", "mountPoints" => Any[], "command" => Any["sleep", "60"], "memory" => 128, "environment" => Any[], "vcpus" => 1), "jobDefinitionName" => "sleep60", "jobDefinitionArn" => "arn:aws:batch:us-east-1:012345678910:job-definition/sleep60:1", "type" => "container")])
end

describe_job_queues_1_patch = @patch function describe_job_queues(a...; b...)
    return Dict{String, Any}("jobQueues" => Any[Dict{String, Any}("statusReason" => "JobQueue Healthy", "priority" => 1, "status" => "VALID", "computeEnvironmentOrder" => Any[Dict{String, Any}("order" => 1, "computeEnvironment" => "arn:aws:batch:us-east-1:012345678910:compute-environment/C4OnDemand")], "jobQueueName" => "HighPriority", "jobQueueArn" => "arn:aws:batch:us-east-1:012345678910:job-queue/HighPriority", "state" => "ENABLED")])
end

describe_jobs_1_patch = @patch function describe_jobs(a...; b...)
    return Dict{String, Any}("jobs" => Any[Dict{String, Any}("createdAt" => 1480460782010, "jobId" => "24fa2d7a-64c4-49d2-8b47-f8da4fbde8e9", "parameters" => Dict{String, Any}(), "jobDefinition" => "sleep60", "status" => "SUCCEEDED", "container" => Dict{String, Any}("ulimits" => Any[], "volumes" => Any[], "image" => "busybox", "mountPoints" => Any[], "command" => Any["sleep", "60"], "containerInstanceArn" => "arn:aws:ecs:us-east-1:012345678910:container-instance/5406d7cd-58bd-4b8f-9936-48d7c6b1526c", "memory" => 128, "environment" => Any[], "vcpus" => 1, "exitCode" => 0), "jobQueue" => "arn:aws:batch:us-east-1:012345678910:job-queue/HighPriority", "dependsOn" => Any[], "startedAt" => 1480460816500, "jobName" => "example", "stoppedAt" => 1480460880699)])
end

list_jobs_1_patch = @patch function list_jobs(a...; b...)
    return Dict{String, Any}("jobSummaryList" => Any[Dict{String, Any}("jobId" => "e66ff5fd-a1ff-4640-b1a2-0b0a142f49bb", "jobName" => "example")])
end

list_jobs_2_patch = @patch function list_jobs(a...; b...)
    return Dict{String, Any}("jobSummaryList" => Any[Dict{String, Any}("jobId" => "68f0c163-fbd4-44e6-9fd1-25b14a434786", "jobName" => "example")])
end

list_tags_for_resource_1_patch = @patch function list_tags_for_resource(a...; b...)
    return Dict{String, Any}("tags" => Dict{String, Any}("User" => "JaneDoe", "Department" => "Engineering", "Stage" => "Alpha"))
end

register_job_definition_1_patch = @patch function register_job_definition(a...; b...)
    return Dict{String, Any}("revision" => 1, "jobDefinitionName" => "sleep10", "jobDefinitionArn" => "arn:aws:batch:us-east-1:012345678910:job-definition/sleep10:1")
end

register_job_definition_2_patch = @patch function register_job_definition(a...; b...)
    return Dict{String, Any}("revision" => 1, "jobDefinitionName" => "sleep30", "jobDefinitionArn" => "arn:aws:batch:us-east-1:012345678910:job-definition/sleep30:1")
end

submit_job_1_patch = @patch function submit_job(a...; b...)
    return Dict{String, Any}("jobId" => "876da822-4198-45f2-a252-6cea32512ea8", "jobName" => "example")
end

tag_resource_1_patch = @patch function tag_resource(a...; b...)
    return Dict{String, Any}()
end

terminate_job_1_patch = @patch function terminate_job(a...; b...)
    return Dict{String, Any}()
end

untag_resource_1_patch = @patch function untag_resource(a...; b...)
    return Dict{String, Any}()
end

update_compute_environment_1_patch = @patch function update_compute_environment(a...; b...)
    return Dict{String, Any}("computeEnvironmentName" => "P2OnDemand", "computeEnvironmentArn" => "arn:aws:batch:us-east-1:012345678910:compute-environment/P2OnDemand")
end

update_job_queue_1_patch = @patch function update_job_queue(a...; b...)
    return Dict{String, Any}("jobQueueName" => "GPGPU", "jobQueueArn" => "arn:aws:batch:us-east-1:012345678910:job-queue/GPGPU")
end


text:
email= ENV["NCBI_EMAIL"] #This is an enviroment variable that you need to setup
pmid = 11748933
pmid_list = [24008025, 24170597]

@testset "Export EndNote Citations" begin

    println("-----------------------------------------")
    println("       Export EndNote Citations")

    citation_type="endnote"
    output_file1="11748933.enw"
    output_file2 = "pmid_list.enw"

    export_citation(email, pmid, citation_type, output_file1)
    export_citation(email, pmid_list, citation_type, output_file2)

    #Read back the files and do minimal testing
    ref_lines=[]
    open("./11748933.enw") do f
       ref_lines = readlines(f)
    end

    @test ref_lines[1]== "%0 Journal Article\n"

    ref_lines=[]
    open("./pmid_list.enw") do f
       ref_lines = readlines(f)
    end

    @test ref_lines[1]== "%0 Journal Article\n"
    @test ref_lines[25]== "%0 Journal Article\n"

    # remove temp files
    if isfile(output_file1)
        rm(output_file1)
    end
    if isfile(output_file2)
        rm(output_file2)
    end
end

@testset "Export BibTEX Citation" begin
    println("-----------------------------------------")
    println("       Export BibTEX Citations")
    citation_type="bibtex"
    output_file1="11748933.bibtex"
    output_file2 = "pmid_list.bibtex"

    export_citation(email, pmid, citation_type, output_file1)
    export_citation(email, pmid_list, citation_type, output_file2)

    #Read back the files and do minimal testing
    ref_lines=[]
    open("./11748933.bibtex") do f
       ref_lines = readlines(f)
    end

    @test ref_lines[1]== "@article {PMID:11748933,\n"

    ref_lines=[]
    open("./pmid_list.bibtex") do f
       ref_lines = readlines(f)
    end

    @test ref_lines[1]== "@article {PMID:24008025,\n"
    @test ref_lines[13]== "@article {PMID:24170597,\n"

    # remove temp files
    if isfile(output_file1)
        rm(output_file1)
    end
    if isfile(output_file2)
        rm(output_file2)
    end
end


text:
export compute_body_force_terms1, StaticDomain1, compute_stiffness_matrix_and_internal_force1

@doc raw"""
    StaticDomain1(nodes::Array{Float64}, elements::Array,
        EBC::Array{Int64,1}, g::Array{Float64, 1}, FBC::Array{Int64,1}, f::Array{Float64,1},
        edge_traction_data::Array{Int64,2}=zeros(Int64,0,3)))

Constructs a domain for a scalar state variable in a static equation.

"""
function StaticDomain1(nodes::Array{Float64}, elements::Array,
    EBC::Array{Int64,1}, g::Array{Float64, 1}, FBC::Array{Int64,1}, f::Array{Float64,1},
    edge_traction_data::Array{Int64,2}=zeros(Int64,0,3))
    EBC = reshape(EBC, :, 1)
    FBC = reshape(FBC, :, 1)
    nnodes = size(nodes, 1)
    neles = length(elements)
    ndims = 1

    state = zeros(nnodes * ndims)
    Dstate = zeros(nnodes * ndims)
    LM = Array{Int64}[]
    DOF = Array{Int64}[]
    ID = Int64[]
    neqs = 0
    eq_to_dof = Int64[]
    dof_to_eq = zeros(Bool, nnodes * ndims)
    fext = Float64[]

    npoints = -1
    node_to_point = Int64[]
    history = Dict()

    domain = Domain(nnodes, nodes, neles, elements, ndims, state, Dstate,
    LM, DOF, ID, neqs, eq_to_dof, dof_to_eq,
    EBC, g, FBC, fext, edge_traction_data, 0.0, npoints, node_to_point,
    Int64[], Int64[], Int64[], Float64[],
    Int64[], Int64[], Int64[], Float64[],
    Int64[], Int64[], Int64[], Float64[], history)

    #set fixed(time-independent) Dirichlet boundary conditions
    setConstantDirichletBoundary!(domain, EBC, g)
    #set constant(time-independent) force load boundary conditions
    setConstantNodalForces!(domain, FBC, f)

    assembleSparseMatrixPattern!(domain)

    domain
end

"""
    compute_stiffness_matrix_and_internal_force1(domain::Domain, Δt::Float64 = 0.0)

Assembles the stiffness matrix and internal force.
"""
function compute_stiffness_matrix_and_internal_force1(domain::Domain, Δt::Float64 = 0.0)
    Fint = zeros(Float64, domain.neqs)
    ii = Int64[]; jj = Int64[]; vv = Float64[]
    neles = domain.neles

    for iele  = 1:neles
        element = domain.elements[iele]

        el_nodes = getNodes(element)

        el_eqns = getEqns(domain, iele)

        el_dofs = getDofs(domain, iele)

        el_state  = getState(domain, el_dofs)

        el_Dstate = getDstate(domain, el_dofs)

        fint, stiff  = getStiffAndForce1(element, el_state, el_Dstate, Δt)

        el_eqns_active = el_eqns .>= 1
        Slocal = stiff[el_eqns_active,el_eqns_active]
        Idx = el_eqns[el_eqns_active]
        for i = 1:length(Idx)
            for j = 1:length(Idx)
                push!(ii, Idx[i])
                push!(jj, Idx[j])
                push!(vv, Slocal[i,j])
            end
        end
        Fint[el_eqns[el_eqns_active]] += fint[el_eqns_active]
    end
    Ksparse = sparse(ii, jj, vv, domain.neqs, domain.neqs)

    return Ksparse, Fint
end

@doc raw"""
    compute_body_force_terms1(domain::Domain, f::Function)

Computes the body force
```math
\int_\Omega f \delta v dx
```
"""
function compute_body_force_terms1(domain::Domain, f::Function)

    Fbody = zeros(Float64, domain.neqs)
    neles = domain.neles

    # Loop over the elements in the elementGroup
    for iele  = 1:neles
        element = domain.elements[iele]

        gauss_pts = getGaussPoints(element)
        fvalue = f.(gauss_pts[:,1], gauss_pts[:,2])

        fbody = getBodyForce(element, fvalue)

      # Assemble in the global array
        el_eqns = getEqns(domain, iele)
        el_eqns_active = (el_eqns .>= 1)
        Fbody[el_eqns[el_eqns_active]] += fbody[el_eqns_active]
    end

    return Fbody
end


text:
using Distributions: Normal, Uniform

is_vector(x) = size(x,2) == 1
nrows(x) = size(x,1)
dim_error(message:: String) = throw(DimensionMismatch(message))

"""Sample(X,Y)

Sample holds X array of size n*k (explanatory variables,
independent variables, features or regressor) and Y, a n*1 sized
vector, for observed data (response, dependent variable or regressand).
"""
struct Sample
    X::Array
    Y::Vector
    function Sample(X, Y)
        nrows(X) == nrows(Y) || dim_error("X and Y must have same number of rows")
        is_vector(Y) || dim_error("Y must be a vector: $Y")
        new(X, Y)
    end
end

function add_intercept(X::Array)::Array
    X0 = ones(nrows(X), 1)
    return [X0 X]
end

# TODO: do not add intercept if it is already there.
function add_intercept(sample::Sample)
    return Sample(add_intercept(sample.X), sample.Y)
end

function sample_factory(;x_process, y_process, error_process)
    function sampler(n:: Int)
        X = x_process(n)
        Y = vec(y_process(X) + error_process(X))
        return Sample(X, Y)
    end
end

function normal_sampler(;a, b, β_0, β, sd_e)
    k = length(β)
    du = Uniform(a, b) # this is squares for x only
    dn = Normal(0, sd_e)
    sample_factory(x_process = n -> rand(du, n, k),
                   y_process = X -> β_0 .+ X * β,
                   error_process = X -> rand(dn, nrows(X), 1)
                   )
end


text:
# -*- coding: utf-8 -*-
import HELICS

const h = HELICS

fedinitstring = "--federates=1"
deltat = 0.01

helicsversion = h.helicsGetVersion()

println("PI RECEIVER: Helics version = $helicsversion")

# Create Federate Info object that describes the federate properties */
println("PI RECEIVER: Creating Federate Info")
fedinfo = h.helicsCreateFederateInfo()

# Set Federate name
println("PI RECEIVER: Setting Federate Info Name")
h.helicsFederateInfoSetCoreName(fedinfo, "TestB Federate")

# Set core type from string
println("PI RECEIVER: Setting Federate Info Core Type")
h.helicsFederateInfoSetCoreTypeFromString(fedinfo, "zmq")

# Federate init string
println("PI RECEIVER: Setting Federate Info Init String")
h.helicsFederateInfoSetCoreInitString(fedinfo, fedinitstring)

# Set the message interval (timedelta) for federate. Note that
# HELICS minimum message time interval is 1 ns and by default
# it uses a time delta of 1 second. What is provided to the
# setTimedelta routine is a multiplier for the default timedelta.

# Set one second message interval
println("PI RECEIVER: Setting Federate Info Time Delta")
h.helicsFederateInfoSetTimeProperty(fedinfo, h.HELICS_PROPERTY_TIME_DELTA, deltat)

# Create value federate
println("PI RECEIVER: Creating Value Federate")
vfed = h.helicsCreateValueFederate("TestB Federate", fedinfo)
println("PI RECEIVER: Value federate created")

# Subscribe to PI SENDER's publication
sub = h.helicsFederateRegisterSubscription(vfed, "testA", "")
println("PI RECEIVER: Subscription registered")

h.helicsFederateEnterExecutingMode(vfed)
println("PI RECEIVER: Entering execution mode")

value = 0.0
prevtime = 0

currenttime = -1

while currenttime <= 100
    global currenttime
    currenttime = h.helicsFederateRequestTime(vfed, Float64(100))

    value = h.helicsInputGetString(sub)
    println(
        "PI RECEIVER: Received value = $value at time $currenttime from PI SENDER"
    )
end

h.helicsFederateFinalize(vfed)

h.helicsFederateFree(vfed)
h.helicsCloseLibrary()
println("PI RECEIVER: Federate finalized")


text:
using Test
using Equate
using Statistics, DataFrames

include("testdata.jl")

ACTmathX = ExpandTable(ACTmath.scale, ACTmath.xcount)
ACTmathY = ExpandTable(ACTmath.scale, ACTmath.ycount)

@testset "ExpandTable" begin
    @test ExpandTable([0, 1, 2], [2, 2, 2]) == [0, 0, 1, 1, 2, 2]
    @test expandtable([0, 1, 2], [2, 2, 2]) == [0, 0, 1, 1, 2, 2]
end

@testset "DataMoments" begin
    @test round(mean(ACTmathX); digits = 4) == 19.8524
    @test round(std(ACTmathX; corrected = false); digits = 4) == 8.2116
    @test round(mean(ACTmathY); digits = 4) == 18.9798
    @test round(std(ACTmathY; corrected = false); digits = 4) == 8.9393
end

sgtabX = freqtab(ACTmathX; scale = 0:1:40)
sgtabY = freqtab(ACTmathY; scale = 0:1:40)

neattabX = freqtab(KBneatX.total, KBneatX.anchor; scaleX = 0:1:36, scaleV = 0:1:12)
neattabY = freqtab(KBneatY.total, KBneatY.anchor; scaleX = 0:1:36, scaleV = 0:1:12)

# Dummy table for NEAT equipercentile method
dummytabX = Equate.NEATFreqTab(
    DataFrame(scale = [0, 1, 2, 3, 4, 5], freq = [.1, .15, .25, .25, .15, .10] .* 100, prob = [.1, .15, .25, .25, .15, .10], cumprob = cumsum([.1, .15, .25, .25, .15, .10]) ),
    DataFrame(scale = [0, 1, 2, 3], freq = [.2, .4, .2, .2] .* 100, prob = [.2, .4, .2, .2], cumprob = cumsum([.2, .4, .2, .2])),
    ExpandTable([0, 1, 2, 3, 4, 5], Int64.([.1, .15, .25, .25, .15, .10] .* 100)),
    ExpandTable([0, 1, 2, 3], Int64.([.2, .4, .2, .2] .* 100)),
    1,
    1,
    [.04 .04 .02 .00
     .04 .08 .02 .01
     .06 .12 .05 .02
     .03 .12 .05 .05
     .02 .03 .04 .06
     .01 .01 .02 .06] .* 100,
    (;dummyX = "dummy"),
    (;dummyV = "dummy")
)

dummytabY = Equate.NEATFreqTab(
    DataFrame(scale = [0, 1, 2, 3, 4, 5], freq = [.08, .2, .22, .25, .15, .1] .* 100, prob = [.08, .2, .22, .25, .15, .1], cumprob = cumsum([.08, .2, .22, .25, .15, .1]) ),
    DataFrame(scale = [0, 1, 2, 3], freq = [.2, .2, .4, .2] .* 100, prob = [.2, .2, .4, .2], cumprob = cumsum([.2, .2, .4, .2])),
    ExpandTable([0, 1, 2, 3, 4, 5], Int64.([.08, .2, .22, .25, .15, .1] .* 100)),
    ExpandTable([0, 1, 2, 3], Int64.([.2, .2, .4, .2] .* 100)),
    1,
    1,
    [.04 .03 .01 .00
     .07 .05 .07 .01
     .03 .05 .12 .02
     .03 .04 .13 .05
     .02 .02 .05 .06
     .01 .01 .02 .06] .* 100,
    (;dummyX = "dummy"),
    (;dummyV = "dummy")
)


text:
using DelimitedFiles
using LinearAlgebra

function convolve_cyclic(signal::Array{T, 1},
                         filter::Array{T, 1}) where {T <: Number}

    # output size will be the size of sign
    output_size = max(length(signal), length(filter))

    # convolutional output
    out = Array{Float64,1}(undef,output_size)
    sum = 0

    for i = 1:output_size
        for j = 1:output_size
            if (mod1(i-j, output_size) <= length(filter))
                sum += signal[mod1(j, output_size)] * filter[mod1(i-j, output_size)]
            end
        end

        out[i] = sum
        sum = 0

    end

    return out
end

function convolve_linear(signal::Array{T, 1}, filter::Array{T, 1},
                         output_size) where {T <: Number}

    # convolutional output
    out = Array{Float64,1}(undef, output_size)
    sum = 0

    for i = 1:output_size
        for j = max(1, i-length(filter)):i
            if j <= length(signal) && i-j+1 <= length(filter)
                sum += signal[j] * filter[i-j+1]
            end
        end

        out[i] = sum
        sum = 0
    end

    return out
end

function main()

    # sawtooth functions for x and y
    x = [float(i)/200 for i = 1:200]
    y = [float(i)/200 for i = 1:200]

    # Normalization is not strictly necessary, but good practice
    normalize!(x)
    normalize!(y)

    # full convolution, output will be the size of x + y
    full_linear_output = convolve_linear(x, y, length(x) + length(y))

    # simple boundaries
    simple_linear_output = convolve_linear(x, y, length(x))

    # cyclic convolution
    cyclic_output = convolve_cyclic(x, y)

    # outputting convolutions to different files for plotting in external code
    writedlm("full_linear.dat", full_linear_output)
    writedlm("simple_linear.dat", simple_linear_output)
    writedlm("cyclic.dat", cyclic_output)

end


text:
# Nonlasing -> Lasing: take the mth nonlasing mode and use it to set the guess information
# about the mth lasing mode.
function Base.push!(lsol::LasingSol, m::Integer, nlsol::NonlasingSol, msgprefix::String="  ")
    printstyled(msgprefix * "Pushing nonlasing mode $m with ωₙₗ[$m] = $(nlsol.ω[m]) to lasing solution.\n", color=:blue)
    lsol.ω[m] = real(nlsol.ω[m])
    lsol.a²[m] = 0
    iₐ = lsol.iₐ[m] = nlsol.iₐ[m]
    ψ = lsol.ψ[m] .= nlsol.ψ[m]

    # Below, we don't test if ψ[iₐ] == 1, because the normalization of ψ was done by complex
    # division.  For a complex scalar z, z/z may not be exactly zero in floating-point
    # arithmetic, because z/w is basically evaluated as z*conj(w) / abs(w)^2, during which
    # rounding that could make z / z different from 1 occurs.
    @assert ψ[iₐ] ≈ 1  # make sure ψₘ is already normalized

    M = length(lsol)
    lsol.activated[m] = true
    lsol.active[m] = true
    lsol.m_active = (1:M)[lsol.active]

    return nothing
end

function Base.pop!(lsol::LasingSol, m::Integer, msgprefix::String="  ")
    printstyled(msgprefix * "Popping lasing mode $m from lasing solution.\n", color=:blue)
    lsol.a²[m] = 0  # indicate this mode is nonlasing
    lsol.ψ[m] .= 0  # good for compressing data when writing in file

    M = length(lsol)
    lsol.active[m] = false
    lsol.m_active = (1:M)[lsol.active]

    return nothing
end

# Lasing -> Nonlasing: take the mth lasing mode and use it to set the guess information
# about the mth nonlasing mode.
function Base.push!(nlsol::NonlasingSol, m::Integer, lsol::LasingSol, msgprefix::String="  ")
    printstyled(msgprefix * "Pushing lasing mode $m with ωₗ[$m] = $(lsol.ω[m]) to nonlasing solution.\n", color=:blue)
    nlsol.ω[m] = lsol.ω[m]  # real scalar
    iₐ = nlsol.iₐ[m] = lsol.iₐ[m]  # scalar
    ψ = nlsol.ψ[m] .= lsol.ψ[m]  # vector

    # Below, we don't test if ψ[iₐ] == 1, because the normalization of ψ was done by complex
    # division.  For a complex scalar z, z/z may not be exactly zero in floating-point
    # arithmetic, because z/w is basically evaluated as z*conj(w) / abs(w)^2, during which
    # rounding that could make z / z different from 1 occurs.
    @assert ψ[iₐ] ≈ 1  # make sure ψₘ is already normalized

    M = length(nlsol)
    nlsol.activated[m] = true
    nlsol.active[m] = true
    nlsol.m_active = (1:M)[nlsol.active]

    return nothing
end

function Base.pop!(nlsol::NonlasingSol, m::Integer, msgprefix::String="  ")
    printstyled(msgprefix * "Popping nonlasing mode $m from nonlasing solution.\n", color=:blue)
    nlsol.ω[m] = real(nlsol.ω[m])  # indicate this mode is lasing
    nlsol.ψ[m] .= 0  # good for compressing data when writing in file

    M = length(nlsol)
    nlsol.active[m] = false
    nlsol.m_active = (1:M)[nlsol.active]

    return nothing
end

# Prepare LasingSol by turning on the nonlasing mode with the largest positive imaginary
# part, if any.  See Sec. III.D of the SALT paper.
# Return the index of the mode to turn on; return 0 otherwise.
function turnon!(lsol::LasingSol, nlsol::NonlasingSol, msgprefix::String="  ")
    m = argmax(imag, nlsol.ω, nlsol.m_active)  # m==0 if nlsol.m_active is empty
    m≠0 && imag(nlsol.ω[m])>0 || return 0  # consider imag(ω) = 0 as nonlasing in order to keep lasing equations minimal

    # Now m ≠ 0 and imag(nlsol.ω[m]) > 0.
    # If the mode m was shut down just now, do not turn it on; otherwise, turn on the mode.
    printstyled(msgprefix * "Turning on mode $m where ωₙₗ[m=1:$(length(nlsol))] = [$(split(string(nlsol.ω), ('[',']'))[2])]...\n", color=:blue)  # 17 is to skip header "Complex{Float64}"
    if nlsol.activated[m]  # mode m is just shut down
        nlsol.ω[m] = real(nlsol.ω[m])
        printstyled(msgprefix * "Don't, because mode $m was shut down just now.\n", color=:red)
        return 0
    else
        push!(lsol, m, nlsol)
        pop!(nlsol, m)
        printstyled(msgprefix * "Done.\n", color=:blue)
        return m
    end
end

# Prepare NonlasingSol by shutting down the lasing mode with negative amplitude, if any.
# Return the index of the mode to shut down; return 0 otherwise.
function shutdown!(lsol::LasingSol, nlsol::NonlasingSol, msgprefix::String="  ")
    m = argmin(identity, lsol.a², lsol.m_active)  # m==0 if lsol.m_active is empty
    m≠0 && lsol.a²[m]≤0 || return 0  # consider a² = 0 as nonlasing in order to keep lasing equations minimal

    # Now m ≠ 0 and lsol.a²[m] ≤ 0.
    # If the mode m was turned on just now, do not shut it down; otherwise, shut down the mode.
    printstyled(msgprefix * "Shutting down mode $m where aₗ²[m=1:$(length(lsol))] = $(lsol.a²)...\n", color=:blue)
    if lsol.activated[m]  # mode m is just turned on
        lsol.a²[m] = 0
        printstyled(msgprefix * "Don't, because mode $m was turned on just now.\n", color=:red)
        return 0
    else
        push!(nlsol, m, lsol)
        pop!(lsol, m)
        printstyled(msgprefix * "Done.\n", color=:blue)
        return m
    end
end

# Once switching has occurred, lasing and nonlasing modes are re-calculated, and then they
# should pass this test.
function check_conflict(lsol, nlsol)
    M = length(lsol)
    length(nlsol) == M ||
        @error "length(lsol) = $M and length(nlsol) = $(length(nlsol)) should be same."

    for m = 1:M
        xor(lsol.active[m], nlsol.active[m]) ||
            @error "Mode $m should be either lasing or nonlasing, exclusively."
    end

    return nothing
end

# Return true if the mode m is simultaneously lasing and nonlasing (i.e., at threshold).
#
# Even in such a case, the mode m is indicated only lasing or nonlasing, not both. (In other
# words, only one of lsol.active[m] and nlsol.active[m] is true.)  For example, if the mode
# m is turned on but its lasing amplitude is still zero, this mode is physically at
# threshold and therefore simultaneously lasing and nonlasing.  However, we have
# lsol.active[m] == true and nlsol.active[m] == false, because technically it is turned on.
# Similarly, if the mode is shut down but its eigenfrequency has a zero imaginary part, this
# mode is again physically at threshold and therefor simultaneously lasing and nonlasing,
# but we have lsol.active[m] == false and nlsol.active[m] == true because technically this
# mode is shut down.
islnl(lsol, nlsol, m) = lsol.a²[m]==0 && imag(nlsol.ω[m])==0


text:
#=
function Base.append!(t::MetidaTable, t2::MetidaTable)
    if !(names(t) ⊆ names(t2)) error("Names for t not in t2") end
    for n in names(t)
        append!(t.table[n], t2.table[n])
    end
    t
end
=#
function MetidaBase.metida_table_(obj::DataSet{T}) where T <: PKSubject
    idset  = Set(keys(first(obj).id))
    if length(obj) > 1
        for i = 2:length(obj)
            union!(idset,  Set(keys(obj[i].id)))
        end
    end
    mt1 = metida_table_((fill(getid(obj, 1, c), length(obj[1])) for c in idset)...; names = idset)
    mt2 = metida_table_(deepcopy(obj[1].time), deepcopy(obj[1].obs); names = [:time, :obs])
    mtm = merge(mt1, mt2)
    if length(obj) > 1
        for i = 2:length(obj)
            mt1 = metida_table_((fill(getid(obj, i, c), length(obj[i])) for c in idset)...; names = idset)
            mt2 = metida_table_(obj[i].time, obj[i].obs; names = [:time, :obs])
            amtm = merge(mt1, mt2)
            for n in keys(mtm)
                append!(mtm[n], amtm[n])
            end
        end
    end
    mtm
end


text:
function DriftModifier(filter::Function, drift::Function)
    function f(data, event)
        elements = filter(data)

        if length(elements) > 0
            df = @view data[elements, :]
            drift(df, event)
        end
        return nothing
    end

    AlterDataModifier(f)
end

sigmoid(x; c1 = 1.0, c2 = 0.0) = 1 / (1 + ℯ ^ (-c1 * (x - c2)))

function ClassDriftModifier(column::Symbol, value::T, drift::Function) where T <: Number
    return DriftModifier((data) -> findall(r-> r == value, data[:, column]), drift)
end

function IncrementalDriftModifier(vetor::Dict{Symbol, T}, filter::Function; c1 = 1.0, c2 = 0.0)::Modifiers where T <: Number
    modifiers = EasyStream.Modifier[]
    for (column, value) in vetor
        drift = DriftModifier(filter, (data, event) ->
                                        data[:, column] = data[:, column] .+ value .* sigmoid(event.time; c1 = c1, c2 = c2))
        push!(modifiers, drift)
    end

    return Modifiers(modifiers)
end

function SuddenDriftModifier(vetor::Dict{Symbol, T}, filter::Function; c2 = 0.0)::Modifiers where T <: Number
    return IncrementalDriftModifier(vetor, filter; c1 = 1.0, c2 = c2)
end


text:
module TestSort

using Random: shuffle
using Test
using ThreadsX

@testset for basesize in 1:8
    @testset for alg in [ThreadsX.MergeSort, ThreadsX.QuickSort, ThreadsX.StableQuickSort]
        @test ThreadsX.sort!(shuffle([1:1000;]); alg = alg, basesize = basesize) == 1:1000
        @test ThreadsX.sort!(
            shuffle([1:1000;]);
            alg = alg,
            basesize = basesize,
            by = inv,
        ) == 1000:-1:1
        @test ThreadsX.sort!([1:1000;]; alg = alg, basesize = basesize) == 1:1000
    end
end

@testset "stable sort" begin
    @testset for alg in [ThreadsX.MergeSort, ThreadsX.StableQuickSort]
        @test ThreadsX.sort(1:45; alg = alg, basesize = 25, by = _ -> 1) == 1:45
        @test ThreadsX.sort(1:1000; alg = alg, basesize = 200, by = _ -> 1) == 1:1000
    end
end

@testset "inplace" begin
    @testset "default" begin
        xs = [1:10;]
        @test ThreadsX.sort!(xs) === xs
        @test issorted(xs)
    end
    @testset for alg in [
        ThreadsX.MergeSort,
        ThreadsX.QuickSort,
        ThreadsX.StableQuickSort,
        MergeSort,
        QuickSort,
    ]
        xs = [1:10;]
        @test ThreadsX.sort!(xs; alg = alg) === xs
        @test issorted(xs)
    end
    @testset for alg in [ThreadsX.MergeSort, ThreadsX.QuickSort, ThreadsX.StableQuickSort]
        xs = [1:10;]
        @test sort!(xs; alg = alg) === xs
        @test issorted(xs)
    end
end

randnans(n) = reinterpret(Float64, [rand(UInt64) | 0x7ff8000000000000 for i in 1:n])

function randn_with_nans(n, p)
    v = randn(n)
    x = findall(rand(n) .< p)
    v[x] = randnans(length(x))
    return v
end

# Taken from "advanced sorting" testset.  This was the only test that
# failed when `fpsort!` was replaced with plain `sort!`.
@testset "NaN" begin
    @testset for n in [0:10; 100; 101; 1000; 1001]
        v = randn_with_nans(n, 0.1)
        @testset for alg in [
                ThreadsX.MergeSort,
                ThreadsX.QuickSort,
                ThreadsX.StableQuickSort,
            ],
            rev in [false, true],
            basesize in 1:8
            # test float sorting with NaNs
            s = ThreadsX.sort(v, alg = alg, rev = rev, basesize = basesize)
            @test issorted(s, rev = rev)
            @test reinterpret(UInt64, v[isnan.(v)]) == reinterpret(UInt64, s[isnan.(s)])
        end
    end
end

@testset "UI" begin
    @test ThreadsX.sort(1:10) == 1:10
    @test ThreadsX.sort(1:10; alg = MergeSort) == 1:10
    @test ThreadsX.sort(1:10; alg = QuickSort) == 1:10
end

end  # module


text:
function EWS(stochasticprogram::StochasticProgram,
             structure::VerticalStructure{2,1,Tuple{SP}}) where SP <: DistributedScenarioProblems
    partial_ews = Vector{Float64}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_ews[i] = remotecall_fetch(
                w,
                scenarioproblems(stochasticprogram)[w-1],
                stochasticprogram.generator[:stage_1],
                stochasticprogram.generator[:stage_2],
                stage_parameters(stochasticprogram, 1),
                stage_parameters(stochasticprogram, 2),
                sub_optimizer(stochasticprogram)) do (sp,gen_one,gen_two,one_params,two_params,opt)
                    scenarioproblems = fetch(sp)
                    num_scenarios(scenarioproblems) == 0 && return zero(T)
                    subproblems = map(scenarios(scenarioproblems)) do scenario
                        ws = _WS(gen_one,
                                 gen_two,
                                 one_params,
                                 two_params,
                                 scenario,
                                 Decisions(),
                                 opt)
                        return ws
                    end
                    return outcome_mean(subproblems, probability.(subproblems))
                end
        end
    end
    return sum(partial_ews)
end

function statistical_EWS(stochasticprogram::StochasticProgram,
                         structure::VerticalStructure{2,1,Tuple{SP}}) where SP <: DistributedScenarioProblems
    partial_welfords = Vector{Tuple{Float64,Float64,Float64,Int}}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_ews[i] = remotecall_fetch(
                w,
                scenarioproblems(stochasticprogram)[w-1],
                stochasticprogram.generator[:stage_1],
                stochasticprogram.generator[:stage_2],
                stage_parameters(stochasticprogram, 1),
                stage_parameters(stochasticprogram, 2),
                sub_optimizer(stochasticprogram)) do (sp,gen_one,gen_two,one_params,two_params,opt)
                    scenarioproblems = fetch(sp)
                    num_scenarios(scenarioproblems) == 0 && return 0.0, 0.0, 0.0, 0
                    ws_models = map(scenarios(scenarioproblems)) do scenario
                        ws = _WS(gen_one,
                                 gen_two,
                                 one_params,
                                 two_params,
                                 scenario,
                                 Decisions(),
                                 opt)
                        return ws
                    end
                    return welford(ws_models, probability.(scenarios(scenarioproblems)))
                end
        end
    end
    𝔼WS, σ², _ = reduce(aggregate_welford, partial_welfords)
    return 𝔼WS, sqrt(σ²)
end


text:
using Fjage

export UnetSource

const RxBasebandSignalNtf = MessageClass(@__MODULE__, "org.arl.unet.bb.RxBasebandSignalNtf")

struct UnetSource
  gw::Gateway
  bb::AgentID
  fs::Float32
  blksize::Int
  nchannels::Int
  streaming::Ref{Bool}
end

function UnetSource(host="127.0.0.1", port=1100; blksize=4096, bb=nothing)
  Fjage.registermessages()
  gw = Gateway(host, port)
  try
    if bb === nothing
      bb = agentforservice(gw, "org.arl.unet.Services.BASEBAND")
    else
      bb = agent(gw, bb)
    end
    subscribe(gw, topic(bb))
    blksize === nothing || (bb.pbsblk = blksize)
    blksize = bb.pbsblk
    bb.pbscnt = 1
    ntf = receive(gw, RxBasebandSignalNtf, 5000)
    if ntf === nothing
      bb.pbscnt = 0
      throw(ErrorException("No data from host"))
    end
    UnetSource(gw, bb, Float32(ntf.fs), blksize, ntf.channels, Ref(false))
  catch ex
    close(gw)
    rethrow(ex)
  end
end

connect(::UnetSource) = true

function Base.close(src::UnetSource)
  stopinputstream(src)
  close(src.gw)
end

inputframerate(src::UnetSource) = src.fs
inputblocksize(src::UnetSource) = src.blksize
inputchannels(src::UnetSource) = src.nchannels

function startinputstream(src::UnetSource, callback)
  src.streaming[] && return
  src.bb.pbscnt = -1
  src.streaming[] = true
  @async begin
    t0 = time()
    while src.streaming[]
      ntf = receive(src.gw, RxBasebandSignalNtf, 1000)
      if ntf !== nothing && iszero(ntf.fc)
        t = time() - t0
        try
          data = transpose(reshape(ntf.signal, src.nchannels, src.blksize))
          callback(round(Int, 1000000t), data)
        catch ex
          @error "$ex"
        end
      end
    end
  end
end

function stopinputstream(src::UnetSource)
  try
    req = ParameterReq(recipient=src.bb)
    set!(req, "pbscnt", 0)
    send(src.gw, req)
  catch ex
    @error "$ex"
  end
  src.streaming[] = false
end


text:
sz = (4,5)
filename = tempname()
#filename = "/tmp/test-6.nc"
# The mode "c" stands for creating a new file (clobber)

NCDatasets.Dataset(filename,"c") do ds

    # define the dimension "lon" and "lat"
    #NCDatasets.defDim(ds,"lon",sz[1])
    #NCDatasets.defDim(ds,"lat",Inf)

    ds.dim["lon"] = sz[1]
    ds.dim["lat"] = Inf

    # variables
    for T in [UInt8,Int8,UInt16,Int16,UInt32,Int32,UInt64,Int64,Float32,Float64]
    #for T in [Float32]
        local data
        data = zeros(T,sz)

        v = NCDatasets.defVar(ds,"var-$T",T,("lon","lat"))

        for j = 1:sz[2]
            data[:,j] .= T(j)
            v[:,j] = T(j)
        end

        @test all(v[:,:] == data)

    end
end


text:
module Controllers

using TrajectoryOptimization
using TrajectoryOptimization.Dynamics
const TO = TrajectoryOptimization

using StaticArrays
using LinearAlgebra
using ForwardDiff

include("../controllers/rbstate.jl")
include("../controllers/tracking_control.jl")

end


text:
using .Flux

import Base: convert
import IncrementalInference: getSample

export MixtureFluxPose2Pose2, PackedMixtureFluxPose2Pose2

struct MixtureFluxPose2Pose2{F <: FunctorInferenceType} <: AbstractRelativeRoots
  Zij::F
  # delta time between variables
  DT::Ref{Float64}
end

mutable struct PackedMixtureFluxPose2Pose2 <: PackedInferenceType
  packedZij::PackedMixture
  DT::Float64
  specialSampler::String
end

function IIF.getSample(cfo::CalcFactor{<:MixtureFluxPose2Pose2}, N::Int=1)
  #
  nfb = cfo.factor
  fmd = cfo.metadata

  # only calculate once after default value of 0
  if nfb.DT[] == 0
    # cache the time difference estimate
    nfb.DT[] = (getTimestamp(fmd.fullvariables[2]) - getTimestamp(fmd.fullvariables[1])).value * 1e-3
  end

  cf_ = CalcFactor( cfo.factor.Zij, cfo.metadata, 0, length(cfo._legacyMeas), cfo._legacyMeas, cfo._legacyParams)

  smpl = getSample(cf_, N)[1]

  return (smpl, cf_.factor.labels)

end

function MixtureFluxPose2Pose2( fluxmodels::AbstractVector,
                                inputDims::NTuple,
                                data::AbstractArray,
                                outputDims::NTuple,
                                otherComponents::AbstractVector{<:SamplableBelief}=[MvNormal(zeros(3),diagm(ones(3)));],
                                diversity=[0.5;0.5];
                                DT::Real=0,
                                shuffle::Bool=true,
                                serializeHollow::Bool=false )
  #
  mix = MixtureFluxModels(Pose2Pose2,
                          fluxmodels,
                          inputDims,
                          data,
                          outputDims,
                          otherComponents,
                          diversity,
                          shuffle=shuffle,
                          serializeHollow=serializeHollow)
  #
  MixtureFluxPose2Pose2(mix, Base.RefValue{Float64}(DT))
end

function Base.convert(::Union{Type{<:PackedInferenceType},Type{<:PackedMixtureFluxPose2Pose2}},
                      obj::MixtureFluxPose2Pose2 )
  #
  toFnc = typeof(obj.specialSampler)
  fncName = string(toFnc.name.module)*"."*string(toFnc.name.name)
  PackedMixtureFluxPose2Pose2(convert(PackedMixture,obj.Zij),
                              obj.DT[],
                              fncName)
end

function Base.convert(::Union{Type{<:FunctorInferenceType},Type{<:MixtureFluxPose2Pose2}},
                      obj::PackedMixtureFluxPose2Pose2 )
  #
  sFnc = split(obj.specialSampler, '.') .|> Symbol
  specFnc = getfield(getfield(Main, sFnc[1]), sFnc[2])
  PackedMixtureFluxPose2Pose2(convert(Mixture,obj.packedZij),
                                      obj.DT,
                                      specFnc)
end

function calcVelocityInterPose2!( nfb::MixtureFluxPose2Pose2,
                                  iPts::AbstractVector{<:Real},
                                  jPts::AbstractVector{<:Real} )
  #
  joyVelData = nfb.Zij.components.fluxnn.data
  # DXY[1:2,i] .= TransformUtils.R(iPts[3,i])'*DXY[1:2,i]
  joyVelData[1,3:4] .= jPts[1:2]
  joyVelData[1,3:4] .-= iPts[1:2]
  joyVelData[1,3:4] ./= nfb.DT[] # convert to velocity
  joyVelData[1,3:4] .= TransformUtils.R(iPts[3])'*joyVelData[1,3:4]
  # just set zero if something is wrong
  if isnan(joyVelData[1,3]) || isinf(abs(joyVelData[1,3])) || isnan(joyVelData[1,4]) || isinf(abs(joyVelData[1,4]))
    joyVelData[1,3:4] .= 0.0
  end
  for i in 2:size(joyVelData,1)
    joyVelData[i,3:4] .= joyVelData[1,3:4]
  end
  return nothing
end

# function calcVelocityInterPose2!( nfb::MixtureFluxPose2Pose2,
#                                   iPts::AbstractMatrix{<:Real},
#                                   jPts::AbstractMatrix{<:Real}  )
#   #
#   @assert size(jPts,2) == size(iPts,2) "calcVelocityInterPose2! can currently only evaluate equal population size variables"

#   # calculate an average velocity component
#   if nfb.DT[] == 0
#     # DXY = (@view jPts[1:2,:]) - (@view jPts[1:2,:])
#     # rotate delta position from world to local iX frame
#     for i in 1:size(iPts,2)
#       calcVelocityInterPose2!(nfb, iPts, jPts, i)
#     end
#   else
#     nfb.joyVelData[:,3:4] .= 0.0
#   end
#   nothing
# end

function (cfo::CalcFactor{<:MixtureFluxPose2Pose2})(meas1, meas2, Xi, Xj)
  #
  userdata = cfo.metadata
  nfb = cfo.factor
  fmd = cfo.metadata

  # if, use prediction sample
  if meas2 == 2
    # get live velocity estimate for each sample (nfb.joyVelData[:,3:4])
    calcVelocityInterPose2!(nfb, Xi, Xj)
    # predict odom for this sample from a specific prediction model
    # NOTE, already integrated as part of Mixture
    # meas[1][1:2,idx] = nfb.allPredModels[meas[3][idx]](nfb.joyVelData)
  end

  # calculate the error for that measurement sample as Pose2Pose2
  #TODO
  cfZij = CalcFactor( cfo.factor.Zij.mechanics, cfo.metadata, 0, length(cfo._legacyMeas), cfo._legacyMeas, cfo._legacyParams)

  return  cfZij(meas1, Xi, Xj)

end

## ============================================================================================
## Legacy support:  FluxModelsPose2Pose2
## ============================================================================================

export FluxModelsPose2Pose2

# Convenience function to help call the right constuctor
FluxModelsPose2Pose2( allModels::Vector{P},
                      jvd::D,
                      naiveModel::M,
                      naiveFrac::Real=0.5,
                      DT::Real=0.0,
                      shuffle::Bool=true ) where {P, M <: SamplableBelief, D <: AbstractMatrix} = MixtureFluxPose2Pose2(
                                        allModels,
                                        (25,4),
                                        jvd,
                                        (3,),
                                        [naiveModel;],
                                        [1-naiveFrac;naiveFrac],
                                        DT=DT,
                                        shuffle=shuffle )
#

# struct MixtureFlux{N,F<:FunctorInferenceType,S,T}
#   mixture::Mixture{N,F,S,T}
#   # special keyword field name used to invoke 'specialSampler' logic
#   specialSampler::Function
# end

# # to be overrided by user via more specialized dispatch
# # FIXME as part of
# function sampleMixtureFlux( nfb::MixtureFluxPose2Pose2,
#                             N::Int,
#                             fmd::FactorMetadata,
#                             legacy... )
#   #
#   getSample(nfb.mixture, N)
# end

# # const _IIFListTypes = Union{<:AbstractVector, <:Tuple, <:NTuple}

# function MixtureFlux( F_::FunctorInferenceType,
#                       compList::_IIFListTypes,
#                       diversity::Union{<:AbstractVector, <:NTuple, <:DiscreteNonParametric})
#   #
#   mix = Mixture(F_, compList, diversity)
#   _populateMixture(mix{N,F,S,T}) where {N,F,S,T} = MixtureFlux{N,F,S,T}(mix, sampleMixtureFlux)
#   return _populateMixture(mix)
# end

# MixtureFlux(::Type{F}, w...;kw...) where F <: FunctorInferenceType = MixtureFlux(F(LA.I), w...;kw...)

#