text:
function encode_num(num::Integer)
    if num == 0
        return [0x00]
    end
    abs_num = abs(num)
    negative = num < 0
    result = UInt8[]
    while abs_num != 0
        push!(result, abs_num & 0xff)
        abs_num >>= 8
    end
    # if the top bit is set,
    # for negative numbers we ensure that the top bit is set
    # for positive numbers we ensure that the top bit is not set
    if result[end] & 0x80 != 0
        if negative
            push!(result, 0x80)
        else
            push!(result, 0x00)
        end
    elseif negative
        result[end] |= 0x80
    end
    return result
end

function decode_num(element::Vector{UInt8})
    if isempty(element)
        return 0
    end
    if length(element) > div(Sys.WORD_SIZE, 8)
        T = BigInt
    else
        T = Int
    end
    result = 0
    # reverse for big endian
    big_endian = reverse(element)
    # top bit being 1 means it's negative
    if big_endian[1] & 0x80 != 0
        negative = true
        result = T(big_endian[1] & 0x7f)
    else
        negative = false
        result = T(big_endian[1])
    end
    for c in big_endian[2:end]
        result <<= 8
        result += c
    end
    if negative
        return -result
    else
        return result
    end
end

StackType = Vector{Vector{UInt8}}
const short_stack_string = "Stack does not contain enough items"

function op_0(stack::StackType)
    push!(stack, encode_num(0))
    return true
end

function op_1negate(stack::StackType)
    push!(stack, encode_num(-1))
    return true
end

function op_1(stack::StackType)
    push!(stack, encode_num(1))
    return true
end

function op_2(stack::StackType)
    push!(stack, encode_num(2))
    return true
end

function op_3(stack::StackType)
    push!(stack, encode_num(3))
    return true
end

function op_4(stack::StackType)
    push!(stack, encode_num(4))
    return true
end

function op_5(stack::StackType)
    push!(stack, encode_num(5))
    return true
end

function op_6(stack::StackType)
    push!(stack, encode_num(6))
    return true
end

function op_7(stack::StackType)
    push!(stack, encode_num(7))
    return true
end

function op_8(stack::StackType)
    push!(stack, encode_num(8))
    return true
end

function op_9(stack::StackType)
    push!(stack, encode_num(9))
    return true
end

function op_10(stack::StackType)
    push!(stack, encode_num(10))
    return true
end

function op_11(stack::StackType)
    push!(stack, encode_num(11))
    return true
end

function op_12(stack::StackType)
    push!(stack, encode_num(12))
    return true
end

function op_13(stack::StackType)
    push!(stack, encode_num(13))
    return true
end

function op_14(stack::StackType)
    push!(stack, encode_num(14))
    return true
end

function op_15(stack::StackType)
    push!(stack, encode_num(15))
    return true
end

function op_16(stack::StackType)
    push!(stack, encode_num(16))
    return true
end

function op_nop(stack::StackType)
    return true
end

# TODO verify actual behavior
function op_if(stack::StackType, items::Vector{UInt8})
    @assert length(stack) ≥ 1 short_stack_string
    # go through and re-make the items array based on the top stack element
    true_items = UInt8[]
    false_items = UInt8[]
    current_array = true_items
    found = false
    num_endifs_needed = 1
    while length(items) > 0
        item = popfirst!(items)
        if item in [99, 100]
            # nested if, we have to go another endif
            num_endifs_needed += 1
            push!(current_array, item)
        elseif num_endifs_needed == 1 && item == 103
            current_array = false_items
        elseif item == 104
            if num_endifs_needed == 1
                found = true
                break
            else
                num_endifs_needed -= 1
                push!(current_array, item)
            end
        else
            push!(current_array, item)
        end
    end
    if !found
        return false
    end
    element = pop!(stack)
    if decode_num(element) == 0
        prepend!(items, false_items)
    else
        prepend!(items, true_items)
    end
    return true
end

function op_notif(stack::StackType, items::Vector{UInt8})
    @assert length(stack) ≥ 1 short_stack_string
    # go through and re-make the items array based on the top stack element
    true_items = UInt8[]
    false_items = UInt8[]
    current_array = true_items
    found = False
    num_endifs_needed = 1
    while length(items) > 0
        item = popfirst!(items)
        if item in [99, 100]
            # nested if, we have to go another endif
            num_endifs_needed += 1
            push!(current_array, item)
        elseif num_endifs_needed == 1 && item == 103
            current_array = false_items
        elseif item == 104
            if num_endifs_needed == 1
                found = true
                break
            else
                num_endifs_needed -= 1
                push!(current_array, item)
            end
        else
            push!(current_array, item)
        end
    end
    if !found
        return false
    end
    element = pop!(stack)
    if decode_num(element) == 0
        prepend!(items, true_items)
    else
        prepend!(items, false_items)
        return true
    end
end

function op_verify(stack::StackType)
    @assert length(stack) ≥ 1 short_stack_string
    element = pop!(stack)
    if decode_num(element) == 0
        return false
    end
    return true
end

function op_return(stack::StackType)
    return false
end

function op_toaltstack(stack::StackType, altstack::StackType)
    @assert length(stack) ≥ 1 short_stack_string
    push!(altstack, pop!(stack))
    return true
end

function op_fromaltstack(stack::StackType, altstack::StackType)
    @assert length(stack) ≥ 1 short_stack_string
    push!(stack, pop!(altstack))
    return true
end

function op_2drop(stack::StackType)
    @assert length(stack) ≥ 2 short_stack_string
    pop!(stack)
    pop!(stack)
    return true
end

function op_2dup(stack::StackType)
    @assert length(stack) ≥ 2 short_stack_string
    append!(stack, stack[end-1:end])
    return true
end

function op_3dup(stack::StackType)
    @assert length(stack) ≥ 3 short_stack_string
    append!(stack, stack[end-2:end])
    return true
end

function op_2over(stack::StackType)
    @assert length(stack) ≥ 4 short_stack_string
    append!(stack, stack[end-3:end-2])
    return true
end

function op_2rot(stack::StackType)
    if length(stack) < 6
        return false
    end
    append!(stack, stack[end-5:end-4])
    return true
end

function op_2swap(stack::StackType)
    @assert length(stack) ≥ 4 short_stack_string
    stack[end-3:end] = append!(stack[end-1:end],stack[end-3:end-2])
    return true
end

function op_ifdup(stack::StackType)
    @assert length(stack) ≥ 1 short_stack_string
    if decode_num(stack[end]) != 0
        push!(stack, stack[end])
    end
    return true
end

function op_depth(stack::StackType)
    push!(stack, encode_num(length(stack)))
    return true
end

function op_drop(stack::StackType)
    @assert length(stack) ≥ 1 short_stack_string
    pop!(stack)
    return true
end

function op_dup(stack::StackType)
    @assert length(stack) ≥ 1 short_stack_string
    push!(stack, stack[end])
    return true
end

function op_nip(stack::StackType)
    @assert length(stack) ≥ 2 short_stack_string
    splice!(stack,length(stack)-1)
    return true
end

function op_over(stack::StackType)
    @assert length(stack) ≥ 2 short_stack_string
    push!(stack, stack[end-1])
    return true
end

function op_pick(stack::StackType)
    @assert length(stack) ≥ 1 short_stack_string
    n = decode_num(pop!(stack))
    if length(stack) < n
        return false
    end
    push!(stack, stack[end-n+1])
    return true
end

function op_roll(stack::StackType)
    @assert length(stack) ≥ 1 short_stack_string
    n = decode_num(pop!(stack))
    if length(stack) < n
        return false
    end
    if n == 0
        return true
    end
    element = splice!(stack, length(stack)-n+1)
    push!(stack, element)
    return true
end

function op_rot(stack::StackType)
    @assert length(stack) ≥ 3 short_stack_string
    stack[end-2:end] = circshift(stack[end-2:end],-1)
    return true
end

function op_swap(stack::StackType)
    @assert length(stack) ≥ 2 short_stack_string
    element = splice!(stack, length(stack) - 1)
    push!(stack, element)
    return true
end

function op_tuck(stack::StackType)
    @assert length(stack) ≥ 2 short_stack_string
    element = stack[end]
    splice!(stack, length(stack)-1:length(stack)-2, [element])
    return true
end

function op_size(stack::StackType)
    @assert length(stack) ≥ 1 short_stack_string
    push!(stack, encode_num(length(stack[end])))
    return true
end

function op_equal(stack::StackType)
    @assert length(stack) ≥ 2 short_stack_string
    element1 = pop!(stack)
    element2 = pop!(stack)
    if element1 == element2
        push!(stack, encode_num(1))
    else
        push!(stack, encode_num(0))
    end
    return true
end

function op_equalverify(stack::StackType)
    return op_equal(stack) && op_verify(stack)

end

function op_1add(stack::StackType)
    @assert length(stack) ≥ 1 short_stack_string
    element = decode_num(pop!(stack))
    push!(stack, encode_num(element + 1))
    return true
end
function op_1sub(stack::StackType)
    @assert length(stack) ≥ 1 short_stack_string
    element = decode_num(pop!(stack))
    push!(stack, encode_num(element - 1))
    return true
end

function op_negate(stack::StackType)
    @assert length(stack) ≥ 1 short_stack_string
    element = decode_num(pop!(stack))
    push!(stack, encode_num(- element))
    return true
end

function op_abs(stack::StackType)
    @assert length(stack) ≥ 1 short_stack_string
    element = decode_num(pop!(stack))
    push!(stack, encode_num(abs(element)))
    return true
end

function op_not(stack::StackType)
    @assert length(stack) ≥ 1 short_stack_string
    element = decode_num(pop!(stack))
    if element == 0
        push!(stack, encode_num(1))
    else
        push!(stack, encode_num(0))
    end
    return true
end

function op_0notequal(stack::StackType)
    @assert length(stack) ≥ 1 short_stack_string
    element = pop!(stack)
    if decode_num((element)) == 0
        push!(stack, encode_num(0))
    else
        push!(stack, encode_num(1))
    end
    return true
end

function op_add(stack::StackType)
    @assert length(stack) ≥ 2 short_stack_string
    element1 = decode_num(pop!(stack))
    element2 = decode_num(pop!(stack))
    push!(stack, encode_num(element2 + element1))
    return true
end

function op_sub(stack::StackType)
    @assert length(stack) ≥ 2 short_stack_string
    element1 = decode_num(pop!(stack))
    element2 = decode_num(pop!(stack))
    push!(stack, encode_num(element2 - element1))
    return true
end

function op_booland(stack::StackType)
    @assert length(stack) ≥ 2 short_stack_string
    element1 = decode_num(pop!(stack))
    element2 = decode_num(pop!(stack))
    if element1 == 1 && element2 == 1
        push!(stack, encode_num(1))
    else
        push!(stack, encode_num(0))
    end
    return true
end

function op_boolor(stack::StackType)
    @assert length(stack) ≥ 2 short_stack_string
    element1 = decode_num(pop!(stack))
    element2 = decode_num(pop!(stack))
    if element1 == 1 || element2 == 1
        push!(stack, encode_num(1))
    else
        push!(stack, encode_num(0))
    end
    return true
end

function op_numequal(stack::StackType)
    @assert length(stack) ≥ 2 short_stack_string
    element1 = decode_num(pop!(stack))
    element2 = decode_num(pop!(stack))
    if element1 == element2
        push!(stack, encode_num(1))
    else
        push!(stack, encode_num(0))
    end
    return true
end

function op_numequalverify(stack::StackType)
    return op_numequal(stack) && op_verify(stack)
end

function op_numnotequal(stack::StackType)
    @assert length(stack) ≥ 2 short_stack_string
    element1 = decode_num(pop!(stack))
    element2 = decode_num(pop!(stack))
    if element1 != element2
        push!(stack, encode_num(1))
    else
        push!(stack, encode_num(0))
    end
    return true
end

function op_lessthan(stack::StackType)
    @assert length(stack) ≥ 2 short_stack_string
    element1 = decode_num(pop!(stack))
    element2 = decode_num(pop!(stack))
    if element2 < element1
        push!(stack, encode_num(1))
    else
        push!(stack, encode_num(0))
    end
    return true
end

function op_greaterthan(stack::StackType)
    @assert length(stack) ≥ 2 short_stack_string
    element1 = decode_num(pop!(stack))
    element2 = decode_num(pop!(stack))
    if element2 > element1
        push!(stack, encode_num(1))
    else
        push!(stack, encode_num(0))
    end
    return true
end

function op_lessthanorequal(stack::StackType)
    @assert length(stack) ≥ 2 short_stack_string
    element1 = decode_num(pop!(stack))
    element2 = decode_num(pop!(stack))
    if element2 <= element1
        push!(stack, encode_num(1))
    else
        push!(stack, encode_num(0))
    end
    return true
end

function op_greaterthanorequal(stack::StackType)
    @assert length(stack) ≥ 2 short_stack_string
    element1 = decode_num(pop!(stack))
    element2 = decode_num(pop!(stack))
    if element2 >= element1
        push!(stack, encode_num(1))
    else
        push!(stack, encode_num(0))
    end
    return true
end

function op_min(stack::StackType)
    @assert length(stack) ≥ 2 short_stack_string
    element1 = decode_num(pop!(stack))
    element2 = decode_num(pop!(stack))
    if element1 < element2
        push!(stack, encode_num(element1))
    else
        push!(stack, encode_num(element2))
    end
    return true
end

function op_max(stack::StackType)
    @assert length(stack) ≥ 2 short_stack_string
    element1 = decode_num(pop!(stack))
    element2 = decode_num(pop!(stack))
    if element1 > element2
        push!(stack, encode_num(element1))
    else
        push!(stack, encode_num(element2))
    end
    return true
end

function op_within(stack::StackType)
    @assert length(stack) ≥ 3 short_stack_string
    maxlimit = decode_num(pop!(stack))
    minlimit = decode_num(pop!(stack))
    element = decode_num(pop!(stack))
    if element >= minlimit && element < maxlimit
        push!(stack, encode_num(1))
    else
        push!(stack, encode_num(0))
    end
    return true
end

"""
Return RIPEMD160 hash of top item
"""
function op_ripemd160(stack::StackType)
    @assert length(stack) ≥ 1 short_stack_string
    element = pop!(stack)
    push!(stack, ripemd160(element))
    return true
end

"""
Return SHA1 hash of top item
"""
function op_sha1(stack::StackType)
    @assert length(stack) ≥ 1 short_stack_string
    element = pop!(stack)
    push!(stack, sha1(element))
    return true
end

"""
Return SHA256 hash of top item
"""
function op_sha256(stack::StackType)
    @assert length(stack) ≥ 1 short_stack_string
    element = pop!(stack)
    push!(stack, sha256(element))
    return true
end

"""
Return RIPEMD160(SHA256(x)) hash of top item
"""
function op_hash160(stack::StackType)
    @assert length(stack) ≥ 1 short_stack_string
    element = pop!(stack)
    push!(stack, ripemd160(sha256(element)))
    return true
end

"""
Return SHA256(SHA256(x)) hash of top item
"""
function op_hash256(stack::StackType)
    @assert length(stack) ≥ 1 short_stack_string
    element = pop!(stack)
    push!(stack, sha256(sha256(element)))
    return true
end

"""
Pop a public key and signature and validate the signature for the transaction's
hashed data, return TRUE if matching
"""
function op_checksig(stack::StackType, z::Integer)
    @assert length(stack) ≥ 2 short_stack_string
    sec_pubkey = pop!(stack)
    der_signature = pop!(stack)[1:end-1]
    point = Secp256k1.Point(sec_pubkey)
    sig = Signature(der_signature)
    if ECDSA.verify(point, z, sig)
        push!(stack, encode_num(1))
    else
        push!(stack, encode_num(0))
    end
    return true

end

"""
Same as CHECKSIG, then OP_VERIFY to halt if not TRUE
"""
function op_checksigverify(stack::StackType, z::Integer)
    return op_checksig(stack, z) && op_verify(stack)
end

"""
Run CHECKSIG for each pair of signature and public key provided. All must match.
Bug in implementation pops an extra value, prefix with OP_NOP as workaround
"""
function op_checkmultisig(stack::StackType, z::Integer)
    @assert length(stack) ≥ 1 short_stack_string
    n = decode_num(pop!(stack))
    if length(stack) < (n + 1)
        return false
    end
    sec_pubkeys = []
    for sec in 1:n
        push!(sec_pubkeys, pop!(stack))
    end
    m = decode_num(pop!(stack))
    if length(stack) < (m + 1)
        return false
    end
    der_signatures = []
    for der in 1:m
        # signature is assumed to be using SIGHASH_ALL
        push!(der_signatures, pop!(stack)[1:end-1])
    end
    # OP_CHECKMULTISIG bug
    pop!(stack)
    try
        points = Secp256k1.Point[]
        for sec in sec_pubkeys
            push!(points, Secp256k1.Point(sec))
        end
        sigs = Signature[]
        for der in der_signatures
            push!(sigs, Signature(der))
        end
        for sig in sigs
            if length(points) == 0
                append!(stack, encode_num(0))
                return true
            end
            while length(points) > 0
                point = popfirst!(points)
                if ECDSA.verify(point, z, sig)
                    break
                end
                println("not valid sig")
            end
        end
        push!(stack, encode_num(1))
    catch
        return false
    end
    return true
end

"""
Same as CHECKMULTISIG, then OP_VERIFY to halt if not TRUE
"""
function op_checkmultisigverify(stack::StackType, z::Integer)
    return op_checkmultisig(stack, z) && op_verify(stack)
end

"""
Marks transaction as invalid if the top stack item is greater than the
transaction's nLockTime field, otherwise script eveluation continues as though
an OP_NOP was executed. Transaction is also invalid if 1. the stack is empty;
or 2. the top stack item is negative; or 3. the top stack item is greater than
or equal to 500000000 while the transaction's nLockTime field is less than
500000000, or vice versa; or 4. the input's nSequence field is equal to
0xffffffff. The precise semantics are described in BIP-65
"""
function op_checklocktimeverify(stack::StackType, locktime::Integer, sequence::Integer)
    if sequence == 0xffffffff || length(stack) < 1
         return false
    end
    element = decode_num(stack[end])
    if element < 0
         return false
    elseif element < 500000000 && locktime > 500000000
         return false
    elseif locktime < element
         return false
    else
        return true
    end
end

"""
Marks transaction as invalid if the relative lock time of the input (enforced by
BIP-68 with nSequence) is not equal to or longer than the value of the top stack
item. The precise semantics are described in BIP-112
"""
function op_checksequenceverify(stack::StackType, version::Integer, sequence::Integer)
    error("not implemented")
end

const OP_CODE_FUNCTIONS = Dict([
    (0,  op_0),
    (79,  op_1negate),
    (81,  op_1),
    (82,  op_2),
    (83,  op_3),
    (84,  op_4),
    (85,  op_5),
    (86,  op_6),
    (87,  op_7),
    (88,  op_8),
    (89,  op_9),
    (90,  op_10),
    (91,  op_11),
    (92,  op_12),
    (93,  op_13),
    (94,  op_14),
    (95,  op_15),
    (96,  op_16),
    (97,  op_nop),
    (99,  op_if),
    (100,  op_notif),
    (105,  op_verify),
    (106,  op_return),
    (107,  op_toaltstack),
    (108,  op_fromaltstack),
    (109,  op_2drop),
    (110,  op_2dup),
    (111,  op_3dup),
    (112,  op_2over),
    (113,  op_2rot),
    (114,  op_2swap),
    (115,  op_ifdup),
    (116,  op_depth),
    (117,  op_drop),
    (118,  op_dup),
    (119,  op_nip),
    (120,  op_over),
    (121,  op_pick),
    (122,  op_roll),
    (123,  op_rot),
    (124,  op_swap),
    (125,  op_tuck),
    (130,  op_size),
    (135,  op_equal),
    (136,  op_equalverify),
    (139,  op_1add),
    (140,  op_1sub),
    (143,  op_negate),
    (144,  op_abs),
    (145,  op_not),
    (146,  op_0notequal),
    (147,  op_add),
    (148,  op_sub),
    (154,  op_booland),
    (155,  op_boolor),
    (156,  op_numequal),
    (157,  op_numequalverify),
    (158,  op_numnotequal),
    (159,  op_lessthan),
    (160,  op_greaterthan),
    (161,  op_lessthanorequal),
    (162,  op_greaterthanorequal),
    (163,  op_min),
    (164,  op_max),
    (165,  op_within),
    (166,  op_ripemd160),
    (167,  op_sha1),
    (168,  op_sha256),
    (169,  op_hash160),
    (170,  op_hash256),
    (172,  op_checksig),
    (173,  op_checksigverify),
    (174,  op_checkmultisig),
    (175,  op_checkmultisigverify),
    (176,  op_nop),
    (177,  op_checklocktimeverify),
    (178,  op_checksequenceverify),
    (179,  op_nop),
    (180,  op_nop),
    (181,  op_nop),
    (182,  op_nop),
    (183,  op_nop),
    (184,  op_nop),
    (185,  op_nop)])

const OP_CODE_NAMES = Dict([
    (0, "OP_0"),
    (76, "OP_PUSHDATA1"),
    (77, "OP_PUSHDATA2"),
    (78, "OP_PUSHDATA4"),
    (79, "OP_1NEGATE"),
    (81, "OP_1"),
    (82, "OP_2"),
    (83, "OP_3"),
    (84, "OP_4"),
    (85, "OP_5"),
    (86, "OP_6"),
    (87, "OP_7"),
    (88, "OP_8"),
    (89, "OP_9"),
    (90, "OP_10"),
    (91, "OP_11"),
    (92, "OP_12"),
    (93, "OP_13"),
    (94, "OP_14"),
    (95, "OP_15"),
    (96, "OP_16"),
    (97, "OP_NOP"),
    (99, "OP_IF"),
    (100, "OP_NOTIF"),
    (103, "OP_ELSE"),
    (104, "OP_ENDIF"),
    (105, "OP_VERIFY"),
    (106, "OP_RETURN"),
    (107, "OP_TOALTSTACK"),
    (108, "OP_FROMALTSTACK"),
    (109, "OP_2DROP"),
    (110, "OP_2DUP"),
    (111, "OP_3DUP"),
    (112, "OP_2OVER"),
    (113, "OP_2ROT"),
    (114, "OP_2SWAP"),
    (115, "OP_IFDUP"),
    (116, "OP_DEPTH"),
    (117, "OP_DROP"),
    (118, "OP_DUP"),
    (119, "OP_NIP"),
    (120, "OP_OVER"),
    (121, "OP_PICK"),
    (122, "OP_ROLL"),
    (123, "OP_ROT"),
    (124, "OP_SWAP"),
    (125, "OP_TUCK"),
    (130, "OP_SIZE"),
    (135, "OP_EQUAL"),
    (136, "OP_EQUALVERIFY"),
    (139, "OP_1ADD"),
    (140, "OP_1SUB"),
    (143, "OP_NEGATE"),
    (144, "OP_ABS"),
    (145, "OP_NOT"),
    (146, "OP_0NOTEQUAL"),
    (147, "OP_ADD"),
    (148, "OP_SUB"),
    (154, "OP_BOOLAND"),
    (155, "OP_BOOLOR"),
    (156, "OP_NUMEQUAL"),
    (157, "OP_NUMEQUALVERIFY"),
    (158, "OP_NUMNOTEQUAL"),
    (159, "OP_LESSTHAN"),
    (160, "OP_GREATERTHAN"),
    (161, "OP_LESSTHANOREQUAL"),
    (162, "OP_GREATERTHANOREQUAL"),
    (163, "OP_MIN"),
    (164, "OP_MAX"),
    (165, "OP_WITHIN"),
    (166, "OP_RIPEMD160"),
    (167, "OP_SHA1"),
    (168, "OP_SHA256"),
    (169, "OP_HASH160"),
    (170, "OP_HASH256"),
    (171, "OP_CODESEPARATOR"),
    (172, "OP_CHECKSIG"),
    (173, "OP_CHECKSIGVERIFY"),
    (174, "OP_CHECKMULTISIG"),
    (175, "OP_CHECKMULTISIGVERIFY"),
    (176, "OP_NOP1"),
    (177, "OP_CHECKLOCKTIMEVERIFY"),
    (178, "OP_CHECKSEQUENCEVERIFY"),
    (179, "OP_NOP4"),
    (180, "OP_NOP5"),
    (181, "OP_NOP6"),
    (182, "OP_NOP7"),
    (183, "OP_NOP8"),
    (184, "OP_NOP9"),
    (185, "OP_NOP10")])


text:
@testset "create_geometry.jl" begin

end

ext: jl


text:
abstract type division_type end
struct Sizer <: division_type end
struct Adder <: division_type end
struct Timer <: division_type end
struct Sizer_Timer <: division_type end
struct Adder_Timer <: division_type end

##### calculate probability of cell division
@inline calc_division(::Sizer, Sz, iS, t, stp, reg, stp2, reg2, P) = P * (tanh(stp*(Sz - reg)) + 1.0)

@inline calc_division(::Adder, Sz, iS, t, stp, reg, stp2, reg2, P) = P * (tanh(stp*(Sz - iS - reg)) + 1.0)

@inline calc_division(::Timer, Sz, iS, t, stp, reg, stp2, reg2, P) = P * (tanh(stp*(t%86400/3600 - reg)) + 1.0)

@inline calc_division(::Sizer_Timer, Sz, iS, t, stp, reg, stp2, reg2, P) =
                        P * (tanh(stp*(Sz - reg)) + 1.0) * (tanh(stp2*(t%86400/3600 - reg2)) + 1.0)

@inline calc_division(::Adder_Timer, Sz, iS, t, stp, reg, stp2, reg2, P) =
                        P * (tanh(stp*(Sz - iS - reg)) + 1.0) * (tanh(stp2*(t%86400/3600 - reg2)) + 1.0)

@inline function divide_type(dvid_type)
    if dvid_type == 1
        return Sizer()
    elseif dvid_type == 2
        return Adder()
    elseif dvid_type == 3
        return Timer()
    elseif dvid_type == 4
        return Sizer_Timer()
    elseif dvid_type == 5
        return Adder_Timer()
    else
        throw(ArgumentError("Wrong cell division type, must be in 1 to 5"))
    end
end

@kernel function calc_dvid_kernel!(plank, dvid_type, p, t)
    i = @index(Global)
    @inbounds plank.dvid[i] = calc_division(dvid_type, plank.Sz[i], plank.iS[i], t,
                                            p.dvid_stp, p.dvid_reg, p.dvid_stp2, p.dvid_reg2, p.dvid_P) *
                                            isless(2.0 * p.Cquota * p.Nsuper, plank.Bm[i]) * plank.ac[i]
end
function calc_dvid!(plank, dvid_type, p, t, arch)
    kernel! = calc_dvid_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, dvid_type, p, t)
    wait(device(arch), event)
    return nothing
end

@kernel function calc_MM_dvid_kernel!(plank, p)
    i = @index(Global)
    @inbounds plank.dvid[i] = p.dvid_P * isless(2.0, plank.DNA[i]/(p.C_DNA*p.Nsuper))
end
function calc_MM_dvid!(plank, p, arch)
    kernel! = calc_MM_dvid_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, p)
    wait(device(arch), event)
    return nothing
end

##### calculate the probability of grazing
##### quadratic grazing
@kernel function calc_graz_quadratic_kernel!(plank, nuts, P)
    i = @index(Global)
    @inbounds plank.graz[i] = nuts.pop[i] * P * plank.ac[i]
end
function calc_graz_quadratic!(plank, nuts, P, arch)
    kernel! = calc_graz_quadratic_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, nuts, P)
    wait(device(arch), event)
    return nothing
end

##### calculate the probability of mortality
@kernel function calc_mort_kernel!(plank, p)
    i = @index(Global)
    @inbounds plank.mort[i] = p.mort_P * (tanh(6.0*(p.mort_reg - plank.Sz[i])) + 1.0) * plank.ac[i]
end
function calc_mort!(plank, p, arch)
    kernel! = calc_mort_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, p)
    wait(device(arch), event)
    return nothing
end

@kernel function calc_MM_mort_kernel!(plank, p)
    i = @index(Global)
    @inbounds plank.mort[i] = p.mort_P * (tanh(1.0*(plank.age[i] - p.mort_reg)) + 1.0) * plank.ac[i]
end
function calc_MM_mort!(plank, p, arch)
    kernel! = calc_MM_mort_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, p)
    wait(device(arch), event)
    return nothing
end

##### generate the random results from probabilities of grazing, mortality and cell division
function get_probability!(plank, rnd, ΔT, arch)
    ##### generate random numbers (0,1)
    rand!(rng_type(arch), rnd.x)
    rand!(rng_type(arch), rnd.y)
    rand!(rng_type(arch), rnd.z)

    ##### compare the random number with the given probability (per time step or per hour whichever is shorter)
    ##### return 1 if random number is smaller
    @inbounds plank.graz .= isless.(rnd.x, plank.graz .* ΔT .* min(10,3600÷ΔT)) .* plank.ac
    @inbounds plank.mort .= isless.(rnd.y, plank.mort .* ΔT .* min(10,3600÷ΔT)) .* plank.ac
    @inbounds plank.dvid .= isless.(rnd.z, plank.dvid .* ΔT .* min(10,3600÷ΔT)) .* plank.ac
    return nothing
end


text:
using DelimitedFiles, Match
using aoc2019.computer

input = joinpath(@__DIR__, "input")
p = computer.load_program(input)

stdin, stdout = computer.io()
task = @async computer.run(copy(p), stdin, stdout)
while !istaskdone(task); yield(); end
close(stdout)

bits = [x for x in stdout]
screen = Dict{Pair, Int}()

while length(bits) > 0
    tile = pop!(bits)
    y = pop!(bits)
    x = pop!(bits)
    screen[Pair(x, y)] = tile
end

p1 = length(filter(x -> x == 2, collect(values(screen))))
@assert p1 == 335

stdin, stdout = computer.io()
program = copy(p)
## freeplay mode
program[1] = 2
task = @async computer.run(program, stdin, stdout)
sleep(1)

function play(task, stdin, stdout)
    screen = Dict{Pair, Int}()
    score = 0
    while true
        i = 1
        x, y, tile = nothing, nothing, nothing
        while isready(stdout)
            val = take!(stdout)
            @match mod(i, 1:3) begin
                1 => begin x = val end
                2 => begin y = val end
                3 => begin tile = val end
            end
            if x == -1 && y == 0 && tile !== nothing
                score = tile
                x, y, tile = nothing, nothing, nothing
            elseif x !== nothing && y !== nothing && tile !== nothing
                screen[Pair(x, y)] = tile
                x, y, tile = nothing, nothing, nothing
            end
            i += 1
        end
        ball = filter(bit -> bit.second == 4, screen)
        paddle = filter(bit -> bit.second == 3, screen)
        (bx, by) = first(collect(keys(ball)))
        (px, py) = first(collect(keys(paddle)))
        if bx > px
            put!(stdin, 1)
        elseif bx < px
            put!(stdin, -1)
        else
            put!(stdin, 0)
        end
        if istaskdone(task)
            break
        else
            yield()
        end
    end
    return score
end

p2 = play(task, stdin, stdout)
@assert p2 == 15706

print("-----------------------------------------------------------------------\n")
print("care package -- part one\n    total blocks : $p1\n")
print("care package -- part two\n    final score  : $p2\n")
print("-----------------------------------------------------------------------\n")


text:
module WeaveUtil

using Weave

function getJMD(to::IOStream, from::IOStream)
	r = r"input\((.+)\)"
	for i in eachline(from)
		if (m = match(r, i)) != nothing
			open(m.captures[1], "r") do subfile
				getJMD(to, subfile)
			end
		else
			write(to, i, '\n')
		end
	end
end

export buildJMD
function buildJMD(inputpath::String;
		doctype = "pandoc2pdf", docname = split(inputpath, '.')[end - 1] * ".pdf")
	open(inputpath, "r") do from
		open("_tmp.jmd", "w") do to
			getJMD(to, from)
		end
	end

	weave("_tmp.jmd"; doctype = doctype)

	mv("_tmp.pdf", docname; force = true)
	rm("_tmp.jmd")
end

function getTex(topath::AbstractString, frompath::AbstractString)
	path::String = splitdir(frompath)[1]
	r1::Regex = r"(.*)\\input\{(.+)\}"
	r2::Regex = r"[^\\]%"
	open(topath, "w") do to
		function rec(from::IOStream)
			for i in eachline(from)
				if (m = match(r1, i)) != nothing &&
						!startswith(m.captures[1], "%") &&
						match(r2, m.captures[1]) == nothing
					open(io->rec(io), joinpath(path, m.captures[2]), "r")
				else
					write(to, i, '\n')
				end
			end
		end
		open(io->rec(io), frompath)
	end
end

setflag(doc::AbstractString, flg::AbstractString) =
		open(io->write(io, '\0'), doc * flg, "w")
#

export donetx
donetx(doc) = setflag(doc, "")

export wvtx
function wvtx(doc::AbstractString)
	getTex(doc * "temp.tex", doc * ".tex")

	weave(doc * "temp.tex", out_path = doc * ".out.tex")
	rm(doc * "temp.tex")
	donetx(doc)
end

fileext(path::AbstractString) = splitext(path)[2]

export clntx
function clntx(doc::AbstractString)
	path = splitdir(doc)[1]
	blacklist = [
		".out.tex",
		".log",
		".xdv",
		".fls",
		".synctex.gz",
		".fdb_latexmk",
		".aux"
	]

	for s in blacklist
		rm(doc * s)
	end
end

end


text:
export ft,ift, rft, irft
export ffts, ffts!, iffts, rffts, irffts
export fft2d,ifft2d, rfft2d, irfft2d
export ft2d,ift2d, rft2d, irft2d
export ffts2d, ffts2d!, iffts2d, rffts2d, irffts2d

"""
    optional_collect(a)

Only collects certain arrays, for a pure `Array` there is no collect
and it returns simply `a`.
"""
# collect
optional_collect(a::AbstractArray) = collect(a)
# no need to collect
optional_collect(a::Array) = a

# for CircShiftedArray we only need collect if shifts is non-zero
function optional_collect(csa::ShiftedArrays.CircShiftedArray)
    if all(iszero.(csa.shifts))
        return optional_collect(parent(csa))
    else
        return collect(csa)
    end
end

"""
    ffts(A [, dims])

Result is semantically equivalent to `fftshift(fft(A, dims), dims)`
However, the shift is done with `ShiftedArrays` and therefore doesn't allocate memory.

See also: [`ft`](@ref ift), [`ift`](@ref ift), [`rft`](@ref rft), [`irft`](@ref irft),
          [`ffts`](@ref ffts),  [`iffts`](@ref iffts),  [`ffts!`](@ref ffts!), [`rffts`](@ref rffts), [`irffts`](@ref irffts),
"""
function ffts(mat::AbstractArray{T, N}, dims=ntuple(identity, Val(N))) where {T, N}
    return fftshift_view(fft(mat, dims), dims)
end

"""
    ffts!(A [, dims])

Result is semantically equivalent to `fftshift(fft!(A, dims), dims)`.
`A` is in-place modified.
However, the shift is done with `ShiftedArrays` and therefore doesn't allocate memory.

See also: [`ft`](@ref ift), [`ift`](@ref ift), [`rft`](@ref rft), [`irft`](@ref irft),
          [`ffts`](@ref ffts),  [`iffts`](@ref iffts),  [`ffts!`](@ref ffts!), [`rffts`](@ref rffts), [`irffts`](@ref irffts),
"""
function ffts!(mat::AbstractArray{T, N}, dims=ntuple(identity, Val(N))) where {T, N}
    return fftshift_view(fft!(mat, dims), dims)
end

"""
    iffts(A [, dims])

Result is semantically equivalent to `ifft(ifftshift(A, dims), dims)`.
`A` is in-place modified.
However, the shift is done with `ShiftedArrays` and therefore doesn't allocate memory.

See also: [`ft`](@ref ift), [`ift`](@ref ift), [`rft`](@ref rft), [`irft`](@ref irft),
          [`ffts`](@ref ffts),  [`iffts`](@ref iffts),  [`ffts!`](@ref ffts!), [`rffts`](@ref rffts), [`irffts`](@ref irffts),
"""
function iffts(mat::AbstractArray{T, N}, dims=ntuple(identity, Val(N))) where {T, N}
    return ifft(optional_collect(ifftshift_view(mat, dims)), dims)
end

"""
    rffts(A [, dims])

Calculates a `rfft(A, dims)` and then shift the frequencies to the center.
`dims[1]` is not shifted, because there is no negative and positive frequency.
The shift is done with `ShiftedArrays` and therefore doesn't allocate memory.

See also: [`ft`](@ref ift), [`ift`](@ref ift), [`rft`](@ref rft), [`irft`](@ref irft),
          [`ffts`](@ref ffts),  [`iffts`](@ref iffts),  [`ffts!`](@ref ffts!), [`rffts`](@ref rffts), [`irffts`](@ref irffts),
"""
function rffts(mat::AbstractArray{T, N}, dims=ntuple(identity, Val(N))) where {T, N}
    rfftshift_view(rfft(mat, dims), dims);
end

"""
    irffts(A, d, [, dims])

Calculates a `irfft(A, d, dims)` and then shift the frequencies back to the corner.
`dims[1]` is not shifted, because there is no negative and positive frequency.
The shift is done with `ShiftedArrays` and therefore doesn't allocate memory.

See also: [`ft`](@ref ift), [`ift`](@ref ift), [`rft`](@ref rft), [`irft`](@ref irft),
          [`ffts`](@ref ffts),  [`iffts`](@ref iffts),  [`ffts!`](@ref ffts!), [`rffts`](@ref rffts), [`irffts`](@ref irffts),
"""
function irffts(mat::AbstractArray{T, N}, d::Int, dims=ntuple(identity, Val(N))) where {T, N}
    irfft(optional_collect(irfftshift_view(mat, dims)), d, dims)
end

"""
    ft(A [, dims])

Digital Fourier-transformation centered in both spaces.
The result is semantically equivalent to `fftshift(fft(ifftshift(A, dims), dims), dims)`
This is a digital Fourier transformation with both coordinate systems in real and Fourier-space being
centered at position CtrFT == size÷2+1

The following identities are true:
```jldoctest
julia> sz = (5,5)
(5, 5)

julia> ft(ones(sz)) ≈ prod(sz) .* δ(sz)
true

julia> ft(δ(sz)) ≈ ones(sz)
true
```

See also: [`ft`](@ref ift), [`ift`](@ref ift), [`rft`](@ref rft), [`irft`](@ref irft),
          [`ffts`](@ref ffts),  [`iffts`](@ref iffts),  [`ffts!`](@ref ffts!), [`rffts`](@ref rffts), [`irffts`](@ref irffts),

"""
function ft(mat::AbstractArray{T, N}, dims=ntuple(identity, Val(N))) where {T, N}
    return fftshift_view(fft(optional_collect(ifftshift_view(mat, dims)), dims), dims)
end

"""
    ift(A [, dims])

Digital inverse Fourier-transformation centered in both spaces.
The result is semantically equivalent to `fftshift(ifft(ifftshift(A, dims), dims), dims)`
This is a digital Fourier transformation with both coordinate systems in real and Fourier-space being
centered at position CtrFT == size÷2+1

The following identities are true:

```jldoctest
julia> sz = (5,6,7)
(5, 6, 7)

julia> ift(ones(sz)) ≈ δ(sz)
true

julia> ift(δ(sz)) ≈ ones(sz) ./ prod(sz)
true
```

See also: [`ft`](@ref ift), [`ift`](@ref ift), [`rft`](@ref rft), [`irft`](@ref irft),
          [`ffts`](@ref ffts),  [`iffts`](@ref iffts),  [`ffts!`](@ref ffts!), [`rffts`](@ref rffts), [`irffts`](@ref irffts),

"""
function ift(mat::AbstractArray{T, N}, dims=ntuple(identity, Val(N))) where {T, N}
    return fftshift_view(ifft(optional_collect(ifftshift_view(mat, dims)), dims), dims)  # is faster than the exp mulitplication
end

"""
    rft(A [, dims])

Digital real-valued Fourier-transformation centered in both spaces.
The result is semantically equivalent to `fftshift(rfft(ifftshift(A, dims), dims), dims)`
This is a digital Fourier transformation with the coordinate systems in real space centered at CtrFT == size÷2+1
and in (half) Fourier-space being centered at CtrRFT == setindex(size÷2 +1,1,1).

The following identities are true:
```jldoctest
julia> sz = (6,6)
(6, 6)

julia> rft(δ(sz)) ≈ ones(rft_size(sz))
true
```

See also: [`ft`](@ref ift), [`ift`](@ref ift), [`rft`](@ref rft), [`irft`](@ref irft),
          [`ffts`](@ref ffts),  [`iffts`](@ref iffts),  [`ffts!`](@ref ffts!), [`rffts`](@ref rffts), [`irffts`](@ref irffts),

"""
function rft(mat::AbstractArray{T, N}, dims=ntuple(identity, Val(N))) where {T, N}
    return rfftshift_view(rfft(optional_collect(ifftshift_view(mat, dims)), dims), dims);
end

"""
    irft(A, d, [, dims])

Digital real-valued inverse Fourier-transformation centered in both spaces.
The result is semantically equivalent to `fftshift(irfft(ifftshift(A, dims), dims), dims)`
This is a digital Fourier transformation with the coordinate systems in real space centered at CtrFT == size÷2+1
and in (half) Fourier-space being centered at CtrRFT == setindex(size÷2 +1,1,1).
Note that the size `d` of the first transform direction [1] is a required argument.

The following identities are true:
```jldoctest
julia> sz = (6,6)
(6, 6)

julia> irft(ones(rft_size(sz)),sz[1]) ≈ δ(sz)
true
```

See also: [`ft`](@ref ift), [`ift`](@ref ift), [`rft`](@ref rft), [`irft`](@ref irft),
          [`ffts`](@ref ffts),  [`iffts`](@ref iffts),  [`ffts!`](@ref ffts!), [`rffts`](@ref rffts), [`irffts`](@ref irffts),

"""
function irft(mat::AbstractArray{T, N}, d::Int, dims=ntuple(identity, Val(N))) where {T, N}
    fftshift_view(irfft(optional_collect(irfftshift_view(optional_collect(mat), dims)), d, dims), dims);
end

## Short-hand versions of the ft functions
"""
    ft2d(mat::AbstractArray{T, N}) where {T, N}

Only over `dims=(1,2)`.
"""
function ft2d(mat::AbstractArray{T, N}) where {T, N}
    ft(mat,(1,2))
end
"""
    ift2d(mat::AbstractArray{T, N}) where {T, N}

Only over `dims=(1,2)`.
"""
function ift2d(mat::AbstractArray{T, N}) where {T, N}
    ift(mat,(1,2))
end
"""
    rft2d(mat::AbstractArray{T, N}) where {T, N}

Only over `dims=(1,2)`.
"""
function rft2d(mat::AbstractArray{T, N}) where {T, N}
    rft(mat,(1,2))
end

"""
    irft2d(mat::AbstractArray{T, N}, d) where {T, N}

Only over `dims=(1,2)`.
"""
function irft2d(mat::AbstractArray{T, N}, d::Int) where {T, N}
    irft(mat,d,(1,2))
end

"""
    ft2d(mat::AbstractArray{T, N}) where {T, N}

Only over `dims=(1,2)`.
"""
function ffts2d(mat::AbstractArray{T, N}) where {T, N}
    ffts(mat,(1,2))
end

"""
    fft2ds!(mat::AbstractArray{T, N}) where {T, N}

Only over `dims=(1,2)`.
"""
function ffts2d!(mat::AbstractArray{T, N}) where {T, N}
    ffts!(mat,(1,2))
end

"""
    iffts2d(mat::AbstractArray{T, N}) where {T, N}

Only over `dims=(1,2)`.
"""
function iffts2d(mat::AbstractArray{T, N}) where {T, N}
    iffts(mat,(1,2))
end
"""
    rffts2d(mat::AbstractArray{T, N}) where {T, N}

Only over `dims=(1,2)`.
"""
function rffts2d(mat::AbstractArray{T, N}) where {T, N}
    rffts(mat,(1,2))
end

"""
    riffts2d(mat::AbstractArray{T, N}, d) where {T, N}

Only over `dims=(1,2)`.
"""
function irffts2d(mat::AbstractArray{T, N}, d::Int) where {T, N}
    irffts(mat,d,(1,2))
end

## Short-hand versions of the fft functions
"""
    fft2d(mat::AbstractArray{T, N}) where {T, N}

Only over `dims=(1,2)`.
"""
function fft2d(mat::AbstractArray{T, N}) where {T, N}
    fft(mat,(1,2))
end
"""
    ifft2d(mat::AbstractArray{T, N}) where {T, N}

Only over `dims=(1,2)`.
"""
function ifft2d(mat::AbstractArray{T, N}) where {T, N}
    ifft(mat,(1,2))
end

"""
    rfft2d(mat::AbstractArray{T, N}) where {T, N}

Only over `dims=(1,2)`.
"""
function rfft2d(mat::AbstractArray{T, N}) where {T, N}
    rfft(mat,(1,2))
end

"""
    irfft2d(mat::AbstractArray{T, N}, d) where {T, N}

Only over `dims=(1,2)`.
"""
function irfft2d(mat::AbstractArray{T, N}, d::Int) where {T, N}
    irfft(mat,d,(1,2))
end


text:
struct TERN <: DataOneV2 end

base_url(::TERN) = "https://dataone.tern.org.au/mn/v2/object/aekos.org.au/collection/nsw.gov.au/nsw_atlas/"

function get_urls(repo::TERN, page)
     links = eachmatch(sel"distribution online url", page.root)
     [text_only(links)]
end

function desc_(repo::TERN, mainpage)
    desc_ele = eachmatch(sel"description", mainpage.root)
    text_only(first(desc_ele))
end

function authors_(repo::TERN, mainpage)
    author_ele = eachmatch(sel"organizationName", mainpage.root)
    [text_only(i) for i in author_ele]
end

function pub_date(repo::TERN, mainpage)
    date_ele = eachmatch(sel"temporalCoverage rangeOfDates endDate", mainpage.root)
    Dates.Date(text_only(first(date_ele)), "yyyy-mm-dd")
end

function data_fullname(::TERN, mainpage)
    text_only(first(eachmatch(sel"project title", mainpage.root)))
end


text:
using ImageTransformations: center
using Compat

"""
    mask_circle!(::AbstractMatrix, npix; fill=0)
    mask_circle!(::AbstractArray, npix; fill=0)

In-place version of [`mask_circle`](@ref)
"""
function mask_circle!(arr::AbstractMatrix, npix; fill = 0)
    y, x = axes(arr)
    yc, xc = center(arr)
    d = @. sqrt((x' - xc)^2 + (y - yc)^2)
    @. arr[d < npix] = fill
    return arr
end

function mask_circle!(cube::AbstractArray{T, 3}, npix; fill=0) where T
    @inbounds for i in axes(cube, 1)
        slice = @view cube[i, :, :]
        mask_circle!(slice, npix; fill=fill)
    end
    return cube
end

"""
    mask_circle(::AbstractMatrix, npix; fill=0)
    mask_circle(::AbstractArray, npix; fill=0)

Mask the inner-circle of an image with radius `npix` with value `fill`. Note that the input type must be compatible with the fill value's type. If the input is a cube it will mask each frame individually.

# See Also
[`mask_circle!`](@ref)
"""
mask_circle(arr::AbstractMatrix, npix; fill = 0) = mask_circle!(deepcopy(arr), npix, fill = fill)
mask_circle(cube::AbstractArray{T, 3}, npix; fill = 0) where T = mask_circle!(deepcopy(cube), npix, fill=fill)

"""
    mask_annulus!(::AbstractMatrix, npix_in, npix_out; fill=NaN)

In-place version of [`mask_annulus`](@ref)
"""
function mask_annulus!(arr::AbstractMatrix, npix_in, npix_out; fill = 0)
    y, x = axes(arr)
    yc, xc = center(arr)
    d = @. sqrt((x' - xc)^2 + (y - yc)^2)
    @. arr[npix_in ≤ d < npix_out] = fill
    return arr
end

"""
    mask_annulus(::AbstractMatrix, npix_in, npix_out; fill=0)

Mask an annular region of an image with inner-radius `npix_in`, outer-radius `npix_out` with value `fill`.
Note that the input type must be compatible with the fill value's type.

# See Also
[`mask_annulus!`](@ref)
"""
mask_annulus(arr::AbstractMatrix, npix_in, npix_out; fill = 0) = mask_annulus!(deepcopy(arr), npix_in, npix_out, fill = fill)

"""
    get_annulus_segments(data, inner_radius, width, [nsegments];
                         theta_init=0, scale_factor=1, mode=:index)

Returns indices or values in segments of a centered annulus.

### Modes
* `:mask` - returns a positive boolean mask for indexing
* `:value` - returns the data values indexed with the boolean mask
* `:apply` - returns the input data weighted by the boolean mask
"""
function get_annulus_segments(data, inner_radius, width, nsegments; theta_init=0, scale_factor=1, mode=:mask)
    cy, cx = center(data)
    ys, xs = axes(data)
    dist = @. sqrt((xs' - cx)^2 + (ys - cy)^2)
    ϕ_stride = deg2rad(ceil(Int, 360 / nsegments))
    ϕ = @. atan(ys - cy, xs - cx)
    ϕ_rot = mod2pi.(ϕ)
    outer_radius = @. inner_radius + scale_factor * width

    masks = similar(data, Matrix{Bool}, nsegments)
    for i in 1:nsegments
        ϕ_start = deg2rad(theta_init) + (i - 1) * ϕ_stride
        ϕ_end = ϕ_start + ϕ_stride

        if ϕ_start < 2π && ϕ_end > 2π
            mask = @. (dist ≥ inner_radius) & (dist < outer_radius) & (ϕ_rot ≥ ϕ_start) & (ϕ_rot ≤ 2π) | (dist ≥ inner_radius) & (dist < outer_radius) & (ϕ_rot ≥ 0) & (ϕ_rot < ϕ_end - 2π)
        elseif ϕ_start ≥ 2π && ϕ_end > 2π
            mask = @. (dist ≥ inner_radius) & (dist < outer_radius) & (ϕ_rot ≥ ϕ_start - 2π) & (ϕ_rot < ϕ_end - 2π)
        else
            mask = @. (dist ≥ inner_radius) & (dist < outer_radius) & (ϕ_rot ≥ ϕ_start) & (ϕ_rot < ϕ_end)
        end
        @inbounds masks[i] = mask
    end

    return _convert_mask.(Val(mode), masks, (data,))
end

@compat get_annulus_segments(data, inner_radius, width; theta_init=0, scale_factor=1, mode=:mask) = get_annulus_segments(data, inner_radius, width, 1; theta_init=theta_init, scale_factor=scale_factor, mode=mode) |> only

_convert_mask(::Val{:mask}, mask, data) = mask
_convert_mask(::Val{:value}, mask, data) = data[mask]
_convert_mask(::Val{:apply}, mask, data) = data .* mask


text:
function lineal_S(M = 3)
    N = M + 1
    # S
    S = zeros(M,N)
    for i in 1:M
        S[i,i] = 1.0
        S[i,i+1] = -1.0
    end
    return S
end

function simple_lineal_MetNet(M)
    net = Dict()
    net[:S] = lineal_S(M)
    M,N = size(net[:S])
    net[:c] = zeros(N)
    net[:b] = zeros(M)
    net[:lb] = zeros(N)
    net[:ub] = ones(N)

    return MetNet(net; reshape = false)
end


text:
using DataFrames, CSV, Statistics, StatsBase

df = CSV.read("./preprocessing/input/land-change.csv")

# remove last (empty) column
select!(df, Not(:Model))
select!(df, Not(:Region))

df = df[.!ismissing.(df.Scenario),:]

df = stack(df, 4:size(df)[2])

rename!(df, [:year, :value, :scenario, :variable, :unit])

df[df.variable .== "Land Cover|Cropland",:value] -= df[df.variable .== "Land Cover|Cropland|Energy Crops",:value]

df[df.variable .== "Land Cover|Cropland",:variable] .= "Food and Feed Crops"
df[df.variable .== "Land Cover|Cropland|Energy Crops",:variable] .= "Bioenergy Crops"
df[df.variable .== "Land Cover|Pasture",:variable] .= "Pasture"
df[df.variable .== "Land Cover|Forest",:variable] .= "Forest"
df[df.variable .== "Land Cover|Other Natural Land",:variable] .= "Other Natural Land"

df = unstack(df, :year, :value)

CSV.write("./src/assets/data/land-change.csv", df)

max_stars_repo_name: SensesProject/land-use

max_issues_repo_name: SensesProject/land-use

max_forks_repo_name: SensesProject/land-use


text:
module Rain

using ..Ahorn, Maple

placements = Maple.Rain

function Ahorn.canFgBg(effect::Maple.Rain)
    return true, true
end

end


text:
module DataArraysBenchmark
using DataArrays, Benchmark

# seed rng for more consistent timings
Random.seed!(1776)

const TEST_NAMES = [
    "Vector",
    "DataVector No missing",
    "DataVector Half missing",
    "Matrix",
    "DataMatrix No missing",
    "DataMatrix Half missing"
]

function make_test_types(genfunc, sz)
    mat = genfunc(abs2(sz))
    na = shuffle!([trues(ifloor(abs2(sz)/2)), falses(iceil(abs2(sz)/2))])
    (
        mat,
        DataArray(mat),
        DataArray(mat, na),
        reshape(mat, sz, sz),
        DataArray(reshape(mat, sz, sz)),
        DataArray(reshape(mat, sz, sz), reshape(na, sz, sz))
    )
end

make_bools(x::NTuple{N}) where {N} = convert(Array{Bool, N}, bitrand(x...))
make_bools(x::Integer...) = make_bools(x)

macro perf(fn, replications, idx...)
    name = fn.args[1]
    idx = isempty(idx) ? (1:length(TEST_NAMES)) : idx[1]
    quote
        println($name)
        gc_disable()
        df = compare([let i=i; ()->$fn; end for i=$idx], $replications)
        gc_enable()
        gc()
        df[:Function] = TEST_NAMES[$idx]
        df[:Relative] = df[:Average]./df[1, :Average]
        println(df)
    end
end

const Float1 = make_test_types(rand, 1000)
const Float2 = make_test_types(rand, 1000)
const Bool1 = make_test_types(make_bools, 1000)
const Bool2 = make_test_types(make_bools, 1000)

# Unary operators
@perf -Float1[i] 100
@perf transpose(Float1[i]) 100
@perf abs(Float1[i]) 100
@perf sin(Float1[i]) 25
@perf round(Float1[i]) 25
@perf any(Bool1[i]) 10000
@perf all(Bool1[i]) 10000

# Binary operators
@perf isequal(Float1[i], Float2[i]) 10000
@perf .==(Float1[i], Float2[i]) 100
@perf +(Float1[i], Float2[i]) 100
@perf .+(Float1[i], Float2[i]) 100
@perf .*(Float1[i], Float2[i]) 100
@perf ./(Float1[i], Float2[i]) 50
@perf *(Float1[i], Float2[i]) 10 div(length(Float1), 2)+1:length(Float1)
@perf Bool1[i] & Bool2[i] 100
@perf Bool1[i] | Bool2[i] 100
@perf Bool1[i] $ Bool2[i] 100

# Vector operators
@perf diff(Float1[i]) 50 1:div(length(Float1), 2)
@perf cumsum(Float1[i]) 50 1:div(length(Float1), 2)
end


text:
@testset "SAVE" begin

  mktempdir() do path
    @testset "array" begin
      file = joinpath(path, "arr.bin")
      x0 = randn(())
      x1 = randn((5,))
      x2 = randn((3,3))
      x3 = randn((4,3,2))
      x4 = randn((5,4,3,2))
      Torch.THsave(file, Dict("x0"=>x0, "x1"=>x1,
                    "x2"=>x2, "x3"=>x3, "x4"=>x4))
      save_arr = libtorch.load(file)
      @test save_arr["x0"].numpy() ≈ x0
      @test save_arr["x1"].numpy() ≈ x1
      @test save_arr["x2"].numpy() ≈ x2
      @test save_arr["x3"].numpy() ≈ x3
      @test save_arr["x4"].numpy() ≈ x4
    end

    x = randn(10, 10)
    @testset "slice" begin
      file = joinpath(path, "slice.bin")
      x_slice = @strided x[2:5, 3:8]
      Torch.THsave(file, [x, x_slice])
      test_save_slice = libtorch.load(file)
      save_x, save_slice = test_save_slice
      @test save_x.numpy() ≈ x
      @test save_slice.numpy() ≈ x_slice
    end

    @testset "stride" begin
      file = joinpath(path, "stride.bin")
      x_stride = @strided x[2:3:9, 1:2:8]
      Torch.THsave(file, [x, x_stride])
      test_save_stride = libtorch.load(file)
      save_x, save_stride = test_save_stride
      @test save_x.numpy() ≈ x
      @test save_stride.numpy() ≈ x_stride
    end

    @testset "reshape" begin
      file = joinpath(path, "reshape.bin")
      x_reshape = @strided reshape(x, (2,5,5,2))
      Torch.THsave(file, [x, x_reshape])
      test_save_reshape = libtorch.load(file)
      save_x, save_reshape = test_save_reshape
      @test save_x.numpy() ≈ x
      @test save_reshape.numpy() ≈ x_reshape
    end

    @testset "offset" begin
      file = joinpath(path, "offset.bin")
      x_offset = @strided reshape(reshape(x, 100)[6:end], (5, 19))
      Torch.THsave(file, [x, x_offset])
      test_save_offset = libtorch.load(file)
      save_x, save_offset = test_save_offset
      @test save_x.numpy() ≈ x
      @test save_offset.numpy() ≈ x_offset
    end

    @testset "mutation" begin
      file = joinpath(path, "mutation.bin")
      x1 = copy(x)
      x2 = @strided x1[6:end, :]
      Torch.THsave(file, [x1, x2])
      test_save_mutate = libtorch.load(file)
      save_x1, save_x2 = test_save_mutate
      set!(save_x2, 0, 0)
      x2[1, :] .= 0
      @test (get(save_x1, 5) == 0).all().numpy()[]
      @test save_x1.numpy() ≈ x1
      @test save_x2.numpy() ≈ x2
    end

    @testset "transpose" begin
      file = joinpath(path, "transpose.bin")
      x_transpose = @strided permutedims(x, (2,1))
      Torch.THsave(file, [x, x_transpose])
      test_save_transpose = libtorch.load(file)
      save_x, save_transpose = test_save_transpose
      @test save_x.numpy() ≈ x
      @test save_transpose.numpy() ≈ x_transpose
    end

  end

end


text:
## Exercise 4-2
## Write a function called square that takes a parameter named t, which is a turtle. It should use the turtle to draw a square.
using ThinkJulia

println("Ans: ")

function square(turtle::Turtle, distance::Int = 100)
    @svg begin
        for i in 1:3
            forward(turtle, distance)
            turn(turtle, -90)
        end

        forward(turtle, distance)
    end
end

turtle = Turtle()
square(turtle)

println("End.")


text:
"""
The method for two vectors is like R's findinterval.
"""
searchsortedfirst2(v::AbstractVector, x::AbstractVector) = searchsortedfirst2(v, x, 1, length(x))
function searchsortedfirst2(v::AbstractVector, x::AbstractVector, lo::Int, hi::Int)
    indices = similar(x)
    n = hi = hi + 1
    @inbounds for (i,query) in enumerate(x)
        # unsorted x, restart left side
        if lo < 1 || query <= v[lo]
            lo = 0
        end
        # cast out exponentially to get hi to the right of query
        jump = 1
        while hi <= n && query > v[hi]
            lo = hi
            hi = hi + jump
            jump = jump * 2
        end
        hi > n ? n : hi
        # binary search for the exact bin
        while lo < hi-1
            m = (lo+hi)>>>1
            if query > v[m]
                lo = m
            else
                hi = m
            end
        end
        indices[i] = hi
        lo = hi - 1
    end
    return(indices)
end


text:
using DataFrames, IndexedTables, CSV

##############################################################################
##
## defVars()
##
##############################################################################

"""
    defVars(vars, df;<keyword arguments>)

Set getData() and setData() type function for each variable binding to the IndexedTable that it returns.

# Arguments
* `vars`:  an array of variables for which to build the get/set functions;
* `df`: the dataframe from which to load the inital values of the variables. The DataFrame must be in long mode with one column for each dimension;
* `tableName` (def "df"): the name of the variable pointing to the IndexedTable returned by this function (a string);
* `varNameCol` (def "varName"): the name of the dataframe column containing the variable name
* `valueCol` (def "value"): the name of the dataframe column storing the value
* `retValue` (def missing): what the `var_()`` functions should return if nothing is found (default to missing)
* `debug` (def false): returns a touple with the getData() and setData() expressions instead of actual parsing and evaluating them

# Notes:
* This function autogenerate two type of functions:
    * getData-type of function in the form of var1_(dim1=key1,dim2=key2,...)
    * setData-type of function in the form of var1!(value,dim1=key1,dim2=key2,..)
* If debug is passed, the function returns a touple with (getData(), setData()) expressions. These can then be parsed and evaluated as needed.

# Examples
```julia
julia> tableData = defVars(["supply","cons","exp","imp","transfCoeff","tranfCost"], data, tableName="tableData",varNameCol="variable",valueCol="value")
```
"""
function defVars(vars, df; tableName="table", varNameCol="varName", valueCol="value", retValue=missing, debug=false)

    colNames = names(df)
    varColPos = findall(x -> x == Symbol(varNameCol), colNames )
    valueColPos = findall(x -> x == Symbol(valueCol), colNames )
    colNamesString = [String(c) for c in colNames]

    deleteat!(colNames , varColPos)
    deleteat!(colNames , findall(x -> x == Symbol(valueCol), colNames )) # index changed bec of previous delete

    # creating empty table with types from the DataFrame
    typeVarCol     = eltype(df[Symbol(varNameCol)])
    typeDimCols    = []
    for dim in colNames
        nmissing = length(findall(x -> ismissing(x), df[dim]))
        if nmissing == 0
            push!(typeDimCols,eltype(df[dim]))
        else
            push!(typeDimCols,Union{eltype(df[dim]),Missing})
        end
    end
    #typeDimCols    = [eltype(df[dim]) for dim in colNames]
    typeValueCol   =  length(findall(x -> ismissing(x), df[Symbol(valueCol)]))>0 ? Union{eltype(df[Symbol(valueCol)]),missing} : eltype(df[Symbol(valueCol)])
    typeVarDimCols = vcat(typeVarCol,typeDimCols)
    dimValues = [Array{T,1}() for T in typeVarDimCols]
    t = IndexedTables.NDSparse(dimValues..., names=vcat(Symbol(varNameCol),colNames), Array{typeValueCol,1}())

    # filling the table with data
    for r in eachrow(df)
        rv = Any[]
        for d in vcat(Symbol(varNameCol),colNames)
            push!(rv,r[d])
        end
        t[rv...] = r[Symbol(valueCol)]
    end

    # Creating get/set functions for each variable
    dimNames = ["$(colNames[i]), " for i in 1:(length(colNames)-1)]
    push!(dimNames, "$(colNames[length(colNames)])")
    dimNamesWithmissing = ["$(colNames[i]) = missing, " for i in 1:(length(colNames)-1)]
    push!(dimNamesWithmissing, "$(colNames[length(colNames)]) = missing")
    expr1 = ""
    expr2 = ""
    for var in vars
        # Get value
        expr1 *=  "\"\"\"Return the value of $(var) under the dimensions  $(dimNames...).\"\"\" "  # documentation string
        expr1  *= "function $(var)_( $(dimNamesWithmissing...)  );"
        expr1  *= "    try;"
        expr1  *= "        return $(tableName)[\"$var\",$(dimNames...)];"
        expr1  *= "    catch  e;"
        expr1  *= "        if isa(e, KeyError);"
        expr1  *= "          return $(retValue);"
        expr1  *= "        end;"
        expr1  *= "        rethrow(e);"
        expr1  *= "    end;"
        expr1  *= "end; "
        # Set value
        expr2  *=  "\"\"\"Set the value of $(var) equal to v under the dimensions $(dimNames...) (either updating existing value(s) or creating a new record).\"\"\" "
        expr2  *= "function $(var)!(v, $(dimNamesWithmissing...)  );"
        expr2  *= "    $(tableName)[\"$var\",$(dimNames...)] = v;"
        expr2  *= "    return v;"
        expr2  *= "end; "
    end
    if debug return (expr1,expr2) end
    pexpr1 = Meta.parse(expr1)
    eval(pexpr1)
    pexpr2 = Meta.parse(expr2)
    eval(pexpr2)
    return t
end

# ##############################################################################
# ##
# ## defVarsDf()
# ##
# ##############################################################################
#
# """
#     defVarsDf(vars, df;<keyword arguments>)
#
# Set getData() and setData() type function for each variable binding to the given dataframe.
#
# Use this function if you want the data-backend of your model being a DataFrame.
# This is slower than `defVars()`` that instead returns a faster IndexedTable and bind get/setData() to that.
# To use the get/setData() functions generated by this function you will need the `DataFramesMeta` package.
#
# # Arguments
# * `vars`:  an array of variables for which to build the get/set functions;
# * `df`: the dataframe to bind the variables with. The DataFrame must be in long mode with one column for each dimension;
# * `dfName` (def "df"): the name of the variable pointing to the DataFrame (this is normally just a string version of the previous parameter);
# * `varNameCol` (def "varName"): the name of the dataframe column containing the variable name
# * `valueCol` (def "value"): the name of the dataframe column storing the value
# * `retValue` (def missing): what the `var_()`` functions should return if nothing is found (default to missing)
# * `debug` (def false): returns a touple with the getData() and setData() expressions instead of actual parsing and evaluating them
#
# # Notes:
# * This function autogenerate two type of functions:
#     * getData-type of function in the form of var1_(dim1=key1,dim2=key2,...)
#     * setData-type of function in the form of var1!(value,dim1=key1,dim2=key2,..)
# * If debug is passed, the function returns a touple with (getData(), setData()) expressions. These can then be parsed and evaluated as needed.
#
#
# # Examples
# ```julia
# julia> defVarsDf(["supply","cons","exp","imp","transfCoeff","tranfCost"], data, dfName="data",varNameCol="variable",valueCol="value")
# ```
# """
# function defVarsDf(vars, df; dfName="df", varNameCol="varName", valueCol="value", retValue=missing, debug=false)
#     colNames = names(df)
#     varColPos = find(x -> x == Symbol(varNameCol), colNames )
#     valueColPos = find(x -> x == Symbol(valueCol), colNames )
#     deleteat!(colNames , varColPos)
#     deleteat!(colNames , find(x -> x == Symbol(valueCol), colNames )) # index changed bec of previous delete
#     dimNames = ["$(colNames[i]), " for i in 1:(length(colNames)-1)]
#     push!(dimNames, "$(colNames[length(colNames)])")
#     dimNamesWithmissing = ["$(colNames[i]) = missing, " for i in 1:(length(colNames)-1)]
#     push!(dimNamesWithmissing, "$(colNames[length(colNames)]) = missing")
#     expr1 = ""
#     expr2 = ""
#     for var in vars
#         # Get value
#         expr1 *=  "\"\"\"Return the value of $(var) under the dimensions  $(dimNames...).\"\"\""  # documentation string
#         expr1  *= "function $(var)_( $(dimNamesWithmissing...)  );"
#         expr1  *= "out = @where($(dfName), :$(varNameCol) .== \"$(var)\", "
#         for (i,c) in enumerate(colNames)
#         expr1  *= "isequal.(:$(c),$(c) ), "
#         end
#         expr1  *= ");"
#         expr1  *= "  if size(out)[1] > 0;"
#         expr1  *= "    return out[end,:$(valueCol)];"
#         expr1  *= "  else;"
#         expr1  *= "    return $(retValue);"
#         expr1  *= "  end;"
#         expr1  *= "end;"
#         # Set value
#         expr2  *=  "\"\"\"Set the value of $(var) equal to v under the dimensions $(dimNames...) (either updating existing value(s) or creating a new record).\"\"\""
#         expr2  *= "function $(var)!(v, $(dimNamesWithmissing...)  );"
#         expr2  *= "dfFilter = "
#         expr2  *= "($(dfName)[:$(varNameCol)] .== \"$var\") "
#         if length(colNames) > 0
#             expr2  *= " .& "
#         end
#         for (i,c) in enumerate(colNames)
#             expr2  *= "isequal.($(dfName)[:$(c)],$(c) ) "
#             if i < length(colNames)
#                 expr2 *=   " .& "
#             else
#                 expr2 *= ";"
#             end
#         end
#         expr2   *= "if any(dfFilter);"
#         expr2   *=  " $(dfName)[dfFilter, :$(valueCol)] = v;"
#         expr2   *= "else;"
#         outNames = []
#         for (i,n) in enumerate(names(df))
#             if i == valueColPos[1]
#                 push!(outNames,"v, ")
#             elseif i == varColPos[1]
#                 push!(outNames,"\"$(var)\", ")
#             else
#                 push!(outNames,"$(n), ")
#             end
#         end
#         expr2   *= " push!($(dfName), [ $(outNames...)  ]);"
#         expr2   *= "end; return v; end;"
#     end
#     if debug return (expr1,expr2) end
#     pexpr1 = parse(expr1)
#     eval(pexpr1)
#     pexpr2 = parse(expr2)
#     eval(pexpr2)
#     return nothing
# end

##############################################################################
##
## meq()
##
##############################################################################

"""
   meq(exp)

Macro to expand functions like `f(dim1,dim2,..) = value` or `t[d1 in dim1, d2 in dim2, dfix,..] = value`

With this macro it is possible to write either:

```
@meq par1!(d1 in DIM1, d2 in DIM2, dfix3) =  par2_(d1,d2)+par3_(d1,d2)
```

and obtain
```
[par1!( par2_(d1,d2)+par3_(d1,d2)   ,d1,d2,dfix3) for d1 in DIM1, d2 in DIM2]
```

or

```
@meq par1[d1 in DIM1, d2 in DIM2, dfix3] =  par2(d1,d2)+par3(d1,d2)
```

and obtain
```
[par1[d1,d2,dfix3] =  par2(d1,d2)+par3(d1,d2) for d1 in DIM1, d2 in DIM2]
```

That is, it is possible to write equations in a concise and readable way

"""
macro meq(eq) # works without MacroTools
    lhs                   = eq.args[1]
    rhs                   = eq.args[2]
    lhs_par               = lhs.args[1]
    lhs_dims              = lhs.args[2:end]
    loop_counters         = [d.args[2] for d in lhs_dims if typeof(d) == Expr]
    loop_sets             = [d.args[2] for d in lhs_dims if typeof(d) == Expr]
    loop_wholeElements    = []
    lhs_dims_placeholders = []
    for d in lhs_dims
        if typeof(d) == Expr && d.args[1] == :in
            push!(lhs_dims_placeholders,d.args[2])
            push!(loop_wholeElements, :($(d.args[2]) = $(d.args[3])))
        else
            push!(lhs_dims_placeholders,d)
        end
    end
    if (lhs.head == :ref)      # lhs is an array
        ret = Expr(:comprehension, :($lhs_par[$(lhs_dims_placeholders...)] = $(rhs)),loop_wholeElements...)
    elseif (lhs.head == :call) # lhs is a function call
        ret = Expr(:comprehension, :($lhs_par($(rhs),$(lhs_dims_placeholders...))),loop_wholeElements...)
    else
        error("Didn't understand the Left Hand Side.")
    end
    #show(ret)
    return ret
end

# # needs MacroTools
# macro meqold(ex)
#    #return dump(ex)
#    @capture(ex, par_(dims__) = rhs_)
#    loopElements = []
#    dimsPlaceholders = []
#    for d in dims
#        @capture(d, di_ in DIMi_) || (push!(dimsPlaceholders, d); continue)
#        # push!(loopElements, x)
#        push!(loopElements, :($di = $DIMi))
#        push!(dimsPlaceholders, di)
#    end
#    ret = Expr(:comprehension, :($par($(rhs),$(dimsPlaceholders...))), loopElements...)
#    #show(ret)
#    return ret
# end


text:
import Base: eltype

# definition of the `Status` struct.
# TODO: this bakes in a lot of assumptions about ordering and padding
if !@isdefined(Status)
    let
        struct_fields = Any[]
        empty_vals = Cint[]
        offset = 0
        i = 0

        while offset < MPI_Status_Source_offset
            push!(struct_fields, :($(Symbol(:_pad, i+=1))::Cint))
            push!(empty_vals, zero(Cint))
            offset += 4
        end
        @assert offset == MPI_Status_Source_offset
        push!(struct_fields, :(source::Cint))
        push!(empty_vals, MPI_ANY_SOURCE)
        offset += 4

        while offset < MPI_Status_Tag_offset
            push!(struct_fields, :($(Symbol(:_pad, i+=1))::Cint))
            push!(empty_vals, zero(Cint))
            offset += 4
        end
        @assert offset == MPI_Status_Tag_offset
        push!(struct_fields, :(tag::Cint))
        push!(empty_vals, MPI_ANY_TAG)
        offset += 4

        while offset < MPI_Status_Error_offset
            push!(struct_fields, :($(Symbol(:_pad, i+=1))::Cint))
            push!(empty_vals, zero(Cint))
            offset += 4
        end
        @assert offset == MPI_Status_Error_offset
        push!(struct_fields, :(error::Cint))
        push!(empty_vals, MPI_SUCCESS)
        offset += 4

        while offset < MPI_Status_size
            push!(struct_fields, :($(Symbol(:_pad, i+=1))::Cint))
            push!(empty_vals, zero(Cint))
            offset += 4
        end

        @eval begin
            struct Status
                $(struct_fields...)
            end
            const STATUS_EMPTY = Status($(empty_vals...))
        end
    end

    @assert fieldoffset(Status, findfirst(isequal(:source), fieldnames(Status))) == MPI_Status_Source_offset
    @assert fieldoffset(Status, findfirst(isequal(:tag), fieldnames(Status))) == MPI_Status_Tag_offset
    @assert fieldoffset(Status, findfirst(isequal(:error), fieldnames(Status))) == MPI_Status_Error_offset
    @assert sizeof(Status) == MPI_Status_size
end

"""
    MPI.Status

The status of an MPI receive communication. It has 3 accessible fields

- `source`: source of the received message
- `tag`: tag of the received message
- `error`: error code. This is only set if a function returns multiple statuses.

Additionally, the accessor function [`MPI.Get_count`](@ref) can be used to determine the
number of entries received.
"""
Status

Get_source(status::Status) = Int(status.source)
Get_tag(status::Status) = Int(status.tag)
Get_error(status::Status) = Int(status.error)

"""
    MPI.Request

An MPI Request object, representing a non-blocking communication. This also contains a
reference to the buffer used in the communication to ensure it isn't garbage-collected
during communication.

The status of a Request can be checked by the [`Wait!`](@ref) and [`Test!`](@ref) functions
or their multiple-request variants, which will deallocate the request once it is
determined to be complete. Alternatively, it will be deallocated at finalization, meaning
that it is safe to ignore the request objects if the status of the communication can be
checked by other means.

See also [`Cancel!`](@ref).
"""
@mpi_handle Request MPI_Request buffer

const REQUEST_NULL = _Request(MPI_REQUEST_NULL, nothing)
Request() = Request(REQUEST_NULL.val, nothing)
isnull(req::Request) = req.val == REQUEST_NULL.val

function free(req::Request)
    if !isnull(req) && !MPI.Finalized()
        @mpichk ccall(:MPI_Request_free, Cint, (Ptr{MPI_Request},), req)
        req.buffer = nothing
    end
    return nothing
end

"""
    status = Probe(src::Integer, tag::Integer, comm::Comm)

Blocks until there is a message that can be received matching `src`, `tag` and
`comm`. Returns the corresponding [`Status`](@ref) object.

# External links
$(_doc_external("MPI_Probe"))
"""
function Probe(src::Integer, tag::Integer, comm::Comm)
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    @mpichk ccall(:MPI_Probe, Cint,
          (Cint, Cint, MPI_Comm, Ptr{Status}),
          src, tag, comm, stat_ref)
    stat_ref[]
end

"""
    ismessage, (status|nothing) = Iprobe(src::Integer, tag::Integer, comm::Comm)

Checks if there is a message that can be received matching `src`, `tag` and `comm`. If so,
returns a tuple `true` and a [`Status`](@ref) object, otherwise returns a tuple `false, nothing`.

# External links
$(_doc_external("MPI_Iprobe"))
"""
function Iprobe(src::Integer, tag::Integer, comm::Comm)
    flag = Ref{Cint}()
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    @mpichk ccall(:MPI_Iprobe, Cint,
          (Cint, Cint, MPI_Comm, Ptr{Cint}, Ptr{Status}),
          src, tag, comm, flag, stat_ref)
    if flag[] == 0
        return false, nothing
    end
    true, stat_ref[]
end

"""
    MPI.Get_count(status::Status, T)

The number of entries received. `T` should match the argument provided by the receive call that set the status variable.

If the number of entries received exceeds the limits of the count parameter, then it returns `MPI_UNDEFINED`.

# External links
$(_doc_external("MPI_Get_count"))
"""
function Get_count(stat::Status, datatype::Datatype)
    count = Ref{Cint}()
    @mpichk ccall(:MPI_Get_count, Cint,
                  (Ptr{Status}, MPI_Datatype, Ptr{Cint}),
                  Ref(stat), datatype, count)
    Int(count[])
end
Get_count(stat::Status, ::Type{T}) where {T} = Get_count(stat, Datatype(T))

"""
    Send(buf, dest::Integer, tag::Integer, comm::Comm)

Perform a blocking send from the buffer `buf` to MPI rank `dest` of communicator `comm`
using the message tag `tag`.

# External links
$(_doc_external("MPI_Send"))
"""
function Send(buf::Buffer, dest::Integer, tag::Integer, comm::Comm)
    # int MPI_Send(const void* buf, int count, MPI_Datatype datatype, int dest,
    #              int tag, MPI_Comm comm)
    @mpichk ccall(:MPI_Send, Cint,
          (MPIPtr, Cint, MPI_Datatype, Cint, Cint, MPI_Comm),
          buf.data, buf.count, buf.datatype, dest, tag, comm)
    return nothing
end
Send(arr::Union{Ref,AbstractArray}, dest::Integer, tag::Integer, comm::Comm) =
    Send(Buffer(arr), dest, tag, comm)

"""
    Send(obj::T, dest::Integer, tag::Integer, comm::Comm) where T

Complete a blocking send of `obj` to MPI rank `dest` of communicator `comm`
using with the message tag `tag`.
"""
function Send(obj::T, dest::Integer, tag::Integer, comm::Comm) where T
    buf = Ref{T}(obj)
    Send(buf, dest, tag, comm)
end

"""
    send(obj, dest::Integer, tag::Integer, comm::Comm)

Complete a blocking send using a serialized version of `obj` to MPI rank
`dest` of communicator `comm` using with the message tag `tag`.
"""
function send(obj, dest::Integer, tag::Integer, comm::Comm)
    buf = MPI.serialize(obj)
    Send(buf, dest, tag, comm)
end

"""
    Isend(data, dest::Integer, tag::Integer, comm::Comm)

Starts a nonblocking send of `data` to MPI rank `dest` of communicator `comm` using with
the message tag `tag`.

`data` can be a [`Buffer`](@ref), or any object for which [`Buffer_send`](@ref) is defined.

Returns the [`Request`](@ref) object for the nonblocking send.

# External links
$(_doc_external("MPI_Isend"))
"""
function Isend(buf::Buffer, dest::Integer, tag::Integer, comm::Comm)
    req = Request()
    # int MPI_Isend(const void* buf, int count, MPI_Datatype datatype, int dest,
    #               int tag, MPI_Comm comm, MPI_Request *request)
    @mpichk ccall(:MPI_Isend, Cint,
          (MPIPtr, Cint, MPI_Datatype, Cint, Cint, MPI_Comm, Ptr{MPI_Request}),
                  buf.data, buf.count, buf.datatype, dest, tag, comm, req)
    req.buffer = buf
    finalizer(free, req)
    return req
end
Isend(data, dest::Integer, tag::Integer, comm::Comm) =
    Isend(Buffer_send(data), dest, tag, comm)

"""
    isend(obj, dest::Integer, tag::Integer, comm::Comm)

Starts a nonblocking send of using a serialized version of `obj` to MPI rank
`dest` of communicator `comm` using with the message tag `tag`.

Returns the commication `Request` for the nonblocking send.
"""
function isend(obj, dest::Integer, tag::Integer, comm::Comm)
    buf = MPI.serialize(obj)
    Isend(buf, dest, tag, comm)
end

"""
    Recv!(data, src::Integer, tag::Integer, comm::Comm)

Completes a blocking receive into the buffer `data` from MPI rank `src` of communicator
`comm` using with the message tag `tag`.

`data` can be a [`Buffer`](@ref), or any object for which `Buffer(data)` is defined.

Returns the [`Status`](@ref) of the receive.

# See also
- [`Recv`](@ref)
- [`recv`](@ref)

# External links
$(_doc_external("MPI_Recv"))
"""
function Recv!(buf::Buffer, src::Integer, tag::Integer, comm::Comm)
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    # int MPI_Recv(void* buf, int count, MPI_Datatype datatype, int source,
    #              int tag, MPI_Comm comm, MPI_Status *status)
    @mpichk ccall(:MPI_Recv, Cint,
                  (MPIPtr, Cint, MPI_Datatype, Cint, Cint, MPI_Comm, Ptr{Status}),
                  buf.data, buf.count, buf.datatype, src, tag, comm, stat_ref)
    return stat_ref[]
end
Recv!(buf, src::Integer, tag::Integer, comm::Comm) =
    Recv!(Buffer(buf), src, tag, comm)

"""
    Recv(::Type{T}, src::Integer, tag::Integer, comm::Comm)

Completes a blocking receive of an object of type `T` from MPI rank `src` of communicator
`comm` using with the message tag `tag`.

Returns a tuple of the object of type `T` and the [`Status`](@ref) of the receive.

# See also
- [`Recv!`](@ref)
- [`recv`](@ref)

# External links
$(_doc_external("MPI_Recv"))
"""
function Recv(::Type{T}, src::Integer, tag::Integer, comm::Comm) where T
    buf = Ref{T}()
    stat = Recv!(buf, src, tag, comm)
    (buf[], stat)
end

"""
    recv(src::Integer, tag::Integer, comm::Comm)

Completes a blocking receive of a serialized object from MPI rank `src` of communicator
`comm` using with the message tag `tag`.

Returns the deserialized object and the [`Status`](@ref) of the receive.
"""
function recv(src::Integer, tag::Integer, comm::Comm)
    stat = Probe(src, tag, comm)
    count = Get_count(stat, UInt8)
    buf = Array{UInt8}(undef, count)
    stat = Recv!(buf, Get_source(stat), Get_tag(stat), comm)
    (MPI.deserialize(buf), stat)
end

"""
    Irecv!(data, src::Integer, tag::Integer, comm::Comm)

Starts a nonblocking receive into the buffer `data` from MPI rank `src` of communicator
`comm` using with the message tag `tag`.

`data` can be a [`Buffer`](@ref), or any object for which `Buffer(data)` is defined.

Returns the [`Request`](@ref) for the nonblocking receive.

# External links
$(_doc_external("MPI_Irecv"))
"""
function Irecv!(buf::Buffer, src::Integer, tag::Integer, comm::Comm)
    req = Request()
    # int MPI_Irecv(void* buf, int count, MPI_Datatype datatype, int source,
    #               int tag, MPI_Comm comm, MPI_Request *request)
    @mpichk ccall(:MPI_Irecv, Cint,
                  (MPIPtr, Cint, MPI_Datatype, Cint, Cint, MPI_Comm, Ptr{MPI_Request}),
                  buf.data, buf.count, buf.datatype, src, tag, comm, req)
    req.buffer = buf
    finalizer(free, req)
    return req
end
Irecv!(data, src::Integer, tag::Integer, comm::Comm) =
    Irecv!(Buffer(data), src, tag, comm)

function irecv(src::Integer, tag::Integer, comm::Comm)
    (flag, stat) = Iprobe(src, tag, comm)
    if !flag
        return (false, nothing, nothing)
    end
    count = Get_count(stat, UInt8)
    buf = Array{UInt8}(undef, count)
    stat = Recv!(buf, Get_source(stat), Get_tag(stat), comm)
    (true, MPI.deserialize(buf), stat)
end

"""
    Sendrecv!(sendbuf, dest::Integer, sendtag::Integer,
              recvbuf, source::Integer, recvtag::Integer,
              comm::Comm)

Complete a blocking send-receive operation over the MPI communicator `comm`. Send
`sendcount` elements of type `sendtype` from `sendbuf` to the MPI rank `dest` using message
tag `tag`, and receive `recvcount` elements of type `recvtype` from MPI rank `source` into
the buffer `recvbuf` using message tag `tag`. Return a [`Status`](@ref) object.

If not provided, `sendtype`/`recvtype` and `sendcount`/`recvcount` are derived from the
element type and length of `sendbuf`/`recvbuf`, respectively.

# External links
$(_doc_external("MPI_Sendrecv"))
"""
function Sendrecv!(sendbuf::Buffer, dest::Integer, sendtag::Integer,
                   recvbuf::Buffer, source::Integer, recvtag::Integer,
                   comm::Comm)
    # int MPI_Sendrecv(const void *sendbuf, int sendcount, MPI_Datatype sendtype, int dest,   int sendtag,
    #                        void *recvbuf, int recvcount, MPI_Datatype recvtype, int source, int recvtag,
    #                    MPI_Comm comm, MPI_Status *status)
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    @mpichk ccall(:MPI_Sendrecv, Cint,
                  (MPIPtr, Cint, MPI_Datatype, Cint, Cint,
                   MPIPtr, Cint, MPI_Datatype, Cint, Cint,
                   MPI_Comm, Ptr{Status}),
                  sendbuf.data, sendbuf.count, sendbuf.datatype, dest, sendtag,
                  recvbuf.data, recvbuf.count, recvbuf.datatype, source, recvtag,
                  comm, stat_ref)
    return stat_ref[]
end
Sendrecv!(sendbuf, dest::Integer, sendtag::Integer, recvbuf, source::Integer, recvtag::Integer, comm::Comm) =
    Sendrecv!(Buffer(sendbuf), dest, sendtag, Buffer(recvbuf), source, recvtag, comm)

"""
    status = Wait!(req::Request)

Block until the request `req` is complete and deallocated. Returns the [`Status`](@ref) of the request.

# External links
$(_doc_external("MPI_Wait"))
"""
function Wait!(req::Request)
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    alreadynull = isnull(req)
    # int MPI_Wait(MPI_Request *request, MPI_Status *status)
    @mpichk ccall(:MPI_Wait, Cint,
                  (Ptr{MPI_Request}, Ptr{Status}),
                  req, stat_ref)
    if !alreadynull
        req.buffer = nothing
    end
    stat
end

"""
    (flag, status) = Test!(req::Request)

Check if the request `req` is complete. If so, the request is deallocated and `flag` is returned `true` and `status` as the [`Status`](@ref) of the request. Otherwise `flag` is returned `false` and `status` is `nothing`.

# External links
$(_doc_external("MPI_Test"))
"""
function Test!(req::Request)
    flag = Ref{Cint}()
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    alreadynull = isnull(req)
    # int MPI_Test(MPI_Request *request, int *flag, MPI_Status *status)
    @mpichk ccall(:MPI_Test, Cint,
                  (Ptr{MPI_Request}, Ptr{Cint}, Ptr{Status}),
                  req, flag, stat_ref)
    if flag[] == 0
        return (false, nothing)
    end
    if !alreadynull
        req.buffer = nothing
    end
    (true, stat_ref[])
end

"""
    statuses = Waitall!(reqs::Vector{Request})

Block until all active requests in the array `reqs` are complete. Returns an array of the
[`Status`](@ref) objects corresponding to each request.

# External links
$(_doc_external("MPI_Waitall"))
"""
function Waitall!(reqs::Vector{Request})
    count = length(reqs)
    reqvals = [reqs[i].val for i in 1:count]
    stats = fill(STATUS_EMPTY, count)
    # int MPI_Waitall(int count, MPI_Request array_of_requests[],
    #                 MPI_Status array_of_statuses[])
    @mpichk ccall(:MPI_Waitall, Cint,
                  (Cint, Ptr{MPI_Request}, Ptr{Status}),
                  count, reqvals, stats)
    for i in 1:count
        req = reqs[i]
        if !isnull(req)
            req.val = reqvals[i]
            req.buffer = nothing
        end
    end
    return stats
end

"""
    (flag, statuses) = Testall!(reqs::Vector{Request})

Check if all active requests in the array `reqs` are complete. If so, the requests are
deallocated and `flag` is returned as `true` and `statuses` is an array of the
[`Status`](@ref) objects corresponding to each request is returned. Otherwise no requests
are modified a tuple of `false, nothing` is returned.

# External links
$(_doc_external("MPI_Testall"))
"""
function Testall!(reqs::Vector{Request})
    count = length(reqs)
    reqvals = [reqs[i].val for i in 1:count]
    flag = Ref{Cint}()
    stats = fill(STATUS_EMPTY, count)
    # int MPI_Testall(int count, MPI_Request array_of_requests[], int *flag,
    #                 MPI_Status array_of_statuses[])
    @mpichk ccall(:MPI_Testall, Cint,
                  (Cint, Ptr{MPI_Request}, Ptr{Cint}, Ptr{Status}),
                  count, reqvals, flag, stats)
    if flag[] == 0
        return (false, nothing)
    end
    for i in 1:count
        req = reqs[i]
        if !isnull(req)
            req.val = reqvals[i]
            req.buffer = nothing
        end
    end
    (true, stats)
end

"""
    (index, status) = Waitany!(reqs::Vector{Request})

Blocks until one of the requests in the array `reqs` is complete: if more than one is
complete, one is chosen arbitrarily. The request is deallocated and a tuple of the index
of the completed request and its [`Status`](@ref) is returned. If there are no active
requests, then `index` is returned as `0`.

# External links
$(_doc_external("MPI_Waitany"))
"""
function Waitany!(reqs::Vector{Request})
    count = length(reqs)
    reqvals = [reqs[i].val for i in 1:count]
    ind = Ref{Cint}()
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    # int MPI_Waitany(int count, MPI_Request array_of_requests[], int *index,
    #                 MPI_Status *status)
    @mpichk ccall(:MPI_Waitany, Cint,
                  (Cint, Ptr{MPI_Request}, Ptr{Cint}, Ptr{Status}),
                  count, reqvals, ind, stat_ref)
    if ind[] == MPI_UNDEFINED
        return (0, stat_ref[])
    end
    i = Int(ind[]) + 1
    req = reqs[i]
    if !isnull(req)
        req.val = reqvals[i]
        req.buffer = nothing
    end
    (i, stat_ref[])
end

"""
    (flag, index, status) = Testany!(reqs::Vector{Request})

Check if any one of the requests in the array `reqs` is complete.

If one or more requests are complete, then one is chosen arbitrarily, deallocated and `flag` is
returned as `true`, along with the index and the [`Status`](@ref) of the request.

Otherwise, if there are no complete requests, then `index` is returned as `0`, `status` as
`nothing`, and `flag` as `true` if there are no active requests and `false` otherwise.

# External links
$(_doc_external("MPI_Testany"))
"""
function Testany!(reqs::Vector{Request})
    count = length(reqs)
    reqvals = [req.val for req in reqs]
    ind = Ref{Cint}()
    flag = Ref{Cint}()
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    # int MPI_Testany(int count, MPI_Request array_of_requests[], int *index,
    #                 int *flag, MPI_Status *status)
    @mpichk ccall(:MPI_Testany, Cint,
                  (Cint, Ptr{MPI_Request}, Ptr{Cint}, Ptr{Cint}, Ptr{Status}),
                  count, reqvals, ind, flag, stat_ref)
    if flag[] == 0
        return (false, 0, nothing)
    elseif ind[] == MPI_UNDEFINED
        return (true, 0, nothing)
    end
    i = Int(ind[]) + 1
    req = reqs[i]
    if !isnull(req)
        req.val = reqvals[i]
        req.buffer = nothing
    end
    (true, i, stat_ref[])
end

"""
    (indices, statuses) = Waitsome!(reqs::Vector{Request})

Block until at least one of the active requests in the array `reqs` is complete. The
completed requests are deallocated, and a tuple of their indices in `reqs` and their
corresponding [`Status`](@ref) objects are returned. If there are no active requests, then
the function returns immediately and `indices` and `statuses` are empty.

# External links
$(_doc_external("MPI_Waitsome"))
"""
function Waitsome!(reqs::Vector{Request})
    count = length(reqs)
    reqvals = [reqs[i].val for i in 1:count]
    outcnt = Ref{Cint}()
    inds = Array{Cint}(undef, count)
    stats = fill(STATUS_EMPTY, count)
    # int MPI_Waitsome(int incount, MPI_Request array_of_requests[],
    #                  int *outcount, int array_of_indices[],
    #                  MPI_Status array_of_statuses[])
    @mpichk ccall(:MPI_Waitsome, Cint,
          (Cint, Ptr{MPI_Request}, Ptr{Cint}, Ptr{Cint}, Ptr{Status}),
          count, reqvals, outcnt, inds, stats)
    outcount = Int(outcnt[])
    # This can happen if there were no valid requests
    if outcount == MPI_UNDEFINED
        outcount = 0
    end
    indices = Array{Int}(undef, outcount)
    for j in 1:outcount
        i = Int(inds[j])+1
        req = reqs[i]
        if !isnull(req)
            req.val = reqvals[i]
            req.buffer = nothing
        end
        indices[j] = i
    end
    resize!(stats, outcount)
    (indices, stats)
end

"""
    (indices, statuses) = Testsome!(reqs::Vector{Request})

Similar to [`Waitsome!`](@ref) except that it returns immediately: if no operations have
completed then `indices` and `statuses` will be empty.

# External links
$(_doc_external("MPI_Testsome"))
"""
function Testsome!(reqs::Vector{Request})
    count = length(reqs)
    reqvals = [req.val for req in reqs]
    outcnt = Ref{Cint}()
    inds = Array{Cint}(undef, count)
    stats = fill(STATUS_EMPTY, count)
    # int MPI_Testsome(int incount, MPI_Request array_of_requests[],
    #                  int *outcount, int array_of_indices[],
    #                  MPI_Status array_of_statuses[])
    @mpichk ccall(:MPI_Testsome, Cint,
                  (Cint, Ptr{MPI_Request}, Ptr{Cint}, Ptr{Cint}, Ptr{Status}),
                  count, reqvals, outcnt, inds, stats)
    outcount = outcnt[]
    # This can happen if there were no valid requests
    if outcount == MPI_UNDEFINED
        outcount = 0
    end
    indices = Array{Int}(undef, outcount)
    for j in 1:outcount
        i = Int(inds[j])+1
        req = reqs[i]
        if !isnull(req)
            req.val = reqvals[i]
            req.buffer = nothing
        end
        indices[j] = i
    end
    resize!(stats, outcount)
    (indices, stats)
end

"""
    Cancel!(req::Request)

Marks a pending [`Irecv!`](@ref) operation for cancellation (cancelling a [`Isend`](@ref),
while supported in some implementations, is deprecated as of MPI 3.1). Note that the
request is not deallocated, and can still be queried using the test or wait functions.

# External links
$(_doc_external("MPI_Cancel"))
"""
function Cancel!(req::Request)
    # int MPI_Cancel(MPI_Request *request)
    @mpichk ccall(:MPI_Cancel, Cint, (Ptr{MPI_Request},), req)
    nothing
end