text:
function clip(
        A::Array{Union{Missing, T}, 2} where T <: Number;
        x::Int64,
        y::Int64,
        distance::Int64
    )

    sizes = size(A)

    xlower = Int64(max(x - distance, 1))
    xupper = Int64(min(x + distance, sizes[2]))
    ylower = Int64(max(y - distance, 1))
    yupper = Int64(min(y + distance, sizes[1]))

    A_sub = A[ylower:yupper, xlower:xupper]

    dim1 = size(A_sub)[1]
    dim2 = size(A_sub)[2]

    new_x = min(distance + 1, x)
    new_y = min(distance + 1, y)

    dist_array = [sqrt((j - new_x)^2 + (i - new_y)^2) for i = 1:dim1, j = 1:dim2]

    A_sub[dist_array .> distance] .= missing

    A_sub
end

function get_targets(
        source_array::Array{Union{T, Missing}, 2} where T <: Number,
        arguments::Dict{String, Int64},
        precision::DataType
    )

    block_size = arguments["block_size"]
    block_radius = arguments["block_radius"]
    nrows = arguments["nrows"]
    ncols = arguments["ncols"]

    start = block_radius + 1

    xs = [start:block_size:(ncols - block_radius);]
    ys = [start:block_size:(nrows - block_radius);]

    ground_points = zeros(precision, (length(xs)*length(ys), 2))

    let
        c = 1
        for i = 1:length(xs)
            for j = 1:length(ys)
                ground_points[c, 1] = xs[i]
                ground_points[c, 2] = ys[j]
                c += 1
            end
        end
    end

    # create column ground_points[, 3] to hold source strengths for each target
    ground_points = cat(ground_points,
                        zeros(precision, size(ground_points)[1], 1);
                        dims = 2
                    )

    # populate ground_points[, 3] with sum of sources in block of size
    # block_size centered on each target
    for i = 1:size(ground_points)[1]
        xlower = Int64(ground_points[i, 1] - block_radius)
        xupper = min(Int64(ground_points[i, 1] + block_radius), ncols)
        ylower = Int64(ground_points[i, 2] - block_radius)
        yupper = min(Int64(ground_points[i, 2] + block_radius), nrows)

        # source_array is Array{Union{T, Missing}, 2} but should have no
        # missings, so not having skipmissings be okay
        ground_points[i, 3] = sum(skipmissing(source_array[ylower:yupper, xlower:xupper]))
    end

    # get rid of ground_points with strength equal to 0
    targets = ground_points[ground_points[:, 3] .> 0, 1:3]
    targets
end

# x and y defined by targets object. Ultimately the for loop will be done by
# iterating through rows of targets object
function get_source(
        source_array::Array{Union{Missing, T}, 2} where T <: Number,
        arguments::Dict{String, Int64},
        os_flags::OmniscapeFlags,
        condition_layers::ConditionLayers,
        conditions::Conditions,
        target::Target
    )

    block_radius = arguments["block_radius"]
    radius = arguments["radius"]
    buffer = arguments["buffer"]
    nrows = arguments["nrows"]
    ncols = arguments["ncols"]

    source_subset = clip(
        source_array,
        x = target.x_coord,
        y = target.y_coord,
        distance = radius
    )

    # Append missing if buffer > 0
    if buffer > 0
        ### Columns
        nrow_sub = size(source_subset)[1]
        left_col_num = max(0, min(buffer, target.x_coord - radius - 1))
        right_col_num = max(0, min(buffer, ncols - (target.x_coord + radius)))

        # Add left columns
        if left_col_num > 0
            source_subset = hcat(fill(missing, (nrow_sub, left_col_num)),
                                 source_subset)
        end
        # Add right columns
        if right_col_num > 0
            source_subset = hcat(source_subset,
                                 fill(missing, (nrow_sub, right_col_num)))
        end

        ### Rows
        ncol_sub = size(source_subset)[2]
        top_row_num = max(0, min(buffer, target.y_coord - radius - 1))
        bottom_row_num = max(0, min(buffer, nrows - (target.y_coord + radius)))

        # Add top rows
        if top_row_num > 0
            source_subset = vcat(fill(missing, (top_row_num, ncol_sub)),
                                 source_subset)
        end
        #Add bottom rows
        if bottom_row_num > 0
            source_subset = vcat(source_subset,
                                 fill(missing, (bottom_row_num, ncol_sub)))
        end
    end

    # Replace nodata vals with 0s
    source_subset[ismissing.(source_subset)] .= 0.0

    # Set any sources inside target block to 0
    xlower_sub = (radius + buffer  + 1) - block_radius
    xupper_sub = min((radius + buffer  + 1)  + block_radius, ncols)
    ylower_sub = (radius + buffer + 1)  - block_radius
    yupper_sub = min((radius + buffer  + 1)  + block_radius, nrows)

    source_subset[ylower_sub:yupper_sub, xlower_sub:xupper_sub] .= 0

    # allocate total current equal to target "strength", divide among sources
    # according to their source strengths
    source_sum = sum(source_subset[coalesce.(source_subset .> 0, false)])
    source_subset[source_subset .> 0] .=
        (source_subset[coalesce.(source_subset .> 0, false)] * target.amps) / source_sum

    if os_flags.conditional

        xlower_buffered = Int64(max(target.x_coord - radius - buffer, 1))
        xupper_buffered = Int64(min(target.x_coord + radius + buffer, ncols))
        ylower_buffered = Int64(max(target.y_coord - radius - buffer, 1))
        yupper_buffered = Int64(min(target.y_coord + radius + buffer, nrows))
        xlower = target.x_coord - block_radius
        xupper = min(target.x_coord + block_radius, ncols)
        ylower = target.y_coord - block_radius
        yupper = min(target.y_coord + block_radius, nrows)

        source_target_match!(source_subset,
                             arguments["n_conditions"],
                             condition_layers,
                             conditions,
                             ylower,
                             yupper,
                             xlower,
                             xupper,
                             ylower_buffered,
                             yupper_buffered,
                             xlower_buffered,
                             xupper_buffered
                             )
    end
    source_subset
end

function source_target_match!(source_subset::Array{Union{T, Missing}, 2} where T <: Number,
                              n_conditions::Int64,
                              condition_layers::ConditionLayers,
                              conditions::Conditions,
                              ylower::Int64,
                              yupper::Int64,
                              xlower::Int64,
                              xupper::Int64,
                              ylower_buffered::Int64,
                              yupper_buffered::Int64,
                              xlower_buffered::Int64,
                              xupper_buffered::Int64
                              )
    condition1_present = condition_layers.condition1_present
    condition1_future = condition_layers.condition1_future
    condition2_present = condition_layers.condition2_present
    condition2_future = condition_layers.condition2_future

    comparison1 = conditions.comparison1
    comparison2 = conditions.comparison2
    condition1_lower = conditions.condition1_lower
    condition1_upper = conditions.condition1_upper
    condition2_lower = conditions.condition2_lower
    condition2_upper = conditions.condition2_upper

    con1_present_subset = condition1_present[ylower_buffered:yupper_buffered,
                                             xlower_buffered:xupper_buffered]

    if comparison1 == "within"
      value1 = median(skipmissing(condition1_future[ylower:yupper, xlower:xupper]))
      source_subset[coalesce.(((con1_present_subset .- value1) .> condition1_upper) .|
          ((con1_present_subset .- value1) .< condition1_lower), false)] .= 0
    elseif comparison1 == "equals"
      value1 = mode(skipmissing(condition1_future[ylower:yupper, xlower:xupper]))
      source_subset[coalesce.(con1_present_subset .!= value1, false)] .= 0
    end

    if n_conditions == 2
      con2_present_subset = condition2_present[ylower_buffered:yupper_buffered,
                                               xlower_buffered:xupper_buffered]
      if comparison2 == "within"
          value2 = median(skipmissing(condition2_future[ylower:yupper, xlower:xupper]))
          source_subset[coalesce.(((con2_present_subset .- value2) .> condition2_upper) .|
              ((con2_present_subset .- value2) .< condition2_lower), false)] .= 0
      elseif comparison2 == "equals"
          value2 = mode(skipmissing(condition2_future[ylower:yupper, xlower:xupper]))
          source_subset[coalesce.(con2_present_subset .!= value2, false)] .= 0
      end
    end
end

function get_ground(arguments::Dict{String, Int64},
                    precision::DataType,
                    target::Target)
    radius = arguments["radius"]
    buffer = arguments["buffer"]
    distance = radius + buffer

    nrows = arguments["nrows"]
    ncols = arguments["ncols"]

    xlower = Int64(max(target.x_coord - radius - buffer, 1))
    xupper = Int64(min(target.x_coord + radius + buffer, ncols))
    ylower = Int64(max(target.y_coord - radius - buffer, 1))
    yupper = Int64(min(target.y_coord + radius + buffer, nrows))

    size_x = xupper - xlower + 1
    size_y = yupper - ylower + 1

    ground = fill(convert(precision, 0.0),
                  size_y,
                  size_x)

    new_x = min(distance + 1, target.x_coord)
    new_y = min(distance + 1, target.y_coord)

    ground[new_y, new_x] = Inf

    ground
end

function get_conductance(
        resistance::Array{Union{T, Missing}, 2} where T <: Number,
        arguments::Dict{String, Int64},
        target::Target,
        os_flags::OmniscapeFlags
    )

    radius = arguments["radius"]
    buffer = arguments["buffer"]

    resistance_clipped = clip(resistance, x = target.x_coord, y = target.y_coord, distance = radius + buffer)

    if os_flags.resistance_is_conductance
        conductance = resistance_clipped
    else
        conductance = 1 ./ resistance_clipped
    end

    convert(typeof(resistance_clipped), conductance)
end

function calculate_current(
        conductance::Array{Union{T, Missing}, 2} where T <: Number,
        source::Array{Union{T, Missing}, 2} where T <: Number,
        ground::Array{T, 2} where T <: Number,
        cs_flags::Circuitscape.RasterFlags,
        cs_cfg::Dict{String, String},
        T::DataType
    )
    V = Int64

    # Replace missings with -9999, then convert to Array{T, 2}
    # prior to circuitscape
    conductance[ismissing.(conductance)] .= -9999
    source[ismissing.(source)] .= -9999
    conductance = convert(Array{T, 2}, conductance)
    source = convert(Array{T, 2}, source)

    # get raster data
    cellmap = conductance
    polymap = Matrix{V}(undef, 0, 0)
    source_map = source
    ground_map = ground
    points_rc = (V[], V[], V[])
    strengths = Matrix{T}(undef, 0, 0)

    included_pairs = Circuitscape.IncludeExcludePairs(:undef,
                                                      V[],
                                                      Matrix{V}(undef,0,0))

    # This is just to satisfy type requirements, most of it not used
    hbmeta = Circuitscape.RasterMeta(size(cellmap)[2],
                                     size(cellmap)[1],
                                     0.,
                                     0.,
                                     1.,
                                     -9999.,
                                     Array{Float64, 1}(undef, 1),
                                     "")

    rasterdata = Circuitscape.RasterData(cellmap,
                                         polymap,
                                         source_map,
                                         ground_map,
                                         points_rc,
                                         strengths,
                                         included_pairs,
                                         hbmeta)

    # Generate advanced data
    data = Circuitscape.compute_advanced_data(rasterdata, cs_flags, cs_cfg)

    G = data.G
    nodemap = data.nodemap
    polymap = data.polymap
    hbmeta = data.hbmeta
    sources = data.sources
    grounds = data.grounds
    finitegrounds = data.finite_grounds
    cc = data.cc
    check_node = data.check_node
    source_map = data.source_map # Need it for one to all mode
    cellmap = data.cellmap

    f_local = Vector{eltype(G)}()
    voltages = Vector{eltype(G)}()
    outcurr = Circuitscape.alloc_map(hbmeta)

    for c in cc
        if check_node != -1 && !(check_node in c)
            continue
        end

        # a_local = laplacian(G[c, c])
        a_local = G[c,c]
        s_local = sources[c]
        g_local = grounds[c]

        if sum(s_local) == 0 || sum(g_local) == 0
            continue
        end

        if finitegrounds != [-9999.]
            f_local = finitegrounds[c]
        else
            f_local = finitegrounds
        end

        voltages = Circuitscape.multiple_solver(cs_cfg,
                                                data.solver,
                                                a_local,
                                                s_local,
                                                g_local,
                                                f_local)

        local_nodemap = Circuitscape.construct_local_node_map(nodemap,
                                                              c,
                                                              polymap)

        Circuitscape.accum_currents!(outcurr,
                                     voltages,
                                     cs_cfg,
                                     a_local,
                                     voltages,
                                     f_local,
                                     local_nodemap,
                                     hbmeta)
    end

    outcurr
end

function solve_target!(
        target::Target,
        int_arguments::Dict{String, Int64},
        source_strength::Array{Union{Missing, T}, 2} where T <: Number,
        resistance::Array{Union{Missing, T}, 2} where T <: Number,
        os_flags::OmniscapeFlags,
        cs_cfg::Dict{String, String},
        cs_flags::Circuitscape.RasterFlags,
        condition_layers::ConditionLayers,
        conditions::Conditions,
        correction_array::Array{T, 2} where T <: Number,
        cum_currmap::Array{T, 3} where T <: Number,
        fp_cum_currmap::Array{T, 3}  where T <: Number,
        precision::DataType
    )

    ## get source
    source = get_source(source_strength,
                        int_arguments,
                        os_flags,
                        condition_layers,
                        conditions,
                        target)

    ## get ground
    ground = get_ground(int_arguments,
                        precision,
                        target)

    ## get conductances for Omniscape
    conductance = get_conductance(resistance,
                                  int_arguments,
                                  target,
                                  os_flags)

    grid_size = size(source)

    ## Run circuitscape
    curr = calculate_current(conductance,
                             source,
                             ground,
                             cs_flags,
                             cs_cfg,
                             precision)

    ## If normalize = True, calculate null map and normalize
    if os_flags.compute_flow_potential
        null_conductance = convert(Array{Union{precision, Missing}, 2}, fill(1, grid_size))

        flow_potential = calculate_current(null_conductance,
                                           source,
                                           ground,
                                           cs_flags,
                                           cs_cfg,
                                           precision)
    end

    if os_flags.correct_artifacts && !(int_arguments["block_size"] == 1)
        correction_array2 = deepcopy(correction_array)
        lowerxcut = 1
        upperxcut = size(correction_array, 2)
        lowerycut = 1
        upperycut = size(correction_array, 1)

        if target.x_coord > int_arguments["ncols"] - (int_arguments["radius"] + int_arguments["buffer"])
            upperxcut = upperxcut - (upperxcut - grid_size[2])
        end

        if target.x_coord < int_arguments["radius"] + int_arguments["buffer"] + 1
            lowerxcut = upperxcut - grid_size[2] + 1
        end

        if target.y_coord > int_arguments["nrows"] - (int_arguments["radius"] + int_arguments["buffer"])
            upperycut = upperycut - (upperycut - grid_size[1])
        end

        if target.y_coord < int_arguments["radius"] + int_arguments["buffer"] + 1
            lowerycut = upperycut - grid_size[1] + 1
        end

        correction_array2 = correction_array[lowerycut:upperycut,
                                             lowerxcut:upperxcut]

        curr = curr .* correction_array2

        if os_flags.compute_flow_potential
            flow_potential = flow_potential .* correction_array2
        end
    end

    ## Accumulate values
    xlower = max(target.x_coord - int_arguments["radius"] - int_arguments["buffer"], 1)
    xupper = min(target.x_coord + int_arguments["radius"] + int_arguments["buffer"],
                 int_arguments["ncols"])
    ylower = max(target.y_coord - int_arguments["radius"] - int_arguments["buffer"], 1)
    yupper = min(target.y_coord + int_arguments["radius"] + int_arguments["buffer"],
                 int_arguments["nrows"])

    cum_currmap[ylower:yupper, xlower:xupper, threadid()] .=
        cum_currmap[ylower:yupper, xlower:xupper, threadid()] .+ curr

    if os_flags.compute_flow_potential
        fp_cum_currmap[ylower:yupper, xlower:xupper, threadid()] .=
            fp_cum_currmap[ylower:yupper, xlower:xupper, threadid()] .+ flow_potential
    end

end

function calc_correction(
        arguments::Dict{String, Int64},
        os_flags::OmniscapeFlags,
        cs_cfg::Dict{String, String},
        cs_flags::Circuitscape.RasterFlags,
        condition_layers::ConditionLayers,
        conditions::Conditions,
        precision::DataType
    )
    buffer = arguments["buffer"]
    # This may not apply seamlessly in the case (if I add the option) that source strengths
    # are not adjusted by target weight, but stay the same according to their
    # original values. Something to keep in mind...

    temp_source = convert(Array{Union{precision, Missing}, 2},
                          fill(1.0,
                               arguments["radius"] * 2 + buffer * 2 + 1,
                               arguments["radius"] * 2 + buffer * 2 + 1))

    source_null = clip(temp_source,
                       x = arguments["radius"] + buffer + 1,
                       y = arguments["radius"] + buffer + 1,
                       distance = arguments["radius"])

    # Append NoData (-9999) if buffer > 0
    if buffer > 0
        column_dims = (size(source_null)[1], buffer)
        # Add columns
        source_null = hcat(fill(missing, column_dims),
                                source_null,
                                fill(missing, column_dims))

        row_dims = (buffer, size(source_null)[2])
        # Add rows
        source_null = vcat(fill(missing, row_dims),
                                source_null,
                                fill(missing, row_dims))
    end
    n_sources = sum(source_null[(!).(ismissing.(source_null))])

    source_null[ismissing.(source_null)] .= 0.0
    source_null[arguments["radius"] + arguments["buffer"] + 1,
                arguments["radius"] + arguments["buffer"] + 1] = 0.0
    source_null[source_null .!= 0.0] .= 1 / (n_sources - 1)

    target = Target((arguments["radius"] + arguments["buffer"] + 1),
                    (arguments["radius"] + arguments["buffer"] + 1),
                    float(arguments["block_size"] ^ 2))

    source_blocked = get_source(temp_source,
                                arguments,
                                os_flags,
                                condition_layers,
                                conditions,
                                target)

    conductance = clip(temp_source,
                       x = arguments["radius"] + arguments["buffer"] + 1,
                       y = arguments["radius"] + arguments["buffer"] + 1,
                       distance = arguments["radius"] + arguments["buffer"])

    ground = fill(convert(precision, 0.0),
                  size(source_null))

    ground[arguments["radius"] + arguments["buffer"] + 1,
           arguments["radius"] + arguments["buffer"] + 1] = Inf

    block_null_current = calculate_current(conductance,
                                           source_blocked,
                                           ground,
                                           cs_flags,
                                           cs_cfg,
                                           precision)

    null_current =  calculate_current(conductance,
                                      source_null,
                                      ground,
                                      cs_flags,
                                      cs_cfg,
                                      precision)
    null_current_total = fill(convert(precision, 0.),
                              arguments["radius"] * 2 + arguments["buffer"] * 2 + arguments["block_size"],
                              arguments["radius"] * 2 + arguments["buffer"] * 2 + arguments["block_size"])

    for i in 1:arguments["block_size"]
        for j in 1:arguments["block_size"]
            null_current_total[i:(i + arguments["radius"] * 2 + arguments["buffer"] * 2),
                               j:(j + arguments["radius"] * 2 + arguments["buffer"] * 2)] += null_current
        end
    end

    null_current_total = null_current_total[(arguments["block_radius"] + 1):(size(null_current, 1) + arguments["block_radius"]),
                                            (arguments["block_radius"] + 1):(size(null_current, 2) + arguments["block_radius"])]

    null_current_total[block_null_current .== 0.] .= 0

    null_current_total[null_current_total .== 0.0] .= 1.0
    block_null_current[block_null_current .== 0.0] .= 1.0

    correction = null_current_total ./ block_null_current

    correction
end

function get_omniscape_flags(cfg::Dict{String, String})
    OmniscapeFlags(
        cfg["calc_flow_potential"] in TRUELIST,
        cfg["calc_normalized_current"] in TRUELIST,
        cfg["calc_flow_potential"] in TRUELIST || cfg["calc_normalized_current"] in TRUELIST,
        cfg["write_raw_currmap"] in TRUELIST,
        cfg["parallelize"] in TRUELIST,
        cfg["correct_artifacts"] in TRUELIST,
        cfg["source_from_resistance"] in TRUELIST,
        cfg["conditional"] in TRUELIST,
        cfg["mask_nodata"] in TRUELIST,
        cfg["resistance_is_conductance"] in TRUELIST,
        cfg["write_as_tif"] in TRUELIST,
        cfg["allow_different_projections"] in TRUELIST,
        cfg["reclassify_resistance"] in TRUELIST,
        cfg["write_reclassified_resistance"] in TRUELIST
    )
end

# Calculate the source layer using resistance surface and arguments from cfg
function source_from_resistance(resistance::Array{Union{T, Missing}, 2} where T <: Number,
                                cfg::Dict{String, String},
                                reclass_table::Array{Union{T, Missing}, 2} where T <: Number)
    full_cfg = init_cfg()
    update_cfg!(full_cfg, cfg)
    r_cutoff = parse(Float64, full_cfg["r_cutoff"])
    precision = full_cfg["precision"] in SINGLE ? Float32 : Float64
    reclassify = full_cfg["reclassify_resistance"] in TRUELIST

    if reclassify
        resistance_for_source = deepcopy(resistance)
        reclassify_resistance!(resistance_for_source, reclass_table)
    else
        resistance_for_source = resistance
    end

    source_strength = deepcopy(resistance_for_source)

    if full_cfg["resistance_is_conductance"] ∉ TRUELIST
        source_strength = Array{Union{precision, Missing}, 2}(1 ./ source_strength)
    end
    source_strength[coalesce.(source_strength .< 1/r_cutoff, true)] .= 0.0 # handles replacing NoData with 0 as well

    source_strength
end

function reclassify_resistance!(resistance::Array{Union{T, Missing}, 2} where T <: Number,
                                reclass_table::Array{Union{T, Missing}, 2} where T <: Number)
    resistance_old = deepcopy(resistance)
    for i in 1:(size(reclass_table)[1])
        resistance[coalesce.(resistance_old .== reclass_table[i, 1], false)] .= reclass_table[i, 2]
    end
    resistance_old = nothing # remove from memory
end

function convert_and_fill_missing(A::Array{T, 2} where T <: Number,
                                  precision::DataType)
    A = convert(Array{Union{precision, Missing}, 2}, A)
    A[A .== -9999] .= missing

    A
end

function arrays_equal(A::Array{Union{T, Missing}, 2} where T <: Number,
                      B::Array{Union{T, Missing}, 2} where T <: Number)
    # Check that non-missing entries are equal
    A[ismissing.(A)] .= -9999
    B[ismissing.(B)] .= -9999

    isapprox(Array{Float64, 2}(A), Array{Float64, 2}(B); rtol = 1e-6)
end


text:
"""
This module provides the [`getconnection`](@ref) function with support for:
- Opening TCP and SSL connections.
- Reusing connections for multiple Request/Response Messages,
- Pipelining Request/Response Messages. i.e. allowing a new Request to be
  sent before previous Responses have been read.

This module defines a [`Connection`](@ref)
struct to manage pipelining and connection reuse and a
[`Transaction`](@ref)`<: IO` struct to manage a single
pipelined request. Methods are provided for `eof`, `readavailable`,
`unsafe_write` and `close`.
This allows the `Transaction` object to act as a proxy for the
`TCPSocket` or `SSLContext` that it wraps.

The [`POOL`](@ref) is used to manage connection pooling. Connections
are identified by their host, port, pipeline limit, whether they require
ssl verification, and whether they are a client or server connection.
If a subsequent request matches these properties of a previous connection
and limits are respected (reuse limit, idle timeout), and it wasn't otherwise
remotely closed, a connection will be reused. Transactions pipeline their
requests and responses concurrently on a Connection by calling `startwrite`
and `closewrite`, with corresponding `startread` and `closeread`.
"""
module ConnectionPool

export Connection, Transaction,
       getconnection, getrawstream, inactiveseconds

using ..IOExtras, ..Sockets

import ..@debug, ..@debugshow, ..DEBUG_LEVEL, ..taskid
import ..@require, ..precondition_error, ..@ensure, ..postcondition_error
using MbedTLS: SSLConfig, SSLContext, setup!, associate!, hostname!, handshake!

const default_connection_limit = 8
const default_pipeline_limit = 16
const nolimit = typemax(Int)

# certain operations, like locking Channels and Conditions
# is only supported in >= 1.3
macro v1_3(expr, elses=nothing)
    esc(quote
        @static if VERSION >= v"1.3"
            $expr
        else
            $elses
        end
    end)
end

@static if VERSION >= v"1.3"
    const Cond = Threads.Condition
else
    const Cond = Condition
end

"""
    Connection{T <: IO}

A `TCPSocket` or `SSLContext` connection to a HTTP `host` and `port`.

Fields:
- `host::String`
- `port::String`, exactly as specified in the URI (i.e. may be empty).
- `pipeline_limit`, number of requests to send before waiting for responses.
- `idle_timeout`, No. of seconds to maintain connection after last transaction.
- `peerport`, remote TCP port number (used for debug messages).
- `localport`, local TCP port number (used for debug messages).
- `io::T`, the `TCPSocket` or `SSLContext.
- `clientconnection::Bool`, whether the Connection was created from client code (as opposed to server code)
- `buffer::IOBuffer`, left over bytes read from the connection after
   the end of a response header (or chunksize). These bytes are usually
   part of the response body.
- `sequence`, number of most recent `Transaction`.
- `writecount`, number of Messages that have been written, protected by `writelock`
- `writelock`, lock writecount and writebusy, and signal that `writecount` was incremented.
- `writebusy`, whether a Transaction currently holds the Connection write lock, protected by `writelock`
- `readcount`, number of Messages that have been read, protected by `readlock`
- `readlock`, lock readcount and readbusy, and signal that `readcount` was incremented.
- `readbusy`, whether a Transaction currently holds the Connection read lock, protectecd by `readlock`
- `timestamp`, time data was last received.
"""
mutable struct Connection{T <: IO}
    host::String
    port::String
    pipeline_limit::Int
    idle_timeout::Int
    require_ssl_verification::Bool
    peerport::UInt16 # debug only
    localport::UInt16 # debug only
    io::T
    clientconnection::Bool
    buffer::IOBuffer
    sequence::Threads.Atomic{Int}
    writecount::Int
    writelock::Cond # protects the writecount and writebusy fields, notifies on closewrite
    writebusy::Bool
    readcount::Int
    readlock::Cond # protects the readcount and readbusy fields, notifies on closeread
    readbusy::Bool
    timestamp::Float64
    closelock::ReentrantLock
    closed::Bool
end

"""
    hashconn

Used for "hashing" a Connection object on just the key properties necessary for determining
connection re-useability. That is, when a new request calls `getconnection`, we take the
request parameters of what socket type, the host and port, any pipeline_limit and if ssl
verification is required, and if an existing Connection was already created with the exact
same parameters, we can re-use it (as long as it's not already being used, obviously).
"""
function hashconn end

hashconn(x::Connection{T}) where {T} = hashconn(T, x.host, x.port, x.pipeline_limit, x.require_ssl_verification, x.clientconnection)
hashconn(T, host, port, pipeline_limit, require_ssl_verification, client) = hash(T, hash(host, hash(port, hash(pipeline_limit, hash(require_ssl_verification, hash(client, UInt(0)))))))

"""
    Transaction

A single pipelined HTTP Request/Response transaction.

Fields:
 - `c`, the shared [`Connection`](@ref) used for this `Transaction`.
 - `sequence::Int`, identifies this `Transaction` among the others that share `c`.
 - `writebusy::Bool`, whether this Transaction holds its parent Connection write lock, protected by c.writelock
 - `readbusy::Bool`, whether this Transaction holds its parent Connection read lock, protected by c.readlock
"""
mutable struct Transaction{T <: IO} <: IO
    c::Connection{T}
    sequence::Int
    writebusy::Bool
    readbusy::Bool
end

Connection(host::AbstractString, port::AbstractString,
           pipeline_limit::Int, idle_timeout::Int,
           require_ssl_verification::Bool, io::T, client=true) where T <: IO =
    Connection{T}(host, port,
                  pipeline_limit, idle_timeout,
                  require_ssl_verification,
                  peerport(io), localport(io),
                  io, client, PipeBuffer(), Threads.Atomic{Int}(0),
                  0, Cond(), false,
                  0, Cond(), false,
                  time(), ReentrantLock(), false)

Connection(io; require_ssl_verification::Bool=true) =
    Connection("", "", default_pipeline_limit, 0, require_ssl_verification, io, false)

Transaction(c::Connection{T}) where T <: IO =
    Transaction{T}(c, (Threads.atomic_add!(c.sequence, 1)), false, false)

function client_transaction(c)
    t = Transaction(c)
    startwrite(t)
    return t
end

getrawstream(t::Transaction) = t.c.io

inactiveseconds(t::Transaction) = inactiveseconds(t.c)

function inactiveseconds(c::Connection)::Float64
    return time() - c.timestamp
end

Base.unsafe_write(t::Transaction, p::Ptr{UInt8}, n::UInt) =
    unsafe_write(t.c.io, p, n)

Base.isopen(c::Connection) = isopen(c.io)

Base.isopen(t::Transaction) = isopen(t.c) &&
                              readcount(t.c) <= t.sequence &&
                              writecount(t.c) <= t.sequence

writebusy(c::Connection) = @v1_3 lock(() -> c.writebusy, c.writelock) c.writebusy
writecount(c::Connection) = @v1_3 lock(() -> c.writecount, c.writelock) c.writecount
readbusy(c::Connection) = @v1_3 lock(() -> c.readbusy, c.readlock) c.readbusy
readcount(c::Connection) = @v1_3 lock(() -> c.readcount, c.readlock) c.readcount

writebusy(t::Transaction) = @v1_3 lock(() -> t.writebusy, t.c.writelock) t.writebusy
readbusy(t::Transaction) = @v1_3 lock(() -> t.readbusy, t.c.readlock) t.readbusy

"""
Is `c` currently in use or expecting a response to request already sent?
"""
isbusy(c::Connection) = isopen(c) && (writebusy(c) || readbusy(c) ||
                                      writecount(c) > readcount(c))

function Base.eof(t::Transaction)
    @require isreadable(t) || !isopen(t)
    if bytesavailable(t) > 0
        return false
    end             ;@debug 4 "eof(::Transaction) -> eof($(typeof(t.c.io))): $t"
    return eof(t.c.io)
end

Base.bytesavailable(t::Transaction) = bytesavailable(t.c)
Base.bytesavailable(c::Connection) = bytesavailable(c.buffer) +
                                     bytesavailable(c.io)

Base.isreadable(t::Transaction) = readbusy(t)
Base.iswritable(t::Transaction) = writebusy(t)

function Base.read(t::Transaction, nb::Int)
    nb = min(nb, bytesavailable(t))
    bytes = Base.StringVector(nb)
    unsafe_read(t, pointer(bytes), nb)
    return bytes
end

function Base.read(t::Transaction, ::Type{UInt8})
    if bytesavailable(t.c.buffer) == 0
        read_to_buffer(t)
    end
    return read(t.c.buffer, UInt8)
end

function Base.unsafe_read(t::Transaction, p::Ptr{UInt8}, n::UInt)
    l = bytesavailable(t.c.buffer)
    if l > 0
        nb = min(l,n)
        unsafe_read(t.c.buffer, p, nb)
        p += nb
        n -= nb
        @debug 4 "↩  read $nb-bytes from buffer."
        t.c.timestamp = time()
    end
    if n > 0
        unsafe_read(t.c.io, p, n)
        @debug 4 "⬅  read $n-bytes from $(typeof(t.c.io))"
        t.c.timestamp = time()
    end
    return nothing
end

function read_to_buffer(t::Transaction, sizehint=4096)
    buf = t.c.buffer

    # Reset the buffer if it is empty.
    if bytesavailable(buf) == 0
        buf.size = 0
        buf.ptr = 1
    end

    # Wait for data.
    if eof(t.c.io)
        throw(EOFError())
    end

    # Read from stream into buffer.
    n = min(sizehint, bytesavailable(t.c.io))
    buf = t.c.buffer
    Base.ensureroom(buf, n)
    unsafe_read(t.c.io, pointer(buf.data, buf.size + 1), n)
    buf.size += n
end

"""
Read until `find_delimiter(bytes)` returns non-zero.
Return view of bytes up to the delimiter.
"""
function Base.readuntil(t::Transaction, f::Function #=Vector{UInt8} -> Int=#,
                                        sizehint=4096)::ByteView
    buf = t.c.buffer
    if bytesavailable(buf) == 0
        read_to_buffer(t, sizehint)
    end
    while (bytes = readuntil(buf, f)) === nobytes
        read_to_buffer(t, sizehint)
    end
    return bytes
end

"""
    startwrite(::Transaction)

Wait for prior pending writes to complete.
"""
function IOExtras.startwrite(t::Transaction)
    @require !iswritable(t)

    @v1_3 lock(t.c.writelock)
    try
        while writecount(t.c) != t.sequence
            @debug 1 "⏳  Wait write: $t"
            wait(t.c.writelock)
        end
        t.writebusy = true
        t.c.writebusy = true
        @ensure iswritable(t)
        @debug 2 "👁  Start write:$t"
    finally
        @v1_3 unlock(t.c.writelock)
    end

    return
end

"""
    closewrite(::Transaction)

Signal that an entire Request Message has been written to the `Transaction`.
"""
function IOExtras.closewrite(t::Transaction)
    @require iswritable(t)

    @v1_3 lock(t.c.writelock)
    try
        t.writebusy = false
        t.c.writecount += 1          ;@debug 2 "🗣  Write done: $t"
        t.c.writebusy = false
        notify(t.c.writelock)
        @ensure !iswritable(t)
    finally
        @v1_3 unlock(t.c.writelock)
    end
    release(t.c)

    return
end

"""
    startread(::Transaction)

Wait for prior pending reads to complete.
"""
function IOExtras.startread(t::Transaction)
    @require !isreadable(t)

    t.c.timestamp = time()
    @v1_3 lock(t.c.readlock)
    try
        while readcount(t.c) != t.sequence
            @debug 1 "⏳  Wait read: $t"
            wait(t.c.readlock)
        end
        t.readbusy = true
        t.c.readbusy = true
        @debug 2 "👁  Start read: $t"
        @ensure isreadable(t)
    finally
        @v1_3 unlock(t.c.readlock)
    end

    return
end

"""
    closeread(::Transaction)

Signal that an entire Response Message has been read from the `Transaction`.

Increment `readcount` and wake up tasks waiting in `startread`.
"""
function IOExtras.closeread(t::Transaction)
    @require isreadable(t)

    @v1_3 lock(t.c.readlock)
    try
        t.readbusy = false
        t.c.readcount += 1         ;@debug 2 "✉  Read done:  $t"
        t.c.readbusy = false
        notify(t.c.readlock)
        @ensure !isreadable(t)
        if !isbusy(t.c)
            @async monitor_idle_connection(t.c)
        end
    finally
        @v1_3 unlock(t.c.readlock)
    end
    release(t.c)

    return
end

"""
Wait for `c` to receive data or reach EOF.
Close `c` on EOF or if response data arrives when no request was sent.
"""
function monitor_idle_connection(c::Connection)
    if eof(c.io)                                  ;@debug 2 "💀  Closed:     $c"
        close(c.io)
    elseif !isbusy(c)                             ;@debug 1 "😈  Idle RX!!:  $c"
        close(c.io)
    end
end

function monitor_idle_connection(c::Connection{SSLContext})
    # MbedTLS.jl monitors idle connections for TLS close_notify messages.
    # https://github.com/JuliaWeb/MbedTLS.jl/pull/145
end

Base.wait_close(t::Transaction) = Base.wait_close(tcpsocket(t.c.io))

function Base.close(t::Transaction)
    close(t.c)
    if iswritable(t)
        closewrite(t)
    end
    if isreadable(t)
        closeread(t)
    end
    return
end

function Base.close(c::Connection)
    close(c.io)
    if bytesavailable(c) > 0
        purge(c)
    end
    return
end

"""
    purge(::Connection)

Remove unread data from a `Connection`.
"""
function purge(c::Connection)
    @require !isopen(c.io)
    while !eof(c.io)
        readavailable(c.io)
    end
    c.buffer.size = 0
    c.buffer.ptr = 1
    @ensure bytesavailable(c) == 0
end

"""
    closeall()

Close all connections in `pool`.
"""
function closeall()
    lock(POOL.lock) do
        for pod in values(POOL.conns)
            @v1_3 lock(pod.conns)
            while isready(pod.conns)
                close(take!(pod.conns))
            end
            pod.numactive = 0
            @v1_3 unlock(pod.conns)
        end
    end
    return
end

mutable struct Pod
    conns::Channel{Connection}
    numactive::Int
end

Pod() = Pod(Channel{Connection}(Inf), 0)

function decr!(pod::Pod)
    @v1_3 @assert islocked(pod.conns.cond_take)
    pod.numactive -= 1
    return
end

function incr!(pod::Pod)
    @v1_3 @assert islocked(pod.conns.cond_take)
    pod.numactive += 1
    return
end

function release(c::Connection)
    h = hashconn(c)
    if haskey(POOL.conns, h)
        pod = getpod(POOL, h)
        @debug 2 "returning connection to pod: $c"
        put!(pod.conns, c)
    end
    return
end

# "release" a Connection, without returning to pod for re-use
# used for https proxy tunnel upgrades which shouldn't be reused
function kill!(c::Connection)
    h = hashconn(c)
    if haskey(POOL.conns, h)
        pod = getpod(POOL, h)
        @v1_3 lock(pod.conns)
        try
            decr!(pod)
        finally
            @v1_3 unlock(pod.conns)
        end
    end
    return
end

struct Pool
    lock::ReentrantLock
    conns::Dict{UInt, Pod}
end

const POOL = Pool(ReentrantLock(), Dict{UInt, Pod}())

function getpod(pool::Pool, x)
    lock(pool.lock) do
        get!(() -> Pod(), pool.conns, x)
    end
end

"""
    getconnection(type, host, port) -> Connection

Find a reusable `Connection` in the `pool`,
or create a new `Connection` if required.
"""
function getconnection(::Type{Transaction{T}},
                       host::AbstractString,
                       port::AbstractString;
                       connection_limit::Int=default_connection_limit,
                       pipeline_limit::Int=default_pipeline_limit,
                       idle_timeout::Int=0,
                       reuse_limit::Int=nolimit,
                       require_ssl_verification::Bool=true,
                       kw...)::Transaction{T} where T <: IO
    pod = getpod(POOL, hashconn(T, host, port, pipeline_limit, require_ssl_verification, true))
    @v1_3 lock(pod.conns)
    try
        while isready(pod.conns)
            conn = take!(pod.conns)
            if isvalid(pod, conn, reuse_limit, pipeline_limit)
                # this is a reuseable connection, so use it
                @debug 2 "1 reusing connection: $conn"
                return client_transaction(conn)
            end
        end
        # If there are not too many connections to this host:port,
        # create a new connection...
        if pod.numactive < connection_limit
            return newconnection(pod, T, host, port, pipeline_limit,
                require_ssl_verification, idle_timeout; kw...)
        end
        # wait for a Connection to be released
        while true
            conn = take!(pod.conns)
            if isvalid(pod, conn, reuse_limit, pipeline_limit)
                # this is a reuseable connection, so use it
                @debug 2 "2 reusing connection: $conn"
                return client_transaction(conn)
            elseif pod.numactive < connection_limit
                return newconnection(pod, T, host, port, pipeline_limit,
                    require_ssl_verification, idle_timeout; kw...)
            end
        end
    finally
        @v1_3 unlock(pod.conns)
    end
end

function isvalid(pod, conn, reuse_limit, pipeline_limit)
    # Close connections that have reached the reuse limit...
    if reuse_limit != nolimit
        if readcount(conn) >= reuse_limit && !readbusy(conn)
            @debug 2 "💀 overuse:         $conn"
            close(conn.io)
        end
    end
    # Close connections that have reached the timeout limit...
    if conn.idle_timeout > 0
        if !isbusy(conn) && inactiveseconds(conn) > conn.idle_timeout
            @debug 2 "💀 idle timeout:    $conn"
            close(conn.io)
        end
    end
    # For closed connections, we decrease active count in pod, and "continue"
    # which effectively drops the connection
    if !isopen(conn.io)
        close(conn)
        lock(conn.closelock) do
            if !conn.closed
                conn.closed = true
                decr!(pod)
            end
        end
        return false
    end
    # If we've hit our pipeline_limit, can't use this one, but don't close
    if (writecount(conn) - readcount(conn)) >= pipeline_limit + 1
        return false
    end

    return !writebusy(conn)
end

function newconnection(pod, T, host, port, pipeline_limit, require_ssl_verification, idle_timeout; kw...)
    io = getconnection(T, host, port;
            require_ssl_verification=require_ssl_verification, kw...)
    c = Connection(host, port, pipeline_limit, idle_timeout, require_ssl_verification, io)
    incr!(pod)
    @debug 1 "🔗  New:            $c"
    return client_transaction(c)
end

function keepalive!(tcp)
    @debug 2 "setting keepalive on tcp socket"
    err = ccall(:uv_tcp_keepalive, Cint, (Ptr{Nothing}, Cint, Cuint),
                                          tcp.handle, 1, 1)
    err != 0 && error("error setting keepalive on socket")
    return
end

struct ConnectTimeout <: Exception
    host
    port
end

function getconnection(::Type{TCPSocket},
                       host::AbstractString,
                       port::AbstractString;
                       keepalive::Bool=false,
                       connect_timeout::Int=0,
                       kw...)::TCPSocket

    p::UInt = isempty(port) ? UInt(80) : parse(UInt, port)

    @debug 2 "TCP connect: $host:$p..."

    if connect_timeout == 0
        tcp = Sockets.connect(host == "localhost" ? ip"127.0.0.1" : Sockets.getalladdrinfo(host)[1], p)
        keepalive && keepalive!(tcp)
        return tcp
    end

    tcp = Sockets.TCPSocket()
    Sockets.bind(tcp, Sockets.getaddrinfo(host), p)

    timeout = Ref{Bool}(false)
    @async begin
        sleep(connect_timeout)
        if tcp.status == Base.StatusConnecting
            timeout[] = true
            tcp.status = Base.StatusClosing
            ccall(:jl_forceclose_uv, Nothing, (Ptr{Nothing},), tcp.handle)
            #close(tcp)
        end
    end
    try
        Base.wait_connected(tcp)
    catch e
        if timeout[]
            throw(ConnectTimeout(host, port))
        end
        rethrow(e)
    end

    keepalive && keepalive!(tcp)
    return tcp
end

const nosslconfig = SSLConfig()
default_sslconfig = nothing
noverify_sslconfig = nothing

function global_sslconfig(require_ssl_verification::Bool)::SSLConfig
    global default_sslconfig
    global noverify_sslconfig
    if default_sslconfig === nothing
        default_sslconfig = SSLConfig(true)
        noverify_sslconfig = SSLConfig(false)
    end
    return require_ssl_verification ? default_sslconfig : noverify_sslconfig
end

function getconnection(::Type{SSLContext},
                       host::AbstractString,
                       port::AbstractString;
                       kw...)::SSLContext

    port = isempty(port) ? "443" : port
    @debug 2 "SSL connect: $host:$port..."
    tcp = getconnection(TCPSocket, host, port; kw...)
    return sslconnection(tcp, host; kw...)
end

function sslconnection(tcp::TCPSocket, host::AbstractString;
                       require_ssl_verification::Bool=true,
                       sslconfig::SSLConfig=nosslconfig,
                       kw...)::SSLContext

    if sslconfig === nosslconfig
        sslconfig = global_sslconfig(require_ssl_verification)
    end

    io = SSLContext()
    setup!(io, sslconfig)
    associate!(io, tcp)
    hostname!(io, host)
    handshake!(io)
    return io
end

function sslupgrade(t::Transaction{TCPSocket},
                    host::AbstractString;
                    require_ssl_verification::Bool=true,
                    kw...)::Transaction{SSLContext}
    tls = sslconnection(t.c.io, host;
                        require_ssl_verification=require_ssl_verification,
                        kw...)
    c = Connection(tls; require_ssl_verification=require_ssl_verification)
    kill!(t.c)
    return client_transaction(c)
end

function Base.show(io::IO, c::Connection)
    nwaiting = applicable(tcpsocket, c.io) ? bytesavailable(tcpsocket(c.io)) : 0
    print(
        io,
        tcpstatus(c), " ",
        lpad(writecount(c), 3),"↑", writebusy(c) ? "🔒  " : "   ",
        lpad(readcount(c), 3), "↓", readbusy(c) ? "🔒 " : "  ",
        "$(lpad(round(Int, time() - c.timestamp), 3))s ",
        c.host, ":",
        c.port != "" ? c.port : Int(c.peerport), ":", Int(c.localport),
        " ≣", c.pipeline_limit,
        bytesavailable(c.buffer) > 0 ?
            " $(bytesavailable(c.buffer))-byte excess" : "",
        nwaiting > 0 ? " $nwaiting bytes waiting" : "",
        DEBUG_LEVEL[] > 1 && applicable(tcpsocket, c.io) ?
            " $(Base._fd(tcpsocket(c.io)))" : "")
end

Base.show(io::IO, t::Transaction) = print(io, "T$(rpad(t.sequence,2)) ", t.c)

function tcpstatus(c::Connection)
    if !applicable(tcpsocket, c.io)
        return ""
    end
    s = Base.uv_status_string(tcpsocket(c.io))
        if s == "connecting" return "🔜🔗"
    elseif s == "open"       return "🔗 "
    elseif s == "active"     return "🔁 "
    elseif s == "paused"     return "⏸ "
    elseif s == "closing"    return "🔜💀"
    elseif s == "closed"     return "💀 "
    else
        return s
    end
end

end # module ConnectionPool

ext: jl


text:
module superoperators

export SuperOperator, DenseSuperOperator, SparseSuperOperator,
        spre, spost, liouvillian, exp

import Base: ==, *, /, +, -
import ..bases
import ..states
import SparseArrays: sparse

using ..bases, ..operators, ..operators_dense, ..operators_sparse
using SparseArrays

"""
Base class for all super operator classes.

Super operators are bijective mappings from operators given in one specific
basis to operators, possibly given in respect to another, different basis.
To embed super operators in an algebraic framework they are defined with a
left hand basis `basis_l` and a right hand basis `basis_r` where each of
them again consists of a left and right hand basis.
```math
A_{bl_1,bl_2} = S_{(bl_1,bl_2) ↔ (br_1,br_2)} B_{br_1,br_2}
\\\\
A_{br_1,br_2} = B_{bl_1,bl_2} S_{(bl_1,bl_2) ↔ (br_1,br_2)}
```
"""
abstract type SuperOperator{B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis}} end

"""
    DenseSuperOperator(b1[, b2, data])

SuperOperator stored as dense matrix.
"""
mutable struct DenseSuperOperator{B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis},T<:Matrix{ComplexF64}} <: SuperOperator{B1,B2}
    basis_l::B1
    basis_r::B2
    data::T
    function DenseSuperOperator{BL,BR,T}(basis_l::BL, basis_r::BR, data::T) where {BL<:Tuple{Basis,Basis},BR<:Tuple{Basis,Basis},T<:Matrix{ComplexF64}}
        if length(basis_l[1])*length(basis_l[2]) != size(data, 1) ||
           length(basis_r[1])*length(basis_r[2]) != size(data, 2)
            throw(DimensionMismatch())
        end
        new(basis_l, basis_r, data)
    end
end
DenseSuperOperator{BL,BR}(basis_l::BL, basis_r::BR, data::T) where {BL<:Tuple{Basis,Basis},BR<:Tuple{Basis,Basis},T<:Matrix{ComplexF64}} = DenseSuperOperator{BL,BR,T}(basis_l, basis_r, data)
DenseSuperOperator(basis_l::BL, basis_r::BR, data::T) where {BL<:Tuple{Basis,Basis},BR<:Tuple{Basis,Basis},T<:Matrix{ComplexF64}} = DenseSuperOperator{BL,BR,T}(basis_l, basis_r, data)

function DenseSuperOperator(basis_l::Tuple{Basis, Basis}, basis_r::Tuple{Basis, Basis})
    Nl = length(basis_l[1])*length(basis_l[2])
    Nr = length(basis_r[1])*length(basis_r[2])
    data = zeros(ComplexF64, Nl, Nr)
    DenseSuperOperator(basis_l, basis_r, data)
end

"""
    SparseSuperOperator(b1[, b2, data])

SuperOperator stored as sparse matrix.
"""
mutable struct SparseSuperOperator{B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis},T<:SparseMatrixCSC{ComplexF64,Int}} <: SuperOperator{B1,B2}
    basis_l::B1
    basis_r::B2
    data::T
    function SparseSuperOperator{B1,B2,T}(basis_l::B1, basis_r::B2, data::T) where {B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis},T<:SparseMatrixCSC{ComplexF64,Int}}
        if length(basis_l[1])*length(basis_l[2]) != size(data, 1) ||
           length(basis_r[1])*length(basis_r[2]) != size(data, 2)
            throw(DimensionMismatch())
        end
        new(basis_l, basis_r, data)
    end
end
SparseSuperOperator(basis_l::B1, basis_r::B2, data::T) where {B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis},T<:SparseMatrixCSC{ComplexF64,Int}} = SparseSuperOperator{B1,B2,T}(basis_l, basis_r, data)

function SparseSuperOperator(basis_l::Tuple{Basis, Basis}, basis_r::Tuple{Basis, Basis})
    Nl = length(basis_l[1])*length(basis_l[2])
    Nr = length(basis_r[1])*length(basis_r[2])
    data = spzeros(ComplexF64, Nl, Nr)
    SparseSuperOperator(basis_l, basis_r, data)
end

SuperOperator(basis_l, basis_r, data::SparseMatrixCSC{ComplexF64, Int}) = SparseSuperOperator(basis_l, basis_r, data)
SuperOperator(basis_l, basis_r, data::Matrix{ComplexF64}) = DenseSuperOperator(basis_l, basis_r, data)
SuperOperator{B1,B2}(basis_l::B1, basis_r::B2, data::SparseMatrixCSC{ComplexF64, Int}) where {B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis}} = SparseSuperOperator(basis_l, basis_r, data)
SuperOperator{B1,B2}(basis_l::B1, basis_r::B2, data::Matrix{ComplexF64}) where {B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis}} = DenseSuperOperator(basis_l, basis_r, data)

Base.copy(a::T) where {T<:SuperOperator} = T(a.basis_l, a.basis_r, copy(a.data))

operators.dense(a::SparseSuperOperator) = DenseSuperOperator(a.basis_l, a.basis_r, Matrix(a.data))
operators.dense(a::DenseSuperOperator) = copy(a)

sparse(a::DenseSuperOperator) = SparseSuperOperator(a.basis_l, a.basis_r, sparse(a.data))
sparse(a::SparseSuperOperator) = copy(a)

==(a::T, b::T) where {T<:SuperOperator} = a.data == b.data
==(a::SuperOperator, b::SuperOperator) = false

Base.length(a::SuperOperator) = length(a.basis_l[1])*length(a.basis_l[2])*length(a.basis_r[1])*length(a.basis_r[2])
bases.samebases(a::SuperOperator, b::SuperOperator) = samebases(a.basis_l[1], b.basis_l[1]) && samebases(a.basis_l[2], b.basis_l[2]) &&
                                                      samebases(a.basis_r[1], b.basis_r[1]) && samebases(a.basis_r[2], b.basis_r[2])
bases.multiplicable(a::SuperOperator, b::SuperOperator) = multiplicable(a.basis_r[1], b.basis_l[1]) && multiplicable(a.basis_r[2], b.basis_l[2])
bases.multiplicable(a::SuperOperator, b::AbstractOperator) = multiplicable(a.basis_r[1], b.basis_l) && multiplicable(a.basis_r[2], b.basis_r)

# Arithmetic operations
function *(a::SuperOperator{B1,B2}, b::DenseOperator{BL,BR}) where {BL<:Basis,BR<:Basis,B1<:Tuple{Basis,Basis},B2<:Tuple{BL,BR}}
    data = a.data*reshape(b.data, length(b.data))
    return DenseOperator(a.basis_l[1], a.basis_l[2], reshape(data, length(a.basis_l[1]), length(a.basis_l[2])))
end
function *(a::SuperOperator{B1,B2}, b::SparseOperator{BL,BR}) where {BL<:Basis,BR<:Basis,B1<:Tuple{Basis,Basis},B2<:Tuple{BL,BR}}
    data = a.data*reshape(b.data, length(b.data))
    return SparseOperator(a.basis_l[1], a.basis_l[2], reshape(data, length(a.basis_l[1]), length(a.basis_l[2])))
end

function *(a::SuperOperator{B1,B2}, b::SuperOperator{B2,B3}) where {B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis},B3<:Tuple{Basis,Basis}}
    return SuperOperator{B1,B3}(a.basis_l, b.basis_r, a.data*b.data)
end

*(a::T, b::Number) where T<:SuperOperator = T(a.basis_l, a.basis_r, a.data*b)
*(a::Number, b::T) where T<:SuperOperator = T(b.basis_l, b.basis_r, a*b.data)

/(a::T, b::Number) where T<:SuperOperator = T(a.basis_l, a.basis_r, a.data/b)

+(a::SuperOperator{B1,B2}, b::SuperOperator{B1,B2}) where {B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis}} = SuperOperator{B1,B2}(a.basis_l, a.basis_r, a.data+b.data)
+(a::SuperOperator, b::SuperOperator) = throw(bases.IncompatibleBases())

-(a::SuperOperator{B1,B2}, b::SuperOperator{B1,B2}) where {B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis}} = SuperOperator{B1,B2}(a.basis_l, a.basis_r, a.data-b.data)
-(a::SuperOperator) = SuperOperator(a.basis_l, a.basis_r, -a.data)
-(a::SuperOperator, b::SuperOperator) = throw(bases.IncompatibleBases())

"""
    spre(op)

Create a super-operator equivalent for right side operator multiplication.

For operators ``A``, ``B`` the relation

```math
    \\mathrm{spre}(A) B = A B
```

holds. `op` can be a dense or a sparse operator.
"""
spre(op::AbstractOperator) = SuperOperator((op.basis_l, op.basis_l), (op.basis_r, op.basis_r), tensor(op, identityoperator(op)).data)

"""
Create a super-operator equivalent for left side operator multiplication.

For operators ``A``, ``B`` the relation

```math
    \\mathrm{spost}(A) B = B A
```

holds. `op` can be a dense or a sparse operator.
"""
spost(op::AbstractOperator) = SuperOperator((op.basis_r, op.basis_r), (op.basis_l, op.basis_l), kron(permutedims(op.data), identityoperator(op).data))

function _check_input(H::AbstractOperator{B1,B2}, J::Vector, Jdagger::Vector, rates::Union{Vector{Float64}, Matrix{Float64}}) where {B1<:Basis,B2<:Basis}
    # TODO: replace basis checks by dispatch; make type of J known
    for j=J
        @assert isa(j, AbstractOperator{B1,B2})
    end
    for j=Jdagger
        @assert isa(j, AbstractOperator{B1,B2})
    end
    @assert length(J)==length(Jdagger)
    if isa(rates, Matrix{Float64})
        @assert size(rates, 1) == size(rates, 2) == length(J)
    elseif isa(rates, Vector{Float64})
        @assert length(rates) == length(J)
    end
end

"""
    liouvillian(H, J; rates, Jdagger)

Create a super-operator equivalent to the master equation so that ``\\dot ρ = S ρ``.

The super-operator ``S`` is defined by

```math
S ρ = -\\frac{i}{ħ} [H, ρ] + \\sum_i J_i ρ J_i^† - \\frac{1}{2} J_i^† J_i ρ - \\frac{1}{2} ρ J_i^† J_i
```

# Arguments
* `H`: Hamiltonian.
* `J`: Vector containing the jump operators.
* `rates`: Vector or matrix specifying the coefficients for the jump operators.
* `Jdagger`: Vector containing the hermitian conjugates of the jump operators. If they
             are not given they are calculated automatically.
"""
function liouvillian(H::AbstractOperator, J::Vector;
            rates::Union{Vector{Float64}, Matrix{Float64}}=ones(Float64, length(J)),
            Jdagger::Vector=dagger.(J))
    _check_input(H, J, Jdagger, rates)
    L = spre(-1im*H) + spost(1im*H)
    if isa(rates, Matrix{Float64})
        for i=1:length(J), j=1:length(J)
            jdagger_j = rates[i,j]/2*Jdagger[j]*J[i]
            L -= spre(jdagger_j) + spost(jdagger_j)
            L += spre(rates[i,j]*J[i]) * spost(Jdagger[j])
        end
    elseif isa(rates, Vector{Float64})
        for i=1:length(J)
            jdagger_j = rates[i]/2*Jdagger[i]*J[i]
            L -= spre(jdagger_j) + spost(jdagger_j)
            L += spre(rates[i]*J[i]) * spost(Jdagger[i])
        end
    end
    return L
end

"""
    exp(op::DenseSuperOperator)

Operator exponential which can for example used to calculate time evolutions.
"""
Base.exp(op::DenseSuperOperator) = DenseSuperOperator(op.basis_l, op.basis_r, exp(op.data))

# Array-like functions
Base.size(A::SuperOperator) = size(A.data)
@inline Base.axes(A::SuperOperator) = axes(A.data)
Base.ndims(A::SuperOperator) = 2
Base.ndims(::Type{<:SuperOperator}) = 2

# Broadcasting
Base.broadcastable(A::SuperOperator) = A

# Custom broadcasting styles
abstract type SuperOperatorStyle{BL<:Tuple{Basis,Basis},BR<:Tuple{Basis,Basis}} <: Broadcast.BroadcastStyle end
struct DenseSuperOperatorStyle{BL<:Tuple{Basis,Basis},BR<:Tuple{Basis,Basis}} <: SuperOperatorStyle{BL,BR} end
struct SparseSuperOperatorStyle{BL<:Tuple{Basis,Basis},BR<:Tuple{Basis,Basis}} <: SuperOperatorStyle{BL,BR} end

# Style precedence rules
Broadcast.BroadcastStyle(::Type{<:DenseSuperOperator{BL,BR}}) where {BL<:Tuple{Basis,Basis},BR<:Tuple{Basis,Basis}} = DenseSuperOperatorStyle{BL,BR}()
Broadcast.BroadcastStyle(::Type{<:SparseSuperOperator{BL,BR}}) where {BL<:Tuple{Basis,Basis},BR<:Tuple{Basis,Basis}} = SparseSuperOperatorStyle{BL,BR}()
Broadcast.BroadcastStyle(::DenseSuperOperatorStyle{B1,B2}, ::SparseSuperOperatorStyle{B1,B2}) where {B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis}} = DenseSuperOperatorStyle{B1,B2}()
Broadcast.BroadcastStyle(::DenseSuperOperatorStyle{B1,B2}, ::DenseSuperOperatorStyle{B3,B4}) where {B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis},B3<:Tuple{Basis,Basis},B4<:Tuple{Basis,Basis}} = throw(bases.IncompatibleBases())
Broadcast.BroadcastStyle(::SparseSuperOperatorStyle{B1,B2}, ::SparseSuperOperatorStyle{B3,B4}) where {B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis},B3<:Tuple{Basis,Basis},B4<:Tuple{Basis,Basis}} = throw(bases.IncompatibleBases())
Broadcast.BroadcastStyle(::DenseSuperOperatorStyle{B1,B2}, ::SparseSuperOperatorStyle{B3,B4}) where {B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis},B3<:Tuple{Basis,Basis},B4<:Tuple{Basis,Basis}} = throw(bases.IncompatibleBases())

# Out-of-place broadcasting
@inline function Base.copy(bc::Broadcast.Broadcasted{Style,Axes,F,Args}) where {BL<:Tuple{Basis,Basis},BR<:Tuple{Basis,Basis},Style<:DenseSuperOperatorStyle{BL,BR},Axes,F,Args<:Tuple}
    bcf = Broadcast.flatten(bc)
    bl,br = states.find_basis(bcf.args)
    bc_ = Broadcasted_restrict_f(bcf.f, bcf.args, axes(bcf))
    # TODO: remove convert
    return DenseSuperOperator{BL,BR}(bl, br, convert(Matrix{ComplexF64}, copy(bc_)))
end
@inline function Base.copy(bc::Broadcast.Broadcasted{Style,Axes,F,Args}) where {BL<:Tuple{Basis,Basis},BR<:Tuple{Basis,Basis},Style<:SparseSuperOperatorStyle{BL,BR},Axes,F,Args<:Tuple}
    bcf = Broadcast.flatten(bc)
    bl,br = states.find_basis(bcf.args)
    bc_ = Broadcasted_restrict_f(bcf.f, bcf.args, axes(bcf))
    return SuperOperator{BL,BR}(bl, br, copy(bc_))
end
states.find_basis(a::SuperOperator, rest) = (a.basis_l, a.basis_r)

const BasicMathFunc = Union{typeof(+),typeof(-),typeof(*)}
function Broadcasted_restrict_f(f::BasicMathFunc, args::Tuple{Vararg{<:SuperOperator}}, axes)
    args_ = Tuple(a.data for a=args)
    return Broadcast.Broadcasted(f, args_, axes)
end
function Broadcasted_restrict_f(f, args::Tuple{Vararg{<:SuperOperator}}, axes)
    throw(error("Cannot broadcast function `$f` on type `$(eltype(args))`"))
end

# In-place broadcasting
@inline function Base.copyto!(dest::SuperOperator{BL,BR}, bc::Broadcast.Broadcasted{Style,Axes,F,Args}) where {BL<:Tuple{Basis,Basis},BR<:Tuple{Basis,Basis},Style<:SuperOperatorStyle{BL,BR},Axes,F,Args}
    axes(dest) == axes(bc) || Base.Broadcast.throwdm(axes(dest), axes(bc))
    # Performance optimization: broadcast!(identity, dest, A) is equivalent to copyto!(dest, A) if indices match
    if bc.f === identity && isa(bc.args, Tuple{<:SuperOperator{BL,BR}}) # only a single input argument to broadcast!
        A = bc.args[1]
        if axes(dest) == axes(A)
            return copyto!(dest, A)
        end
    end
    # Get the underlying data fields of operators and broadcast them as arrays
    bcf = Broadcast.flatten(bc)
    bc_ = Broadcasted_restrict_f(bcf.f, bcf.args, axes(bcf))
    copyto!(dest.data, bc_)
    return dest
end
@inline Base.copyto!(A::SuperOperator{BL,BR},B::SuperOperator{BL,BR}) where {BL<:Tuple{Basis,Basis},BR<:Tuple{Basis,Basis}} = (copyto!(A.data,B.data); A)
@inline function Base.copyto!(dest::SuperOperator{B1,B2}, bc::Broadcast.Broadcasted{Style,Axes,F,Args}) where {
        B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis},B3<:Tuple{Basis,Basis},
        B4<:Tuple{Basis,Basis},Style<:SuperOperatorStyle{B3,B4},Axes,F,Args
        }
    throw(bases.IncompatibleBases())
end

end # module


text:
# TODO: Generic implementation of this
function rand(::Type{<:Normal}, pri::NormalInverseChisq, X)
    if length(X) > 0
        pri = posterior_canon(pri, suffstats(Normal, X))
    end
    μ, σ2 = rand(pri)
    Normal(μ, sqrt(σ2))
end

# abstract type ConjguateModel end

# struct NormalNIX <: ConjguateModel
#     distn::Normal
#     prior::NormalInverseChisq
# end

# struct MixtureModel{T} <: ConjguateModel
#     components::Vector{T}
#     weights::Vector{Float64}
#     prior::Dirichlet
# end

# # Rename resample to rand ?
# function resample(m::NormalNIX, X)
#     distn = rand(typeof(distn), m.prior, X)
#     NormalNIX(distn, prior)
# end

# function resample(m::MixtureModel{T}, X) where T
#     # TODO: Argcheck X is vector of vector
#     components = map(i -> resample(m.components[i], X[i]))
#     weights = rand()
# end

max_stars_repo_name: maxmouchet/Statistikk

max_issues_repo_name: maxmouchet/Statistikk

max_forks_repo_name: maxmouchet/Statistikk


text:
function get_mutated_instruction(instruction)
    choices = [1, 2, 4, 8, 16]
    choice = rand(choices) # determine which bit to flip
    return instruction ⊻ choice # XOR
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
# gantt chart

# VegaLite

using VegaLite

@vlplot(
    :bar,
    data={
        values=[
            {task="A",start=1,stop=3},
            {task="B",start=3,stop=8},
            {task="C",start=8,stop=10}
        ]
    },
    y="task:o",
    x="start:q",
    x2="stop:q"
)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
## ----------------------------------------------------------------------------
# utils
const _CLEAR_BIBTEX_ENTRY_STRIP_SET = ['{', '}', '\n', ' ', ',', ';']
const _alphanum_str = "a-zA-Z0-9"
const _punct_str = "\\,\\-\\;\\:\\."
const _reg1 = Regex("[$_punct_str]{2}")
const _reg2 = Regex("\\s{2}")
const _reg3 = Regex("[^$(_alphanum_str)$(_punct_str)/\\s]")

function _clear_bibtex_entry(str::AbstractString)

    str = strip(str, _CLEAR_BIBTEX_ENTRY_STRIP_SET)
    str = replace(str, "\n" => "")

    _till_unchanged() do
        str = replace(str, _reg1 => ",")
        str = replace(str, _reg2 => " ")
        for c in ",;:."
            str = replace(str, " $c" => c)
        end
        str = replace(str, _reg3 => "")
        return str
    end
    return str
end

## ----------------------------------------------------------------------------
function _till_unchanged(f!::Function; n = 100)
    ret = f!()
    ret0 = ret
    for i in 1:n
        ret = f!()
        (ret0 == ret) && return ret
    end
    return ret
end

## ----------------------------------------------------------------------------
function _push_csv!(col, str::AbstractString)
    vals = split(str, ",", keepempty = false)
    for vali in vals
        push!(col, strip(vali))
    end
    return col
end

function _push_csv!(col, str::AbstractString, strs::AbstractString...)
    _push_csv!(col, str)
    for stri in strs
        _push_csv!(col, stri)
    end
    return col
end

## ----------------------------------------------------------------------------
# zip
# function zip(srcfile, destfile = string(srcfile, ".zip"))
#     zdir = ZipFile.Writer(srcfile)
#     for (root, dirs, files) in walkdir(destfile)
#         for file in files
#             filepath = joinpath(root, file)
#             f = open(filepath, "r")
#             content = read(f, String)
#             close(f)
#             zf = ZipFile.addfile(zdir, basename(filepath));
#             write(zf, content)
#         end
#     end
#     close(zdir)
# end

function zip(srcfile, destfile = string(srcfile, ".tar.gz"))
    # p = run(`zip -r $(destfile) $(srcfile)`; wait = false)
    p = run(`tar -zcvf $(destfile) $(srcfile)`; wait = false)
    wait(p)
    return destfile
end

function unzip(srcfile, destfile = replace(srcfile, r".zip$" => ""))
    fileFullPath = isabspath(srcfile) ?  srcfile : joinpath(pwd(),srcfile)
    basePath = dirname(fileFullPath)
    outPath = (destfile == "" ? basePath : (isabspath(destfile) ? destfile : joinpath(pwd(),destfile)))
    isdir(outPath) ? "" : mkdir(outPath)
    zarchive = ZipFile.Reader(fileFullPath)
    for f in zarchive.files
        fullFilePath = joinpath(outPath,f.name)
        if (endswith(f.name,"/") || endswith(f.name,"\\"))
            mkdir(fullFilePath)
        else
            write(fullFilePath, read(f))
        end
    end
    close(zarchive)
end

## ----------------------------------------------------------------------------
function _dirsize(dir)
    size = 0
    for (root, dirs, files) in walkdir(dir)
        for file in files
            filepath = joinpath(root, file)
            size += filesize(filepath)
        end
    end
    return size
end

_filesize(path) = isdir(path) ? _dirsize(path) : filesize(path)

## ----------------------------------------------------------------------------
function _doi_to_url(doi::AbstractString)
    if startswith(doi, "http")
        return doi
    elseif startswith(doi, "doi.org")
        return string("https://", doi)
    elseif startswith(doi, r"\d{1,10}\.\d{3,20}")
        return string("https://doi.org/", doi)
    end
    return ""
end

function camel_alphanumeric(str)
    strs = split(str, r"[^a-zA-Z0-9]"; keepempty = false)
    join(uppercasefirst.(strs), " ")
end

function _if_nonempty(f::Function, str, dflt = "")
    isempty(str) ? dflt : f(str)
end
_if_nonempty(str, val, dflt = "") = _if_nonempty((x) -> val, str, dflt)

## ----------------------------------------------------------------------------
function _hash_file(file::AbstractString)
    h = hash("")
    !isfile(file) && return h

    open(file) do io
        for byte in read(io)
            h = hash(byte, h)
        end
    end
    return h
end


text:
using FluxGAN
using Test

@testset "FluxGAN.jl" begin
    # Write your tests here.
end

max_stars_repo_name: aarontrowbridge/GAN

max_issues_repo_name: aarontrowbridge/GAN

max_forks_repo_name: aarontrowbridge/GAN


text:
__precompile__()
module DigitSetSudoku

    using DigitSets

    type SudokuPuzzle
        squares::Array{Int8, 4}
    end

    SudokuPuzzle(spec::AbstractString) = parse(SudokuPuzzle, spec)

    function Base.parse(::Type{SudokuPuzzle}, spec::AbstractString)
        squares = Int8[]
        sizehint!(squares, 81)

        for c in spec
            if '0' <= c <= '9'
                push!(squares, parse(Int, c))
            elseif c == '.'
                push!(squares, 0)
            end
        end

        if length(squares) != 81
            error(@sprintf(
                "Expected 81 characters from [\\.0-9] but saw %d.",
                length(squares)
            ))
        end

        SudokuPuzzle(reshape(squares, (3, 3, 3, 3)))
    end

    function Base.show(io::IO, puzzle::SudokuPuzzle)
        for band in 1:3
            for row in 1:3
                for stack in 1:3
                    for col in 1:3
                        digit = puzzle.squares[col, stack, row, band]
                        print(io, digit == 0 ? "." : digit)
                        (stack == 3 && col == 3) || print(io, " ")
                    end
                    stack < 3 && print(io, "| ")
                end
                println(io, "")
            end
            band < 3 && println(io, "------+-------+------")
        end
    end

    type SudokuBoard
        squares::Array{DigitSet, 4}
    end

    SudokuBoard() = SudokuBoard(
        [
            DigitSet(1:9)
            for col in 1:3, stack in 1:3, row in 1:3, band in 1:3
        ]
    )

    const squareindices = reshape(1:81, 3, 3, 3, 3)
    const units = [
        (
            collect(squareindices[1:3, stack, 1:3, band]),
            collect(squareindices[1:3, 1:3, row, band]),
            collect(squareindices[col, stack, 1:3, 1:3])
        )
        for col in 1:3, stack in 1:3, row in 1:3, band in 1:3
    ]

    function solve(puzzle::SudokuPuzzle)
        board = SudokuBoard()
        assigngivens!(board, puzzle) || error("Inconsistent board")
        search!(board) || error("Inconsistent board")
        board
    end

    function assigngivens!(board::SudokuBoard, puzzle::SudokuPuzzle)
        for i in eachindex(puzzle.squares)
            digit = puzzle.squares[i]
            digit == 0 && continue
            assign!(board, DigitSet(digit), i) || return false
        end
        true
    end

    SudokuBoard(puzzle::SudokuPuzzle) = solve(puzzle)

    function assign!(board::SudokuBoard, ds::DigitSet, i)
        ds == board.squares[i] && return true
        length(ds) == 0 && return false

        board.squares[i] = ds

        # Strategies
        nakedsingle!(board, ds, i) || return false
        hiddensingle!(board, i) || return false

        true
    end

    function nakedsingle!(board::SudokuBoard, ds::DigitSet, i)
        length(ds) == 1 || return true
        for u in units[i]
            nakedsingleunit!(board, ds, i, u) || return false
        end
        true
    end

    function nakedsingleunit!(board::SudokuBoard, ds::DigitSet, i, unit)
        for j in unit
            i == j && continue
            difference = setdiff(board.squares[j]::DigitSet, ds)
            assign!(board, difference, j) || return false
        end
        true
    end

    function hiddensingle!(board::SudokuBoard, i)
        for u in units[i]
            hiddensingleunit!(board, u) || return false
        end
        true
    end

    function hiddensingleunit!(board::SudokuBoard, unit)
        singles = lonelydigits(board, unit)
        length(singles) > 0 || return true
        for i in unit
            intersection = intersect(singles, board.squares[i])
            # A board is inconsistent if we ever find a square that contains
            # more than one lonely digit for the same unit
            length(intersection) <= 1 || return false
            length(intersection) == 1 || continue
            assign!(board, intersection, i) || return false
        end
        true
    end

    function lonelydigits(board::SudokuBoard, unit)
        # Digits that occur only once in a unit
        singles = DigitSet()
        # Digits that occur multiple times
        multiples = DigitSet()
        for i in unit
            square = board.squares[i]
            # digits that occurred once before and also occur in square
            # have now occurred multiple times
            multiples = union(multiples, intersect(square, singles))
            # Add digits in square to singles and remove known multiples
            singles = setdiff(union(singles, square), multiples)
        end
        singles
    end

    function search!(board::SudokuBoard)
        i = searchindex(board)
        # Board is solved. Could adjust this to push board onto
        # an accumulator if we wanted to find all solutions.
        i == 0 && return true
        ds = board.squares[i]
        for digit in ds
            tmp = copy(board.squares)
            status = assign!(board, DigitSet(digit), i)
            status = status && search!(board)
            status && return true
            board.squares = tmp
        end
        false
    end

    function searchindex(board::SudokuBoard)
        out = 0
        minlen = 10
        for i in eachindex(board.squares)
            len = length(board.squares[i])
            if 1 < len < minlen
                minlen = len
                out = i
            end
        end
        out
    end

    function Base.show(io::IO, board::SudokuBoard)
        width = 1 + maximum(map(length, board.squares))
        dashes = repeat("-", 3*width)
        line = string(dashes, "+", dashes, "-", "+", dashes)
        for band in 1:3
            for row in 1:3
                for stack in 1:3
                    for col in 1:3
                        square = board.squares[col, stack, row, band]
                        for digit in square
                            print(io, digit)
                        end
                        print(io, repeat(" ", width - length(square) - 1))
                        (stack == 3 && col == 3) || print(io, " ")
                    end
                    stack < 3 && print(io, "| ")
                end
                println(io, "")
            end
            band < 3 && println(io, line)
        end
    end

    export SudokuBoard, SudokuPuzzle, solve

end # module


text:
for N in [Float64, Float32, Rational{Int}]
    c1 = N[1, 2]
    c2 = N[3, 4, 5]
    Z1 = Zonotope(c1, N[1 0; 0 1])
    Z2 = Zonotope(c2, N[1 0 0; 0 1 0; 0 0 1])
    H1 = Hyperrectangle(c1, N[1, 1])
    H2 = Hyperrectangle(c2, N[1, 1, 1])
    B1 = BallInf(c1, N(1))
    B2 = BallInf(c2, N(1))

    Z1Z2 = cartesian_product(Z1, Z2)
    Z1H2 = cartesian_product(Z1, H2)
    H1Z2 = cartesian_product(H1, Z2)
    H1H2 = cartesian_product(H1, H2)
    H1B2 = cartesian_product(H1, B2)
    B1B2 = cartesian_product(B1, B2)
    for X in [Z1Z2, Z1H2, H1Z2]
        @test X isa Zonotope{N}
    end
    for X in [H1H2, H1B2, B1B2]
        @test X isa Hyperrectangle{N}
    end
    if test_suite_polyhedra || N <: AbstractFloat
        for X in [Z1Z2, Z1H2, H1Z2, H1H2, H1B2]
            @test isequivalent(X, B1B2)
        end
    end

    # Universe is handled correctly (#2735)
    P = HalfSpace(N[1, 2], N(3))
    U = Universe{N}(2)
    @test isequivalent(cartesian_product(P, U), HalfSpace(N[1, 2, 0, 0], N(3)))
    @test isequivalent(cartesian_product(U, P), HalfSpace(N[0, 0, 1, 2], N(3)))
    @test cartesian_product(U, U) == Universe{N}(4)
end


text:
using JLD

include("network_library.jl")

function n_rand_in_range(A::Array,N::Int64)
    # generate an array of length N within the range specified by A
    # length 1 means parameter is fixed, length 2 means random within bounds
    @assert length(A) == 1 || length(A) == 2
    eltype = typeof(A).parameters[1]
    @assert eltype == Float64 || eltype == Int64

    if length(A) == 1
        # length 1 indicates no range, meaning all are identical
        ary = ones(N)*A[1]

        if eltype == Int64
            return convert(Vector{Int64},ary)
        else
            return ary
        end
    else
       # length 2 means we have the bounds
       if eltype == Int64
           return A[1] + rand(A[1]:A[2],N)
       else
           A_range = diff(A)[1]
           return A[1] + rand(N)*A_range
       end
   end
end

function get_exit_node(g::MetaGraphs.MetaGraph)
    e_nodes = g.gprops[:exit_nodes]
    if length(e_nodes) > 1
        error("something is wrong, should only be length 1 in this experiment")
    else
        return e_nodes[1]
    end
end

function populate_net_dict(i,training_set_size,n,exit_rwd,caught_rwd,sensor_rwd,seed_list,degree,mcts_its,mcts_depth,mcts_e,discount_fact,transition_prob,fname,net_type)
    problem_dict = Dict()
        #  try
            println("Making network $i of $training_set_size")

            if net_type == :random
                mthd = rand([:watts_strogatz,:expected_degree,:erdos_n_e,:static_scale_free],1)[1]
                g = rand_network(n,exit_rwd=exit_rwd,caught_rwd=caught_rwd,
                                 sensor_rwd=sensor_rwd,net_seed=seed_list,approx_E=2*n,
                                 exit_nodes=[8],target_mean_degree=degree,method=mthd)
            elseif net_type == :original
                mthd = :original
                g = original_roadnet(exit_rwd=exit_rwd,caught_rwd=caught_rwd,sensor_rwd=sensor_rwd)
            else
                println("no other net_types implemented at this time")
            end

            evader_start = 1
            pursuer_start = 4
            exit_loc = get_exit_node(g)
            #  display_network(g,evader_locs=[evader_start],pursuer_locs=[pursuer_start],fname="logs/net$i")

            problem_dict[i] = Dict(:graph=>g,:mcts_its=>mcts_its,:mcts_depth=>mcts_depth,
                                   :mcts_e=>mcts_e,:net_seed=>seed_list,:num_nodes=>n,
                                   :exit_rwd=>exit_rwd,:caught_rwd=>caught_rwd,
                                   :discount=>discount_fact,:trans_prob=>transition_prob,
                                   :sensor_rwd=>sensor_rwd,:target_degree=>degree,
                                   :evader_start=>evader_start,:pursuer_start=>pursuer_start,
                                   :exit_loc=>exit_loc,:net_gen_method=>mthd)
    return problem_dict
end

function make_nets(net::Dict)
    net_dict = net[:nets]
    net_type = net_dict[:net_type]
    discount = net_dict[:d_bnds]
    t_bnds = net_dict[:t_bnds]
    mcts_e = net_dict[:e_bnds]
    mcts_d = net_dict[:solver_depth]
    mcts_its_bnds = net_dict[:mcts_its_bnds]
    n_bnds = net_dict[:n_bnds]
    deg_bnds = net_dict[:deg_bnds]
    caught_bnds = net_dict[:caught_rwd_bnds]
    sense_bnds = net_dict[:sense_rwd_bnds]
    ext_bnds = net_dict[:exit_rwd_bnds]
    fldr = net[:fldr]
    fname = net[:fname]
    num_nets = net_dict[:num_nets]
    seed = net_dict[:seed]

    fn = joinpath(fldr,fname*".jld")

    make_nets(num_nets,fname=fn,exit_rwd_bounds=ext_bnds,
              sensor_rwd_bounds=sense_bnds,caught_rwd_bounds=caught_bnds,degree_bounds=deg_bnds,
              n_bounds=n_bnds,mcts_its_bounds=mcts_its_bnds,mcts_depth_bounds=mcts_d,
              mcts_e_bounds=mcts_e,trans_prob_bounds=t_bnds,
              discount_fact_bounds=discount,net_type=net_type,random_seed=seed)
end
function make_nets(training_set_size::Int64;exit_rwd_bounds::Array{Float64}=[2000.],
                   sensor_rwd_bounds::Array{Float64}=[-100.],caught_rwd_bounds::Array{Float64}=[-2000.],
                   degree_bounds::Array{Float64}=[3.0,7.0],n_bounds::Array{Int64}=[15,75],
                   mcts_its_bounds::Array{Int64}=[100],mcts_depth_bounds::Array{Int64}=[5],
                   mcts_e_bounds::Array{Float64}=[1.],trans_prob_bounds::Array{Float64}=[1.],
                   discount_fact_bounds::Array{Float64}=[0.9],net_type::Symbol=:original,
                   random_seed::Int64=1,fname::String="logs/generated_training_nets.jld")

    srand(random_seed)
    seed_list = collect(1:training_set_size)

    # POMDP Problem Properties
    exit_rwd_ary = n_rand_in_range(exit_rwd_bounds,training_set_size)
    sensor_rwd_ary = n_rand_in_range(sensor_rwd_bounds,training_set_size)
    caught_rwd_ary = n_rand_in_range(caught_rwd_bounds,training_set_size)
    transition_prob_ary = n_rand_in_range(trans_prob_bounds,training_set_size)
    discount_fact_ary = n_rand_in_range(discount_fact_bounds,training_set_size)

    # Solver Properties
    mcts_its_ary = n_rand_in_range(mcts_its_bounds,training_set_size)
    mcts_depth_ary = n_rand_in_range(mcts_depth_bounds,training_set_size)
    mcts_e_ary = n_rand_in_range(mcts_e_bounds,training_set_size)

    # Network Properties
    degree_ary = n_rand_in_range(degree_bounds,training_set_size)
    n_ary = n_rand_in_range(n_bounds,training_set_size)

    @save "$fname" training_set_size seed_list exit_rwd_bounds exit_rwd_ary sensor_rwd_bounds sensor_rwd_ary caught_rwd_bounds caught_rwd_ary degree_bounds degree_ary n_bounds n_ary

    println(exit_rwd_ary,sensor_rwd_ary,caught_rwd_ary,degree_ary,n_ary)

    # store problems in a dictionary
    pd_array = pmap(populate_net_dict,collect(1:training_set_size),repmat([training_set_size],training_set_size,1),n_ary,exit_rwd_ary,caught_rwd_ary,sensor_rwd_ary,seed_list,degree_ary,mcts_its_ary,mcts_depth_ary,mcts_e_ary,discount_fact_ary,transition_prob_ary,repmat([fname],training_set_size,1),repmat([net_type],training_set_size,1))
    # pmap gives an array, we want a Dict, so put the entries into a dict
    # ther must be a better way to do this, but I don't have time to figure it out
    problem_dict = Dict()
    for i=1:length(pd_array)
        #  display(pd_array[i][i])
        # HACK, HACKETY HACK HACK HACK
        problem_dict[i] = pd_array[i][i]
    end
    #  display(problem_dict)
    #  display(problem_dict[1])
    #  display(problem_dict[2])
    #  error()

    println("writing data to $fname")
    jldopen("$fname", "w") do file
        addrequire(file, MetaGraphs)
        write(file,"problem_dict", problem_dict)
    end
end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
include("agent.jl")
include("dyna_agent.jl")


text:
module Models
#computations
using Statistics, LinearAlgebra, SpecialFunctions, Distributions, Printf, Distributed, DataStructures
import ..Utils: log_message, log_progress, PROJECT_ROOT
import ..Data: Experiment, CellLine, Drug, Gene, Protein, RNASeqCall, ExomeSeq, RPPA, Pathway
import ..Data: DataView, KeyType, ViewType, ValueViewType, NAViewType
import ..Data: get_experiment, get_measurement_value, normalize_data_views, filter_data_views, unfilter_data_views
using ..Utils, ..Data
include("distributional_parameters.jl")
include("types.jl")
include("utils.jl")
include("bmtmkl/bmtmkl.jl")
# <======== additional models need to be included here
include("run_models.jl")
include("io.jl")
include("pipeline.jl")

end


text:
# utility functions and macros, which are independent of `FixedPoint`
bitwidth(T::Type) = 8sizeof(T)

widen1(T::Type)         = T # fallback
widen1(::Type{Int8})    = Int16
widen1(::Type{UInt8})   = UInt16
widen1(::Type{Int16})   = Int32
widen1(::Type{UInt16})  = UInt32
widen1(::Type{Int32})   = Int64
widen1(::Type{UInt32})  = UInt64
widen1(::Type{Int64})   = Int128
widen1(::Type{UInt64})  = UInt128
widen1(x::Integer) = x % widen1(typeof(x))

signedtype(::Type{T}) where {T <: Integer} = typeof(signed(zero(T)))

const ShortInts = Union{Int8, UInt8, Int16, UInt16}
const LongInts = Union{Int64, UInt64, Int128, UInt128, BigInt}

const ShorterThanInt = Int === Int32 ? ShortInts : Union{ShortInts, Int32, UInt32}
const NotBiggerThanInt = Union{ShorterThanInt, Int, UInt}
const NotBiggerThanInt64 = Union{ShortInts, Int32, UInt32, Int64, UInt64}
const SShorterThanInt = typeintersect(ShorterThanInt, Signed)
const UShorterThanInt = typeintersect(ShorterThanInt, Unsigned)

macro f32(x::Float64) # just for hexadecimal floating-point literals
    :(Float32($x))
end
macro exp2(n)
     :(_exp2(Val($(esc(n)))))
end
_exp2(::Val{N}) where {N} = exp2(N)

# these are defined in julia/float.jl or julia/math.jl, but not exported
significand_bits(::Type{Float32}) = 23
significand_bits(::Type{Float64}) = 52
exponent_bias(::Type{Float32}) = 127
exponent_bias(::Type{Float64}) = 1023

_unsafe_trunc(::Type{T}, x::Integer) where {T} = x % T
_unsafe_trunc(::Type{T}, x) where {T} = unsafe_trunc(T, x)
# issue #202, #211
_unsafe_trunc(::Type{T}, x::BigFloat) where {T <: Integer} = trunc(BigInt, x) % T

if !signbit(signed(unsafe_trunc(UInt, -12.345)))
    # a workaround for ARM (issue #134)
    function _unsafe_trunc(::Type{T}, x::AbstractFloat) where {T <: Integer}
        if T === UInt32
            copysign(unsafe_trunc(T, abs(x)), x)
        else
            unsafe_trunc(T, unsafe_trunc(signedtype(T), x))
        end
    end
end

wrapper(@nospecialize(T)) = Base.typename(T).wrapper


text:
using WriteVTK
using Test: @test
using SHA: sha1

include("extent.jl")

const tests = [
    "polyhedron_cube.jl",
    "multiblock.jl",
    "rectilinear.jl",
    "imagedata.jl",
    "structured.jl",
    "unstructured.jl",
    "polydata.jl",
    "bezier.jl",
    "pvdCollection.jl",
    "array.jl",
    "pvtk_grid.jl",
]

# Only toggle to generate new checksums, if new tests are added.
const OVERWRITE_CHECKSUMS = false
const checksums_file = joinpath(dirname(@__FILE__), "checksums.sha1")
const checksum_list = read(checksums_file, String)

if OVERWRITE_CHECKSUMS
    csio = open(checksums_file, "w")
end

const EXECDIR = "output"
mkpath(EXECDIR)
cd(EXECDIR)

# Run the test scripts.
for test in tests
    println("TEST (first run): ", test)
    outfiles = evalfile(test)::Vector{String}

    # Check that the generated files match the stored checksums.
    for file in outfiles
        sha_str = bytes2hex(open(sha1, file)) * "  $file\n"
        if OVERWRITE_CHECKSUMS
            write(csio, sha_str)
        else
            # Returns `nothing` if string is not found.
            @info "Verifying $file"
            @test findfirst(sha_str, checksum_list) !== nothing
        end
    end
    println()
end

OVERWRITE_CHECKSUMS && close(csio)

println("="^60, "\n")

# Run the tests again, just to measure the time and allocations once all the
# functions have already been compiled.
for test in tests
    println("TEST (second run): ", test)
    outfiles = evalfile(test)::Vector{String}
    println()
end


text:
function next_permutation{Container}(arr::Container)

	len=length(arr)
	i=len
	while i>1 && arr[i-1]>=arr[i]
		i-=1
	end
	if i<=1
		return false
	end
	begin
		j=len
		while arr[j]<=arr[i-1]
			j-=1
		end
		arr[i-1],arr[j]=arr[j],arr[i-1]
	end
	begin
	    j=len
	    while i<j
	    	arr[i],arr[j]=arr[j],arr[i]
	    	i+=1
	    	j-=1
	    end
	end

	    return true

	end

"""
# test
   a=collect(1:5)
   while next_permutation(a)
       println(a)
   end

"""