text:
#=
Implementation of the KC2 pipeline.
=#

# There's a lot of twiddling in dealing with some of the bitmasks for some
# instructions. This is a helpful function to avoid writing out a lot of this
# manually.
@inline isbitset(x::Integer, i) = (x & (1 << i)) != 0
@inline setbit(x::Integer, i) = x | (1 << i)
isbitset(x::Integer, ::Val{N}) where N = (x & (1 << N)) != 0

# Basically clock the core.
function update!(core::AsapCore)

    # First, get a record of all the stall signals that are active this cycle.
    stalls = stall_check(core)
    # Look at Stage 4 to see if a misprediction happened. Must do this before
    # next PC is calculated.
    mispredict = checkbranch(core)

    # Evaluate stage 0 - get the next-state signals to avoid changing the
    # state of "core" until the end of the cycle.
    stage0_nextstate = pipeline_stage0(core, stalls.stall_01, mispredict)

    # Update pipe stages 1, 2, and 3.
    # These are mutating functions and will change their start-of-stage pipeline
    # entries.
    #
    # At the end, we will have to shift these entries in the pipeline according
    # to the set stall signals.
    if !stalls.stall_01
        stage1_next = pipeline_stage1(
            core,
            stage0_nextstate.pc,
            stage0_nextstate.next_pc_unbranched,
        )
    end

    if !stalls.stall_234
        stage2_next = pipeline_stage2(
            core,
            stalls.nop_2,
            mispredict,
            stage0_nextstate.next_pc_unbranched,
        )

        stage3_next = pipeline_stage3(core, mispredict)

        # Stage 4 sets new ALU Flags. Since these flags are used in Stage 5 to
        # evaluate conditional execution, store them here and update the flags
        # at the end of this function.
        stage4_next, stage4_nextstate = pipeline_stage4(core, mispredict)
    end

    if !stalls.stall_5678

        # Evaluate Stages 5, 6, and 7. These again are mutating functiong with
        # no returned next state.
        stage5_next = pipeline_stage5(core, stalls.nop_5)
        stage6_next = pipeline_stage6(core)
        stage7_next = pipeline_stage7(core)

        # Do Stage 8 write-back. Note that this stage never stalls, but after it
        # executes, it turns its instruction into a NOP which avoids further
        # writebacks.
        pipeline_stage8(core)
    end

    # Stage 0 update
    core.pc             = stage0_nextstate.pc
    core.repeat_count   = stage0_nextstate.repeat_count
    core.repeat_start   = stage0_nextstate.repeat_start
    core.repeat_end     = stage0_nextstate.repeat_end
    core.return_address = stage0_nextstate.return_address

    # Shuffle the pipeline - start at the end and move forward.
    pipeline = core.pipeline
    if !stalls.stall_01
        pipeline.stage1 = stage1_next
    end

    if !stalls.stall_234
        pipeline.stage2 = stage2_next
        pipeline.stage3 = stage3_next
        pipeline.stage4 = stage4_next
        core.aluflags = stage4_nextstate
    end

    if !stalls.stall_5678
        pipeline.stage5 = stage5_next
        pipeline.stage6 = stage6_next
        pipeline.stage7 = stage7_next
    end

    # Clock the read side of the input fifos.
    for fifo in core.fifos
        readupdate!(fifo)
    end

    # Clock the write side of the output fifos.
    for fifo in values(core.outputs)
        writeupdate!(fifo)
    end

    return nothing
end

#-------------------------------------------------------------------------------
#                          Stage 0 - PC Control
#-------------------------------------------------------------------------------

# There's a circular dependency between the PC update and updates of various
# signals down the pipeline. We don't want this stage to mutate the state of
# the Core and for down-stream pipeline stages to see those changes until the
# next call to "update". So, have pipeline_stage0 return a struct with the
# next state for items in the core. This can then be updated later to mutate
# "core" without those changes being visible to other pipestages this cycle.
struct Stage0NextState
    # Next state of the program counter.
    pc :: Int
    # Hardware looper
    repeat_count :: Int
    repeat_start :: Int
    repeat_end   :: Int
    # New return address
    return_address :: Int

    # Data that's not part of the state but needs to be passed to Stage 1 to
    # create new instructions linked to the correct return address
    next_pc_unbranched :: Int
end

# NOTE: Need to think about how to best structure this for testing.
function pipeline_stage0(core::AsapCore, stall::Bool, mispredict)
    # Initialize the next state values to the current state of the core.
    pc_next = core.pc
    repeat_count_next = core.repeat_count
    repeat_start_next = core.repeat_start
    repeat_end_next   = core.repeat_end
    return_address_next = core.return_address

    # Grab stage 4 to determine if a write is occuring to the return_address
    # register
    stage4 = core.pipeline.stage4

    ## Setup some preliminary variables to use elsewhere
    repeat_loop = core.pc == core.repeat_end && core.repeat_count > 1

    # Compute the next PC if no branch happens
    next_pc_unbranched = (!stall && repeat_loop) ? core.repeat_start : core.pc + 1

    # --------------------- #
    # Update repeat counter #
    # --------------------- #

    # If a past branch was mispredicted, restore the old repeat count because the
    # new one is potentially out-dated
    if mispredict
        repeat_count_next = core.pipeline.stage3.old_repeat_count

    # Check if stage 3 is a RPT instruction. If so, grab the value of the repeat
    # counter and store it.
    elseif core.pipeline.stage3.instruction.op == RPT
        repeat_count_next = Int(core.pipeline.stage3.src1_value)

    # Check if PC is equal to the current repeat end block and there are pending
    # repeats. If so, decrement the repeat counter.
    elseif !stall && repeat_loop
        repeat_count_next = core.repeat_count - 1
    end

    # ---------------------------------- #
    # Update repeat start and end blocks #
    # ---------------------------------- #
    if core.pipeline.stage3.instruction.op == RPT
        repeat_start_next = repeat_start(core.pipeline.stage3.instruction)
        repeat_end_next = repeat_end(core.pipeline.stage3.instruction)
    end

    # --------------------- #
    # Update return address #
    # --------------------- #

    # Four cases:
    # 1. Mispredicted jump - update to jump's next instruction
    # 2. Stage 5 write to return address
    # 3. Mispredicted non-jump. Restore the old return address.
    # 4. S1 has a new predicted taken branch.
    if mispredict && core.pipeline.stage3.instruction.jump
        return_address_next = core.pipeline.stage3.old_alternate_pc

    # If explicitly writing to return_address in Stage 5
    elseif stage4.instruction.dest == RET
        # TODO: Do a bounds check to make sure this stays in bounds.
        return_address_next = Int(stage4.result)

    # Check if rolling back a non-jump instruction
    elseif mispredict && !core.pipeline.stage3.instruction.jump
        return_address_next = core.pipeline.stage3.old_return_address

    # Check if stage 1 has a predicted-taken jump. If so, save the return
    # address as the next untaken PC.
    elseif !stall && core.pipeline.stage1.instruction.op == BRL && core.pipeline.stage1.instruction.jump
        return_address_next = next_pc_unbranched
    end

    # ---------------- #
    # Update PC Conter #
    # ---------------e #

    # Rollback from mispredicted branches
    if mispredict
        # If mispredicted branch was a return, go to its return address
        if core.pipeline.stage3.instruction.isreturn
            pc_next = core.pipeline.stage3.old_return_address

        # If this branch was auto-taken, restore to its PC + 1
        elseif core.pipeline.stage3.instruction.op == BRL
            pc_next = core.pipeline.stage3.old_alternate_pc

        # Otherwise, go to its target.
        else
            pc_next = Int(branch_target(core.pipeline.stage3.instruction))
        end

    # Do nothing if stalling this stage.
    elseif stall
        pc_next = core.pc

    # Check for autobranching for the instruction in stage 1.
    elseif core.pipeline.stage1.instruction.op == BRL
        if core.pipeline.stage1.instruction.isreturn
            pc_next = core.return_address
        else
            pc_next = Int(branch_target(core.pipeline.stage1.instruction))
        end
    else
        pc_next = next_pc_unbranched
    end

    # Return state updates to caller
    return Stage0NextState(
        pc_next,
        repeat_count_next,
        repeat_start_next,
        repeat_end_next,
        return_address_next,
        next_pc_unbranched,
    )
end

#-------------------------------------------------------------------------------
#                               Stage 1
#-------------------------------------------------------------------------------
@inline function pipeline_stage1(core::AsapCore, pc, next_pc)

    # TODO: Make this correct.
    if pc <= length(core.program)
        instruction = core.program[pc]
    else
        instruction = InstNOP()
    end
    # Insert this instruction into the pipeline.
    nextstage = PipelineEntry(instruction)

    return nextstage
end

#-------------------------------------------------------------------------------
#                               Stage 2
#-------------------------------------------------------------------------------
@inline function dereference(core::AsapCore, loc::Loc)
    operand = sym(loc)
    index   = ind(loc)
    # If the provided "operand" is accessing an address generater, hardware
    # pointer, or bypass pointer, replace the operand with ":dmem" and the
    # index with the value of the pointer. Otherwise, return `operand` and
    # `index` unchanged.

    # Check for "pointer" reference
    if operand == POINTER
        operand = DMEM
        # Add 1 to "index" to convert from 0-based indexing to 1-based indexing.
        index = Int(core.pointers[index + 1])

    # Check for address generator with no increment.
    elseif operand == AG
        operand = DMEM
        index = Int(read(core.address_generators[index + 1]))

    # Address generator with increment.
    #
    # This just marks the address generator as needing an increment. The actual
    # incremeting operation will happen after all operands in the instruction
    # have had a chance to be dereferenced.
    elseif operand == AG_PI
        # Mark this generator as needing an update.
        #
        # Do this before updating index - otherwise we're gonna have a bad
        # time ...
        mark(core.address_generators[index + 1])

        operand = DMEM
        index = Int(read(core.address_generators[index + 1]))

    elseif operand == POINTER_BYPASS
        operand = DMEM
        # Read from the result of Stage 5 (since pipe stages here are the
        # beginning of stage values, this is the computed result of the last
        # Stage 4)
        index = Int(core.pipeline.stage4.result)
    end

    # Return new operand and index. If the provided operand was not a
    # dereference, we just return them unchanded.
    return Loc(operand, index)
end

# TODO: Early kill for conditional check.
@inline function pipeline_stage2(
        core        :: AsapCore,
        nop         :: Bool,
        mispredict  :: Bool,
        pc_unbranched :: Integer,
    )

    if mispredict
        # Must clean out pending NOPs if mispredicted.
        core.pending_nops = 0
        return PipelineEntry()
    end

    # Grab the input from stage 1
    stage1 = core.pipeline.stage1
    instruction = stage1.instruction

    # Check if we're supposed to insert NOPs. If so, just do that and
    # return.
    if nop
        core.pending_nops -= 1
        nextstate = PipelineEntry()

    else
        # Do pointer dereference checks.
        # If any source or destination is targeting an address generator or hardware
        # pointer, get the DMEM address of that pointer and convert the instruction
        # into a DMEM instruction.

        # Do dereference checks for each operand of the instruction.
        dest = dereference(core, instruction.dest)
        src1 = dereference(core, instruction.src1)
        src2 = dereference(core, instruction.src2)

        # Update all of the address generators that have been read from this cycle.
        for ag in core.address_generators
            increment!(ag)
        end

        # Update the instruction in the pipeline.
        operands = SrcDestCollection(dest, src1, src2)
        new_instruction = set(instruction, operands)
        # Replace the instruction in the pipeline with possibly dereferenced version
        # of the instruction.
        nextstate = set(stage1, INSTRUCTION(new_instruction))

        # Check if the instruction was a branch. If so, record the unbranched
        # PC as well as repeat count and return address.
        if nextstate.instruction.optype == BRANCH_TYPE
            return_values = RETURNS(
                core.return_address,
                pc_unbranched,
                core.repeat_count,
            )
            nextstate = set(nextstate, return_values)
        end

        # Check if this instruction requests any NOPs.
        # If so, modify the the NOP counter.
        if instruction.nops > 0
            core.pending_nops = instruction.nops
        end
    end

    return nextstate
end

#-------------------------------------------------------------------------------
#                               Stage 3
#-------------------------------------------------------------------------------

function stage3_read(core::AsapCore, loc::Loc)
    operand = sym(loc)
    index   = ind(loc)
    # Default the return value to 0. If the operand doesn't perform a read
    # this stage, this is what will be returned.
    #
    # TODO: Think about changing this with a Union{Missing,Int16} when using
    # Julia 0.7+. May be a little more precise about if something has been
    # read or not.
    return_value::Int16 = 0

    # Do a check for immediates.
    if operand == IMMEDIATE
        # For immediates, the index in the instruction is the value of the
        # immediate. Note that this also takes core of extended immediates
        # because I'm treating that very loosely.
        return_value = signed(UInt16(unsigned(index) & 0xFFFF))

    # Check accumulator
    elseif operand == ACC
        # Take the lower 16-bits of the accumulator.
        return_value = signed(UInt16(core.accumulator & 0xFFFF))

    # Check bypass registers 3 and 5
    elseif operand == BYPASS
        if index == 3
            # Get the result from stage 6
            return_value = core.pipeline.stage5.result
        elseif index == 5
            # Get the result from stage 8
            return_value = core.pipeline.stage7.result
        end

    # Check if using the return address.
    elseif operand == RET
        return_value = core.return_address

    # DMEM read - while technically this is routed in Stage 4, we get the
    # correct behavior if we just do the read here.
    elseif operand == DMEM
        # Add 1 to index to convert from index 0 to index 1.
        return_value = core.dmem[index + 1]
    end

    return return_value
end

@inline function pipeline_stage3(core::AsapCore, mispredict)
    # Don't do anything if stalled - skip if instruction is a nop.
    stage2 = core.pipeline.stage2
    instruction = stage2.instruction

    instruction.op == NOP   && return PipelineEntry()

    if mispredict
        nextstate = PipelineEntry()
    else

        # Handle the reading of operands.
        src1_value = stage3_read(core, instruction.src1)
        src2_value = stage3_read(core, instruction.src2)

        nextstate = set(stage2, SRC(src1_value, src2_value))
    end

    return nextstate
end

#-------------------------------------------------------------------------------
#                               Stage 4
#-------------------------------------------------------------------------------
@inline function pipeline_stage4(core::AsapCore, mispredict)
    stage3 = core.pipeline.stage3
    instruction = stage3.instruction
    new_aluflags = core.aluflags

    instruction.op == NOP   && return PipelineEntry(), new_aluflags

    if mispredict
        return PipelineEntry(), new_aluflags
    end

    # Do a stage 4 read. Provide a default value for already read instructions.
    src1_value = stage4_read(core, instruction.src1, stage3.src1_value)
    src2_value = stage4_read(core, instruction.src2, stage3.src2_value)

    nextstate = set(stage3, SRC(src1_value, src2_value))

    # If this is an ALU operation ... do an ALU evaluation.
    if instruction.optype == ALU_TYPE
        result, new_aluflags = stage4_alu(
            instruction.op,
            src1_value,
            src2_value,
            core.aluflags,
            core.condexec
        )
        nextstate = set(nextstate, RESULT(result))
    end

    return nextstate, new_aluflags
end

function stage4_read(core::AsapCore, loc::Loc, default::U) where U
    operand = sym(loc)
    index   = ind(loc)

    # Default the return value
    return_value::U = default

    # Check if this is a fifo read.
    if operand == IBUF
        # TODO: Maybe check that the FIFO is indeed ready to be read here.
        # The stall check ... SHOULD ... take care of this, but if there are
        # bugs going on, it might be helpful.
        fifo = core.fifos[index + 1]
        return_value = read(fifo)

    # Read from FIFO without increment.
    elseif operand == IBUF_NEXT
        fifo = core.fifos[index + 1]
        return_value = peek(fifo)

    # TODO: Packet Routers and memory

    # Read from bypass registers 1, 2, and 4.
    elseif operand == BYPASS
        if index == 1
            return_value = core.pipeline.stage4.result
        elseif index == 2
            return_value = core.pipeline.stage5.result
        elseif index == 4
            return_value = core.pipeline.stage7.result
        end
    end

    return return_value
end

# Helper function.
function asap_add(x::Int32, y::Int32)
    # Julia default to 2's complement arithmetic, so this will return the
    # correct result even if arguments are supposed to be unsigned.
    result_32 = x + y
    result = trunc_int(Int16, result_32)
    cout = isbitset(result_32, 17)
    overflow = cout ⊻ isbitset(result_32, 16)

    return result, cout, overflow
end

function asap_add(::Type{Unsigned}, x::Int16, y::Int16, c::Bool = false)
    asap_add(Int32(unsigned(x)), Int32(unsigned(y)) + Int32(c))
end

function asap_add(::Type{Signed}, x::Int16, y::Int16, c::Bool = false)
    asap_add(Int32(x), Int32(y) + Int32(c))
end

function asap_sub(::Type{Unsigned}, x::Int16, y::Int16, c::Bool = false)
    asap_add(Int32(unsigned(x)), Int32(-y) - Int32(c))
end

function asap_sub(::Type{Signed}, x::Int16, y::Int16, c::Bool = false)
    asap_add(Int32(x), Int32(-y) - Int32(c))
end

function stage4_alu(
        op :: AsapOpcode,
        src1,
        src2,
        flags::ALUFlags,
        cxflags::Vector{CondExec}
    )
    # Unpack the instruction from the pipeline stage to get the instruction
    # op code.
    #instruction = stage.instruction

    # # Get the instruction inputs.
    # src1 = stage.src1_value
    # src2 = stage.src2_value

    # Convert the carry flag and overflow flags to Int16s
    carry_next = flags.carry
    overflow_next = flags.overflow

    # # Begin a big long chain of IF-else cases.
    # op = instruction.op

    # Break instructions into two categories:
    #
    # Addition and subtraction logic, which will share carry detection and
    # saturation at the end, and other more complicated logic.

    # Unsigned Addition
    if op == ADDSU || op == ADDU
        result, carry_next, overflow_next = asap_add(Unsigned, src1, src2)
    elseif op == ADDS || op == ADD
        result, carry_next, overflow_next = asap_add(Signed, src1, src2)
    elseif op == ADDCSU || op == ADDCU
        result, carry_next, overflow_next = asap_add(Unsigned, src1, src2, flags.carry)
    elseif op == ADDCS || op == ADDC
        result, carry_next, overflow_next = asap_add(Signed, src1, src2, flags.carry)

    # Unsigned subtraction
    elseif op == SUBSU || op == SUBU
        result, carry_next, overflow_next = asap_sub(Unsigned, src1, src2)
    elseif op == SUBS || op == SUB
        result, carry_next, overflow_next = asap_sub(Signed, src1, src2)
    elseif op == SUBCSU || op == SUBCU
        result, carry_next, overflow_next = asap_sub(Unsigned, src1, src2, flags.carry)
    elseif op == SUBCS || op == SUBC
        result, carry_next, overflow_next = asap_sub(Signed, src1, src2, flags.carry)
    # Logic operations.
    elseif op == OR
        result = src1 | src2
    elseif op == AND
        result = src1 & src2
    elseif op == XOR
        result = src1 ⊻ src2

    # Move Ops
    elseif op == MOVE || op == MOVI
        result = src1
    elseif op == MOVC
        result = flags.carry ? src1 : src2
    elseif op == MOVZ
        result = flags.zero ? src1 : src2
    elseif op == MOVCX0
        result = cxflags[1].flag ? src1 : src2
    elseif op == MOVCX1
        result = cxflags[2].flag ? src2 : src2

    # # Reduction operators.
    # elseif op == ANDWORD
    #     result = (src1 == Int16(0xFFFF)) ? Int16(1) : Int16(0)
    # elseif op == ORWORD
    #     result = (src1 == 0) ? Int16(0) : Int16(1)
    # elseif op == XORWORD
    #     result = isodd(count_ones(src1)) ? Int16(1) : Int16(0)

    # # Special operators
    # elseif op == BTRV
    #     error("BTRV not implemented yet")
    # elseif op == LSD
    #     result = ((src1 < 0) ? count_ones(src1) - 1 : count_zeros(src1) - 1) |> Int16
    # elseif op == LSDU
    #     result = count_zeros(src1) |> Int16

    # # Shift operations
    # elseif op == SHL
    #     # Since we store numbers as signed, must check if src2 is negative.
    #     temp = (src2 >= 16 || src2 < 0) ? 0 : Int(src1) << src2
    #     carry_next = (temp & 0x10000) != 0
    #     result = trunc_int(Int16, result)
    # elseif op == SHR
    #     # Shift left 1 so we can capture the carry bit.
    #     temp1 = Int(unsigned(src1)) << 1
    #     temp2 = (src2 > 16 || src2 < 0) ? 0 : temp1 >> src2
    #     # Carry the low bit
    #     carry_next = (temp2 & 1) != 0
    #     result = trunc_int(Int16, temp2 >> 1)
    # elseif op == SRA
    #     temp1 = Int(src1) << 1
    #     temp2 = (src2 > 16 || src2 < 0) ? 0 : temp1 >> src2
    #     carry_next = (result & 1) != 0
    #     result = trunc_int(Int16, temp2 >> 1)

    # # Shift left and carry
    # elseif op == SHLC
    #     # Insert carry on the right - preshift by 1
    #     temp1 = Int(src1) << 1 | Int(flags.carry)
    #     temp2 = (src2 >= 16 || src2 < 0) ? 0 : temp1 << src2
    #     # Because of the pre-shift, the new carry is at bit 17
    #     carry_next = (temp2 & 0x20000) != 0
    #     result = trunc_int(Int16, temp2 >> 1)

    # elseif op == SHRC # Shift right and carry
    #     # Insert carry-in on the left.
    #     # Make room for carry out on the right by left-shifting by 1
    #     temp1 = (Int(unsigned(src1)) | (Int(flags.carry) << 16)) << 1
    #     temkp2 = (src2 > 16 || src2 < 0) ? 0 : temp1 >> src2
    #     carry_next = (temp2 & 1) != 0
    #     result = trunc_int(Int16, temp2 >> 1)

    # elseif op == SRAC
    #     temp1 = (Int(src1) | Int(flags.carry) << 16) << 1
    #     temp2 = (src2 > 16 || src2 < 0) ? 0 : temp1 >> src2
    #     carry_next = (temp2 & 1) != 0
    #     result = trunc_int(Int16, temp2 >> 1)

    # Unrecognized instruction
    else
        error("Unrecognized op $op.")
    end

    # Update the zero and negative flats.
    zero_next = result == 0
    negative_next = (result & 0x8000) != 0

    # Return the next carry flags.
    return result, ALUFlags(
        carry_next,
        negative_next,
        overflow_next,
        zero_next,
    )
end

"""
    checkbranch(core::AsapCore, stage)

Return `true` if branch instruction in `stage` was mispredicted. Return `true`
otherwise.
"""
function checkbranch(core::AsapCore)
    # Get stage 4 of the pipeline.
    stage = core.pipeline.stage3

    # If this is not a branch, exit early.
    stage.instruction.op == BR || stage.instruction.op == BRL || return false

    # Determine if this branch is predicted taken or not.
    #
    # Assume that the check for the branch instruction has happened before entry
    # into this function.
    predict_taken = (stage.instruction.op == BRL)

    # Get the branch mask from the instruction.
    aluflags = core.aluflags
    mask = stage.src1_value
    instruction = stage.instruction

    # The bit positions of the mask are as follows:
    # 12 - Negate results
    # 11 - Unconditional (always taken)
    # 10 - Full memory output (TODO)
    #  9 - Full packet output (TODO)
    #  8 - Any obuf direction full
    #  7 - Empty memory input (TODO)
    #  6 - Empty packet input (TODO)
    #  5 - Empty IBUF1
    #  4 - EmptyeIBUF0
    #  3 - Overflow
    #  2 - Carry
    #  1 - Zero
    #  0 - Negative

    # Do a BIT short-circuiting evaulation of everything.
    #
    # This is kind of messy - but not ... TOO .. bad ...
    branch_taken_preinversion =
        # Always bit check
        isbitset(mask, 11) ||
        # ALU Flags
        (isbitset(mask, 0) && aluflags.negative)   ||
        (isbitset(mask, 1) && aluflags.zero)       ||
        (isbitset(mask, 2) && aluflags.carry)      ||
        (isbitset(mask, 3) && aluflags.overflow)   ||
        (isbitset(mask, 4) && isempty(core.fifo[1])) || # NOTE: Converted to index 1
        (isbitset(mask, 5) && isempty(core.fifo[2])) ||
        (isbitset(mask, 8) && checkoutputs(core))

    # If the invert bit is set, do that to switch the branch evaluation.
    branch_taken = isbitset(mask, 12) ?
        ~branch_taken_preinversion :
        branch_taken_preinversion

    # Do a conditional execution check for this branch.
    #
    # If this branch fails the conditional execution check, it's not taken.
    if (instruction.cxflag == CX_TRUE && !core.condexec[instruction.cx_index + 1]) ||
        (instruction.cxflag == CX_FALSE && core.condexec[instruction.cx_index + 1])
        branch_taken = false
    end

    # Check if the branch was mispredicted.
    mispredicted = branch_taken ⊻ predict_taken
    return mispredicted
end

#-------------------------------------------------------------------------------
#                               Stage 5
#-------------------------------------------------------------------------------
@inline function pipeline_stage5(core::AsapCore, nop::Bool)
    # Skip if stalled or ALU. If given NOP, replace stage 5 with a NOP.
    stage4 = core.pipeline.stage4
    instruction = stage4.instruction

    instruction.op == NOP  && return PipelineEntry()

    # If indicated that this instruction should return a NOP - make it so.
    if nop
        return PipelineEntry()
    end

    # Okay, we're going ahead with normal execution. Do the stage-5 saturation
    # if requested.
    if instruction.saturate
        nothing
        #saturate!()
    end

    # Handle writeback to return address.
    # Done in Stage 0
    # if instruction.dest == :ret
    #     core.return_address = stage5.result
    # end

    # ----------------------------- #
    # --- Conditional Execution --- #
    # ----------------------------- #

    # Set flags
    if instruction.cxflag == CX_SET
        # Get the conditional execution unit to use.
        condexec = core.condexec[instruction.cxindex + 1]

        setflag!(condexec, core, stage4)

    # Check if this instruction has a conditional execution flag and if
    # it should continue or not.
    #
    # CXT : Pass if flag is set.
    elseif instruction.cxflag == CX_TRUE
        # If the flag is not set, make this instruction a NOP
        if !core.condexec[instruction.cxindex + 1].flag
            return PipelineEntry()
        end

    # CXF : Pass if flag is NOT set
    elseif instruction.cxflag == CX_FALSE
        # If the flag IS set, make this a NOP
        if core.condexec[instruction.cxindex + 1].flag
            return PipelineEntry()
        end
    end

    return stage4

end

function saturate!(stage::PipelineEntry, aluflags::ALUFlags)
    # Check if this instruction was signed.
    issigned = stage.instruction.signed

    if issigned
        # Check the carry and negative flags to determine how to saturate it.
        #
        # If you think about this really hard, it kind of makes sense.
        #
        # TODO: Think about this more to figure out if it REALLY makes sense
        # or not.
        if (aluflags.carry == false) && (aluflags.negative == true)
            stage.result = typemax(Int16)
        elseif (aluflags.carry == true) && (aluflags.negative == false)
            stage.result = typemin(Int16)
        end
    else
        if aluflags.carry == true
            stage.result = signed(typemax(UInt16))
        end
    end
    return nothing
end

function setflag!(condexec::CondExec, core::AsapCore, stage)
    # Use an integer to store flags.
    numflags = 11
    flags = zero(UInt16)
    mask = condexec.mask

    # The condexec mask bits are as follows:
    # 10 - Full memory output (TODO)
    #  9 - Full packet output (TODO)
    #  8 - Any OBUF Mask full
    #  7 - Empty memory input (TODO)
    #  6 - Empty packet input (TODO)
    #  5 - Empty ibuf 1
    #  4 - Empty ibuf 0
    #  3 - ALU Overflow
    #  2 - ALU Carry
    #  1 - ALU zero
    #  0 - ALU Negative.
    # Just kind of brute force check these things.

    # Check OBUF masks
    if isbitset(mask, 8)
        # Iterate through the selected outputs of the obuf_mask.
        #
        # If any output is full, set bit 8 of the flags to zero and move on.
        for (index, bit) in enumerate(core.obuf_mask)
            bit || continue
            if isfull(core.outputs[index])
                flags = setbit(flags, 8)
                break
            end
        end
    end

    # Check IBUFS
    if isbitset(mask, 5)
        if isempty(core.fifo[2])
            flags = setbit(flags, 5)
        end
    end
    if isbitset(mask, 4)
        if isempty(core.fifo[1])
            flags = setbit(flags, 4)
        end
    end

    # Check ALU Flags.
    aluflags = core.aluflags

    if isbitset(mask, 3)
        if aluflags.overflow
            flags = setbit(flags, 3)
        end
    end
    if isbitset(mask, 2)
        if aluflags.carry
            flags = setbit(flags, 2)
        end
    end
    if isbitset(mask, 1)
        if aluflags.zero
            flags = setbit(flags, 1)
        end
    end
    if isbitset(mask, 0)
        if aluflags.negative
            flags = setbit(flags, 0)
        end
    end

    # Perfrom the reduction and set the flag.
    if condexec.unary_op == :OR
        condexec.flag = count_ones(flags) > 0
    elseif condexec.unary_op == :AND
        condexec.flag = count_ones(flags) == numflags
    else
        condexec.flag = isodd(count_ones(flags))
    end
    return nothing
end

#-------------------------------------------------------------------------------
#                               Stage 6
#-------------------------------------------------------------------------------

# Dummy stage to simulate pipelined multiplier. Actual multiplication happens
# in stage 7.
pipeline_stage6(core::AsapCore) = core.pipeline.stage5

#-------------------------------------------------------------------------------
#                               Stage 7
#-------------------------------------------------------------------------------
multiply(::Type{Signed}, a, b) = Int(signed(a)) * Int(signed(b))
multiply(::Type{Unsigned}, a, b) = UInt(unsigned(a)) * UInt(unsigned(b))

low(x) = signed(UInt16(x & 0xFFFF))
high(x) = signed(UInt16((x >> 16) & 0xFFFF))

@inline function pipeline_stage7(core::AsapCore)
    # Get the start of pipe stage and instruction.
    stage6 = core.pipeline.stage6
    instruction = stage6.instruction

    # Don't do anything if stalled.
    instruction.op == NOP  && return PipelineEntry()

    # Multiply if the instruction is of the correct type.
    if instruction.optype == MAC_TYPE
        result, accumulator = stage4567_multiply(stage6, core.accumulator)

        # Set core acculator and store result.
        core.accumulator = accumulator
        nextstate = set(stage6, RESULT(result))
    else
        nextstate = PipelineEntry()
    end
    return nextstate
end

function stage4567_multiply(stage, accumulator)
    # Another long list of conditionals depending on OP type.

    instruction = stage.instruction
    op = instruction.op
    src1_value = stage.src1_value
    src2_value = stage.src2_value

    result::Int16 = 0
    # Signed Multiply - return low
    if op == MULTL
        result = low(multiply(Signed, src1_value, src2_value))
    # Unsigned multiply - return low
    elseif op == MULTLU
        result = low(multiply(Unsigned, src1_value, src2_value))
    # Signed multiply - return high
    elseif op == MULTH
        result = high(multiply(Signed, src1_value, src2_value))
    # Unsigned multiply - return high
    elseif op == MULTHU
        result = high(multiply(Unsigned, src1_value, src2_value))
    # Signed multiply accumulate. Return low bits of accumulator.
    elseif op == MACL
        # Since the accumulator is signed, just do normal arithmetic with it.
        accumulator += multiply(Signed, src1_value, src2_value)
        result = low(accumulator)
    # Unsigned multiply accumulate. Return low bits.
    elseif op == MACLU
        unsigned_result = multiply(Unsigned, src1_value, src2_value) +
            unsigned(accumulator)

        accumulator = signed(unsigned_result)
        result = low(accumulator)
    # Signed multiply accumulate. Return high bits of accumulator.
    elseif op == MACH
        accumulator += multiply(Signed, src1_value, src2_value)
        result = high(accumulator)
    # Unsigned multiply accumulate. Return high bits.
    elseif op == MACHU
        unsigned_result = multiply(Unsigned, src1_value, src2_value) +
            unsigned(accumulator)

        accumulator = signed(unsigned_result)
        result = high(accumulator)

    # Signed multiply accumulate, wipe accumulator. Return low.
    elseif op == MACCL
        accumulator = multiply(Signed, src1_value, src2_value)
        result = low(accumulator)

    # Unsigned multiply accumulate, wipe accumulator. Return low.
    elseif op == MACCLU
        core.accumulaor = signed(multiply(Unsigned, src1_value, src2_value))
        result = low(accumulator)

    # Signed multiply accumulate, wipe accumulator. Return high.
    elseif op == MACCH
        accumulator = multiply(Signed, src1_value, src2_value)
        result = high(accumulator)

    # Unsigned multiply accumulate, wipe accumulator. Return high.
    elseif op == MACCHU
        core.accumulaor = signed(multiply(Unsigned, src1_value, src2_value))
        result = high(accumulator)

    # --- Accumulator shifts
    # Pick one of 4 shift amounts depending on src1_value
    # 0: >> 1
    # 1: >> 8
    # 2: >> 16
    # 3: << 16
    # Always return low bits of accumulator.

    # Signed shift - preserve sign bits.
    elseif op == ACCSH
        if src1_value == 0
            accumulator >>= 1
        elseif src1_value == 1
            accumulator >>= 8
        elseif src1_value == 2
            accumulator >>= 16
        elseif src1_value == 3
            accumulator <<= 16
        end

    # Unsigned shift. Use the ">>>" operator.
    elseif op == ACCSHU
        if src1_value == 0
            accumulator >>>= 1
        elseif src1_value == 1
            accumulator >>>= 8
        elseif src1_value == 2
            accumulator >>>= 16
        elseif src1_value == 3
            accumulator <<= 16
        end
    else
        error("Unrecognized MAC op: $(op)")
    end

    # Truncate the accumulator to only use 4 bits. Do this by shifting left
    # then arithmetic shifting right to preserve signed bits.
    #
    # NOTE: The behavior of this might not quite be right is an unsigned
    # instruction follows a signed instruction ...
    accumulator = ((accumulator << 64 - 40) >> (64 - 40))

    return result, accumulator
end

#-------------------------------------------------------------------------------
#                               Stage 8
#-------------------------------------------------------------------------------

# Write back stage
@inline function pipeline_stage8(core::AsapCore)
    stage5 = core.pipeline.stage5
    stage7 = core.pipeline.stage7

    # Select stage 8 for write back if is is not a NOP, it's a MAC operation,
    # and it's destination is not ":null"
    if stage7.instruction.op != NOP &&
        sym(stage7.instruction.dest) != NULL &&
        stage7.instruction.optype == MAC_TYPE

        writeback!(core, stage7)

    elseif stage5.instruction.op != NOP &&
            sym(stage5.instruction.dest) != NULL &&
            stage5.instruction.optype != MAC_TYPE

        writeback!(core, stage5)
    end

    return nothing
end

# TODO: Fill this out more completely.
function writeback!(core::AsapCore, stage::PipelineEntry)
    # Unpack some common information
    instruction = stage.instruction
    dest        = sym(instruction.dest)
    index       = ind(instruction.dest)
    result      = stage.result

    # Writeback to DMEM
    if dest == DMEM
        # For now, only handle write to DMEM
        @assert stage.instruction.sw
        core.dmem[index + 1] = stage.result

    # Write to an output buffer
    elseif dest == OUTPUT
        write(core.outputs[index + 1], result)

    # Write to hardware pointers.
    elseif dest == SET_POINTER
        # Write to the pointer in question.
        core.pointers[index + 1] = result

    # Configure AG Start and Stop
    elseif dest == AG_START
        # Get the start and stop values.
        start = unsigned(result & 0xFF)
        stop  = unsigned((result >> 8) & 0xFF)
        set!(core.address_generators[index + 1], start, stop)

    # Configure AG Stride
    elseif dest == AG_STRIDE
        set_stride!(core.address_generators[index + 1], result)

    # Set CX Mask
    elseif dest == CXMASK
        set!(core.condexec[index + 1], result)

    # Set OBUF Mask
    elseif dest == OBUF_MASK
        setobufmask!(core, result)

    end

    return nothing
end


text:
@testset "RootsImage: constructor" begin
  img = RootsImage(-1.0-2.0im, 3.0+4.0im, 4, 4)

  @test size(img.data) == (4, 4)
  @test img.remin ≈ -1.0
  @test img.remax ≈  3.0
  @test img.immin ≈ -2.0
  @test img.immax ≈  4.0
end

@testset "RootsImage: add_root!" begin
  img = RootsImage(-1.0-1.0im, 1.0+1.0im, 2, 2)

  ChristensenSet.add_root!(img, -0.5-0.5im)

  @test img.data ≈ [1.0 0.0; 0.0 0.0]

  ChristensenSet.add_root!(img, 0.5+0.5im)

  @test img.data ≈ [1.0 0.0; 0.0 1.0]

  ChristensenSet.add_root!(img, -0.5+0.5im)

  @test img.data ≈ [1.0 1.0; 0.0 1.0]

  ChristensenSet.add_root!(img, -0.5+0.5im)

  @test img.data ≈ [1.0 2.0; 0.0 1.0]
end

ext: jl


text:
# DE

# Import Packages

using PyCall
using PyPlot
@pyimport numpy as np

pygui(true)

# Drop Wave Function

f(x) = -(1 + cos(12 * x)) / (0.5 * x ^ 2 + 2)

limit = [-5, 5]
x = np.linspace(limit[1], limit[2], 500)
y = f.(x)

plot(x, y)

# Rats

N = 50
iter = 100
AC = 0.5
CR = 0.8

rat = zeros(N, 2)

for ii = 1: N

    global rat

    rat[ii, 1] = rand() * 2 * limit[2] + limit[1]
    rat[ii, 2] = f.(rat[ii, 1])

end

scatter(rat[:, 1], rat[:, 2])

globalBest = zeros(1, 2)
globalBest[1, 1] = rand() * 2 + 0.5
globalBest[1, 2] = f(globalBest[1, 1])

# Iteration

for jj = 1: iter

    for kk = 1: N

        global rat, globalBest

        ratCR = zeros(3, 1)

        while true

            ratCR = Int64.(floor.(rand(3, 1) * 30) .+ 1)

            if ratCR[1] ≠ ratCR[2] ≠ ratCR[3] ≠ kk

                break

            end

        end

        ratNew = zeros(1, 2)

        if rand() > CR

            ratNew[1, 1] = rat[ratCR[1]] + AC * (rat[ratCR[2]] - rat[ratCR[3]])

            if limit[1] ≤ ratNew[1, 1] ≤ limit[2]

                ratNew[1, 2] = f(ratNew[1, 1])

            else

                ratNew[1, 2] = 0

            end

        else

            ratNew[1, :] = rat[kk, :]

        end

        if ratNew[1, 2] ≤ rat[kk, 2]

            rat[kk, :] = ratNew[1, :]

            if rat[kk, 2] ≤ globalBest[1, 2]

                globalBest[1, :] = rat[kk, :]

            end

        end

    end

end

scatter(globalBest[1, 1], globalBest[1, 2])
println("min = $globalBest")

max_stars_repo_name: mizu-bai/Try-Julia

max_issues_repo_name: mizu-bai/Data-Science-Demos

max_forks_repo_name: mizu-bai/Data-Science-Demos


text:
module PolynomialProbabilities
using Distributions

import Base.getindex, Base.append!, Base.isequal, Base.==, Base.length

export PolynomialProbability, PolynomialProbabilityPiece, append!, ==, getindex, length

type PolynomialProbabilityPiece
    t0::Float64
    t1::Float64
    polys::Array{Float64, 2}
    exp::Vector{Float64}

    # Setting all attributes
    PolynomialProbabilityPiece(t0::Float64, t1::Float64, polys::Array{Float64, 2}, exp::Vector{Float64}) = new(t0, t1, polys, exp)
    # Uniform distribution
    PolynomialProbabilityPiece(t0::Float64, t1::Float64, constant::Float64) = new(t0, t1, eye(1) * constant, [0])
end

function ==(pp1::PolynomialProbabilityPiece, pp2::PolynomialProbabilityPiece)
    return pp1.t0 == pp2.t0 && pp1.t1 == pp2.t1 && pp1.polys == pp2.polys && pp1.exp == pp2.exp
end

type PolynomialProbability <: ContinuousUnivariateDistribution
    pieces::Vector{PolynomialProbabilityPiece}

    # Empty distribution
    PolynomialProbability() = new([])
    # Full domain simple distribution
    PolynomialProbability(polys::Array{Float64, 2}, exp::Vector{Float64}) = new([PolynomialProbabilityPiece(-Inf, Inf, polys, exp)])
end

length(pp::PolynomialProbability) = length(pp.pieces)
getindex(pp::PolynomialProbability, ii::Int64) = pp.pieces[ii]

function ==(pp1::PolynomialProbability, pp2::PolynomialProbability)
    if length(pp1) != length(pp2)
        return false
    end

    for ii = 1:length(pp1)
        if pp1[ii] != pp2[ii]
            return false
        end
    end

    return true
end

function append!(pp::PolynomialProbability, t0::Float64, t1::Float64, polys::Array{Float64, 2}, exp::Vector{Float64})
    piece = PolynomialProbabilityPiece(t0, t1, polys, exp)
    push!(pp.pieces, piece)

    pp
end

function append!(pp::PolynomialProbability, t0::Float64, t1::Float64, constant::Float64)
    piece = PolynomialProbabilityPiece(t0, t1, constant)
    push!(pp.pieces, piece)

    pp
end

include("ppsimple.jl")

end # module


text:
# syntax: proto3
using Compat
using ProtoBuf
import ProtoBuf.meta
import Base: hash, isequal, ==

type __enum_NullValue <: ProtoEnum
    NULL_VALUE::Int32
    __enum_NullValue() = new(0)
end #type __enum_NullValue
const NullValue = __enum_NullValue()

type Struct
    fields::Dict{AbstractString,Any} # map entry
    Struct(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type Struct
hash(v::Struct) = ProtoBuf.protohash(v)
isequal(v1::Struct, v2::Struct) = ProtoBuf.protoisequal(v1, v2)
==(v1::Struct, v2::Struct) = ProtoBuf.protoeq(v1, v2)

type ListValue
    values::Array{Any,1}
    ListValue(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type ListValue
hash(v::ListValue) = ProtoBuf.protohash(v)
isequal(v1::ListValue, v2::ListValue) = ProtoBuf.protoisequal(v1, v2)
==(v1::ListValue, v2::ListValue) = ProtoBuf.protoeq(v1, v2)

type Value
    null_value::Int32
    number_value::Float64
    string_value::AbstractString
    bool_value::Bool
    struct_value::Struct
    list_value::ListValue
    Value(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type Value
const __oneofs_Value = Int[1,1,1,1,1,1]
const __oneof_names_Value = [@compat(Symbol("kind"))]
meta(t::Type{Value}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, __oneofs_Value, __oneof_names_Value)
hash(v::Value) = ProtoBuf.protohash(v)
isequal(v1::Value, v2::Value) = ProtoBuf.protoisequal(v1, v2)
==(v1::Value, v2::Value) = ProtoBuf.protoeq(v1, v2)

function meta(t::Type{ListValue})
    haskey(ProtoBuf._metacache, t) && (return ProtoBuf._metacache[t])
    m = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
    m.ordered[1].meta = meta(Value)
    m
end

function meta(t::Type{Struct})
    haskey(ProtoBuf._metacache, t) && (return ProtoBuf._metacache[t])
    m = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
    m.ordered[1].meta = ProtoBuf.mapentry_meta(Dict{AbstractString,Value})
    m
end

export NullValue, Struct, Value, ListValue


text:
"""
$(DocStringExtensions.TYPEDSIGNATURES)

Description

```jldoctest
julia> 1 + 1
2
```
"""
function graph_to_kmers(g)
    kmers = [g.vprops[v][:kmer] for v in Graphs.vertices(g)]
    return kmers
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

Description

```jldoctest
julia> 1 + 1
2
```
"""
function graph_to_edge_sequences(g)
    edges = Set{BioSequences.BigDNAMer{g.gprops[:k]+1}}()
    for edge in Graphs.edges(g)
        src_kmer = g.vprops[edge.src][:kmer]
        dst_kmer = g.vprops[edge.dst][:kmer]
        for orientation in g.eprops[edge][:orientations]
            if orientation.source_orientation
                oriented_src_kmer = src_kmer
            else
                oriented_src_kmer = BioSequences.reverse_complement(src_kmer)
            end
            if orientation.destination_orientation
                oriented_dst_kmer = dst_kmer
            else
                oriented_dst_kmer = BioSequences.reverse_complement(dst_kmer)
            end
            for i in 1:g.gprops[:k]-1
                @assert oriented_src_kmer[i+1] == oriented_dst_kmer[i]
            end
            edge_mer = BioSequences.BigDNAMer((nuc for nuc in oriented_src_kmer)..., last(oriented_dst_kmer))
            push!(edges, BioSequences.canonical(edge_mer))
        end
    end
    return edges
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

Description

```jldoctest
julia> 1 + 1
2
```
"""
function kmer_graph_distances(g1, g2)
    g1_kmers = Set(graph_to_kmers(g1))
    g1_edges = graph_to_edge_sequences(g1)

    g2_kmers = Set(graph_to_kmers(g2))
    g2_edges = graph_to_edge_sequences(g2)

    kmer_distance = 1 - LSHFunctions.jaccard(g1_kmers, g2_kmers)
    edge_distance = 1 - LSHFunctions.jaccard(g1_edges, g2_edges)

    result = (
        kmer_distance = kmer_distance,
        edge_distance = edge_distance
    )
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function set_metadata!(kmer_graph, vertex::V, key, value) where V <: Integer
    if MetaGraphs.has_prop(kmer_graph, vertex, key)
        push!(kmer_graph.vprops[vertex][key], value)
    else
        MetaGraphs.set_prop!(kmer_graph, vertex, key, Set([value]))
    end
    return true
end

function set_metadata!(kmer_graph, edge::E, key, value) where E <: Graphs.Edge
    if MetaGraphs.has_prop(kmer_graph, edge, key)
        current_value = MetaGraphs.get_prop(kmer_graph, edge, key)
        updated_value = push!(current_value, value)
        MetaGraphs.set_prop!(kmer_graph, edge, key, updated_value)
    else
        MetaGraphs.set_prop!(kmer_graph, edge, key, Set([value]))
    end
    return true
end

max_stars_repo_name: snystrom/Mycelia

max_issues_repo_name: snystrom/Mycelia

max_forks_repo_name: cjprybol/Eisenia


text:
"""
    check_method_implementation(interface::Type,
                                func_name,
                                args_funcs::AbstractVector{Function};
                                [ignore_types]::Vector{Type}=Type[],
                                [print_results]::Bool=false
                               )::Bool

Check that a given (interface) function is implemented by all subtypes.

### Input

- `interface`     -- parent interface type
- `func_name`     -- function name
- `args_funcs`    -- list of functions that each map a type to an argument
                     signature tuple
- `ignore_types`  -- (optional, default: `Type[]`) list of types that should be
                     ignored
- `print_results` -- (optional, default: `false`) flag for printing intermediate
                     results

### Output

`true` iff all subtypes implement the given function.

### Notes

It is sufficient that a subinterface implements the function, i.e., it is not
required that every *concrete* type implements the function.

This function can also print all intermediate results to STDOUT.

### Examples

```jldoctest
julia> check_method_implementation(LazySet, σ,
        Function[S -> (AbstractVector{Float64}, S)])
true
```
"""
function check_method_implementation(interface::Type,
                                     func_name,
                                     args_funcs::AbstractVector{Function};
                                     ignore_types::Vector{Type}=Type[],
                                     print_results::Bool=false
                                    )::Bool
    # first collect all base types that are subtypes of this interface
    # NOTE: 'isleaftype' does not work due to type parameters
    base_types = LazySets.subtypes(interface, true)

    # now check all base types
    for subtype in base_types
        if subtype ∈ ignore_types
            if print_results
                println("ignoring type $subtype")
            end
            continue
        end
        found = false
        for args_func in args_funcs
            if hasmethod(func_name, args_func(subtype))
                if print_results
                    println("found implementation of $func_name for $subtype")
                end
                found = true
                break
            end
        end
        if !found
            if print_results
                println("no implementation of $func_name for $subtype")
            end
            return false
        end
    end
    return true
end


text:
using Distributions, Dates, TimeSeries

module FinancialAnalytics

using Distributions, Dates, TimeSeries

export annualized_return, equity,
       greet, swilk, shapiro, ∑matrix

include("returns.jl")
include("methods.jl")

end


text:
# ---
# title: 686. Repeated String Match
# id: problem686
# author: Indigo
# date: 2021-06-29
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/repeated-string-match/description/>
# hidden: true
# ---
#
# Given two strings `a` and `b`, return the minimum number of times you should
# repeat string `a` so that string `b` is a substring of it. If it is impossible
# for `b` to be a substring of `a` after repeating it, return `-1`.
#
# **Notice:**  string `"abc"` repeated 0 times is `""`,  repeated 1 time is
# `"abc"` and repeated 2 times is `"abcabc"`.
#
#
#
# **Example 1:**
#
#
#
#     Input: a = "abcd", b = "cdabcdab"
#     Output: 3
#     Explanation: We return 3 because by repeating a three times "ab **cdabcdab** cd", b is a substring of it.
#
#
# **Example 2:**
#
#
#
#     Input: a = "a", b = "aa"
#     Output: 2
#
#
# **Example 3:**
#
#
#
#     Input: a = "a", b = "a"
#     Output: 1
#
#
# **Example 4:**
#
#
#
#     Input: a = "abc", b = "wxyz"
#     Output: -1
#
#
#
#
# **Constraints:**
#
#   * `1 <= a.length <= 104`
#   * `1 <= b.length <= 104`
#   * `a` and `b` consist of lower-case English letters.
#
#
## @lc code=start
using LeetCode

function repeated_string_match(a::String, b::String)
    minres = (length(b) - 1) ÷ length(a) + 1
    s = a ^ minres
    for i in 0:2
        occursin(b, s) && return minres + i
        s *= a
    end
    return -1
end
## @lc code=end


text:
tanh(x::Float64, ::Type{RoundingMode{:Nearest}}) =
   sinh(x, RoundingMode{:Nearest}) / cosh(x, RoundingMode{:Nearest})

function tanh(x::Float64, ::Type{RoundingMode{:Down}})
   t = sinh(x, RoundingMode{:Nearest})
   c = cosh(x, RoundingMode{:Nearest})
   with_rounding(Float64, RoundDown) do
       t = t/c
   end
   t
end
function tanh(x::Float64, ::Type{RoundingMode{:Up}})
   t = sinh(x, RoundingMode{:Nearest})
   c = cosh(x, RoundingMode{:Nearest})
   with_rounding(Float64, RoundUp) do
       t = t/c
   end
   t
end

for (fn, invfn) in ((:csc,:sin),(:sec,:cos),(:cot,:tan))
    @eval begin
        function ($fn)(x::Float64, ::Type{RoundingMode{:Up}})
            y= $(invfn)(x, RoundingMode{:Nearest})
            with_rounding(Float64, RoundDown) do
               y = 1.0 / y
            end
            y
        end
        function ($fn)(x::Float64, ::Type{RoundingMode{:Down}})
            y= $(invfn)(x, RoundingMode{:Nearest})
            with_rounding(Float64, RoundUp) do
               y = 1.0 / y
            end
            y
        end
        function ($fn)(x::Float64, ::Type{RoundingMode{:Nearest}})
            y = $(invfn)(x, RoundingMode{:Nearest})
            z = $(fn)(y)
            s = sign(y-z)
            if s < 0.0
               with_rounding(Float64, RoundDown) do
                   y = 1.0 / y
               end
            elseif x > 1.0
               with_rounding(Float64, RoundUp) do
                   y = 1.0 / y
               end
            else
               y = 1.0 / y
            end
            y
        end

    end
end


text:
module FractionalDiffEq

using LinearAlgebra, SpecialFunctions, InvertedIndices, MittagLeffler

include("main.jl")
include("matrix.jl")

export FractionalDiffEqAlgorithm
export PECE, MatrixDiscrete

# Export some special equtions
export bagleytorvik

export solve, FDEProblem

end


text:
using OpenQuantumBase, Test
# test suite for Ohmic bath
η = 1e-4
ωc = 8 * pi
β = 1 / 2.23
bath = OpenQuantumBase.OhmicBath(η, ωc, β)
cfun_test = OpenQuantumBase.build_correlation(bath)
γfun_test = OpenQuantumBase.build_spectrum(bath)

@test correlation(0.02, 0.01, bath) == correlation(0.01, bath)
@test cfun_test[1, 1](0.02, 0.01) == correlation(0.01, bath)
@test γ(0.0, bath) == 2 * pi * η / β
@test γfun_test(0.0) == 2 * pi * η / β
@test spectrum(0.0, bath) == 2 * pi * η / β
@test S(0.0, bath) ≈ OpenQuantumBase.lambshift_cpvagk(0.0, (x)->γ(x, bath)) atol = 1e-4 rtol = 1e-4
# the following test is kept as a consistency check
@test OpenQuantumBase.lambshift_cpvagk(0.0, (x)->γ(x, bath)) ≈ -0.0025132734115775254 rtol = 1e-4

η = 1e-4;fc = 4;T = 16
bath = Ohmic(η, fc, T)
τsb, err_τsb = τ_SB(bath)
τb, err_τb = τ_B(bath, 100, τsb)
@test τsb ≈ 284.61181493 atol = 1e-6 rtol = 1e-6
@test τb ≈ 0.07638653 atol = 1e-6 rtol = 1e-6
τc, err_τc = coarse_grain_timescale(bath, 100)
@test τc ≈ sqrt(τsb * τb / 5) atol = 1e-6 rtol = 1e-6

# test suite for CustomBath
cfun = (t) -> exp(-abs(t))
sfun = (ω) -> 2 / (1 + ω^2)
bath = CustomBath(correlation=cfun, spectrum=sfun)
@test correlation(1, bath) ≈ exp(-1)
@test correlation(2, 1, bath) == correlation(1, bath)
@test spectrum(0, bath) ≈ 2
@test γ(0, bath) ≈ 2
@test S(0, bath) == 0

# test suite for ensemble fluctuators
rtn = OpenQuantumBase.SymetricRTN(2.0, 2.0)
@test 4 * exp(-2 * 3) == correlation(3, rtn)
@test 2 * 4 * 2 / (4 + 4) == spectrum(2, rtn)

ensemble_rtn = EnsembleFluctuator([1.0, 2.0], [2.0, 1.0])
@test exp(-2 * 3) + 4 * exp(-3) == correlation(3, ensemble_rtn)
@test 2 * 2 / (9 + 4) + 2 * 4 / (9 + 1) == spectrum(3, ensemble_rtn)

# test suite for HybridOhmic bath
η = 0.01; W = 5; fc = 4; T = 12.5
bath = HybridOhmic(W, η, fc, T)
@test spectrum(0.0, bath) ≈ 1.7045312175373621
@test S(0.0, bath) ≈ OpenQuantumBase.lambshift_cpvagk(0.0, (x)->γ(x, bath)) atol = 1e-4 rtol = 1e-4
# the following test is kept as a consistency check
@test OpenQuantumBase.lambshift_cpvagk(0.0, (x)->γ(x, bath)) ≈ -0.2872777516270734

# test suite for correlated bath
coupling = ConstantCouplings([σ₊, σ₋], unit=:ħ)
γfun= (w) -> w >= 0 ? exp(-w) : exp(0.8w)
cbath = CorrelatedBath(((1, 2), (2, 1)), spectrum=[(w) -> 0 γfun; γfun (w) -> 0])
γm = OpenQuantumBase.build_spectrum(cbath)
@test γm[1, 1](0.0) == 0
@test γm[2, 2](0.0) == 0
@test γm[1, 2](0.5) == exp(-0.5)
@test γm[2, 1](-0.5) == exp(-0.5*0.8)

@test_throws ArgumentError OpenQuantumBase.build_correlation(cbath)
lambfun_1 = OpenQuantumBase.build_lambshift([0.0], false, cbath, nothing)
@test lambfun_1[1,1](0.0) == 0
@test lambfun_1[2,2](0.1) == 0
@test lambfun_1[1,2](0.5) == 0
@test lambfun_1[2,2](1.0) == 0

lambfun_2 = OpenQuantumBase.build_lambshift([], true, cbath, nothing)
lambfun_3 = OpenQuantumBase.build_lambshift(range(-5,5,length=1000), true, cbath, nothing)
lambfun_4 = OpenQuantumBase.build_lambshift([], true, cbath, Dict(:order=>1))
lambfun_5 = OpenQuantumBase.build_lambshift([0.0, 0.01], true, cbath, Dict(:order=>1))
@test isapprox(lambfun_2[1,2](0.0), 0.03551, atol=1e-4)
@test isapprox(lambfun_3[1,2](0.0), 0.03551, atol=1e-4)
@test lambfun_2[1,1](0) == 0
@test lambfun_3[2,2](0) == 0
@test lambfun_4[1,2](0) ≈ 0.08681 atol=1e-4
@test lambfun_5[2, 1](0.01) ≈ 0.12268 atol=1e-4


text:
"""
    _preduceBF!(M, N, Q, Z, L::Union{AbstractMatrix,Missing}, R::Union{AbstractMatrix,Missing};
                fast = true, atol = 0, rtol, roff = 0, coff = 0, rtrail = 0, ctrail = 0,
                withQ = true, withZ = true) -> n, m, p

Reduce the partitioned matrix pencil `M - λN`

              [   *         *         *  ] roff
     M - λN = [   0      M22-λN22     *  ] npp
              [   0         0         *  ] rtrail
                coff       npm     ctrail

to an equivalent basic form `F - λG = Q1'*(M - λN)*Z1` using orthogonal transformation matrices `Q1` and `Z1`
such that the subpencil `M22 - λN22` is transformed into the following standard form

                       | B  | A-λE |
          F22 - λG22 = |----|------| ,
                       | D  |  C   |

where `E` is an `nxn` non-singular matrix, and `A`, `B`, `C`, `D` are `nxn`-, `nxm`-, `pxn`- and `pxm`-dimensional matrices,
respectively. The order `n` of `E` is equal to the numerical rank of `N` determined using the absolute tolerance `atol` and
relative tolerance `rtol`. `M` and `N` are overwritten by `F` and `G`, respectively.

The performed orthogonal or unitary transformations are accumulated in `Q` (i.e., `Q <- Q*Q1`), if `withQ = true`, and
`Z` (i.e., `Z <- Z*Z1`), if `withZ = true`.

The  matrix `L` is overwritten by `Q1'*L` unless `L = missing` and the  matrix `R` is overwritten by `R*Z1` unless `R = missing`.

If `fast = true`, `E` is determined upper triangular using a rank revealing QR-decomposition with column pivoting of `N22`
and `n` is evaluated as the number of nonzero diagonal elements of the R factor, whose magnitudes are greater than
`tol = max(atol,abs(R[1,1])*rtol)`.
If `fast = false`,  `E` is determined diagonal using a rank revealing SVD-decomposition of `N22` and
`n` is evaluated as the number of singular values greater than `tol = max(atol,smax*rtol)`, where `smax`
is the largest singular value.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.
"""
function _preduceBF!(M::AbstractMatrix{T}, N::AbstractMatrix{T},
                     Q::Union{AbstractMatrix{T},Nothing}, Z::Union{AbstractMatrix{T},Nothing},
                     L::Union{AbstractVecOrMat{T},Missing} = missing, R::Union{AbstractMatrix{T},Missing} = missing;
                     atol::Real = zero(real(T)), rtol::Real = (min(size(M)...)*eps(real(float(one(T)))))*iszero(atol),
                     fast::Bool = true, roff::Int = 0, coff::Int = 0, rtrail::Int = 0, ctrail::Int = 0,
                     withQ::Bool = true, withZ::Bool = true) where T <: BlasFloat
   # In interest of performance, no dimensional checks are performed
   mM, nM = size(M)
   npp = mM-rtrail-roff
   npm = nM-ctrail-coff
   # assume the partitioned form
   #
   #              [M11-λN11  M12-λN12  M13-λN13] roff
   #     M - λN = [   0      M22-λN22  M23-λM23] npp
   #              [   0         0      M33-λN33] rtrail
   #                coff       npm      ctrail
   #
   # where M22 and N22 are npp x npm matrices

   # fast return for null dimensions
   (npp == 0 || npm == 0) && (return 0, npm, npp)

   # Step 0: Reduce M22 -λ N22 to the standard form
   i11 = 1:roff
   i22 = roff+1:roff+npp
   j22 = coff+1:coff+npm
   i12 = 1:roff+npp
   j23 = coff+1:nM
   jN23 = coff+npm+1:nM
   T <: Complex ? tran = 'C' : tran = 'T'
   Q === nothing && (withQ = false)
   Z === nothing && (withZ = false)
   if fast
      # compute in-place the QR-decomposition N22*P2 = Q2*[R2;0] with column pivoting
      N22 = view(N,i22,j22)
      _, τ, jpvt = LinearAlgebra.LAPACK.geqp3!(N22)
      tol = max(atol, rtol*abs(N22[1,1]))
      n = count(x -> x > tol, abs.(diag(N22)))
      m = npm-n
      p = npp-n
      i1 = 1:n
      # [M22 M23] <- Q2'*[M22 M23]
      LinearAlgebra.LAPACK.ormqr!('L',tran,N22,τ,view(M,i22,j23))
      # [M12; M22] <- [M12; M22]*P2
      M[i12,j22] = M[i12,j22[jpvt]]
      # N23 <- Q2'*N23
      LinearAlgebra.LAPACK.ormqr!('L',tran,N22,τ,view(N,i22,jN23))
      # N12 <- N12*P2
      N[i11,j22] = N[i11,j22[jpvt]]
      # Q <- Q*Q2
      withQ && LinearAlgebra.LAPACK.ormqr!('R','N',N22,τ,view(Q,:,i22))
      # Z <- Z*P2
      withZ && (Z[:,j22] = Z[:,j22[jpvt]])
      # L <- Q2'*L
      ismissing(L) || LinearAlgebra.LAPACK.ormqr!('L',tran,N22,τ,view(L,i22,:))
      # N22 = [R2;0]
      N22[:,:] = [ triu(N22[i1,:]); zeros(T,p,npm) ]
       # R <- R*P2
      ismissing(R) || (R[:,j22] = R[:,j22[jpvt]])

      # compute in-place the complete orthogonal decomposition N22*Z2*Pc2 = [0 E; 0 0] with E nonsingular and UT
      _, tau = LinearAlgebra.LAPACK.tzrzf!(view(N22,i1,:))
      jp = [coff+n+1:coff+npm; coff+1:coff+n]
      N22r = view(N22,i1,:)
      # [M12; M22] <- [M12; M22]*Z2*Pc2
      LinearAlgebra.LAPACK.ormrz!('R',tran,N22r,tau,view(M,i12,j22))
      M[i12,j22] = M[i12,jp]
      # N12 <- N12*Z2*Pc2
      LinearAlgebra.LAPACK.ormrz!('R',tran,N22r,tau,view(N,i11,j22))
      N[i11,j22] = N[i11,jp]
      if withZ
         LinearAlgebra.LAPACK.ormrz!('R',tran,N22r,tau,view(Z,:,j22))
         Z[:,j22] = Z[:,jp]
      end
      # R <- R*Z2*Pc2
      ismissing(R) || (LinearAlgebra.LAPACK.ormrz!('R',tran,N22r,tau,view(R,:,j22)); R[:,j22] = R[:,jp])
      N22[:,:] = [ zeros(T,n,m) triu(N22[i1,i1]); zeros(T,p,npm) ]
   else
      # compute the complete orthogonal decomposition of N22 using the SVD-decomposition
      U, S, Vt = LinearAlgebra.LAPACK.gesdd!('A',view(N,i22,j22))
      tol = max(atol, rtol*S[1])
      n = count(x -> x > tol, S)
      m = npm-n
      p = npp-n
      jp = [coff+n+1:coff+npm; coff+1:coff+n]
      # [M12; M22] <- [M12; M22]*V*Pc
      M[i12,j22] = M[i12,j22]*Vt'
      M[i12,j22] = M[i12,jp]
      # [M22 M23] <- U'*[M22 M23]
      M[i22,j23] = U'*M[i22,j23]
      # N12 <- N12*V*Pc
      N[i11,j22] = N[i11,j22]*Vt'
      N[i11,j22] = N[i11,jp]
      # N23 <- U'*N23
      N[i22,jN23] = U'*N[i22,jN23]
      # Q <- Q*U
      withQ && (Q[:,i22] = Q[:,i22]*U)
      # Z <- Q*V
      if withZ
         Z[:,j22] = Z[:,j22]*Vt'
         Z[:,j22] = Z[:,jp]
      end
      # L <- U'*L
      ismissing(L) || (L[i22,:] = U'*L[i22,:])
      # R <- R*V
      ismissing(R) || (R[:,j22] = R[:,j22]*Vt'; R[:,j22] = R[:,jp] )
      N[i22,j22] = [ zeros(T,n,m) Diagonal(S[1:n]) ; zeros(T,p,npm) ]
   end
   return n, m, p
end
"""
    _preduce1!(n::Int, m::Int, p::Int, M::AbstractMatrix, N::AbstractMatrix,
               Q::Union{AbstractMatrix,Nothing}, Z::Union{AbstractMatrix,Nothing}, tol,
               L::Union{AbstractMatrix{T},Missing} = missing, R::Union{AbstractMatrix{T},Missing} = missing;
               fast = true, roff = 0, coff = 0, rtrail = 0, ctrail = 0, withQ = true, withZ = true)

Reduce the structured pencil  `M - λN`

     M = [ *  * *  *  ] roff          N = [ *  * *  *  ] roff
         [ 0  B A  *  ] n                 [ 0  0 E  *  ] n
         [ 0  D C  *  ] p                 [ 0  0 0  *  ] p
         [ 0  * *  *  ] rtrail            [ 0  * *  *  ] rtrail
         coff m n ctrail                  coff m n ctrail

with E upper triangular and nonsingular to the following form `M1 - λN1 = Q1'*(M - λN)*Z1` with

    M1 = [ *   *    *    *    *  ] roff       N1 = [ *   *  *   *    *   ] roff
         [ 0   B1   A11 A12   *  ] τ+ρ             [ 0   0  E11 E12  *   ] τ+ρ
         [ 0   0    B2  A22   *  ] n-ρ             [ 0   0  0   E22  *   ] n-ρ
         [ 0   0    D2  C2    *  ] p-τ             [ 0   0  0   0    *   ] p-τ
         [ 0   *    *   *     *  ] rtrail          [ 0   *  *   *    *   ] rtrail
          coff m    ρ  n-ρ ctrail                   coff m  ρ  n-ρ ctrail

where τ = rank D, B1 is full row rank, and E22 is upper triangular and nonsingular.
The performed orthogonal or unitary transformations are accumulated in `Q` (i.e., `Q <- Q*Q1`), if `withQ = true`, and
`Z` (i.e., `Z <- Z*Z1`), if `withZ = true`. The rank decisions use the absolute tolerance `tol` for the nonzero elements of `M`.

The  matrix `L` is overwritten by `Q1'*L` unless `L = missing` and the  matrix `R` is overwritten by `R*Z1` unless `R = missing`.
"""
function _preduce1!(n::Int, m::Int, p::Int, M::AbstractMatrix{T}, N::AbstractMatrix{T},
                    Q::Union{AbstractMatrix{T},Nothing}, Z::Union{AbstractMatrix{T},Nothing}, tol::Real,
                    L::Union{AbstractVecOrMat{T},Missing} = missing, R::Union{AbstractMatrix{T},Missing} = missing;
                    fast::Bool = true, roff::Int = 0, coff::Int = 0, rtrail::Int = 0, ctrail::Int = 0,
                    withQ::Bool = true, withZ::Bool = true) where T <: BlasFloat
   #  Steps 1 and 2: QR- or SVD-decomposition based full column compression of [B; D] with the UT-form preservation of E
   #  Reduce the structured pencil
   #
   #  M = [ *  * *  *  ] roff          N = [ *  * *  *  ] roff
   #      [ 0  B A  *  ] n                 [ 0  0 E  *  ] n
   #      [ 0  D C  *  ] p                 [ 0  0 0  *  ] p
   #      [ 0  * *  *  ] rtrail            [ 0  * *  *  ] rtrail
   #      coff m n ctrail                  coff m n ctrail
   #
   #  with E upper triangular and nonsingular to the following form
   #
   #  M1 = [ *   *    *    *    *  ] roff        N1 = [ *   *  *   *    *   ] roff
   #       [ 0   B1   A11 A12   *  ] τ+ρ              [ 0   0  E11 E12  *   ] τ+ρ
   #       [ 0   0    B2  A22   *  ] n-ρ              [ 0   0  0   E22  *   ] n-ρ
   #       [ 0   0    D2  C2    *  ] p-τ              [ 0   0  0   0    *   ] p-τ
   #       [ 0   *    *   *     *  ] rtrail           [ 0   *  *   *    *   ] rtrail
   #        coff m    ρ  n-ρ ctrail                    coff m  ρ  n-ρ ctrail
   #
   #  where τ = rank D, B1 is full row rank and E22 is upper triangular and nonsingular.

   npm = n+m
   npp = n+p
   ZERO = zero(T)
   mM = roff + npp + rtrail
   nM = coff + npm + ctrail
   # Step 1:
   ia = roff+1:roff+n
   ja = coff+m+1:nM
   jb = coff+1:coff+m
   B = view(M,ia,jb)
   BD = view(M,roff+1:roff+npp,jb)
   if p > 0
      # compress D to [D1 D2;0 0] with D1 invertible
      ic = roff+n+1:roff+npp
      D = view(M,ic,jb)
      CE = view(M,ic,ja)
      EE = view(N,ic,coff+npm+1:nM)
      if fast
         QR = qr!(D, Val(true))
         τ = count(x -> x > tol, abs.(diag(QR.R)))
      else
         τ = count(x -> x > tol, svdvals(D))
         QR = qr!(D, Val(true))
      end
      B[:,:] = B[:,QR.p]
      lmul!(QR.Q',CE)
      lmul!(QR.Q',EE)
      withQ && rmul!(view(Q,:,ic),QR.Q)
      ismissing(L) || lmul!(QR.Q',view(L,ic,:))
      D[:,:] = [ QR.R[1:τ,:]; zeros(T,p-τ,m) ]
      # Step 2:
      k = 1
      for j = coff+1:coff+τ
         for ii = roff+n+k:-1:roff+1+k
            iim1 = ii-1
            G, M[iim1,j] = givens(M[iim1,j],M[ii,j],iim1,ii)
            M[ii,j] = ZERO
            lmul!(G,view(M,:,j+1:nM))
            lmul!(G,view(N,:,ja))  # more efficient computation possible by selecting only the relevant columns
            withQ && rmul!(Q,G')
            ismissing(L) || lmul!(G,L)
         end
         k += 1
      end
   else
      τ = 0
   end
   ρ = _preduce3!(n, m-τ, M, N, Q, Z, tol, L, R, fast = fast, coff = coff+τ, roff = roff+τ, rtrail = rtrail, ctrail = ctrail, withQ = withQ, withZ = withZ)
   if p > 0
      irt = 1:(τ+ρ)
      BD[irt,:] = BD[irt,invperm(QR.p)]
   end
   return τ, ρ
end
"""
    _preduce2!(n::Int, m::Int, p::Int, M::AbstractMatrix, N::AbstractMatrix,
               Q::Union{AbstractMatrix,Nothing}, Z::Union{AbstractMatrix,Nothing}, tol,
               L::Union{AbstractMatrix{T},Missing} = missing, R::Union{AbstractMatrix{T},Missing} = missing;
               fast = true, roff = 0, coff = 0, rtrail = 0, ctrail = 0, withQ = true, withZ = true)

Reduce the structured pencil

     M = [ *  * *   *  ] roff          N = [ *  * *  *   ] roff
         [ *  B A   *  ] n                 [ 0  0 E  *   ] n
         [ *  D C   *  ] p                 [ 0  0 0  *   ] p
         [ 0  0 0   *  ] rtrail            [ 0  0 0  *   ] rtrail
         coff m n ctrail                   coff m n ctrail

with `E` upper triangular and nonsingular to the following form `M1 - λN1 = Q1'*(M - λN)*Z1` with

    M1 = [ *    *   *    *    *  ] roff       N1 = [ *    *   *   *    *   ] roff
         [ *    B1  A11 A12   *  ] n-ρ             [ 0    0  E11 E12   *   ] n-ρ
         [ *    D1  C1  A22   *  ] ρ               [ 0    0   0  E22   *   ] ρ
         [ 0    0   0   C2    *  ] p               [ 0    0   0   0    *   ] p
         [ 0    0   0   0     *  ] rtrail          [ 0    0   0   0    *   ] rtrail
          coff m-τ n-ρ τ+ρ ctrail                   coff m-τ n-ρ τ+ρ ctrail

where `τ = rank D`, `C2` is full column rank and `E11` and `E22` are upper triangular and nonsingular.
The performed orthogonal or unitary transformations are accumulated in `Q` (i.e., `Q <- Q*Q1`), if `withQ = true`, and
`Z` (i.e., `Z <- Z*Z1`), if `withZ = true`. The rank decisions use the absolute tolerance `tol` for the nonzero elements of `M`.

The  matrix `L` is overwritten by `Q1'*L` unless `L = missing` and the  matrix `R` is overwritten by `R*Z1` unless `R = missing`.
"""
function _preduce2!(n::Int, m::Int, p::Int, M::AbstractMatrix{T}, N::AbstractMatrix{T},
                    Q::Union{AbstractMatrix{T},Nothing}, Z::Union{AbstractMatrix{T},Nothing}, tol::Real,
                    L::Union{AbstractVecOrMat{T},Missing} = missing, R::Union{AbstractMatrix{T},Missing} = missing;
                    fast::Bool = true, roff::Int = 0, coff::Int = 0, rtrail::Int = 0, ctrail::Int = 0,
                    withQ::Bool = true, withZ::Bool = true) where T <: BlasFloat
   #  Steps 1 and 2: QR- or SVD-decomposition based full column compression of [D C] with the UT-form preservation of E
   #  Reduce the structured pencil
   #
   #  M = [ *  * *  *  ] roff          N = [ *  * *  *  ] roff
   #      [ *  B A  *  ] n                 [ 0  0 E  *  ] n
   #      [ *  D C  *  ] p                 [ 0  0 0  *  ] p
   #      [ 0  0 0  *  ] rtrail            [ 0  0 0  *  ] rtrail
   #      coff m n ctrail                  coff m n ctrail
   #
   #  with E upper triangular and nonsingular to the following form
   #
   #  M1 = [ *   *   *    *    *  ] roff       N1 = [ *    *   *   *    *   ] roff
   #       [ *   B1  A11 A12   *  ] n-ρ             [ 0    0  E11 E12   *   ] n-ρ
   #       [ *   D1  C1  A22   *  ] ρ               [ 0    0   0  E22   *   ] ρ
   #       [ 0   0   0   C2    *  ] p               [ 0    0   0   0    *   ] p
   #       [ 0   0   0   0     *  ] rtrail          [ 0    0   0   0    *   ] rtrail
   #       coff m-τ n-ρ τ+ρ ctrail                   coff m-τ n-ρ τ+ρ ctrail
   #
   #  where τ = rank D, C2 is full column rank and E11 and E22 are upper triangular and nonsingular.

   npp = n+p
   npm = n+m
   ZERO = zero(T)
   # Step 1:
   mM = roff + npp + rtrail
   nM = coff + npm + ctrail
   ia = 1:roff+n
   ic = roff+n+1:roff+npp
   jc = coff+m+1:coff+npm
   jb = coff+1:coff+m
   jdc = 1:nM
   BE = view(M,ia,jb)
   EE = view(N,1:roff,jb)
   DC = view(M,ic,coff+1:coff+npm)
   C = view(M,ic,jc)
   D = view(M,ic,jb)

   if m > 0
      # compress D' to [D1' D2'; 0 0] = Q'*D'*P with D1 invertible
      ic = roff+n+1:roff+npp
      D = view(M,ic,jb)
      if fast
         QR = qr!(copy(D'), Val(true))
         τ = count(x -> x > tol, abs.(diag(QR.R)))
      else
         # τ = rank(D; atol = tol)
         τ = count(x -> x > tol, svdvals(D))
         QR = qr!(copy(D'), Val(true))
      end
      rmul!(BE,QR.Q)  # BE*Q
      jt = m:-1:1
      BE[:,:] = BE[:,jt]      # BE*Q*P2
      if withZ
         Z1 = view(Z,:,jb)
         rmul!(Z1,QR.Q)
         Z1[:,:] = Z1[:,jt]
      end
      ismissing(R) || ( rmul!(view(R,:,jb),QR.Q); R[:,jb] = R[:,reverse(jb)] )
      rmul!(EE,QR.Q)
      EE[:,:] = EE[:,jt]      # BE*Q*P2
      D[:,:] = [ zeros(p,m-τ)  QR.R[τ:-1:1,p:-1:1]' ]
      C[:,:] = C[QR.p,:]
      C[:,:] = reverse(C,dims=1)
      # Step 2:
      k = 1
      for i = roff+npp:-1:roff+npp+1-τ
         for jj = coff+m+1-k:coff+m+n-k
             jjp1 = jj+1
             G, r = givens(conj(M[i,jjp1]),conj(M[i,jj]),jjp1,jj)
             M[i,jjp1] = conj(r)
             M[i,jj] = ZERO
             rmul!(view(M,1:i-1,jdc),G')
             rmul!(view(N,1:i,jdc),G')
             withZ && rmul!(Z,G')
             ismissing(R) || rmul!(R,G')
         end
         k += 1
      end
   else
      τ = 0
   end
   ρ = _preduce4!(n, m-τ, p-τ, M, N, Q, Z, tol, L, R, fast = fast, coff = coff, roff = roff, rtrail = rtrail+τ, ctrail = ctrail+τ, withQ = withQ, withZ = withZ)
   if m > 0
      jrt = npm-(τ+ρ)+1:npm
      DC[:,jrt] = reverse(DC[:,jrt],dims=1)
      DC[:,jrt] = DC[invperm(QR.p),jrt]
   end
   return τ, ρ
end
"""
    _preduce3!(n::Int, m::Int, M::AbstractMatrix, N::AbstractMatrix,
               Q::Union{AbstractMatrix,Nothing}, Z::Union{AbstractMatrix,Nothing}, tol,
               L::Union{AbstractMatrix{T},Missing} = missing, R::Union{AbstractMatrix{T},Missing} = missing;
               fast = true, roff = 0, coff = 0, rtrail = 0, ctrail = 0,  withQ = true, withZ = true)

Reduce the structured pencil

        [ *  *  *   *  ] roff        [ *  *  *   *  ] roff
    M = [ 0  B  A   *  ] n       N = [ 0  0  E   *  ] n
        [ 0  *  *   *  ] rtrail      [ 0  *  *   *  ] rtrail
        coff m  n ctrail             coff m  n ctrail

with `E` upper triangular and nonsingular to the following form `M1 - λN1 = Q1'*(M - λN)*Z1` with

         [ *  *   *   *     *  ] roff        [ *  *  *   *     *  ] roff
    M1 = [ 0  B1  A11 A12   *  ] ρ      N1 = [ 0  0  E11 E12   *  ] ρ
         [ 0  0   A21 A22   *  ] n-ρ         [ 0  0  0   E22   *  ] n-ρ
         [ 0  *   *   *     *  ] rtrail      [ 0  *  *   *     *  ] rtrail
         coff m   ρ   n-ρ ctrail             coff m  ρ   n-ρ ctrail

where `B1` has full row rank `ρ` and `E11` and `E22` are upper triangular and nonsingular.
The performed orthogonal or unitary transformations are accumulated in `Q` (i.e., `Q <- Q*Q1`), if `withQ = true`, and
`Z` (i.e., `Z <- Z*Z1`), if `withZ = true`. The rank decisions use the absolute tolerance `tol` for the nonzero elements of `M`.

The  matrix `L` is overwritten by `Q1'*L` unless `L = missing` and the  matrix `R` is overwritten by `R*Z1` unless `R = missing`.
"""
function _preduce3!(n::Int, m::Int, M::AbstractMatrix{T}, N::AbstractMatrix{T},
                    Q::Union{AbstractMatrix{T},Nothing}, Z::Union{AbstractMatrix{T},Nothing}, tol::Real,
                    L::Union{AbstractVecOrMat{T},Missing} = missing, R::Union{AbstractMatrix{T},Missing} = missing;
                    fast::Bool = true, roff::Int = 0, coff::Int = 0, rtrail::Int = 0, ctrail::Int = 0,
                    withQ::Bool = true, withZ::Bool = true) where T <: BlasFloat
   #  Step 3: QR- or SVD-decomposition based full row compression of B and UT-form preservation of E
   #  Reduce the structured pencil
   #
   #       [ *  *  *   *  ] roff        [ *  *  *   *  ] roff
   #   M = [ 0  B  A   *  ] n       N = [ 0  0  E   *  ] n
   #       [ 0  *  *   *  ] rtrail      [ 0  *  *   *  ] rtrail
   #       coff m  n ctrail             coff m  n ctrail
   #
   #  with E upper triangular and nonsingular to the following form
   #
   #       [ *  *   *   *     *  ] roff        [ *  *  *   *     *  ] roff
   #  M1 = [ 0  B1  A11 A12   *  ] ρ      N1 = [ 0  0  E11 E12   *  ] ρ
   #       [ 0  0   A21 A22   *  ] n-ρ         [ 0  0  0   E22   *  ] n-ρ
   #       [ 0  *   *   *     *  ] rtrail      [ 0  *  *   *     *  ] rtrail
   #       coff m  n-ρ  ρ   ctrail             coff m  n-ρ ρ   ctrail
   #
   #  where B1 has full row rank and E11 and E22 are upper triangular and nonsingular.
   #
   npm = n+m
   mM = roff + n + rtrail
   nM = coff + npm + ctrail
   ZERO = zero(T)
   ib = roff+1:roff+n
   je = coff+m+1:coff+npm
   B = view(M,ib,coff+1:coff+m)
   E = view(N,ib,je)
   if fast
      ρ = 0
      nrm = similar(real(M),m)
      jp = Vector(1:m)
      nm = min(n,m)
      for j = 1:nm
         for l = j:m
            nrm[l] = norm(B[j:n,l])
         end
         nrmax, ind = findmax(nrm[j:m])
         ind += j-1
         if nrmax < tol
            break
         else
            ρ += 1
         end
         if ind != j
            (jp[j], jp[ind]) = (jp[ind], jp[j])
            (B[:,j],B[:,ind]) = (B[:,ind],B[:,j])
         end
         for ii = n:-1:j+1
             iim1 = ii-1
             G, B[iim1,j] = givens(B[iim1,j],B[ii,j],iim1,ii)
             B[ii,j] = ZERO
             lmul!(G,view(M,ib,coff+j+1:nM))
             lmul!(G,view(N,ib,coff+m+iim1:nM))
             withQ && rmul!(view(Q,:,ib),G')
             ismissing(L) || lmul!(G,view(L,ib,:))
             G, r = givens(conj(E[ii,ii]),conj(E[ii,iim1]),ii,iim1)
             E[ii,ii] = conj(r)
             E[ii,iim1] = ZERO
             rmul!(view(N,1:roff+iim1,je),G')
             withZ && rmul!(view(Z,:,je),G')
             rmul!(view(M,:,je),G')
             ismissing(R) || rmul!(view(R,:,je),G')
         end
      end
      B[:,:] = [ B[1:ρ,invperm(jp)]; zeros(T,n-ρ,m) ]
   else
      if n > m
        for j = 1:m
          for ii = n:-1:j+1
            iim1 = ii-1
            G, B[iim1,j] = givens(B[iim1,j],B[ii,j],iim1,ii)
            B[ii,j] = ZERO
            lmul!(G,view(M,ib,coff+j+1:nM))
            lmul!(G,view(N,ib,coff+m+iim1:nM))
            withQ && rmul!(view(Q,:,ib),G')
            ismissing(L) || lmul!(G,view(L,ib,:))
            G, r = givens(conj(E[ii,ii]),conj(E[ii,iim1]),ii,iim1)
            E[ii,ii] = conj(r)
            E[ii,iim1] = ZERO
            rmul!(view(N,1:roff+iim1,je),G')
            withZ && rmul!(view(Z,:,je),G')
            rmul!(view(M,:,je),G')
            ismissing(R) || rmul!(view(R,:,je),G')
         end
        end
      end
      mn = min(n,m)
      if mn > 0
         ics = 1:mn
         jcs = 1:m
         SVD = svd(B[ics,jcs], full = true)
         ρ = count(x -> x > tol, SVD.S)
         if ρ == mn
            return ρ
         else
            B[ics,jcs] = [ Diagonal(SVD.S[1:ρ])*SVD.Vt[1:ρ,:]; zeros(T,mn-ρ,m) ]
            if ρ == 0
               return ρ
            end
         end
         ibt = roff+1:roff+mn
         jt = coff+m+1:nM
         withQ && (Q[:,ibt] = Q[:,ibt]*SVD.U)
         N[ibt,jt] = SVD.U'*N[ibt,jt]
         M[ibt,jt] = SVD.U'*M[ibt,jt]
         ismissing(L) || (L[ibt,:] = SVD.U'*L[ibt,:])
         tau = similar(N,mn)
         jt1 = coff+m+1:coff+m+mn
         E11 = view(N,ibt,jt1)
         LinearAlgebra.LAPACK.gerqf!(E11,tau)
         T <: Complex ? tran = 'C' : tran = 'T'
         LinearAlgebra.LAPACK.ormrq!('R',tran,E11,tau,view(M,:,jt1))
         withZ && LinearAlgebra.LAPACK.ormrq!('R',tran,E11,tau,view(Z,:,jt1))
         LinearAlgebra.LAPACK.ormrq!('R',tran,E11,tau,view(N,1:roff,jt1))
         ismissing(R) || LinearAlgebra.LAPACK.ormrq!('R',tran,E11,tau,view(R,:,jt1))
         triu!(E11)
      else
         ρ = 0
      end
   end
   return ρ
end
"""
    _preduce4!(n::Int, m::Int, p::Int, M::AbstractMatrix, N::AbstractMatrix,
               Q::Union{AbstractMatrix,Nothing}, Z::Union{AbstractMatrix,Nothing}, tol,
               L::Union{AbstractMatrix{T},Missing} = missing, R::Union{AbstractMatrix{T},Missing} = missing;
               fast = true, roff = 0, coff = 0, rtrail = 0, ctrail = 0,  withQ = true, withZ = true)

Reduce the structured pencil

    M = [ *  *  *   *  ] roff          N = [ *  *  *   *  ] roff
        [ 0  B  A   *  ]  n                [ 0  0  E   *  ]  n
        [ 0  0  C   *  ]  p                [ 0  0  0   *  ]  p
        [ 0  *  *   *  ] rtrail            [ 0  *  *   *  ] rtrail
        coff m  n ctrail                   coff m  n ctrail

with `E` upper triangular and nonsingular to the following form `M1 - λN1 = Q1'*(M - λN)*Z1` with

    M1 = [ *  *   *   *    * ] roff         N1 = [ *  *  *   *    * ] roff
         [ 0  B1  A11 A12  * ]  n-ρ              [ 0  0  E11 E12  * ]  n-ρ
         [ 0  B2  A21 A22  * ]  ρ                [ 0  0  0   E22  * ]  ρ
         [ 0  0   0   C1   * ]  p                [ 0  0  0   0    * ]  p
         [ 0  *   *   *    * ] rtrail            [ 0  *  *   *    * ] rtrail
         coff m  n-ρ  ρ ctrail                   coff m n-ρ  ρ ctrail

where `C1` has full column rank and `E11` and `E22` are upper triangular and nonsingular.
The performed orthogonal or unitary transformations are accumulated in `Q` (i.e., `Q <- Q*Q1`), if `withQ = true`, and
`Z` (i.e., `Z <- Z*Z1`), if `withZ = true`. The rank decisions use the absolute tolerance `tol` for the nonzero elements of `M`.

The  matrix `L` is overwritten by `Q1'*L` unless `L = missing` and the  matrix `R` is overwritten by `R*Z1` unless `R = missing`.
"""
function _preduce4!(n::Int, m::Int, p::Int, M::AbstractMatrix{T},N::AbstractMatrix{T},
                    Q::Union{AbstractMatrix{T},Nothing}, Z::Union{AbstractMatrix{T},Nothing}, tol::Real,
                    L::Union{AbstractVecOrMat{T},Missing} = missing, R::Union{AbstractMatrix{T},Missing} = missing;
                    fast::Bool = true, roff::Int = 0, coff::Int = 0, rtrail::Int = 0, ctrail::Int = 0,
                    withQ::Bool = true, withZ::Bool = true) where T <: BlasFloat
   #  Step 4: QR- or SVD-decomposition based full column compression of C and UT-form preservation of E
   #  Reduce the structured pencil
   #
   #  M = [ *  *  *   *  ] roff          N = [ *  *  *   *  ] roff
   #      [ 0  B  A   *  ]  n                [ 0  0  E   *  ]  n
   #      [ 0  0  C   *  ]  p                [ 0  0  0   *  ]  p
   #      [ 0  *  *   *  ] rtrail            [ 0  *  *   *  ] rtrail
   #      coff m  n ctrail                   coff m  n ctrail
   #
   #  with E upper triangular and nonsingular to the following form
   #
   #  M1 = [ *  *   *   *    * ] roff         N1 = [ *  *  *   *    * ] roff
   #       [ 0  B1  A11 A12  * ]  n-ρ              [ 0  0  E11 E12  * ]  n-ρ
   #       [ 0  B2  A21 A22  * ]  ρ                [ 0  0  0   E22  * ]  ρ
   #       [ 0  0   0   C1   * ]  p                [ 0  0  0   0    * ]  p
   #       [ 0  *   *   *    * ] rtrail            [ 0  *  *   *    * ] rtrail
   #       coff m  n-ρ  ρ ctrail                   coff m n-ρ  ρ ctrail
   #
   #  where C1 has full column rank and E11 and E22 are upper triangular and nonsingular.
   npp = n+p
   npm = n+m
   mM = roff + npp + rtrail
   nM = coff + npm + ctrail
   ZERO = zero(T)
   ie = roff+1:roff+n
   ic = roff+n+1:roff+npp
   jc = coff+m+1:coff+npm
   C = view(M,ic,jc)
   E = view(N,ie,jc)
   AE = view(M,1:roff+n,jc)
   if fast
      ρ = 0
      nrm = similar(real(M),p)
      jp = Vector(1:p)
      np = min(n,p)
      for i = 1:np
         ii = p-i+1
         for l = 1:ii
            nrm[l] = norm(C[l,1:n-i+1])
         end
         nrmax, ind = findmax(nrm[1:ii])
         if nrmax < tol
            break
         else
            ρ += 1
         end
         if ind != ii
            (jp[ii], jp[ind]) = (jp[ind], jp[ii])
            (C[ii,:],C[ind,:]) = (C[ind,:],C[ii,:])
         end
         for jj = 1:n-i
             jjp1 = jj+1
             G, r = givens(conj(C[ii,jjp1]),conj(C[ii,jj]),jjp1,jj)
             C[ii,jjp1] = conj(r)
             C[ii,jj] = ZERO
             rmul!(view(C,1:ii-1,:),G')
             rmul!(AE,G')
             rmul!(view(N,1:roff+jjp1,jc),G')
             withZ && rmul!(view(Z,:,jc),G')
             ismissing(R) || rmul!(view(R,:,jc),G')
             G, E[jj,jj] = givens(E[jj,jj],E[jjp1,jj],jj,jjp1)
             E[jjp1,jj] = ZERO
             lmul!(G,view(N,ie,coff+m+jjp1:nM))
             withQ && rmul!(view(Q,:,ie),G')
             lmul!(G,view(M,ie,coff+1:nM))
             ismissing(L) || lmul!(G,view(L,ie,:))
         end
      end
      C[:,1:n] = [ zeros(T,p,n-ρ)  C[invperm(jp),n-ρ+1:n]];
   else
      if n > p
         for i = 1:p
             ii = p-i+1
             for jj = 1:n-i
                jjp1 = jj+1
                G, r = givens(conj(C[ii,jjp1]),conj(C[ii,jj]),jjp1,jj)
                C[ii,jjp1] = conj(r)
                C[ii,jj] = ZERO
                rmul!(view(C,1:ii-1,:),G')
                rmul!(AE,G')
                rmul!(view(N,1:roff+jjp1,jc),G')
                withZ && rmul!(view(Z,:,jc),G')
                ismissing(R) || rmul!(view(R,:,jc),G')
                G, E[jj,jj] = givens(E[jj,jj],E[jjp1,jj],jj,jjp1)
                E[jjp1,jj] = ZERO
                lmul!(G,view(N,ie,coff+m+jjp1:nM))
                withQ && rmul!(view(Q,:,ie),G')
                lmul!(G,view(M,ie,coff+1:nM))
                ismissing(L) || lmul!(G,view(L,ie,:))
               end
         end
      end
      pn = min(n,p)
      if pn > 0
         ics = 1:p
         jcs = n-pn+1:n
         SVD = svd(C[ics,jcs], full = true)
         ρ = count(x -> x > tol, SVD.S)
         if ρ == pn
            return ρ
         else
            Q1 = reverse(SVD.U,dims=2)
            C[ics,jcs] = [ zeros(T,p,pn-ρ) Q1[:,p-ρ+1:end]*Diagonal(reverse(SVD.S[1:ρ])) ]
            if ρ == 0
               return ρ
            end
         end
         Z1 = reverse(SVD.V,dims=2)
         jt = coff+npm-pn+1:coff+npm
         withZ && (Z[:,jt] = Z[:,jt]*Z1)
         M[1:roff+n,jt] = M[1:roff+n,jt]*Z1
         N[1:roff+n,jt] = N[1:roff+n,jt]*Z1    # more efficient computation possible
         ismissing(R) || (M[:,jt] = M[:,jt]*Z1)
         it = roff+n-pn+1:roff+n
         jt1 = coff+n+m+1:nM
         tau = similar(N,pn)
         E22 = view(N,it,jt)
         LinearAlgebra.LAPACK.geqrf!(E22,tau)
         T <: Complex ? tran = 'C' : tran = 'T'
         LinearAlgebra.LAPACK.ormqr!('L',tran,E22,tau,view(M,it,coff+1:nM))
         withQ && LinearAlgebra.LAPACK.ormqr!('R','N',E22,tau,view(Q,:,it))
         LinearAlgebra.LAPACK.ormqr!('L',tran,E22,tau,view(N,it,jt1))
         ismissing(L) || LinearAlgebra.LAPACK.ormqr!('L',tran,E22,tau,view(L,it,:))
         triu!(E22)
      else
         ρ = 0
      end
   end
   return ρ
end


text:
using ObjectDetector
using Test, PrettyTables
using FileIO, ImageCore

include("prepareImage.jl")
include("maintests.jl")
include("drawBoxes.jl")


text:
"""
Parallel Tempering
"""
module ParTemp

"""
Compute the log probability of swapping chains of different temperatures.

# Example

lpdf_normal(z) = -(log(2*pi) + z*z) / 2
loglike(s) = sum(lpdf_normal.(s))
states = (randn(100), randn(100))  # needs to be the same type
temps = (1.0, 2.0)  # each has to be greater than 1
compute_log_accept_ratio(loglike, states, temps)
compute_accept_ratio(loglike, states, temps)
"""
function compute_log_accept_ratio(loglike::Function, states, temps)
  s1, s2 = states
  t1, t2 = temps
  return min((1/t1 - 1/t2) * (loglike(s2) - loglike(s1)), 0)
end

"""
Compute the probability of swapping chains of different temperatures.

See doc for compute_log_accept_ratio
"""
function compute_accept_ratio(loglike::Function, states, temps)
  return exp(compute_log_accept_ratio(loglike, states, temps))
end

end  # ParTemp


text:
@test ! isprime(2 + 0im)
@test isprime(3 + 0im)


text:
import Base: run

run(agent, env::AbstractEnv, args...) = run(DynamicStyle(env), agent, env, args...)

function run(
    ::Sequential,
    agent::AbstractAgent,
    env::AbstractEnv,
    stop_condition,
    hook::AbstractHook = EmptyHook(),
)

    reset!(env)
    obs = observe(env)
    agent(PRE_EPISODE_STAGE, obs)
    hook(PRE_EPISODE_STAGE, agent, env, obs)
    action = agent(PRE_ACT_STAGE, obs)
    hook(PRE_ACT_STAGE, agent, env, obs, action)

    while true
        env(action)
        obs = observe(env)
        agent(POST_ACT_STAGE, obs)
        hook(POST_ACT_STAGE, agent, env, obs)

        if get_terminal(obs)
            agent(POST_EPISODE_STAGE, obs)  # let the agent see the last observation
            hook(POST_EPISODE_STAGE, agent, env, obs)

            stop_condition(agent, env, obs) && break

            reset!(env)
            obs = observe(env)
            agent(PRE_EPISODE_STAGE, obs)
            hook(PRE_EPISODE_STAGE, agent, env, obs)
            action = agent(PRE_ACT_STAGE, obs)
            hook(PRE_ACT_STAGE, agent, env, obs, action)
        else
            stop_condition(agent, env, obs) && break
            action = agent(PRE_ACT_STAGE, obs)
            hook(PRE_ACT_STAGE, agent, env, obs, action)
        end
    end
    hook
end

function run(
    ::Sequential,
    agent::AbstractAgent,
    env::MultiThreadEnv,
    stop_condition,
    hook::AbstractHook = EmptyHook(),
)

    while true
        reset!(env)
        obs = observe(env)
        action = agent(PRE_ACT_STAGE, obs)
        hook(PRE_ACT_STAGE, agent, env, obs, action)

        env(action)
        obs = observe(env)
        agent(POST_ACT_STAGE, obs)
        hook(POST_ACT_STAGE, agent, env, obs)

        if stop_condition(agent, env, obs)
            agent(PRE_ACT_STAGE, obs)  # let the agent see the last observation
            break
        end
    end
    hook
end

function run(
    ::Sequential,
    agents::Tuple{Vararg{<:AbstractAgent}},
    env::AbstractEnv,
    stop_condition,
    hooks = [EmptyHook() for _ in agents],
)
    reset!(env)
    observations = [observe(env, get_role(agent)) for agent in agents]

    valid_action = rand(get_action_space(env))  # init with a dummy value

    # async here?
    for (agent, obs, hook) in zip(agents, observations, hooks)
        agent(PRE_EPISODE_STAGE, obs)
        hook(PRE_EPISODE_STAGE, agent, env, obs)
        action = agent(PRE_ACT_STAGE, obs)
        hook(PRE_ACT_STAGE, agent, env, obs, action)
        # for Sequential environments, only one action is valid
        if get_current_player(env) == get_role(agent)
            valid_action = action
        end
    end

    while true
        env(valid_action)

        observations = [observe(env, get_role(agent)) for agent in agents]

        for (agent, obs, hook) in zip(agents, observations, hooks)
            agent(POST_ACT_STAGE, obs)
            hook(POST_ACT_STAGE, agent, env, obs)
        end

        if get_terminal(observations[1])
            for (agent, obs, hook) in zip(agents, observations, hooks)
                agent(POST_EPISODE_STAGE, obs)
                hook(POST_EPISODE_STAGE, agent, env, obs)
            end

            stop_condition(agents, env, observations) && break

            reset!(env)

            observations = [observe(env, get_role(agent)) for agent in agents]

            # async here?
            for (agent, obs, hook) in zip(agents, observations, hooks)
                agent(PRE_EPISODE_STAGE, obs)
                hook(PRE_EPISODE_STAGE, agent, env, obs)
                action = agent(PRE_ACT_STAGE, obs)
                hook(PRE_ACT_STAGE, agent, env, obs, action)
                # for Sequential environments, only one action is valid
                if get_current_player(env) == get_role(agent)
                    valid_action = action
                end
            end
        else
            stop_condition(agents, env, observations) && break
            for (agent, obs, hook) in zip(agents, observations, hooks)
                action = agent(PRE_ACT_STAGE, obs)
                hook(PRE_ACT_STAGE, agent, env, obs, action)
                # for Sequential environments, only one action is valid
                if get_current_player(env) == get_role(agent)
                    valid_action = action
                end
            end
        end
    end
    hooks
end


text:
# using Dates

@testset "Are.jl" begin
    filename = "data\\test.are"
    dc = DDR2import.Are.read(filename)
    @test dc["AZ"].nb_point == 48
    @test dc["AZ"].bottom_fl == 0.0
    @test dc["AZ"].top_fl == 0.0
    @test dc["AZ"].surface == 0.0
    @test dc["AZ"].sector_num == 0.0
    @test dc["AZ"].points[5,1] == 45.0
    @test dc["AZ"].points[4,2] == -40.0
    @test dc["AZ"].box[1,1] ≈ 45.00775771902942 atol = 0.00001
    #LAST TEST
    @test dc["EB"].nb_point == 550
    @test dc["EB"].points[550,1] == 51.5
end


text:
macro gdal(args...)
    @assert length(args) > 0
    @assert args[1].head == :(::)
    fhead = (args[1].args[1], GDAL.libgdal)
    returntype = args[1].args[2]
    argtypes = Expr(:tuple, [esc(a.args[2]) for a in args[2:end]]...)
    args = [esc(a.args[1]) for a in args[2:end]]
    return quote ccall($fhead, $returntype, $argtypes, $(args...)) end
end

macro ogrerr(code, message)
    return quote
        if $(esc(code)) != GDAL.OGRERR_NONE
            error($message)
        end
    end
end

macro cplerr(code, message)
    return quote
        if $(esc(code)) != GDAL.CE_None
            error($message)
        end
    end
end

macro cplwarn(code, message)
    return quote
        if $(esc(code)) != GDAL.CE_None
            warn($message)
        end
    end
end

macro cplprogress(progressfunc)
    return quote
        cfunction($(esc(progressfunc)),Cint,(Cdouble,Cstring,Ptr{Void}))
    end
end

"""
Load a `NULL`-terminated list of strings

That is it expects a "StringList", in the sense of the CPL functions, as a
NULL terminated array of strings.
"""
function unsafe_loadstringlist(pstringlist::Ptr{Cstring})
    stringlist = Vector{String}()
    (pstringlist == C_NULL) && return stringlist
    i = 1
    item = unsafe_load(pstringlist, i)
    while Ptr{UInt8}(item) != C_NULL
        push!(stringlist, unsafe_string(item))
        i += 1
        item = unsafe_load(pstringlist, i)
    end
    stringlist
end

"Fetch list of (non-empty) metadata domains. (Since: GDAL 1.11)"
metadatadomainlist(obj) =
    unsafe_loadstringlist(@gdal(GDALGetMetadataDomainList::Ptr{Cstring},
        obj.ptr::Ptr{GDAL.GDALMajorObjectH}
    ))

"Fetch metadata. Note that relatively few formats return any metadata."
metadata(obj; domain::AbstractString = "") =
    unsafe_loadstringlist(@gdal(GDALGetMetadata::Ptr{Cstring},
        obj.ptr::Ptr{GDAL.GDALMajorObjectH},
        domain::Ptr{UInt8}
    ))

"""
Set a configuration option for GDAL/OGR use.

Those options are defined as a (key, value) couple. The value corresponding to a
key can be got later with the `getconfigoption()` method.

### Parameters
* `option`  the key of the option
* `value`   the value of the option, or NULL to clear a setting.

This mechanism is similar to environment variables, but options set with
`setconfigoption()` overrides, for `getconfigoption()` point of view, values
defined in the environment.

If `setconfigoption()` is called several times with the same key, the value
provided during the last call will be used.
"""
setconfigoption(option::AbstractString, value) =
    @gdal(CPLSetConfigOption::Void, option::Cstring, value::Ptr{UInt8})

"""
This function can be used to clear a setting.

Note: it will not unset an existing environment variable; it will
just unset a value previously set by `setconfigoption()`.
"""
clearconfigoption(option::AbstractString) = setconfigoption(option, C_NULL)

"""
Get the value of a configuration option.

The value is the value of a (key, value) option set with `setconfigoption()`.
If the given option was not defined with `setconfigoption()`, it tries to find
it in environment variables.

### Parameters
* `option`  the key of the option to retrieve
* `default` a default value if the key does not match existing defined options

### Returns
the value associated to the key, or the default value if not found.
"""
function getconfigoption(option::AbstractString, default = C_NULL)
    result = @gdal(CPLGetConfigOption::Ptr{UInt8},
        option::Cstring,
        default::Ptr{UInt8}
    )
    return (result == C_NULL) ? "" : unsafe_string(result)
end

"""
Set a configuration option for GDAL/OGR use.

Those options are defined as a (key, value) couple. The value corresponding to a
key can be got later with the `getconfigoption()` method.

### Parameters
* `option`  the key of the option
* `value`   the value of the option

This function sets the configuration option that only applies in the current
thread, as opposed to `setconfigoption()` which sets an option that applies on
all threads.
"""
setthreadconfigoption(option::AbstractString, value) =
    @gdal(CPLSetThreadLocalConfigOption::Void,
        option::Cstring,
        value::Ptr{UInt8}
    )

"""
This function can be used to clear a setting.

Note: it will not unset an existing environment variable; it will
just unset a value previously set by `setthreadconfigoption()`.
"""
clearthreadconfigoption(option::AbstractString) =
    setthreadconfigoption(option, C_NULL)

"Same as `getconfigoption()` but with settings from `setthreadconfigoption()`."
function getthreadconfigoption(option::AbstractString, default = C_NULL)
    result = @gdal(CPLGetThreadLocalConfigOption::Ptr{UInt8},
        option::Cstring,
        default::Ptr{UInt8}
    )
    return (result == C_NULL) ? "" : unsafe_string(result)
end


text:
using ControlSystems, GtkInteract, Immerse

k = slider(1:0.1:10)

cg1 = cairographic(width=200,height=150);

w = window(vbox(k,grow(cg1)))

Reactive.preserve(map(k) do k
  fig1 = Immerse.plot(y=(1:10).^k)
  push!(cg1, fig1)
end)
w