text:
#=
Operations on TS objects
=#

# TODO: increase efficiency running these operations
import Base: ones, zeros, trues, falses, isnan, sum, mean, maximum, minimum,
prod, cumsum, cumprod, diff, all, any, countnz, sign, find, findfirst
importall Base.Operators

# titlecase(str::String) = join([ucfirst(word) for word in split(str, ' ')], ' ')
# titlecase(arr::AbstractArray{String}) = map(titlecase, arr)
islogical(fun::Function) = fun in (<, .<, <=, .<=, .>, >, >=, .>=, ==, .==, !=, .!=, !)
isarithmetic(fun::Function) = fun in (+, .+, -, .-, *, .*, /, ./, %, .%, ^, .^, \, .\)
const logicals = Dict("<" => "LessThan", ".<" => "LessThan", "<=" => "LessThanEqual", ".<=" => "LessThanEqual",
                      ">" => "GreaterThan", ".>" => "GreaterThan", ">=" => "GreaterThanEqual", ".>=" => "GreaterThanEqual",
                      "==" => "Equal", ".==" => "Equal", "!=" => "NotEqual", ".!=" => "NotEqual", "!" => "Not")
const arithmetics = Dict("+" => "Plus", ".+" => "Plus", "-" => "Minus", ".-" => "Minus",
                         "*" => "Times", ".*" => "Times", "/" => "Over", "./" => "Over",
                         "%" => "Mod", ".%" => "Mod", "^" => "Power", ".^" => "Power", "\\" => "Inverse", ".\\" => "Under")

find(x::TS) = find(x.values)
findfirst(x::TS) = findfirst(x.values)
ones(x::TS) = ts(ones(x.values), x.index, x.fields)
zeros(x::TS) = ts(zeros(x.values), x.index, x.fields)
trues(x::TS) = ts(trues(x.values), x.index, x.fields)
falses(x::TS) = ts(falses(x.values), x.index, x.fields)
isnan(x::TS) = ts(isnan(x.values), x.index, x.fields)
countnz(x::TS) = countnz(x.values)
sign(x::TS) = ts(sign(x.values), x.index, x.fields)

# Pass Array operators through to underlying TS values
function  op{V,T}(x::TS{V,T}, y::TS{V,T}, fun::Function; args...)
    if islogical(fun)
        fldstr = "$(string(x.fields[1]))$(logicals[string(fun)])$(string(y.fields[1]))"
    elseif isarithmetic(fun)
        fldstr = "$(string(x.fields[1]))$(arithmetics[string(fun)])$(string(y.fields[1]))"
    else
        fldstr = ucfirst(string(fun))
    end
    fldsym = Symbol(fldstr)
    if x.index == y.index
        return ts(fun(x.values, y.values; args...), x.index, fldsym)
    end
    idx = intersect(x.index, y.index)
    return ts(fun(x[idx].values, y[idx].values; args...), idx, fldsym)
end
op{V,T}(x::TS{V,T}, fun::Function; args...) = ts(fun(x.values), x.index, x.fields; args...)

# Number functions
sum(x::TS) = sum(x.values)
sum(x::TS, dim::Int) = sum(x.values, dim)
sum(f::Function, x::TS) = sum(f, x.values)
mean(x::TS) = mean(x.values)
mean(x::TS, dim::Int) = mean(x, dim)
prod(x::TS) = prod(x.values)
prod(x::TS, dim::Int) = prod(x.values, dim)
maximum(x::TS) = maximum(x.values)
maximum(x::TS, dim::Int) = maximum(x.values, dim)
minimum(x::TS) = minimum(x.values)
minimum(x::TS, dim::Int) = minimum(x.values, dim)
cumsum(x::TS, dim::Int=1) = ts(cumsum(x.values, dim), x.index, x.fields)
cummin(x::TS, dim::Int=1) = ts(cummin(x.values, dim), x.index, x.fields)
cummax(x::TS, dim::Int=1) = ts(cummax(x.values, dim), x.index, x.fields)
cumprod(x::TS, dim::Int=1) = ts(cumprod(x.values, dim), x.index, x.fields)

nans(r::Int, c::Int) = fill(NaN, 1, 2)
nans(dims::Tuple{Int,Int}) = fill(NaN, dims)
function rowdx!{T,N}(dx::AbstractArray{T,N}, x::AbstractArray{T,N}, n::Int, r::Int=size(x,1))
    idx = n > 0 ? (n+1:r) : (1:r+n)
    @inbounds for i=idx
        dx[i,:] = x[i,:] - x[i-n,:]
    end
    nothing
end
function coldx!{T,N}(dx::AbstractArray{T,N}, x::AbstractArray{T,N}, n::Int, c::Int=size(x,2))
    idx = n > 0 ? (n+1:c) : (1:c+n)
    @inbounds for j=idx
        dx[:,j] = x[:,j] - x[:,j-n]
    end
    nothing
end
function diffn{T<:Number,N}(x::Array{T,N}, dim::Int=1, n::Int=1)
    @assert dim == 1 || dim == 2 "Argument `dim` must be 1 (rows) or 2 (columns)."
    @assert abs(n) < size(x,dim) "Argument `n` out of bounds."
    if n == 0
        return x
    end
    dx = zeros(x)
    if dim == 1
        rowdx!(dx, x, n)
    else
        coldx!(dx, x, n)
    end
    return dx
end
function diff{V,T}(x::TS{V,T}, n::Int=1; dim::Int=1, pad::Bool=true, padval::V=zero(V))
    @assert dim == 1 || dim == 2 "Argument dim must be either 1 (rows) or 2 (columns)."
    r = size(x, 1)
    c = size(x, 2)
    dx = diffn(x.values, dim, n)
    if dim == 1
        if pad
            idx = n>0 ? (1:n) : (r+n+1:r)
            dx[idx,:] = padval
            return ts(dx, x.index, x.fields)
        else
            idx = n > 0 ? (n+1:r) : (1:r+n)
            return ts(dx[idx,:], x.index[idx], x.fields)
        end
    else
        if pad
            idx = n > 0 ? (1:c) : (c+1+1:c)
            dx[:,idx] = padval
            return ts(dx, x.index, x.fields[idx])
        else
            idx = n > 0 ? (n+1:c) : (1:c+n)
            return ts(dx[:,idx], x.index, x.fields[idx])
        end
    end
end
function lag{V,T}(x::TS{V,T}, n::Int=1; pad::Bool=true, padval::V=zero(V))
	@assert abs(n) < size(x,1) "Argument `n` out of bounds."
	if n == 0
		return x
	elseif n > 0
		out = zeros(x.values)
		out[n+1:end,:] = x.values[1:end-n,:]
	elseif n < 0
		out = zeros(x.values)
		out[1:end+n,:] = x.values[1-n:end,:]
	end
    r = size(x, 1)
    c = size(x, 2)
    if pad
        idx = n>0 ? (1:n) : (r+n+1:r)
        out[idx,:] = padval
        return ts(out, x.index, x.fields)
    else
        idx = n > 0 ? (n+1:r) : (1:r+n)
        return ts(out[idx,:], x.index[idx], x.fields)
    end
end

# Artithmetic operators
+(x::TS) = ts(+x.values, x.index, x.fields)
-(x::TS) = ts(-x.values, x.index, x.fields)
+(x::TS, y::TS) = op(x, y, +)
-(x::TS, y::TS) = op(x, y, -)
*(x::TS, y::TS) = op(x, y, *)
/(x::TS, y::TS) = op(x, y, /)
.+(x::TS, y::TS) = op(x, y, .+)
.-(x::TS, y::TS) = op(x, y, .-)
.*(x::TS, y::TS) = op(x, y, .*)
./(x::TS, y::TS) = op(x, y, ./)
.\(x::TS, y::TS) = op(x, y, .\)
.^(x::TS, y::TS) = op(x, y, .^)
.%(x::TS, y::TS) = op(x, y, .%)

+(x::TS, y::AbstractArray) = ts(x.values + y, x.index, x.fields)
-(x::TS, y::AbstractArray) = ts(x.values - y, x.index, x.fields)
*(x::TS, y::AbstractArray) = ts(x.values * y, x.index, x.fields)
/(x::TS, y::AbstractArray) = ts(x.values / y, x.index, x.fields)
.+(x::TS, y::AbstractArray) = ts(x.values .+ y, x.index, x.fields)
.-(x::TS, y::AbstractArray) = ts(x.values .- y, x.index, x.fields)
.*(x::TS, y::AbstractArray) = ts(x.values .* y, x.index, x.fields)
./(x::TS, y::AbstractArray) = ts(x.values ./ y, x.index, x.fields)
.\(x::TS, y::AbstractArray) = ts(x.values .\ y, x.index, x.fields)
.^(x::TS, y::AbstractArray) = ts(x.values .^ y, x.index, x.fields)
.%(x::TS, y::AbstractArray) = ts(x.values .% y, x.index, x.fields)
+(y::AbstractArray, x::TS) = x + y
-(y::AbstractArray, x::TS) = x - y
*(y::AbstractArray, x::TS) = x * y
.+(y::AbstractArray, x::TS) = x .+ y
/(y::AbstractArray, x::TS) = x / y
.-(y::AbstractArray, x::TS) = x .- y
.*(y::AbstractArray, x::TS) = x .* y
./(y::AbstractArray, x::TS) = x ./ y
.%(y::AbstractArray, x::TS) = x .% y

+(x::TS, y::Number) = ts(x.values + y, x.index, x.fields)
-(x::TS, y::Number) = ts(x.values - y, x.index, x.fields)
*(x::TS, y::Number) = ts(x.values * y, x.index, x.fields)
/(x::TS, y::Number) = ts(x.values / y, x.index, x.fields)
.+(x::TS, y::Number) = ts(x.values .+ y, x.index, x.fields)
.-(x::TS, y::Number) = ts(x.values .- y, x.index, x.fields)
.*(x::TS, y::Number) = ts(x.values .* y, x.index, x.fields)
./(x::TS, y::Number) = ts(x.values ./ y, x.index, x.fields)
.\(x::TS, y::Number) = ts(x.values .\ y, x.index, x.fields)
.^(x::TS, y::Number) = ts(x.values .^ y, x.index, x.fields)
.%(x::TS, y::Number) = ts(x.values .% y, x.index, x.fields)
+(y::Number, x::TS) = x + y
-(y::Number, x::TS) = x - y
*(y::Number, x::TS) = x * y
/(y::Number, x::TS) = x / y
.+(y::Number, x::TS) = x .+ y
.-(y::Number, x::TS) = x .- y
.*(y::Number, x::TS) = x .* y
./(y::Number, x::TS) = x ./ y
.%(y::Number, x::TS) = x .% y

# Logical operators
all(x::TS) = all(x.values)
any(x::TS) = any(x.values)
!(x::TS) = ts(!x.values, x.index, x.fields)
==(x::TS, y::TS) = x.values == y.values && x.index == y.index
!=(x::TS, y::TS) = !(x == y)

.==(x::TS, y::TS) = op(x, y, .==)
.>(x::TS, y::TS) = op(x, y, .>)
.<(x::TS, y::TS) = op(x, y, .<)
.!=(x::TS, y::TS) = op(x, y, .!=)
.<=(x::TS, y::TS) = op(x, y, .<=)
.>=(x::TS, y::TS) = op(x, y, .>=)

.==(x::TS, y::Number) = ts(x.values .== y, x.index, x.fields)
.>(x::TS, y::Number) = ts(x.values .> y, x.index, x.fields)
.<(x::TS, y::Number) = ts(x.values .< y, x.index, x.fields)
.!=(x::TS, y::Number) = ts(x.values .!= y, x.index, x.fields)
.<=(x::TS, y::Number) = ts(x.values .<= y, x.index, x.fields)
.>=(x::TS, y::Number) = ts(x.values .>= y, x.index, x.fields)
.==(y::Number, x::TS) = ts(x.values .== y, x.index, x.fields)
.>(y::Number, x::TS) = ts(x.values .> y, x.index, x.fields)
.<(y::Number, x::TS) = ts(x.values .< y, x.index, x.fields)
.!=(y::Number, x::TS) = ts(x.values .!= y, x.index, x.fields)
.<=(y::Number, x::TS) = ts(x.values .<= y, x.index, x.fields)
.>=(y::Number, x::TS) = ts(x.values .>= y, x.index, x.fields)

.==(x::TS, y::AbstractArray) = ts(x.values .== y, x.index, x.fields)
.>(x::TS, y::AbstractArray) = ts(x.values .> y, x.index, x.fields)
.<(x::TS, y::AbstractArray) = ts(x.values .< y, x.index, x.fields)
.!=(x::TS, y::AbstractArray) = ts(x.values .!= y, x.index, x.fields)
.<=(x::TS, y::AbstractArray) = ts(x.values .<= y, x.index, x.fields)
.>=(x::TS, y::AbstractArray) = ts(x.values .>= y, x.index, x.fields)
.==(y::AbstractArray, x::TS) = ts(x.values .== y, x.index, x.fields)
.>(y::AbstractArray, x::TS) = ts(x.values .> y, x.index, x.fields)
.<(y::AbstractArray, x::TS) = ts(x.values .< y, x.index, x.fields)
.!=(y::AbstractArray, x::TS) = ts(x.values .!= y, x.index, x.fields)
.<=(y::AbstractArray, x::TS) = ts(x.values .<= y, x.index, x.fields)
.>=(y::AbstractArray, x::TS) = ts(x.values .>= y, x.index, x.fields)

>(x::TS, y::TS) = x .> y
<(x::TS, y::TS) = x .< y
<=(x::TS, y::TS) = x .<= y
>=(x::TS, y::TS) = x .>= y


text:
# Test suite for the ALU, saturator, MAC, and branch evaluator.

mutable struct TestBundle
    src1 :: Vector{Int16}
    src2 :: Vector{Int16}
    cin :: Vector{Bool}
end
TestBundle() = TestBundle(Int16[], Int16[], Bool[])

# @testset "ALU Tests" begin
#     # Set this flag to benchmark ALU Operations
#     benchmark = false
#
#     # This is going to be painful, but necessary. Have to check all of the ALU
#     # type ops
#
#     # Need to set up some preliminary data structures that will be used for
#     # all of the tests.
#
#     # Alias some commands for more convenient testing.
#     CondExec        = AsapSim.CondExec
#     PipelineEntry   = AsapSim.PipelineEntry
#     ALUFlags        = AsapSim.ALUFlags
#     AsapInstruction = AsapSim.AsapInstructionKeyword
#
#     # 2-long vector for conditional execution.
#     condexec = [CondExec(), CondExec()]
#
#     # Template Pipeline Stage.
#     stage = PipelineEntry()
#
#     # Functions for performing various tests
#     function op_tester(
#             ::Type{T},
#             optype,
#             op :: Symbol,
#             test :: TestBundle;
#             benchmark = false
#         ) where T
#         # Instantiate new ALU Flags
#
#         for (src1v, src2v, cin) in zip(test.src1, test.src2, test.cin)
#             # Set the src1_value and src2_value fields
#             # First, convert the source values to 16 bit unsigned numbers
#             # (will error if the conversion is inexact)
#             #
#             # Then reinterpret to Int16 since that is the desired type of the
#             # pipe stage fields.
#             #stage = AsapSim.set(stage, AsapSim.SRC(src1v, src2v))
#
#             # Edit the carry flag.
#             aluflags = ALUFlags(cin, false, false, false)
#             if benchmark
#                 b = @benchmark AsapSim.stage4_alu($op, $src1v, $src2v, $aluflags, $condexec)
#                 println("Average time for $(op): $(mean(b.times))")
#             end
#             result, newflags = AsapSim.stage4_alu(op, src1v, src2v, aluflags, condexec)
#
#             # Next the "op" twice so it works correctly if op is "+" or "-".
#             if T == Unsigned
#                 # Do 2's complement arithmetic on the operands, then convert
#                 # everything to "Unsigned" for prettier printing.
#                 expected_result = reinterpret(UInt16, optype(optype(src1v, src2v), Int16(cin)))
#                 @test reinterpret(UInt16, result) == expected_result
#
#                 # Convert the input operands to UInt64. Compute whether there
#                 # should be a carry out and test for it.
#                 #
#                 # Note that since we converted to UInt64, if subtraction yields
#                 # a carry-out, the following logic should still be correct.
#                 src1_full = UInt(reinterpret(UInt16, src1v))
#                 src2_full = UInt(reinterpret(UInt16, src2v))
#                 full_result = optype(optype(src1_full, src2_full), UInt(cin))
#
#                 # If the full result is greater than the maximum expressible
#                 # value of a unsigned 16 bit number.
#                 expected_carry = full_result > typemax(UInt16)
#                 @test newflags.carry == expected_carry
#
#             elseif T == Signed
#                 expected_result = optype(optype(src1v, src2v), Int16(cin))
#                 @test result == expected_result
#
#                 # Compute if overflow should have happened by performing
#                 # the same operation with full 64-bit numbers and checking if
#                 # the result is expressible in 16 bits.
#                 full_result = optype(optype(Int64(src1v), Int64(src2v)), Int64(cin))
#                 expected_overflow = full_result > typemax(Int16) || full_result < typemin(Int16)
#                 @test newflags.overflow == expected_overflow
#             end
#         end
#     end
#
#
#     test = TestBundle()
#
#     @testset "Arithmetic" begin
#         max_s = typemax(Int16)
#         min_s = typemin(Int16)
#
#         # Lets start going down the list.
#
#         # --- ADDSU / ADDU ---
#         op = :ADDSU
#
#         # Loop through a few operand pairs, making sure to have a mix of numbers
#         # that would be negative 16 bit values is interpreted as signed numbers.
#         # Use a full 64-bit number for comparison.
#         test.src1 = reinterpret.(Int16, [0x0000, 0xFFFF, 0x8001])
#         test.src2 = reinterpret.(Int16, [0x0000, 0xFFFF, 0x1234])
#         test.cin  = [false, false, false]
#         op_tester(Unsigned, +, op, test, benchmark = benchmark)
#
#         # --- ADDCSU / ADDCU ---
#         op = :ADDCSU
#         test.src1 = reinterpret.(Int16, [0xFFFE, 0xFFFF, 0x0000, 0x0001, 0x0001])
#         test.src2 = reinterpret.(Int16, [0x0000, 0x0000, 0x0000, 0xFFFF, 0xFFFE])
#         test.cin  = [true, true, true, false, true]
#         op_tester(Unsigned, +, op, test, benchmark = benchmark)
#
#         # --- ADDS / ADD ---
#         op = :ADD
#         test.src1 = [ -1  , max_s , min_s,    0]
#         test.src2 = [min_s, 3     , max_s, min_s]
#         test.cin  = zeros(Bool, length(test.src1))
#         op_tester(Signed, +, op, test, benchmark = benchmark)
#     end
# end

@testset "Testing Carry Chains" begin
    # Testing subtraction across 5-word elements. This is an operation in the
    # SORT kernal.
    flags = AsapSim.ALUFlags()
    # Don't need condiitonal execution flags.
    condexec = AsapSim.CondExec[]

    result, flags = AsapSim.stage4_alu(AsapSim.SUBU,  Int16(0), Int16(1), flags, condexec)
    # Subtracting 1 from 0 should result in a carry.
    # Carry should propogate through all stages.
    @test flags.carry == true
    @test result == -1

    result, flags = AsapSim.stage4_alu(AsapSim.SUBCU, Int16(0), Int16(0), flags, condexec)
    @test flags.carry == true
    @test result == -1

    result, flags = AsapSim.stage4_alu(AsapSim.SUBCU, Int16(0), Int16(0), flags, condexec)
    @test flags.carry == true
    @test result == -1

    result, flags = AsapSim.stage4_alu(AsapSim.SUBCU, Int16(0), Int16(0), flags, condexec)
    @test flags.carry == true
    @test result == -1

    # Perform a signed subtraction at the end. Carry should still be set.
    result, flags = AsapSim.stage4_alu(AsapSim.SUBC, Int16(0), Int16(0), flags, condexec)
    @test flags.carry == true
    @test result == -1
end

ext: jl


text:
module TestPkg

end


text:
@testset "Java API" begin
    @testset "Unit tests" begin
        @testset "Type hierarchy" begin
            for T in [
                IloConstraint,
                IloAlternative,
                IloIsomorphism,
                IloNoOverlap,
                IloRange,
                IloSpan,
                IloSynchronize,
            ]
                @test T <: Constraint
                @test Vector{T} <: ConstraintArray
                @test T <: Addable
                @test Vector{T} <: AddableArray

                @test T <: IntExpr
                @test Vector{T} <: IntExprArray
            end

            for T in [IloIntVar, IloIntExpr, Constraint]
                @test T <: IntExpr
                @test Vector{T} <: IntExprArray
            end

            for T in [IloIntVar, IloNumVar]
                @test T <: NumVar
                @test Vector{T} <: NumVarArray
            end

            for T in [IntExpr, IloNumVar, IloNumExpr]
                @test T <: NumExpr
                @test Vector{T} <: NumExprArray
            end
        end

        @testset "Model initialisation" begin
            model = JavaCPOModel()
            @test typeof(model) == JavaCPOModel
        end

        @testset "Variable creation" begin
            model = JavaCPOModel()

            # Integer variables
            @test typeof(cpo_java_intvar(model, Int32(5), Int32(10))) ==
                  IloIntVar
            @test typeof(cpo_java_intvar(model, Int32(10), Int32(5))) ==
                  IloIntVar # No check on CPLEXCP side
            @test typeof(cpo_java_intvar(model, Int32(5), Int32(10), "var")) ==
                  IloIntVar

            @test typeof(cpo_java_intvar(model, Int32[5, 10])) == IloIntVar
            @test typeof(cpo_java_intvar(model, Int32[5, 10], "var")) ==
                  IloIntVar

            l = cpo_java_intvararray(model, Int32(20), Int32(5), Int32(10))
            @test typeof(l) == Vector{IloIntVar}
            @test eltype(l) == IloIntVar
            @test length(l) == 20
            l = cpo_java_intvararray(model, Int32(20), Int32[5, 10])
            @test typeof(l) == Vector{IloIntVar}
            @test eltype(l) == IloIntVar
            @test length(l) == 20

            # Numerical variables
            @test typeof(cpo_java_numvar(model, 5.0, 10.0)) == IloNumVar

            l = cpo_java_numvararray(model, Int32(20), 5.0, 10.0)
            @test typeof(l) == Vector{IloNumVar}
            @test eltype(l) == IloNumVar
            @test length(l) == 20

            # Interval variables
            @test typeof(cpo_java_intervalvar(model)) == IloIntervalVar
            @test typeof(cpo_java_intervalvar(model, "var")) == IloIntervalVar
            @test typeof(cpo_java_intervalvar(model, Int32(5))) ==
                  IloIntervalVar
            @test typeof(cpo_java_intervalvar(model, Int32(5), "var")) ==
                  IloIntervalVar
            @test typeof(cpo_java_intervalvar(model, Int32(5), Int32(10))) ==
                  IloIntervalVar

            # Sequence-of-intervals variables
            i1 = cpo_java_intervalvar(model)
            i2 = cpo_java_intervalvar(model)
            @test typeof(cpo_java_intervalsequencevar(model, [i1, i2])) ==
                  IloIntervalSequenceVar
        end
    end

    @testset "Integration tests" begin
        # Based on the provided examples, translated from Java.

        @testset "Color" begin
            model = JavaCPOModel()
            belgium = cpo_java_intvar(model, Int32(0), Int32(3))
            denmark = cpo_java_intvar(model, Int32(0), Int32(3))
            france = cpo_java_intvar(model, Int32(0), Int32(3))
            germany = cpo_java_intvar(model, Int32(0), Int32(3))
            luxembourg = cpo_java_intvar(model, Int32(0), Int32(3))
            netherlands = cpo_java_intvar(model, Int32(0), Int32(3))

            cpo_java_add(model, cpo_java_neq(model, belgium, france))
            cpo_java_add(model, cpo_java_neq(model, belgium, germany))
            cpo_java_add(model, cpo_java_neq(model, belgium, netherlands))
            cpo_java_add(model, cpo_java_neq(model, belgium, luxembourg))
            cpo_java_add(model, cpo_java_neq(model, denmark, germany))
            cpo_java_add(model, cpo_java_neq(model, france, germany))
            cpo_java_add(model, cpo_java_neq(model, france, luxembourg))
            cpo_java_add(model, cpo_java_neq(model, germany, luxembourg))
            cpo_java_add(model, cpo_java_neq(model, germany, netherlands))

            status = cpo_java_solve(model)
            @test status

            @test cpo_java_getvalue(model, belgium) ==
                  cpo_java_getvalue(model, belgium)
            @test cpo_java_getvalue(model, denmark) ==
                  cpo_java_getvalue(model, denmark)
            @test cpo_java_getvalue(model, france) ==
                  cpo_java_getvalue(model, france)
            @test cpo_java_getvalue(model, germany) ==
                  cpo_java_getvalue(model, germany)
            @test cpo_java_getvalue(model, luxembourg) ==
                  cpo_java_getvalue(model, luxembourg)
            @test cpo_java_getvalue(model, netherlands) ==
                  cpo_java_getvalue(model, netherlands)

            @test cpo_java_getvalue(model, belgium) !=
                  cpo_java_getvalue(model, france)
            @test cpo_java_getvalue(model, belgium) !=
                  cpo_java_getvalue(model, germany)
            @test cpo_java_getvalue(model, belgium) !=
                  cpo_java_getvalue(model, netherlands)
            @test cpo_java_getvalue(model, belgium) !=
                  cpo_java_getvalue(model, luxembourg)
            @test cpo_java_getvalue(model, denmark) !=
                  cpo_java_getvalue(model, germany)
            @test cpo_java_getvalue(model, france) !=
                  cpo_java_getvalue(model, germany)
            @test cpo_java_getvalue(model, france) !=
                  cpo_java_getvalue(model, luxembourg)
            @test cpo_java_getvalue(model, germany) !=
                  cpo_java_getvalue(model, luxembourg)
            @test cpo_java_getvalue(model, germany) !=
                  cpo_java_getvalue(model, netherlands)
        end

        @testset "Alloc" begin
            # Model data.
            ncell = Int32(25)
            nfreq = Int32(256)
            nchannel = Int32[
                8,
                6,
                6,
                1,
                4,
                4,
                8,
                8,
                8,
                8,
                4,
                9,
                8,
                4,
                4,
                10,
                8,
                9,
                8,
                4,
                5,
                4,
                8,
                1,
                1,
            ]
            dist = Int32[
                16 1 1 0 0 0 0 0 1 1 1 1 1 2 2 1 1 0 0 0 2 2 1 1 1
                1 16 2 0 0 0 0 0 2 2 1 1 1 2 2 1 1 0 0 0 0 0 0 0 0
                1 2 16 0 0 0 0 0 2 2 1 1 1 2 2 1 1 0 0 0 0 0 0 0 0
                0 0 0 16 2 2 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 1 1
                0 0 0 2 16 2 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 1 1
                0 0 0 2 2 16 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 1 1
                0 0 0 0 0 0 16 2 0 0 1 1 1 0 0 1 1 1 1 2 0 0 0 1 1
                0 0 0 0 0 0 2 16 0 0 1 1 1 0 0 1 1 1 1 2 0 0 0 1 1
                1 2 2 0 0 0 0 0 16 2 2 2 2 2 2 1 1 1 1 1 1 1 0 1 1
                1 2 2 0 0 0 0 0 2 16 2 2 2 2 2 1 1 1 1 1 1 1 0 1 1
                1 1 1 0 0 0 1 1 2 2 16 2 2 2 2 2 2 1 1 2 1 1 0 1 1
                1 1 1 0 0 0 1 1 2 2 2 16 2 2 2 2 2 1 1 2 1 1 0 1 1
                1 1 1 0 0 0 1 1 2 2 2 2 16 2 2 2 2 1 1 2 1 1 0 1 1
                2 2 2 0 0 0 0 0 2 2 2 2 2 16 2 1 1 1 1 1 1 1 1 1 1
                2 2 2 0 0 0 0 0 2 2 2 2 2 2 16 1 1 1 1 1 1 1 1 1 1
                1 1 1 0 0 0 1 1 1 1 2 2 2 1 1 16 2 2 2 1 2 2 1 2 2
                1 1 1 0 0 0 1 1 1 1 2 2 2 1 1 2 16 2 2 1 2 2 1 2 2
                0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 2 2 16 2 2 1 1 0 2 2
                0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 16 2 1 1 0 2 2
                0 0 0 1 1 1 2 2 1 1 2 2 2 1 1 1 1 2 2 16 1 1 0 1 1
                2 0 0 0 0 0 0 0 1 1 1 1 1 1 1 2 2 1 1 1 16 2 1 2 2
                2 0 0 0 0 0 0 0 1 1 1 1 1 1 1 2 2 1 1 1 2 16 1 2 2
                1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 1 1 16 1 1
                1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 1 2 2 1 16 2
                1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 1 2 2 1 2 16
            ]

            function transmitter_idx(cell::Integer, channel::Integer)::Int32
                idx = 0
                c = 1
                while c < cell
                    idx += nchannel[c]
                    c += 1
                end
                return idx + channel
            end

            # Test the helper function for discrepancies due to 1-based indexing.
            @test transmitter_idx(1, 1) == 1
            @test transmitter_idx(ncell, 1) == sum(nchannel)

            # Build the model.
            model = JavaCPOModel()
            ntransmitter = transmitter_idx(ncell, 1)
            freq = cpo_java_intvararray(
                model,
                ntransmitter,
                Int32(1),
                nfreq,
                "freq",
            )

            for cell in 1:ncell
                for channel1 in 1:nchannel[cell]
                    for channel2 in (channel1 + 1):nchannel[cell]
                        diff = cpo_java_diff(
                            model,
                            freq[transmitter_idx(cell, channel1)],
                            freq[transmitter_idx(cell, channel2)],
                        )
                        cstr = cpo_java_ge(model, cpo_java_abs(model, diff), 16)
                        cpo_java_add(model, cstr)
                    end
                end
            end

            for cell1 in 1:ncell
                for cell2 in (cell1 + 1):ncell
                    if dist[cell1, cell2] > 0
                        for channel1 in 1:nchannel[cell1]
                            for channel2 in 1:nchannel[cell2]
                                diff = cpo_java_diff(
                                    model,
                                    freq[transmitter_idx(cell1, channel1)],
                                    freq[transmitter_idx(cell2, channel2)],
                                )
                                cstr = cpo_java_ge(
                                    model,
                                    cpo_java_abs(model, diff),
                                    dist[cell1, cell2],
                                )
                                cpo_java_add(model, cstr)
                            end
                        end
                    end
                end
            end

            obj = cpo_java_countdifferent(model, freq)
            cpo_java_minimize(model, obj)

            status = cpo_java_solve(model)
            @test status
        end

        @testset "Sports" begin
            # Model data.
            n = Int32(10)
            nweeks = Int32(2 * (n - 1))
            ngamesperweek = Int32(n / 2)
            ngames = n * (n - 1)

            function game(h::Int32, a::Int32, n::Int32)
                if a > h
                    return h * (n - 1) + a - 1
                else
                    return h * (n - 1) + a
                end
            end

            # Build the model.
            model = JavaCPOModel()
            games = Dict{Tuple{Int, Int}, IloIntVar}() # nweeks, ngamesperweek
            home = Dict{Tuple{Int, Int}, IloIntVar}() # nweeks, ngamesperweek
            away = Dict{Tuple{Int, Int}, IloIntVar}() # nweeks, ngamesperweek

            for i in 0:(nweeks - 1)
                for j in 0:(ngamesperweek - 1)
                    home[i, j] = cpo_java_intvar(model, Int32(0), Int32(n - 1))
                    away[i, j] = cpo_java_intvar(model, Int32(0), Int32(n - 1))
                    games[i, j] =
                        cpo_java_intvar(model, Int32(0), Int32(ngames - 1))
                end
            end

            # For each play slot, set up correspondance between game id, home team, and away team
            gha = cpo_java_inttable(model, Int32(3))
            for i in Int32(0):Int32(n - 1)
                for j in Int32(0):Int32(n - 1)
                    if i != j
                        cpo_java_inttupleset_addtuple(
                            model,
                            gha,
                            Int32[i, j, game(i, j, n)],
                        )
                    end
                end
            end

            for i in 0:(nweeks - 1)
                for j in 0:(ngamesperweek - 1)
                    vars = [home[i, j], away[i, j], games[i, j]]
                    cpo_java_add(
                        model,
                        cpo_java_allowedassignments(model, vars, gha),
                    )
                end
            end

            # All teams play each week
            for i in 0:(nweeks - 1)
                teamsthisweek = Dict{Int, IloIntVar}() # n
                for j in 0:(ngamesperweek - 1)
                    teamsthisweek[j] = home[i, j]
                    teamsthisweek[ngamesperweek + j] = away[i, j]
                end
                teamsthisweek_vec = [teamsthisweek[i] for i in 0:(n - 1)]
                cpo_java_add(model, cpo_java_alldiff(model, teamsthisweek_vec))
            end

            # Dual representation: for each game id, the play slot is maintained
            weekofgame = Dict{Int, IloIntVar}() # ngames
            allgames = Dict{Int, IloIntVar}() # ngames
            allslots = Dict{Int, IloIntVar}() # ngames

            for i in 0:(ngames - 1)
                weekofgame[i] =
                    cpo_java_intvar(model, Int32(0), Int32(nweeks - 1))
                allslots[i] =
                    cpo_java_intvar(model, Int32(0), Int32(ngames - 1))
            end

            for i in 0:(nweeks - 1)
                for j in 0:(ngamesperweek - 1)
                    allgames[i * ngamesperweek + j] = games[i, j]
                end
            end
            allgames_vec = [allgames[i] for i in 0:(ngames - 1)]
            allslots_vec = [allslots[i] for i in 0:(ngames - 1)]
            cpo_java_add(
                model,
                cpo_java_inverse(model, allgames_vec, allslots_vec),
            )

            for i in 0:(ngames - 1)
                cpo_java_add(
                    model,
                    cpo_java_eq(
                        model,
                        weekofgame[i],
                        cpo_java_div(model, allslots[i], ngamesperweek),
                    ),
                )
            end

            # Two half schedules.  Cannot play the same pair twice in the same half.
            # Plus, impose a minimum number of weeks between two games involving
            # the same teams (up to six weeks)
            mid = Int32(nweeks / 2)
            overlap = Int32(0)
            if n >= 6
                overlap = min(Int32(n / 2), Int32(6))
            end

            for i in Int32(0):Int32(n - 1)
                for j in Int32(i + 1):Int32(n - 1)
                    g1 = game(i, j, n)
                    g2 = game(j, i, n)

                    c = cpo_java_equiv(
                        model,
                        cpo_java_ge(model, weekofgame[g1], mid),
                        cpo_java_lt(model, weekofgame[g2], mid),
                    )
                    cpo_java_add(model, c)

                    if (overlap != 0)
                        c = cpo_java_ge(
                            model,
                            cpo_java_abs(
                                model,
                                cpo_java_diff(
                                    model,
                                    weekofgame[g1],
                                    weekofgame[g2],
                                ),
                            ),
                            overlap,
                        )
                        cpo_java_add(model, c)
                    end
                end
            end

            # Can't have three homes or three aways in a row.
            playhome = Dict{Tuple{Int, Int}, IloIntVar}() # n, nweeks
            for i in Int32(0):Int32(n - 1)
                for j in Int32(0):Int32(nweeks - 1)
                    playhome[i, j] = cpo_java_boolvar(model)
                    home_vec = [home[j, k] for k in 0:(ngamesperweek - 1)]
                    cpo_java_add(
                        model,
                        cpo_java_eq(
                            model,
                            playhome[i, j],
                            cpo_java_count(model, home_vec, i),
                        ),
                    )
                end

                for j in 0:(nweeks - 3 - 1)
                    windowsum =
                        cpo_java_sum(model, [playhome[i, k] for k in j:(j + 2)])
                    cpo_java_add(model, cpo_java_ge(model, windowsum, 1))
                    cpo_java_add(model, cpo_java_le(model, windowsum, 2))
                end
            end

            # If we start the season home, we finish away and vice versa.
            for i in 0:(n - 1)
                cpo_java_add(
                    model,
                    cpo_java_neq(
                        model,
                        playhome[i, 0],
                        playhome[i, nweeks - 1],
                    ),
                )
            end

            # Objective is skipped, not useful for tests. Finds a solution in 0.5s on a powerful machine,
            # no need for more complexity.

            # Each team plays home the same number of times as away
            for i in 0:(n - 1)
                playhomesum = cpo_java_sum(
                    model,
                    [playhome[i, j] for j in 0:(nweeks - 1)],
                )
                cpo_java_add(
                    model,
                    cpo_java_eq(model, playhomesum, Int(nweeks / 2)),
                )
            end

            # Teams are interchangeable. Fix first week. Also breaks reflection symmetry of the whole schedule.
            for i in 0:(ngamesperweek - 1)
                cpo_java_add(model, cpo_java_eq(model, home[0, i], i * 2))
                cpo_java_add(model, cpo_java_eq(model, away[0, i], i * 2 + 1))
            end

            # Order of games in each week is arbitrary. Break symmetry by forcing an order.
            for i in 0:(nweeks - 1)
                for j in 1:(ngamesperweek - 1)
                    cpo_java_add(
                        model,
                        cpo_java_gt(model, games[i, j], games[i, j - 1]),
                    )
                end
            end

            status = cpo_java_solve(model)
            @test status
        end

        @testset "Facility" begin
            # Model data.
            nlocations = 5
            nstores = 8
            capacity = Int32[3, 1, 2, 4, 1]
            fixedcost = Int32[480, 200, 320, 340, 300]
            cost = Int32[
                24 74 31 51 84
                57 54 86 61 68
                57 67 29 91 71
                54 54 65 82 94
                98 81 16 61 27
                13 92 34 94 87
                54 72 41 12 78
                54 64 65 89 89
            ]

            # Build the model.
            model = JavaCPOModel()
            supplier = cpo_java_intvararray(
                model,
                Int32(nstores),
                Int32(0),
                Int32(nlocations - 1),
            )
            open = cpo_java_boolvararray(model, Int32(nlocations))

            for i in 1:nstores
                cpo_java_add(
                    model,
                    cpo_java_eq(
                        model,
                        cpo_java_element(model, open, supplier[i]),
                        1,
                    ),
                )
            end

            for j in Int32(1):Int32(nlocations)
                cpo_java_add(
                    model,
                    cpo_java_le(
                        model,
                        cpo_java_count(model, supplier, j),
                        capacity[j],
                    ),
                )
            end

            obj = cpo_java_scalprod(model, fixedcost, open)
            for i in Int32(1):Int32(nstores)
                obj = cpo_java_sum(
                    model,
                    IntExpr[
                        obj,
                        cpo_java_element(model, cost[i, :], supplier[i]),
                    ],
                )
            end
            cpo_java_add(model, cpo_java_minimize(model, obj))

            status = cpo_java_solve(model)
            @test status
        end
    end
end


text:
module TestLevels
    using Base.Test
    using CategoricalArrays
    using CategoricalArrays: DefaultRefType

    pool = CategoricalPool([2, 1, 3])

    @test isa(levels(pool), Vector{Int})
    @test length(levels(pool)) === 3
    @test levels(pool) == pool.index == [2, 1, 3]
    @test pool.invindex == Dict(1=>2, 2=>1, 3=>3)
    @test pool.order == [1, 2, 3]
    @test pool.valindex == [CategoricalValue(i, pool) for i in 1:3]

    for rep in 1:3
        push!(pool, 4)

        @test isa(pool.index, Vector{Int})
        @test length(pool) === 4
        @test pool.index == [2, 1, 3, 4]
        @test pool.invindex == Dict(1=>2, 2=>1, 3=>3, 4=>4)
        @test pool.order == [1, 2, 3, 4]
        @test pool.levels == [2, 1, 3, 4]
        @test get(pool, 4) === DefaultRefType(4)
        @test pool[4] === CategoricalValue(4, pool)
        @test pool.valindex == [CategoricalValue(i, pool) for i in 1:4]
    end

    for rep in 1:3
        push!(pool, 0)

        @test isa(pool.index, Vector{Int})
        @test length(pool) === 5
        @test levels(pool) == pool.index == [2, 1, 3, 4, 0]
        @test pool.invindex == Dict(1=>2, 2=>1, 3=>3, 4=>4, 0=>5)
        @test pool.order == [1, 2, 3, 4, 5]
        @test pool.levels == [2, 1, 3, 4, 0]
        @test get(pool, 0) === DefaultRefType(5)
        @test pool[5] === CategoricalValue(5, pool)
        @test pool.valindex == [CategoricalValue(i, pool) for i in 1:5]
    end

    for rep in 1:3
        push!(pool, 10, 11)

        @test isa(pool.index, Vector{Int})
        @test length(pool) === 7
        @test levels(pool) == pool.index == [2, 1, 3, 4, 0, 10, 11]
        @test pool.invindex == Dict(1=>2, 2=>1, 3=>3, 4=>4, 0=>5, 10=>6, 11=>7)
        @test pool.order == [1, 2, 3, 4, 5, 6, 7]
        @test pool.levels == [2, 1, 3, 4, 0, 10, 11]
        @test get(pool, 10) === DefaultRefType(6)
        @test get(pool, 11) === DefaultRefType(7)
        @test pool[6] === CategoricalValue(6, pool)
        @test pool[7] === CategoricalValue(7, pool)
        @test pool.valindex == [CategoricalValue(i, pool) for i in 1:7]
    end

    for rep in 1:3
        push!(pool, 12, 13)

        @test isa(pool.index, Vector{Int})
        @test length(pool) === 9
        @test levels(pool) == pool.index == [2, 1, 3, 4, 0, 10, 11, 12, 13]
        @test pool.invindex == Dict(1=>2, 2=>1, 3=>3, 4=>4, 0=>5, 10=>6, 11=>7, 12=>8, 13=>9)
        @test pool.order == [1, 2, 3, 4, 5, 6, 7, 8, 9]
        @test pool.levels == [2, 1, 3, 4, 0, 10, 11, 12, 13]
        @test get(pool, 12) === DefaultRefType(8)
        @test get(pool, 13) === DefaultRefType(9)
        @test pool[8] === CategoricalValue(8, pool)
        @test pool[9] === CategoricalValue(9, pool)
        @test pool.valindex == [CategoricalValue(i, pool) for i in 1:9]
    end

    for rep in 1:3
        delete!(pool, 13)

        @test isa(pool.index, Vector{Int})
        @test length(pool) == 8
        @test levels(pool) == pool.index == [2, 1, 3, 4, 0, 10, 11, 12]
        @test pool.invindex == Dict(1=>2, 2=>1, 3=>3, 4=>4, 0=>5, 10=>6, 11=>7, 12=>8)
        @test pool.order == [1, 2, 3, 4, 5, 6, 7, 8]
        @test pool.levels == [2, 1, 3, 4, 0, 10, 11, 12]
        @test_throws KeyError get(pool, 13)
        @test pool.valindex == [CategoricalValue(i, pool) for i in 1:8]
    end

    for rep in 1:3
        delete!(pool, 12, 11)

        @test isa(pool.index, Vector{Int})
        @test length(pool) == 6
        @test levels(pool) == pool.index == [2, 1, 3, 4, 0, 10]
        @test pool.invindex == Dict(1=>2, 2=>1, 3=>3, 4=>4, 0=>5, 10=>6)
        @test pool.order == [1, 2, 3, 4, 5, 6]
        @test pool.levels == [2, 1, 3, 4, 0, 10]
        @test_throws KeyError get(pool, 11)
        @test_throws KeyError get(pool, 12)
        @test pool.valindex == [CategoricalValue(i, pool) for i in 1:6]
    end

    for rep in 1:3
        delete!(pool, 4)

        @test isa(pool.index, Vector{Int})
        @test length(pool) == 5
        @test levels(pool) == pool.index == [2, 1, 3, 0, 10]
        @test pool.invindex == Dict(1=>2, 2=>1, 3=>3, 0=>4, 10=>5)
        @test pool.order == [1, 2, 3, 4, 5]
        @test pool.levels == [2, 1, 3, 0, 10]
        @test_throws KeyError get(pool, 4)
        @test pool.valindex == [CategoricalValue(i, pool) for i in 1:5]
    end

    @test levels!(pool, [1, 2, 3]) == [1, 2, 3]

    @test isa(pool.index, Vector{Int})
    @test length(pool) == 3
    @test length(pool.valindex) == 3
    @test levels(pool) == pool.index == [1, 2, 3]
    @test pool.invindex == Dict(1=>1, 2=>2, 3=>3)
    @test pool.order == [1, 2, 3]
    @test pool.levels == [1, 2, 3]
    @test get(pool, 1) === DefaultRefType(1)
    @test_throws KeyError get(pool, 0)
    @test_throws KeyError get(pool, 10)
    @test pool.valindex == [CategoricalValue(i, pool) for i in 1:3]

    @test levels!(pool, [1, 2, 4]) == [1, 2, 4]

    @test isa(pool.index, Vector{Int})
    @test length(pool) == 3
    @test length(pool.valindex) == 3
    @test levels(pool) == pool.index == [1, 2, 4]
    @test pool.invindex == Dict(1=>1, 2=>2, 4=>3)
    @test pool.order == [1, 2, 3]
    @test pool.levels == [1, 2, 4]
    @test get(pool, 1) === DefaultRefType(1)
    @test_throws KeyError get(pool, 3)
    @test pool.valindex == [CategoricalValue(i, pool) for i in 1:3]

    @test levels!(pool, [6, 5, 4]) == [6, 5, 4]

    @test isa(pool.index, Vector{Int})
    @test length(pool) == 3
    @test length(pool.valindex) == 3
    @test levels(pool) == pool.index == [6, 5, 4]
    @test pool.invindex == Dict(6=>1, 5=>2, 4=>3)
    @test pool.order == [1, 2, 3]
    @test pool.levels == [6, 5, 4]
    @test get(pool, 5) === DefaultRefType(2)
    @test_throws KeyError get(pool, 3)
    @test pool.valindex == [CategoricalValue(i, pool) for i in 1:3]

    # Changing order while preserving existing levels
    @test levels!(pool, [5, 6, 4]) == [5, 6, 4]

    @test isa(pool.index, Vector{Int})
    @test length(pool) == 3
    @test length(pool.valindex) == 3
    @test levels(pool) == [5, 6, 4]
    @test pool.index == [6, 5, 4]
    @test pool.invindex == Dict(6=>1, 5=>2, 4=>3)
    @test pool.order == [2, 1, 3]
    @test pool.levels == [5, 6, 4]
    @test get(pool, 5) === DefaultRefType(2)
    @test pool.valindex == [CategoricalValue(i, pool) for i in 1:3]

    # Adding levels while preserving existing ones
    @test levels!(pool, [5, 2, 3, 6, 4]) == [5, 2, 3, 6, 4]

    @test isa(pool.index, Vector{Int})
    @test length(pool) == 5
    @test length(pool.valindex) == 5
    @test levels(pool) == [5, 2, 3, 6, 4]
    @test pool.index == [6, 5, 4, 2, 3]
    @test pool.invindex == Dict(6=>1, 5=>2, 4=>3, 2=>4, 3=>5)
    @test pool.order == [4, 1, 5, 2, 3]
    @test pool.levels == [5, 2, 3, 6, 4]
    @test get(pool, 2) === DefaultRefType(4)
    @test get(pool, 3) === DefaultRefType(5)
    @test pool.valindex == [CategoricalValue(i, pool) for i in 1:5]

    for rep in 1:3
        delete!(pool, 6)

        @test isa(pool.index, Vector{Int})
        @test length(pool) == 4
        @test length(pool.valindex) == 4
        @test levels(pool) == [5, 2, 3, 4]
        @test pool.index == [5, 4, 2, 3]
        @test pool.invindex == Dict(5=>1, 4=>2, 2=>3, 3=>4)
        @test pool.order == [1, 4, 2, 3]
        @test pool.levels == [5, 2, 3, 4]
        @test get(pool, 4) === DefaultRefType(2)
        @test_throws KeyError get(pool, 6)
        @test pool.valindex == [CategoricalValue(i, pool) for i in 1:4]
    end
end


text:
module UWDDynam
using Catlab
using Catlab.WiringDiagrams
using Catlab.CategoricalAlgebra
using Catlab.CategoricalAlgebra.FinSets
using Catlab.Theories

using Catlab.WiringDiagrams.UndirectedWiringDiagrams: AbstractUWD
import Catlab.WiringDiagrams: oapply

using OrdinaryDiffEq, DynamicalSystems
import OrdinaryDiffEq: ODEProblem
import DynamicalSystems: DiscreteDynamicalSystem

export AbstractResourceSharer, ContinuousResourceSharer, DiscreteResourceSharer,
euler_approx, nstates, nports, portmap, portfunction,
eval_dynamics, eval_dynamics!, exposed_states, fills, induced_states

using Base.Iterators
import Base: show, eltype

"""

An undirected open dynamical system operating on information of type `T`.

A resource sharer `r` has signature `r.nports`.
"""
abstract type AbstractResourceSharer{T} end

"""

An undirected open continuous system. The dynamics function `f` defines an ODE ``\\dot u(t) = f(u(t),p,t)``.
"""
struct ContinuousResourceSharer{T} <: AbstractResourceSharer{T}
  nports::Int
  nstates::Int
  dynamics::Function
  portmap::Vector{Int64}
end

"""

An undirected open discrete system. The dynamics function `f` defines a discrete update rule ``u_{n+1} = f(u_n, p, t)``.
"""
struct DiscreteResourceSharer{T} <: AbstractResourceSharer{T}
    nports::Int
    nstates::Int
    dynamics::Function
    portmap::Vector{Int64}
end

ContinuousResourceSharer{T}(nstates::Int, dynamics::Function) where T =
    ContinuousResourceSharer{T}(nstates,nstates, dynamics, Vector{Int64}(1:nstates))
DiscreteResourceSharer{T}(nstates::Int, dynamics::Function) where T =
    DiscreteResourceSharer{T}(nstates,nstates, dynamics, Vector{Int64}(1:nstates))

nstates(r::AbstractResourceSharer) = r.nstates
nports(r::AbstractResourceSharer)  = r.nports
portmap(r::AbstractResourceSharer) = r.portmap
portfunction(r::AbstractResourceSharer) = FinFunction(r.portmap, nstates(r))
exposed_states(r::AbstractResourceSharer, u::AbstractVector) = getindex(u, portmap(r))

"""    eval_dynamics(r::AbstractResourceSharer, u::AbstractVector, p, t)

Evaluates the dynamics of the resource sharer `r` at state `u`, parameters `p`, and time `t`.

Omitting `t` and `p` is allowed if the dynamics of `r` does not depend on them.
"""
eval_dynamics(r::AbstractResourceSharer, u::AbstractVector, p, t::Real) = r.dynamics(u, p, t)
eval_dynamics!(du, r::AbstractResourceSharer, u::AbstractVector, p, t::Real) = begin
    du .= eval_dynamics(r, u, p, t)
end
eval_dynamics(r::AbstractResourceSharer, u::AbstractVector) = eval_dynamics(r, u, [], 0)
eval_dynamics(r::AbstractResourceSharer, u::AbstractVector, p) = eval_dynamics(r, u, p, 0)

show(io::IO, vf::ContinuousResourceSharer) = print("ContinuousResourceSharer(ℝ^$(vf.nstates) → ℝ^$(vf.nstates)) with $(vf.nports) exposed ports")
show(io::IO, vf::DiscreteResourceSharer) = print("DiscreteResourceSharer(ℝ^$(vf.nstates) → ℝ^$(vf.nstates)) with $(vf.nports) exposed ports")
eltype(r::AbstractResourceSharer{T}) where T = T

"""    euler_approx(r::ContinuousResourceSharer, h)

Transforms a continuous resource sharer `r` into a discrete resource sharer via Euler's method with step size `h`. If the dynamics of `r` is given by ``\\dot{u}(t) = f(u(t),p,t)`` the the dynamics of the new discrete system is given by the update rule ``u_{n+1} = u_n + h f(u_n, p, t)``.
"""
euler_approx(f::ContinuousResourceSharer{T}, h::Float64) where T = DiscreteResourceSharer{T}(
        nports(f), nstates(f),
        (u, p, t) -> u + h*eval_dynamics(f, u, p, t),
        f.portmap
)

"""    euler_approx(r::ContinuousResourceSharer)

Transforms a continuous resource sharer `r` into a discrete resource sharer via Euler's method where the step size is introduced as a new parameter, the last in the list of parameters.
"""
euler_approx(f::ContinuousResourceSharer{T}) where T = DiscreteResourceSharer{T}(
    nports(f), nstates(f),
    (u, p, t) -> u + p[end]*eval_dynamics(f, u, p[1:end-1], t),
    f.portmap
)

euler_approx(fs::Vector{ContinuousResourceSharer{T}}, args...) where T =
    map(f->euler_approx(f,args...), fs)

euler_approx(fs::AbstractDict{S, ContinuousResourceSharer{T}}, args...) where {S, T} =
    Dict(name => euler_approx(f, args...) for (name, f) in fs)

"""    ODEProblem(r::ContinuousResourceSharer, u0::Vector, tspan)

Constructs an ODEProblem from the vector field defined by `r.dynamics(u,p,t)`.
"""
ODEProblem(r::ContinuousResourceSharer, u0::AbstractVector, tspan::Tuple{Real, Real}, p=nothing) =
    ODEProblem(r.dynamics, u0, tspan, p)

"""    DiscreteDynamicalSystem(r::DiscreteResourceSharer, u0::Vector, p)

Constructs a DiscreteDynamicalSystem from the equation of motion `r.dynamics(u,p,t)`.  Pass `nothing` in place of `p` if your system does not have parameters.
"""
DiscreteDynamicalSystem(r::DiscreteResourceSharer{T}, u0::AbstractVector, p; t0::Int = 0) where T = begin
    if nstates(r) == 1
      DiscreteDynamicalSystem1d(r, u0[1], p; t0 = t0)
    else
      !(T <: AbstractFloat) || error("Cannot construct a DiscreteDynamicalSystem if the type is a float")
      DiscreteDynamicalSystem((u,p,t) -> SVector{nstates(r)}(eval_dynamics(r,u,p,t)), u0, p; t0=t0)
    end
  end

  DiscreteDynamicalSystem(r::DiscreteResourceSharer, u0::Real, p; t0::Int = 0) =
    DiscreteDynamicalSystem1d(r, u0, p; t0 = t0)

  # if the system is 1D then the state must be represented by a number NOT by a 1D array
  DiscreteDynamicalSystem1d(r::DiscreteResourceSharer{T}, u0::Real, p; t0::Int = 0) where T = begin
    nstates(r) == 1 || error("The resource sharer must have exactly 1 state")
    !(T <: AbstractFloat) || error("Cannot construct a DiscreteDynamicalSystem if the type is a float")
    DiscreteDynamicalSystem((u,p,t) -> eval_dynamics(r,[u],p,t)[1], u0, p; t0 = t0)
  end

"""    fills(r::AbstractResourceSharer, d::AbstractUWD, b::Int)

Checks if `r` is of the correct signature to fill box `b` of the undirected wiring diagram `d`.
"""
function fills(r::AbstractResourceSharer, d::AbstractUWD, b::Int)
    b <= nparts(d, :Box) || error("Trying to fill box $b, when $d has fewer than $b boxes")
    return nports(r) == length(incident(d, b, :box))
end

"""     oapply(d::AbstractUWD, rs::Vector)

Implements the operad algebras for undirected composition of dynamical systems given a composition pattern (implemented
by an undirected wiring diagram `d`) and primitive systems (implemented by
a collection of resource sharers `rs`). Returns the composite resource sharer.

Each box of `d` must be filled by a resource sharer of the appropriate type signature.
"""
oapply(d::AbstractUWD, xs::Vector{ResourceSharer}) where {ResourceSharer <: AbstractResourceSharer} =
    oapply(d, xs, induced_states(d, xs))

"""    oapply(d::AbstractUWD, r::AbstractResourceSharer)

A version of `oapply` where each box of `d` is filled with the resource sharer `r`.
"""
oapply(d::AbstractUWD, x::AbstractResourceSharer) =
    oapply(d, collect(repeated(x, nboxes(d))))

"""     oapply(d::AbstractUWD, generators::Dict)

A version of `oapply` where `generators` is a dictionary mapping the name of each box to its corresponding resource sharer.
"""
oapply(d::AbstractUWD, xs::AbstractDict{S,T}) where {S, T <: AbstractResourceSharer} =
    oapply(d, [xs[name] for name in subpart(d, :name)])

function oapply(d::AbstractUWD, xs::Vector{ResourceSharer}, S′::Pushout) where {ResourceSharer <: AbstractResourceSharer}

    S = coproduct((FinSet∘nstates).(xs))
    states(b::Int) = legs(S)[b].func

    v = induced_dynamics(d, xs, legs(S′)[1], states)

    junction_map = legs(S′)[2]
    outer_junction_map = FinFunction(subpart(d, :outer_junction), nparts(d, :Junction))

    return ResourceSharer(
        nparts(d, :OuterPort),
        length(apex(S′)),
        v,
        compose(outer_junction_map, junction_map).func)
end

function induced_states(d::AbstractUWD, xs::Vector{ResourceSharer}) where {ResourceSharer <: AbstractResourceSharer}
    for box in parts(d, :Box)
        fills(xs[box], d, box) || error("$(xs[box]) does not fill box $box")
    end

    S = coproduct((FinSet∘nstates).(xs))
    P = coproduct((FinSet∘nports).(xs))
    total_portfunction = copair([compose( portfunction(xs[i]), legs(S)[i]) for i in 1:length(xs)])

    return pushout(total_portfunction, FinFunction(subpart(d, :junction), nparts(d, :Junction)))
end

function induced_dynamics(d::AbstractUWD, xs::Vector{ContinuousResourceSharer{T}}, state_map::FinFunction, states::Function) where T

    function v(u′::AbstractVector, p, t::Real)
      u = getindex(u′,  state_map.func)
      du = zero(u)
      # apply dynamics
      for b in parts(d, :Box)
        eval_dynamics!(view(du, states(b)), xs[b], view(u, states(b)), p, t)
      end
      # add along junctions
      du′ = [sum(Array{T}(view(du, preimage(state_map, i)))) for i in codom(state_map)]
      return du′
    end

end

function induced_dynamics(d::AbstractUWD, xs::Vector{DiscreteResourceSharer{T}}, state_map::FinFunction, states::Function) where T
    function v(u′::AbstractVector, p, t::Real)
        u0 = getindex(u′,  state_map.func)
        u1 = zero(u0)
        # apply dynamics
        for b in parts(d, :Box)
          eval_dynamics!(view(u1, states(b)), xs[b], view(u0, states(b)), p, t)
        end
        Δu = u1 - u0
        # add along junctions
        return u′+ [sum(Array{T}(view(Δu, preimage(state_map, i)))) for i in codom(state_map)]
    end
end

end #module


text:
```python
import sympy
```

```python
import tensorflow as tf
```

    2021-07-24 12:11:31.418536: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcudart.so.11.0'; dlerror: libcudart.so.11.0: cannot open shared object file: No such file or directory
    2021-07-24 12:11:31.418577: I tensorflow/stream_executor/cuda/cudart_stub.cc:29] Ignore above cudart dlerror if you do not have a GPU set up on your machine.

```python
%timeit list(sympy.primerange(1,100))
```

    64.9 µs ± 4.19 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

```python
%timeit [p for p in sympy.primerange(1,100)]
```

    62.4 µs ± 892 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

help(sympy.primerange)

```python
list(sympy.primerange(1,12))
```

    [2, 3, 5, 7, 11]

Since the paper seems to
- not only want the prime numbers
- but also `1`

we'd better start from `0` instead of from `1` (as `Chia-Yen Chiang` did).

```python
list(sympy.primerange(0,13))
```

    [2, 3, 5, 7, 11]

No!? But of course, `sympy` was right: `1` is not a prime.

Let's see how `sympy` think when it comes to negative prime numbers.

```python
list(sympy.primerange(-13,13))
```

    [2, 3, 5, 7, 11]

```python
list(sympy.primerange(2,2))
```

    []

```python
list(sympy.primerange(2, 3.14))
```

    [2, 3]

```python
def kernel_size_generator(k):
    """
    This is just a simple wrapper around sympy.primerange()
    which yields 1 and all the prime numbers up until < k
    """
    if not isinstance(k, int) and not isinstance(k, float):
        raise ValueError("input arg k should be either an int or a float")
    if k < 1:
        raise ValueError("input arg k should be >= 1")
    if k == 1:
        yield 1
    else:
        yield 1
        for p in sympy.primerange(2,k):
            yield p
```

```python
list(kernel_size_generator(13))
```

    [1, 2, 3, 5, 7, 11]

```python
list(kernel_size_generator(3.14))
```

    [1, 2, 3]

```python
list(kernel_size_generator(-10))
```

```python
list(kernel_size_generator("10"))
```

```python

```

ext: ipynb