text:
abstract type Data end

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(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},
        conditional::Bool,
        condition1_present::Array{Union{Missing, T}, 2} where T <: Number,
        condition1_future::Array{Union{Missing, T}, 2} where T <: Number,
        condition2_present::Array{Union{Missing, T}, 2} where T <: Number,
        condition2_future::Array{Union{Missing, T}, 2} where T <: Number,
        comparison1::String,
        comparison2::String,
        condition1_lower::Number,
        condition1_upper::Number,
        condition2_lower::Number,
        condition2_upper::Number,
        precision::DataType;
        x::Int64,
        y::Int64,
        strength::Number
    )

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

    source_subset = clip(source_array,
                         x = x,
                         y = y,
                         distance = radius)

    # Append missing if buffer > 0
    if buffer > 0
        ### Columns
        nrow_sub = size(source_subset)[1]
        left_col_num = max(0, min(buffer, x - radius - 1))
        right_col_num = max(0, min(buffer, ncols - (x + 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, y - radius - 1))
        bottom_row_num = max(0, min(buffer, nrows - (y + 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)] * strength) / source_sum

    if conditional
        xlower_buffered = Int64(max(x - radius - buffer, 1))
        xupper_buffered = Int64(min(x + radius + buffer, ncols))
        ylower_buffered = Int64(max(y - radius - buffer, 1))
        yupper_buffered = Int64(min(y + radius + buffer, nrows))
        xlower = x - block_radius
        xupper = min(x + block_radius, ncols)
        ylower = y - block_radius
        yupper = min(y + block_radius, nrows)

        source_target_match!(source_subset,
                             arguments["n_conditions"],
                             condition1_present,
                             condition1_future,
                             condition2_present,
                             condition2_future,
                             comparison1,
                             comparison2,
                             condition1_lower,
                             condition1_upper,
                             condition2_lower,
                             condition2_upper,
                             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,
                              condition1_present::Array{Union{T, Missing}, 2} where T <: Number,
                              condition1_future::Array{Union{T, Missing}, 2} where T <: Number,
                              condition2_present::Array{Union{T, Missing}, 2} where T <: Number,
                              condition2_future::Array{Union{T, Missing}, 2} where T <: Number,
                              comparison1::String,
                              comparison2::String,
                              condition1_lower::Number,
                              condition1_upper::Number,
                              condition2_lower::Number,
                              condition2_upper::Number,
                              ylower::Int64,
                              yupper::Int64,
                              xlower::Int64,
                              xupper::Int64,
                              ylower_buffered::Int64,
                              yupper_buffered::Int64,
                              xlower_buffered::Int64,
                              xupper_buffered::Int64
                              )
    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;
        x::Int64,
        y::Int64,

    )
    radius = arguments["radius"]
    buffer = arguments["buffer"]
    distance = radius + buffer

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

    xlower = Int64(max(x - radius - buffer, 1))
    xupper = Int64(min(x + radius + buffer, ncols))
    ylower = Int64(max(y - radius - buffer, 1))
    yupper = Int64(min(y + 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, x)
    new_y = min(distance + 1, y)

    ground[new_y, new_x] = Inf

    ground
end

function get_conductance(
        resistance::Array{Union{T, Missing}, 2} where T <: Number,
        arguments::Dict{String, Int64},
        x::Int64,
        y::Int64,
        resistance_file_is_conductance::Bool
    )

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

    resistance_clipped = clip(resistance,
                              x = x,
                              y = y,
                              distance = radius + buffer)

    if resistance_file_is_conductance
        conductance = resistance_clipped
    else
        conductance = 1 ./ resistance_clipped
    end

    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,
        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.RasData(cellmap,
                                      polymap,
                                      source_map,
                                      ground_map,
                                      points_rc,
                                      strengths,
                                      included_pairs,
                                      hbmeta)

    # Generate advanced data
    data = Circuitscape.compute_advanced_data(rasterdata, 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
    src = data.src
    check_node = data.check_node
    source_map = data.source_map # Need it for one to all mode
    cellmap = data.cellmap

    # Flags
    is_raster = flags.is_raster
    is_alltoone = flags.is_alltoone
    is_onetoall = flags.is_onetoall
    write_v_maps = flags.outputflags.write_volt_maps
    write_c_maps = flags.outputflags.write_cur_maps
    write_cum_cur_map_only = flags.outputflags.write_cum_cur_map_only

    volt = zeros(eltype(G), size(nodemap))
    ind = findall(x -> x != 0, nodemap)
    f_local = Vector{eltype(G)}()
    solver_called = false
    voltages = Vector{eltype(G)}()
    outvolt = Circuitscape.alloc_map(hbmeta)
    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,
                                                a_local,
                                                s_local,
                                                g_local,
                                                f_local)

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

        solver_called = true

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

    outcurr
end

function solve_target!(
        i::Int64,
        n_targets::Int64,
        int_arguments::Dict{String, Int64},
        targets::Array{T, 2} where T <: Number,
        source_strength::Array{Union{Missing, T}, 2} where T <: Number,
        resistance::Array{Union{Missing, T}, 2} where T <: Number,
        cs_cfg::Dict{String, String},
        flags::Circuitscape.RasterFlags,
        o::Circuitscape.OutputFlags,
        calc_flow_potential::Bool,
        correct_artifacts::Bool,
        conditional::Bool,
        condition1_present::Array{Union{Missing, T}, 2} where T <: Number,
        condition1_future::Array{Union{Missing, T}, 2} where T <: Number,
        condition2_present::Array{Union{Missing, T}, 2} where T <: Number,
        condition2_future::Array{Union{Missing, T}, 2} where T <: Number,
        comparison1::String,
        comparison2::String,
        condition1_lower::Number,
        condition1_upper::Number,
        condition2_lower::Number,
        condition2_upper::Number,
        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,
        resistance_file_is_conductance::Bool
    )
    ## get source
    x_coord = Int64(targets[i, 1])
    y_coord = Int64(targets[i, 2])
    source = get_source(source_strength,
                        int_arguments,
                        conditional,
                        condition1_present,
                        condition1_future,
                        condition2_present,
                        condition2_future,
                        comparison1,
                        comparison2,
                        condition1_lower,
                        condition1_upper,
                        condition2_lower,
                        condition2_upper,
                        precision,
                        x = x_coord,
                        y = y_coord,
                        strength = float(targets[i, 3]))

    ## get ground
    ground = get_ground(int_arguments,
                        precision,
                        x = x_coord,
                        y = y_coord)

    ## get conductances for Omniscape
    conductance = get_conductance(resistance,
                                  int_arguments,
                                  x_coord,
                                  y_coord,
                                  resistance_file_is_conductance)

    grid_size = size(source)

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

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

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

    if correct_artifacts
        correction_array2 = deepcopy(correction_array)
        lowerxcut = 1
        upperxcut = size(correction_array, 2)
        lowerycut = 1
        upperycut = size(correction_array, 1)

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

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

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

        if 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 calc_flow_potential
            flow_potential = flow_potential .* correction_array2
        end
    end

    ## Accumulate values
    xlower = max(x_coord - int_arguments["radius"] - int_arguments["buffer"], 1)
    xupper = min(x_coord + int_arguments["radius"] + int_arguments["buffer"],
                 int_arguments["ncols"])
    ylower = max(y_coord - int_arguments["radius"] - int_arguments["buffer"], 1)
    yupper = min(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 calc_flow_potential == true
        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},
        cs_cfg::Dict{String, String},
        flags::Circuitscape.RasterFlags,
        o,
        conditional::Bool,
        condition1_present::Array{Union{T, Missing}, 2} where T <: Number,
        condition1_future::Array{Union{T, Missing}, 2} where T <: Number,
        condition2_present::Array{Union{T, Missing}, 2}where T <: Number,
        condition2_future::Array{Union{T, Missing}, 2} where T <: Number,
        comparison1::String,
        comparison2::String,
        condition1_lower::Number,
        condition1_upper::Number,
        condition2_lower::Number,
        condition2_upper::Number,
        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)

    source_blocked = get_source(temp_source,
                              arguments,
                              conditional,
                              condition1_present,
                              condition1_future,
                              condition2_present,
                              condition2_future,
                              comparison1,
                              comparison2,
                              condition1_lower,
                              condition1_upper,
                              condition2_lower,
                              condition2_upper,
                              precision,
                              x = (arguments["radius"] + arguments["buffer"] + 1),
                              y = (arguments["radius"] + arguments["buffer"] + 1),
                              strength = float(arguments["block_size"] ^ 2))

    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,
                                           flags,
                                           cs_cfg,
                                           precision)

    null_current =  calculate_current(conductance,
                                      source_null,
                                      ground,
                                      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


text:
@noinline ThrowBoundsError(A, i) = (println("A of length $(length(A))."); throw(BoundsError(A, i)))

Base.IndexStyle(::Type{<:AbstractStrideArray}) = IndexCartesian()
Base.IndexStyle(::Type{<:AbstractStrideVector{<:Any,<:Any,<:Any,1}}) = IndexLinear()
@generated function Base.IndexStyle(::Type{<:AbstractStrideArray{S,D,T,N,C,B,R}}) where {S,D,T,N,C,B,R}
    # if is column major || is a transposed contiguous vector
    if all(D) && ((isone(C) && R === ntuple(identity, Val(N))) || (C === 2 && R === (2,1) && S <: Tuple{One,Integer}))
        :(IndexLinear())
    else
        :(IndexCartesian())
    end
end

Base.@propagate_inbounds Base.getindex(A::AbstractStrideVector, i::Int, j::Int) = A[i]
@inline function Base.getindex(A::PtrArray{S,D,T,K}, i::Vararg{Integer,K}) where {S,D,T,K}
    @boundscheck checkbounds(A, i...)
    vload(stridedpointer(A), i)
end
@inline function Base.getindex(A::AbstractStrideArray{S,D,T,K}, i::Vararg{Integer,K}) where {S,D,T,K}
    b = preserve_buffer(A)
    P = PtrArray(A)
    GC.@preserve b begin
        @boundscheck checkbounds(P, i...)
        vload(stridedpointer(P), i)
    end
end
@inline function Base.setindex!(A::PtrArray{S,D,T,K}, v, i::Vararg{Integer,K}) where {S,D,T,K}
    @boundscheck checkbounds(A, i...)
    vstore!(stridedpointer(A), v, i)
    v
end
@inline function Base.setindex!(A::AbstractStrideArray{S,D,T,K}, v, i::Vararg{Integer,K}) where {S,D,T,K}
    b = preserve_buffer(A)
    P = PtrArray(A)
    GC.@preserve b begin
        @boundscheck checkbounds(P, i...)
        vstore!(stridedpointer(P), v, i)
    end
    v
end
@inline function Base.getindex(A::PtrArray, i::Integer)
    @boundscheck checkbounds(A, i)
    vload(stridedpointer(A), (i - one(i),))
end
@inline function Base.getindex(A::AbstractStrideArray, i::Integer)
    b = preserve_buffer(A)
    P = PtrArray(A)
    GC.@preserve b begin
        @boundscheck checkbounds(P, i)
        vload(stridedpointer(P), (i - one(i),))
    end
end
@inline function Base.setindex!(A::PtrArray, v, i::Integer)
    @boundscheck checkbounds(A, i)
    vstore!(stridedpointer(A), v, (i - one(i),))
    v
end
@inline function Base.setindex!(A::AbstractStrideArray, v, i::Integer)
    b = preserve_buffer(A)
    P = PtrArray(A)
    GC.@preserve b begin
        @boundscheck checkbounds(P, i)
        vstore!(stridedpointer(P), v, (i - one(i),))
    end
    v
end

ext: jl


text:
module Lib

export PersistentOID, ObjectPool

# Wrapper for libpmemobj
using Libdl

# constants
const _DEPS = joinpath(@__DIR__, "..", "deps")
const _LIB = joinpath(_DEPS, "usr", "lib")
const _INCLUDE = joinpath(_DEPS, "usr", "include")

const libpmem = joinpath(_LIB, "libpmem.so")
const libpmemobj = joinpath(_LIB, "libpmemobj.so")

# dependencies
function __init__()
    global libpmem
    global libpmemobj
    Libdl.dlopen(libpmem, Libdl.RTLD_GLOBAL)
    Libdl.dlopen(libpmemobj, Libdl.RTLD_GLOBAL)
end

_typenum(::Type{T}) where {T} = hash(T)

# Interface into the low level C calls.
#
# Implemented Functions
#
# pmemobj_open
# pmemobj_close
# pmemobj_create
# pmemobj_check
#
# pmemobj_root
#
# pmemobj_tx_begin
# pmemobj_tx_end
# pmemobj_tx_stage
# pmemobj_tx_commit
# pmemobj_tx_abort
#
# pmemobj_tx_alloc
#
# List of calls that probably have to get made at some point:
#
# pmemobj_tx_zalloc
# pmemobj_type_num
# pmemobj_tx_realloc
# pmemobj_tx_zrealloc
# pmemobj_tx_free
# pmemobj_direct
# pmemobj_tx_add_range_direct
# pmemobj_first
# pmemobj_next
#
#
### Constants
#
# PMEMOBJ_MIN_POOL
# PMEMOBJ_MAX_ALLOC_SIZE

function check_version(major, minor)
    result = ccall((:pmemobj_check_version, libpmemobj), Cstring, (Cint, Cint), major, minor)
    #result = @cxx pmemobj_check_version(major, minor)
    # Returned result is null if we're good to go.
    # Otherwise, it's a Cstring and we have an error
    if result != C_NULL
        error(unsafe_string(result))
    end
    return nothing
end

#####
##### Pool creation and handling
#####

geterrno() = unsafe_load(cglobal((:errno, libpmemobj), Int32))

# Wrapper for
struct ObjectPool
    ptr::Ptr{Nothing}
end
isnull(pool::ObjectPool) = pool.ptr == C_NULL

@enum TX_Param::Int32 begin
    TX_PARAM_NONE
    #TX_PARAM_MUTEX
    #TX_PARAM_RWLOCK
    #TX_PARAM_CB
end

@enum TX_Stage::Int32 begin
    TX_STAGE_NONE
    TX_STAGE_WORK
    TX_STAGE_ONCOMMIT
    TX_STAGE_ONABORT
    TX_STAGE_FINALLY
end

"""
    Persistent Object ID

Implementation Notes
--------------------
Matches the layout of

```c
typedef struct pmemoid {
    uint64_t pool_uuid_lo;
    uint64_t off;
} PMEMoid;
```
in `pmdk/src/include/libpmemobj/base.h`, line 108-111. Since this type is composed of
bits-type objects, it can be passed back and forth by value to the C code.
"""
struct PersistentOID{T}
    pool_uuid_lo::UInt64
    off::UInt64
end
isnull(oid::PersistentOID) = iszero(oid.pool_uuid_lo) && iszero(oid.off)

#####
##### pmemobj_open
#####

# These are are funcions defined in:
# http://pmem.io/pmdk/manpages/linux/v1.5/libpmemobj/pmemobj_open.3

"""
    close(pool::ObjectPool)

Close the memory pool indicated by `pool` and delete the memory pool handle. The object
store itself lives on in the file that contains it and may be reopened at a later time using
[`open`](@ref).
"""
function close(pool::ObjectPool)
    ccall(
        (:pmemobj_close, libpmemobj),
        Cvoid,
        (Ptr{Nothing},),
        pool.ptr
    )
end

"""
    create(path, poolsize; layout = "", mode = UInt32(0o666) -> ObjectPool

Create a transactional object store with the given total `poolsize` in bytes. Argument
`path` specifies the name of the memory pool file to be created. Argument `layout` specifies
the application's layout type in the form of a string. The layout name is not interpreted by
`libpmemobj`, but may be used as a check when [`pmemobj_open`](@ref`) is called.
"""
function create(path, poolsize; layout = "", mode = UInt32(0o0666))
    ptr = ccall(
        (:pmemobj_create, libpmemobj),
        Ptr{Nothing},
        (Cstring, Cstring, Csize_t, Base.Cmode_t),
        path, layout, poolsize, mode
    )
    return ObjectPool(ptr)
end

"""
    open(path; layout = "") -> ObjectPool

Open an existing object store memory pool. Similar to [`create`](@ref), `path` must identify
either an existing object memory pool file.

Return an handle `ObjectPool`. If allocation fails, returns a null `ObjectPool` which can
be tested with `isnull`.
"""
function open(path; layout = "")
    ptr = ccall(
        (:pmemobj_open, libpmemobj),
        Ptr{Nothing},
        (Cstring, Cstring),
        path, layout,
    )
    return ObjectPool(ptr)
end

"""
    check(path; layout = "") -> Bool

Perform a consistency check of the file indicated by `path`. Return `true` if the file
passes the check.
"""
function check(path; layout = "")
    val = ccall(
        (:pmemobj_check, libpmemobj),
        Cint,
        (Cstring, Cstring),
        path, layout
    )
    # According to the man page, a return value of zero indicates that something
    # went wrong
    return val != 0
end

#####
##### oid_is_null
#####

# http://pmem.io/pmdk/manpages/linux/v1.5/libpmemobj/oid_is_null.3

"""
    direct(oid::PersistentOID{T}) -> Ptr{T}

Return a pointer to the persistent memory object with the handle `oid`.
"""
function direct(oid::PersistentOID{T}) where {T}
    ccall((:pmemobj_direct, libpmemobj), Ptr{T}, (PersistentOID{T},), oid)
end

"""
    oid(ptr::Ptr{T}) -> PersistentOID{T}

Return a [`PersistentOID{T}`](@ref) handle to the object pointed to by `ptr`.
"""
oid(ptr::Ptr{T}) where {T} = ccall(
    (:pmemobj_oid, libpmemobj),
    PersistentOID{T},
    (Ptr{Cvoid},),
    ptr
)

"""
    pool(ptr::Ptr) -> ObjectPool

Return a [`ObjectPool`](@ref) handle to the pool containing the object containt the addresses
pointed to by `ptr`.
"""
function pool(ptr::Ptr)
    pool = ccall(
        (:pmemobj_pool_by_ptr, libpmemobj),
        Ptr{Cvoid},
        (Ptr{Cvoid},),
        ptr
    )
    return ObjectPool(pool)
end

"""
    pool(oid::PersistentOid} -> ObjectPool

Return an [`ObjectPool`](@ref) handle to the pool containing the object with handle `oid`.
"""
function pool(oid::PersistentOID{T}) where {T}
    pool = ccall(
        (:pmemobj_pool_by_oid, libpmemobj),
        Ptr{Cvoid},
        (PersistentOID{T},),
        oid
    )
    return ObjectPool(pool)
end

#####
##### pmemobj_root
#####

# http://pmem.io/pmdk/manpages/linux/v1.5/libpmemobj/pmemobj_root.3

"""
    root(pool, size, [::Type{T}]) -> PersistentOID{T}

Create or resize the root object for the persistent memory pool `pool`.  Return a
[`PersistentOID{T}`](@ref) handle to the allocated root object. If `T` is not provided,
default to `Nothing`.

*   If this is the first call to `root`, the requested size is greater than zero and the
    root object does not exist, it is implicitly allocated in a thread-safe manner.

*   If the requested size is larger than the current size, the root object is automatically
    resized. In such case, the old data is preserved and the extra space is zeroed.

*   If the requested size is equal to zero, the root object is not allocated.
"""
function root(pool::ObjectPool, size, ::Type{T} = Nothing) where {T}
    oid = ccall(
        (:pmemobj_root, libpmemobj),
        PersistentOID{T},
        (Ptr{Cvoid}, Csize_t),
        pool.ptr, size
    )
    return oid
end

"""
    root_size(pool::ObjectPool) -> Int

Return the size in bytes of the root object of `pool`.
"""
function root_size(pool::ObjectPool)
    size = ccall(
        (:pmemobj_root_size, libpmemobj),
        Csize_t,
        (Ptr{Cvoid},),
        pool.ptr
    )
    return Int(size)
end

#####
##### pmemobj_tx_begin
#####

# TODO: Last argument is likely wrong.
# Need to deal with the "enums" directly
function tx_begin(pop::ObjectPool)
    ret = ccall(
        (:pmemobj_tx_begin, libpmemobj),
        Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}, Cint),
        pop.ptr, C_NULL, TX_PARAM_NONE
    )

    ret == 0 || throw(error("Failed to begin transaction"))
    return ret
end

tx_abort(errnum::Integer) = ccall((:pmemobj_tx_abort, libpmemobj), Cvoid, (Cint,), errnum)
tx_commit() = ccall((:pmemobj_tx_commit, libpmemobj), Cvoid, ())
tx_end() = ccall((:pmemobj_tx_end, libpmemobj), Cint, ())
tx_process() = ccall((:pmemobj_tx_process, libpmemobj), Cvoid, ())

function tx_stage()
    stage = ccall((:pmemobj_tx_stage, libpmemobj), Cint, ())
    return TX_Stage(stage)
end

#####
##### pmemobj_tx_alloc
#####

# http://pmem.io/pmdk/manpages/linux/v1.5/libpmemobj/pmemobj_tx_alloc.3
#
function tx_alloc(size, ::Type{T} = Nothing, typenum = _typenum(T)) where {T}
    pmemoid = ccall(
        (:pmemobj_tx_alloc, libpmemobj),
        PersistentOID{T},
        (Csize_t, Culonglong),
        size, typenum
    )
    return pmemoid
end

function tx_zalloc(size, ::Type{T} = Nothing, typenum = _typenum(T)) where {T}
    pmemoid = ccall(
        (:pmemobj_tx_zalloc, libpmemobj),
        PersistentOID{T},
        (Csize_t, Culonglong),
        size, typenum
    )
    return pmemoid
end

function tx_free(oid::PersistentOID{T}) where {T}
    ret = ccall(
        (:pmemobj_tx_free, libpmemobj),
        Cint,
        (PersistentOID{T},),
        oid
    )
    ret != 0 && throw(error("Unable to free object: $oid"))
    return nothing
end

function tx_strdup(str, ::Type{T} = String) where {T}
    typenum = _typenum(T)
    oid = ccall(
        (:pmemobj_tx_strdup, libpmemobj),
        PersistentOID{T},
        (Cstring, Culonglong),
        str, typenum
    )
    return oid
end

#####
##### pmemobj_tx_add_range
#####

# http://pmem.io/pmdk/manpages/linux/v1.5/libpmemobj/pmemobj_tx_add_range.3
function add_range(oid::PersistentOID{T}, off, size) where {T}
    ret = ccall(
        (:pmemobj_tx_add_range, libpmemobj),
        Cint,
        (PersistentOID{T}, Culonglong, Csize_t),
        oid, off, size
    )
    ret != 0 && throw(error("Adding failed: error: $ret"))
    return nothing
end

function add_range_direct(ptr::Ptr, size)
    ret = ccall(
        (:pmemobj_tx_add_range_direct, libpmemobj),
        Cint,
        (Ptr{Cvoid}, Csize_t),
        ptr, size
    )
    ret != 0 && throw(error("Add range direct failed: error: $ret"))
    return nothing
end

#####
##### pmemobj_memcpy_persist
#####

# http://pmem.io/pmdk/manpages/linux/v1.5/libpmemobj/pmemobj_memcpy_persist.3
function persist(pool::ObjectPool, addr, size)
    ccall(
        (:pmemobj_persist, libpmemobj),
        Cvoid,
        (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
        pool.ptr, addr, size
    )
end

#####
##### pmemobj_alloc
#####

# http://pmem.io/pmdk/manpages/linux/v1.5/libpmemobj/pmemobj_alloc.3

function _constructor(pool_pointer::Ptr{Nothing}, obj_pointer::Ptr{T}, arg_pointer::Ptr{T}) where {T}
    unsafe_store!(obj_pointer, unsafe_load(arg_pointer))
    ccall(
        (:pmemobj_persist, libpmemobj),
        Cvoid,
        (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
        pool_pointer, obj_pointer, sizeof(T)
    )
    return Int32(0)
end

_donothing(args...) = Int32(0)

function alloc(pop::ObjectPool, size, x::T) where {T}
    typenum = _typenum(T)
    oidp = Ref{PersistentOID{T}}()
    f = @cfunction(_donothing, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}))

    ret = ccall(
        (:pmemobj_alloc, libpmemobj),
        Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t, Culonglong, Ptr{Cvoid}, Ptr{Cvoid}),
        pop.ptr, oidp, size, typenum, f, Ref(x)
    )
    ret != 0 && throw(error("Allocation failed"))

    return oidp[]
end

function zalloc(pop::ObjectPool, size, ::Type{T} = Nothing) where {T}
    typenum = _typenum(T)
    oidp = Ref{PersistentOID{T}}()

    ret = ccall(
        (:pmemobj_zalloc, libpmemobj),
        Cint,
        (Ptr{Cvoid}, Ptr{PersistentOID{T}}, Csize_t, Culonglong),
        pop.ptr, oidp, size, typenum
    )
    ret != 0 && throw(error("Allocation failed"))

    return oidp[]
end

function free(oid::PersistentOID{T}) where {T}
    ccall(
        (:pmemobj_free, libpmemobj),
        Cvoid,
        (Ptr{PersistentOID{T}},),
        Ref(oid)
    )

end

function alloc_usable_size(oid::PersistentOID{T}) where {T}
    return ccall(
        (:pmemobj_alloc_usable_size, libpmemobj),
        Csize_t,
        (PersistentOID{T},),
        oid,
    )
end

function strdup(pool::ObjectPool, str, ::Type{T} = UInt8) where {T}
    typenum = _typenum(T)
    oidp = Ref{PersistentOID{T}}()

    ret = ccall(
        (:pmemobj_strdup, libpmemobj),
        Cint,
        (Ptr{Cvoid}, Ptr{PersistentOID{T}}, Cstring, Culonglong),
        pool.ptr, oidp, str, typenum
    )
    ret != 0 && throw(error("Allocation failed"))

    return oidp[]
end

end # module Lib


text:
using TimeSpans
using Documenter

makedocs(modules=[TimeSpans],
         sitename="TimeSpans",
         authors="Beacon Biosignals, Inc.",
         pages=["API Documentation" => "index.md"])

deploydocs(repo="github.com/beacon-biosignals/TimeSpans.jl.git",
           devbranch="main")


text:
# This file holds the definition of the functions pertaining to the MPcondition
#   type.

function Base.show( io::IO, condition::MPcondition )::Nothing

    print( io, condition.attribute, " ", condition.operator, " " )

    if condition.value isa String
        print( io, "\"", condition.value, "\"" )
    elseif condition.value isa Vector{String}
        print( io, "{ \"", join( condition.value, "\", \"" ) ,"\" }" )
    else
        print( io, condition.value )
    end  # if condition.value isa String

    return

end  # show( io, condition )

include( joinpath( privPath, "condition.jl" ) )


text:
# # Point Group Elements

# Search elements by their Seitz name, and show their matrix representations

using LatticeTools

all_elements = Dict()
print("# Reading groups")
for i in 1:32
    print(" $i")
    psym = PointSymmetryDatabase.get(i)
    for (elname, matrep) in zip(psym.element_names, psym.matrix_representations)
        if !haskey(all_elements, elname)
            all_elements[elname] = []
        end
        push!(all_elements[elname], (matrep, psym.hermann_mauguin))
    end
end
println()

println("elements:")
for (k, v) in all_elements
    println("- element_name: $k")
    println("  matrix_representations:")
    for v2 in v
        m = [collect(x) for x in eachrow(v2[1])]
        g = v2[2]
        println("  - { matrix: $(m), group: \"$g\" }")
    end
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
using LinearAlgebra

S(x, A) = dot(x, A, x)
S(x) = S(x, A)

niters = 10^6
n = 10
A = randn(n, n)
x = randn(n);

# %%
@code_warntype S(x, A)

# %%
@code_warntype S(x)

# %%
function f(niters, x, A)
    y = Vector{eltype(x)}(undef, niters)
    for i in 1:niters
        y[i] = S(x, A)
    end
    y
end

@code_warntype f(niters, x, A)

# %%
function f(niters, x)
    y = Vector{eltype(x)}(undef, niters)
    for i in 1:niters
        y[i] = S(x)
    end
    y
end

@code_warntype f(niters, x)

# %%
using BenchmarkTools
@show f(niters, x, A) == f(niters, x)
@btime f(niters, $x, $A)
@btime f(niters, $x);

# %%
function g(niters, x, A, S)
    y = Vector{eltype(x)}(undef, niters)
    for i in 1:niters
        y[i] = S(x, A)
    end
    y
end

@code_warntype g(niters, x, A, S)

# %%
function g(niters, x, S)
    y = Vector{eltype(x)}(undef, niters)
    for i in 1:niters
        y[i] = S(x)
    end
    y
end

@code_warntype g(niters, x, S)

# %%
@show g(niters, x, A, S) == g(niters, x, S)
@btime g(niters, $x, $A, $S)
@btime g(niters, $x, $S);

# %%

max_stars_repo_path: 0018/S(x) = S(x, A) for global A causes memory allocations.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
import Base: show, display, length, collect, indices

type NeighborList{T}
    indices::Vector{Int}
    distances::Vector{T}
    length::Int
    max_length::Int

    NeighborList{T}(max_length::Int) where T = new(Vector{Int}(max_length), Vector{T}(max_length), 0, max_length)
end

length(neighbors::NeighborList) = neighbors.length
indices(neighbors::NeighborList) = neighbors.indices[1:neighbors.length]
distances(neighbors::NeighborList) = neighbors.distances[1:neighbors.length]

function show(io::IO, neighbors::NeighborList)
    for i in 1:neighbors.length
        @printf(io, "%-10d%.4f\n", neighbors.indices[i], neighbors.distances[i])
    end

    for i in 1:neighbors.max_length-neighbors.length
        @printf(io, "*\n")
    end
end

display(neighbors::NeighborList) = show(neighbors)

clear!(neighbors::NeighborList) = neighbors.length = 0

function insert_neighbor!{T}(neighbors::NeighborList{T}, index::Int, distance::T)
    ni = neighbors.indices
    nd = neighbors.distances
    k = neighbors.length

    if k < neighbors.max_length
        neighbors.length += 1
        k += 1
        @inbounds nd[k] = T(Inf)
    end

    j = 1
    @inbounds while distance > nd[j]
        j += 1
    end

    for m in k:-1:j+1
        @inbounds nd[m] = nd[m - 1]
        @inbounds ni[m] = ni[m - 1]
    end

    @inbounds ni[j] = index
    @inbounds nd[j] = distance

    @inbounds τ = k == neighbors.max_length ? nd[k] : T(Inf)

    return τ
end

@inline function add_neighbor!{T}(neighbors::NeighborList{T}, index::Int, distance::T)
    k = neighbors.length + 1
    @inbounds neighbors.indices[k] = index
    @inbounds neighbors.distances[k] = distance
    neighbors.length = k
end

function add_neighbors!{T}(metric, neighbors::NeighborList{T}, node::Node{T}, point_index::Int)
    node.index == 0 && return

    #if node.index != point_index
    #if d > 0.001

    d = metric(node.index, point_index)::T
    add_neighbor!(neighbors, node.index, d)

    node.isleaf && return

    add_neighbors!(metric, neighbors, node.inside, point_index)
    add_neighbors!(metric, neighbors, node.outside, point_index)
end

# Is sqrt(x1) + sqrt(x2) >= sqrt(y)?
@inline function sumsqrtgt{T}(x1::T, x2::T, y::T)
    d = y - x1 - x2
    d <= 0 || d^2 <= 4*x1*x2
end


text:
"""
    getcollection(db_name::String, collection_name::String, client=Mongoc.Client())

Return a `Mongoc.Collection` object corresponding to a MongoDB collection.
"""
function getcollection(db_name::String, collection_name::String, client=Mongoc.Client())
    db = client[db_name]
    return db[collection_name]
end

"""
    isinstance(collection::Mongoc.Collection, instance_name::String)

Check if an instance with "_id" == `instance_name` exists in `collection`.
See also [`getinstance`].
"""
function isinstance(collection::Mongoc.Collection, instance_name::String)
    document = Mongoc.BSON("_id" => instance_name)
    return Mongoc.count_documents(collection, document) >= 1
end

"""
    isdecomposition(collection::Mongoc.Collection, instance_name::String, added_edges::AbstractVector{N})) where N<:AbstractVector{Int}

Check if a decomposition with "_id" = {"instance_name" : `instance_name`, "added_edges" : `added_edges`} exists in `collection`.
See also [`getdecomposition`].
"""
function isdecomposition(collection::Mongoc.Collection, instance_name::String, oid::Union{UUID,AbstractString})
    document = Mongoc.BSON("_id" => Mongoc.BSON("instance_name" => instance_name, "oid" => oid))
    return Mongoc.count_documents(collection, document) >= 1
end

"""
    getinstance(collection::Mongoc.Collection, instance_name::String)

Return a `Mongoc.BSON` corresponding to the document in `collection` with "_id" == `instance_name`. Return `nothing` if document does not exist.
"""
function getinstance(collection::Mongoc.Collection, instance_name::String)
    selector = Mongoc.BSON("_id" => instance_name)
    return Mongoc.find_one(collection, selector)
end

"""
    getdecomposition(collection::Mongoc.Collection, instance_name::String, added_edges::Vector{Vector{Int}})

Return a `Mongoc.BSON` correspong to the document in `collection` width "_id" == { "instance_name" : `instance_name`, "added_edges" : `added_edges`}. Return `nothing` if document does not exist.
"""
function getdecomposition(collection::Mongoc.Collection, instance_name::String, oid::Union{UUID,AbstractString})
    selector = Mongoc.BSON("_id" => Mongoc.BSON("instance_name" => instance_name, "oid" => oid))
    return Mongoc.find_one(collection, selector)
end

"""
    issolved(collection::Mongoc.Collection, instance_name::String, added_edges::Vector{Vector{Int}})

Check if the field `features.solver.solving_time` exists in the specified decomposition.
"""
function issolved(collection::Mongoc.Collection, instance_name::String, oid::Union{UUID,AbstractString})
    document = Mongoc.BSON("_id" => Mongoc.BSON("instance_name" => instance_name, "oid" => oid),
                           "features.solver.solving_time" => Mongoc.BSON("\$exists" => true))
    return Mongoc.count_documents(collection, document) >= 1
end

"""
    getOPFpath_all(collection::Mongoc.Collection)

Return a dict containing the fields "paths.OPF_ctr" and "paths.OPF_mat" for each instance.

# Examples
```julia-repl
julia> getOPFpath_all(instances)
Dict{Any,Any} with 2 entries:
  "case1354pegase" => Dict("mat"=>"path/to/OPF/mat","ctr"=>"path/to/OPF/ctr")
  "case9241pegase" => Dict("mat"=>"path/to/OPF/mat","ctr"=>"path/to/OPF/ctr")
```
"""
function getOPFpath_all(collection::Mongoc.Collection)
    query = Mongoc.BSON("""
        { "paths.OPF_ctr": { "\$exists": true },
          "paths.OPF_mat": { "\$exists": true }
        }""")
    options = Mongoc.BSON("""
        { "projection":
            {
                "features": false,
                "paths.matpower": false,
                "decompositions": false
            }
        }""")
    path_dict = Dict()
    for instance in Mongoc.find(collection, query; options=options)
        path_dict[instance["_id"]] = Dict("ctr" => instance["paths"]["OPF_ctr"],
                                           "mat" => instance["paths"]["OPF_mat"])
    end
    return path_dict
end

"""
    getOPFpath_one(collection::Mongoc.Collection, instance_name::String)

Return a dict containing the path of the specified instance.

```julia-repl
julia> getOPFpath_one(instances, "case118mod")
Dict{String, String} with 2 entries:
  "mat" => "path/to/OPF/ctr"
  "ctr" => "path/to/OPF/ctr"
```
"""
function getOPFpath_one(collection::Mongoc.Collection, instance_name::String)
    query = Mongoc.BSON("_id" => instance_name)
    options = Mongoc.BSON("projection" => Mongoc.BSON("paths" => true))
    path_dict = Mongoc.find_one(collection, query; options=options)["paths"]
    isempty(path_dict) && return
    path_dict = Dict("ctr" => path_dict["OPF_ctr"], "mat" => path_dict["OPF_mat"])
    return path_dict
end

function getpath_one(collection::Mongoc.Collection, instance_name::String)
    query = Mongoc.BSON("_id" => instance_name)
    options = Mongoc.BSON("projection" => Mongoc.BSON("paths" => true))
    path_dict = Mongoc.find_one(collection, query; options=options)["paths"]
    isempty(path_dict) && return
    return path_dict
end

"""
    getcholesky(collection::Mongoc.Collection)

Return an array of BSON document where each BSON document is a decomposition where no edges has been added
(called "cholesky graph" in the context of this experiment).
"""
function getcholesky(collection::Mongoc.Collection)
    query = Mongoc.BSON("nb_added_edges" => 0, "dec_type" => "cholesky")
    return Mongoc.find(collection, query)
end

"""
    get_all_matpower_path(collection::Mongoc.Collection)

Get all the matpower path field in `collection` in the form of an array of `Mongoc.BSON` object.
"""
function getmatpowerpath_all(collection::Mongoc.Collection)
    results = collect(Mongoc.find(collection, Mongoc.BSON(); options=Mongoc.BSON("projection" => Mongoc.BSON("paths.matpower" => true))))
    return [path["paths"]["matpower"] for path in results if haskey(path["paths"], "matpower")]
end

function getdecompositions_index(collection::Mongoc.Collection;
                                 unsolved=[true],
                                 original=[false],
                                 merge=[true, false],
                                 combined=[true, false],
                                 select=nothing)::Vector{Int}
    indexes = []
    for (index, document) in enumerate(collection)
        isoriginal = document["nb_added_edges"] == 0
        isunsolved = !haskey(document["features"], "solver") || !haskey(document["features"]["solver"], "solving_time")
        ismerge = document["nb_added_edges"] == -1
        iscombined = document["nb_added_edges"] == -2
        if (isoriginal in original) && (isunsolved in unsolved) && (ismerge in merge) && (iscombined in combined)
            push!(indexes, index)
        end
    end
    if !isnothing(select)
        decompositions = collect(collection)
        filter!(i -> select(decompositions[i]), indexes)
    end
    return indexes
end

function getcholeskytime(collection::Mongoc.Collection, instance_name::String)
    query = Mongoc.BSON("_id.instance_name" => instance_name, "nb_added_edges" => 0, "dec_type" => "cholesky")
    options = Mongoc.BSON("projection" => Mongoc.BSON("features.solver.solving_time" => true))
    decomposition = Mongoc.find_one(collection, query; options=options)
    if decomposition != nothing && haskey(decomposition["features"], "solver") && haskey(decomposition["features"]["solver"], "solving_time")
        return decomposition["features"]["solver"]["solving_time"]
    end
end

function getcholeskytime_all(collection::Mongoc.Collection)
    query = Mongoc.BSON("nb_added_edges" => 0, "dec_type" => "cholesky")
    options = Mongoc.BSON("projection" => Mongoc.BSON("features.solver.solving_time" => true))
    decompositions = Mongoc.find(collection, query; options=options)
    if decompositions != nothing
        choleskytimes = Dict()
        for d in decompositions
            if haskey(d["features"], "solver")
                choleskytimes[d["_id"]["instance_name"]] = d["features"]["solver"]["solving_time"]
            end
        end
        return choleskytimes
    end
end

function getinstancevalue(collection::Mongoc.Collection, instance_name::String, key::String)
    query = Mongoc.BSON("_id" => instance_name)
    options = Mongoc.BSON("projection" => Mongoc.BSON(key => true))
    instance = Mongoc.find_one(collection, query; options=options)
    isnothing(instance) && return
    keys = split(key, ".")
    for k in keys
        if haskey(instance, k)
            if typeof(instance[k]) <: AbstractDict
                instance = instance[k]
            else
                return instance[k]
            end
        end
    end
end

function getnv(collection::Mongoc.Collection, instance_name::String)
    query = Mongoc.BSON("_id" => instance_name)
    options = Mongoc.BSON("projection" => Mongoc.BSON("features.graph.nv" => true))
    return Mongoc.find_one(collection, query; options=options)["features"]["graph"]["nv"]
end

function getne(collection::Mongoc.Collection, instance_name::String)
    query = Mongoc.BSON("_id" => instance_name)
    options = Mongoc.BSON("projection" => Mongoc.BSON("features.graph.ne" => true))
    return Mongoc.find_one(collection, query; options=options)["features"]["graph"]["ne"]
end

function getkneighbors_indexes(collection::Mongoc.Collection,
                               instance_name::String,
                               oid::Union{UUID,AbstractString},
                               k=3,
                               neighbors_g=getgraphs_dict(collection, x -> x["_id"]["instance_name"] == instance_name && x["features"]["solver"]["cholesky_diff"] <= 0))
    decomposition = getdecomposition(collection, instance_name, oid)
    g = construct_graph_from_cliques(decomposition["cliques"])
    dists = [edge_distance(g, ng[1]) for ng in neighbors_g[instance_name]]
    indexes = sortperm(dists)[begin:k]
    return [n[2] for n in neighbors_g[instance_name][indexes]]
end

function getgraphs_dict(collection::Mongoc.Collection, select=x -> haskey(x["features"], "solver") && x["features"]["solver"]["cholesky_diff"] <= 0)
    graphs = Dict()
    for (i, d) in enumerate(collection)
        if select(d)
            instance_name = d["_id"]["instance_name"]
            gtuple = (construct_graph_from_cliques(d["cliques"]), i)
            if instance_name in keys(graphs)
                push!(graphs[instance_name], gtuple)
            else
                graphs[instance_name] = [gtuple]
            end
        end
    end
    return graphs
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
using Newton
using Test

f1(x)=x^3
f1_prime(x)=3*x^2

f2(x) = sin(x)
f2_prime(x) = cos(x)

f3(x) = log(x)
f3_prime(x) = 1//x

fs=[f1,f2,f3]
f_primes=[f1_prime,f2_prime,f3_prime]
f_nc(x) = 2+ x^2
f_nc_prime(x) = 2*x

@testset "Newton.jl" begin

    for (f,f_prime) in zip(fs,f_primes) # Several test for the root of a known function,given f and analytical f′
        @test newtonroot(f,x₀=1).x≈newtonroot(f,f_prime,x₀=1).x
    end

    for (f,f_prime) in zip(fs,f_primes)
        for i in 1:200:1000 # test to ensure that maxiter is working
            @test  newtonroot(f,f_prime,x₀=1, maxiter=i)===nothing || newtonroot(f,f_prime,x₀=1, maxiter=i).iteration <= i
        end
    end

    for (f,f_prime) in zip(fs,f_primes)
        for i in 1:200:1000 # test to ensure that maxiter is working
            @test  newtonroot(f,f_prime,x₀=1, maxiter=i)===nothing || newtonroot(f,f_prime,x₀=1, maxiter=i).iteration <= i
        end
    end

    @test newtonroot(f_nc,x₀=1) === nothing # Test of non-convergence
    @test newtonroot(f_nc,f_nc_prime,x₀=1)=== nothing # Test of non-convergence

    @test newtonroot(f1,f1_prime,x₀=BigFloat(1)).x |> typeof == BigFloat # BigFloat roots
    @test newtonroot(f1,x₀=BigFloat(1)).x |> typeof == BigFloat # BigFloat roots

    @test newtonroot(f1,x₀=1) == nothing ? true : newtonroot(f1,x₀=1).norm_dif <=1E-5
    @test newtonroot(f1,f1_prime, x₀=1) == nothing ? true : newtonroot(f1,x₀=1).norm_dif <=1E-5

end


text:
function _JPoly(pointlist, color, linewidth, action, close, reversepath)
    sethue(color)
    setline(linewidth)
    poly(pointlist, action; close = close, reversepath = reversepath)
end

"""
    JPoly(pointlist::Vector{Point}; color="black", linewidth=2, action = :stroke, close=true, reversepath=false)

Draw a polygon around points in the pointlist.

# Keywords
- `color` specifies the color of the outline or the fill of it (depends on action)
- `linewidth` linewidth of the outline
- `action` can be `:stroke`, `:fill` or other symbols (check the Luxor documentation for details) (default: :stroke)
- `close` whether the polygon should be closed or not (default: closed)
- `reversepath` can be set to `true` to reverse the path and create a polygon hole
"""
JPoly(
    pointlist::Vector{Point};
    color = "black",
    linewidth = 2,
    action = :stroke,
    close = true,
    reversepath = false,
) =
    (
        args...;
        pointlist = pointlist,
        color = color,
        linewidth = linewidth,
        action = action,
    ) -> _JPoly(pointlist, color, linewidth, action, close, reversepath)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# data stored by the owner of a remote reference
def_rv_channel() = Channel(1)
mutable struct RemoteValue
    c::AbstractChannel
    clientset::BitSet # Set of workerids that have a reference to this channel.
                      # Keeping ids instead of a count aids in cleaning up upon
                      # a worker exit.

    waitingfor::Int   # processor we need to hear from to fill this, or 0

    synctake::Union{ReentrantLock, Nothing}  # A lock used to synchronize the
                      # specific case of a local put! / remote take! on an
                      # unbuffered store. github issue #29932

    function RemoteValue(c)
        c_is_buffered = false
        try
            c_is_buffered = isbuffered(c)
        catch
        end

        if c_is_buffered
            return new(c, BitSet(), 0, nothing)
        else
            return new(c, BitSet(), 0, ReentrantLock())
        end
    end
end

wait(rv::RemoteValue) = wait(rv.c)

# A wrapper type to handle issue #29932 which requires locking / unlocking of
# RemoteValue.synctake outside of lexical scope.
struct SyncTake
    v::Any
    rv::RemoteValue
end

## core messages: do, call, fetch, wait, ref, put! ##
struct RemoteException <: Exception
    pid::Int
    captured::CapturedException
end

"""
    RemoteException(captured)

Exceptions on remote computations are captured and rethrown locally.  A `RemoteException`
wraps the `pid` of the worker and a captured exception. A `CapturedException` captures the
remote exception and a serializable form of the call stack when the exception was raised.
"""
RemoteException(captured) = RemoteException(myid(), captured)
function showerror(io::IO, re::RemoteException)
    (re.pid != myid()) && print(io, "On worker ", re.pid, ":\n")
    showerror(io, re.captured)
end

function run_work_thunk(thunk, print_error)
    local result
    try
        result = thunk()
    catch err
        ce = CapturedException(err, catch_backtrace())
        result = RemoteException(ce)
        print_error && showerror(stderr, ce)
    end
    return result
end
function run_work_thunk(rv::RemoteValue, thunk)
    put!(rv, run_work_thunk(thunk, false))
    nothing
end

function schedule_call(rid, thunk)
    return lock(client_refs) do
        rv = RemoteValue(def_rv_channel())
        (PGRP::ProcessGroup).refs[rid] = rv
        push!(rv.clientset, rid.whence)
        errormonitor(@async run_work_thunk(rv, thunk))
        return rv
    end
end

function deliver_result(sock::IO, msg, oid, value)
    #print("$(myid()) sending result $oid\n")
    if msg === :call_fetch || isa(value, RemoteException)
        val = value
    else
        val = :OK
    end
    try
        send_msg_now(sock, MsgHeader(oid), ResultMsg(val))
    catch e
        # terminate connection in case of serialization error
        # otherwise the reading end would hang
        @error "Fatal error on process $(myid())" exception=e,catch_backtrace()
        wid = worker_id_from_socket(sock)
        close(sock)
        if myid()==1
            rmprocs(wid)
        elseif wid == 1
            exit(1)
        else
            remote_do(rmprocs, 1, wid)
        end
    end
end

## message event handlers ##
function process_messages(r_stream::TCPSocket, w_stream::TCPSocket, incoming::Bool=true)
    errormonitor(@async process_tcp_streams(r_stream, w_stream, incoming))
end

function process_tcp_streams(r_stream::TCPSocket, w_stream::TCPSocket, incoming::Bool)
    Sockets.nagle(r_stream, false)
    Sockets.quickack(r_stream, true)
    wait_connected(r_stream)
    if r_stream != w_stream
        Sockets.nagle(w_stream, false)
        Sockets.quickack(w_stream, true)
        wait_connected(w_stream)
    end
    message_handler_loop(r_stream, w_stream, incoming)
end

"""
    process_messages(r_stream::IO, w_stream::IO, incoming::Bool=true)

Called by cluster managers using custom transports. It should be called when the custom
transport implementation receives the first message from a remote worker. The custom
transport must manage a logical connection to the remote worker and provide two
`IO` objects, one for incoming messages and the other for messages addressed to the
remote worker.
If `incoming` is `true`, the remote peer initiated the connection.
Whichever of the pair initiates the connection sends the cluster cookie and its
Julia version number to perform the authentication handshake.

See also [`cluster_cookie`](@ref).
"""
function process_messages(r_stream::IO, w_stream::IO, incoming::Bool=true)
    errormonitor(@async message_handler_loop(r_stream, w_stream, incoming))
end

function message_handler_loop(r_stream::IO, w_stream::IO, incoming::Bool)
    wpid=0          # the worker r_stream is connected to.
    boundary = similar(MSG_BOUNDARY)
    try
        version = process_hdr(r_stream, incoming)
        serializer = ClusterSerializer(r_stream)

        # The first message will associate wpid with r_stream
        header = deserialize_hdr_raw(r_stream)
        msg = deserialize_msg(serializer)
        handle_msg(msg, header, r_stream, w_stream, version)
        wpid = worker_id_from_socket(r_stream)
        @assert wpid > 0

        readbytes!(r_stream, boundary, length(MSG_BOUNDARY))

        while true
            reset_state(serializer)
            header = deserialize_hdr_raw(r_stream)
            # println("header: ", header)

            try
                msg = invokelatest(deserialize_msg, serializer)
            catch e
                # Deserialization error; discard bytes in stream until boundary found
                boundary_idx = 1
                while true
                    # This may throw an EOF error if the terminal boundary was not written
                    # correctly, triggering the higher-scoped catch block below
                    byte = read(r_stream, UInt8)
                    if byte == MSG_BOUNDARY[boundary_idx]
                        boundary_idx += 1
                        if boundary_idx > length(MSG_BOUNDARY)
                            break
                        end
                    else
                        boundary_idx = 1
                    end
                end

                # remotecalls only rethrow RemoteExceptions. Any other exception is treated as
                # data to be returned. Wrap this exception in a RemoteException.
                remote_err = RemoteException(myid(), CapturedException(e, catch_backtrace()))
                # println("Deserialization error. ", remote_err)
                if !null_id(header.response_oid)
                    ref = lookup_ref(header.response_oid)
                    put!(ref, remote_err)
                end
                if !null_id(header.notify_oid)
                    deliver_result(w_stream, :call_fetch, header.notify_oid, remote_err)
                end
                continue
            end
            readbytes!(r_stream, boundary, length(MSG_BOUNDARY))

            # println("got msg: ", typeof(msg))
            handle_msg(msg, header, r_stream, w_stream, version)
        end
    catch e
        # Check again as it may have been set in a message handler but not propagated to the calling block above
        if wpid < 1
            wpid = worker_id_from_socket(r_stream)
        end

        if wpid < 1
            println(stderr, e, CapturedException(e, catch_backtrace()))
            println(stderr, "Process($(myid())) - Unknown remote, closing connection.")
        elseif !(wpid in map_del_wrkr)
            werr = worker_from_id(wpid)
            oldstate = werr.state
            set_worker_state(werr, W_TERMINATED)

            # If unhandleable error occurred talking to pid 1, exit
            if wpid == 1
                if isopen(w_stream)
                    @error "Fatal error on process $(myid())" exception=e,catch_backtrace()
                end
                exit(1)
            end

            # Will treat any exception as death of node and cleanup
            # since currently we do not have a mechanism for workers to reconnect
            # to each other on unhandled errors
            deregister_worker(wpid)
        end

        close(r_stream)
        close(w_stream)

        if (myid() == 1) && (wpid > 1)
            if oldstate != W_TERMINATING
                println(stderr, "Worker $wpid terminated.")
                rethrow()
            end
        end

        return nothing
    end
end

function process_hdr(s, validate_cookie)
    if validate_cookie
        cookie = read(s, HDR_COOKIE_LEN)
        if length(cookie) < HDR_COOKIE_LEN
            error("Cookie read failed. Connection closed by peer.")
        end

        self_cookie = cluster_cookie()
        for i in 1:HDR_COOKIE_LEN
            if UInt8(self_cookie[i]) != cookie[i]
                error("Process($(myid())) - Invalid connection credentials sent by remote.")
            end
        end
    end

    # When we have incompatible julia versions trying to connect to each other,
    # and can be detected, raise an appropriate error.
    # For now, just return the version.
    version = read(s, HDR_VERSION_LEN)
    if length(version) < HDR_VERSION_LEN
        error("Version read failed. Connection closed by peer.")
    end

    return VersionNumber(strip(String(version)))
end

function handle_msg(msg::CallMsg{:call}, header, r_stream, w_stream, version)
    schedule_call(header.response_oid, ()->msg.f(msg.args...; msg.kwargs...))
end
function handle_msg(msg::CallMsg{:call_fetch}, header, r_stream, w_stream, version)
    errormonitor(@async begin
        v = run_work_thunk(()->msg.f(msg.args...; msg.kwargs...), false)
        if isa(v, SyncTake)
            try
                deliver_result(w_stream, :call_fetch, header.notify_oid, v.v)
            finally
                unlock(v.rv.synctake)
            end
        else
            deliver_result(w_stream, :call_fetch, header.notify_oid, v)
        end
        nothing
    end)
end

function handle_msg(msg::CallWaitMsg, header, r_stream, w_stream, version)
    errormonitor(@async begin
        rv = schedule_call(header.response_oid, ()->msg.f(msg.args...; msg.kwargs...))
        deliver_result(w_stream, :call_wait, header.notify_oid, fetch(rv.c))
        nothing
    end)
end

function handle_msg(msg::RemoteDoMsg, header, r_stream, w_stream, version)
    errormonitor(@async run_work_thunk(()->msg.f(msg.args...; msg.kwargs...), true))
end

function handle_msg(msg::ResultMsg, header, r_stream, w_stream, version)
    put!(lookup_ref(header.response_oid), msg.value)
end

function handle_msg(msg::IdentifySocketMsg, header, r_stream, w_stream, version)
    # register a new peer worker connection
    w = Worker(msg.from_pid, r_stream, w_stream, cluster_manager; version=version)
    send_connection_hdr(w, false)
    send_msg_now(w, MsgHeader(), IdentifySocketAckMsg())
    notify(w.initialized)
end

function handle_msg(msg::IdentifySocketAckMsg, header, r_stream, w_stream, version)
    w = map_sock_wrkr[r_stream]
    w.version = version
end

function handle_msg(msg::JoinPGRPMsg, header, r_stream, w_stream, version)
    LPROC.id = msg.self_pid
    controller = Worker(1, r_stream, w_stream, cluster_manager; version=version)
    notify(controller.initialized)
    register_worker(LPROC)
    topology(msg.topology)

    if !msg.enable_threaded_blas
        Base.disable_library_threading()
    end

    lazy = msg.lazy
    PGRP.lazy = lazy

    @sync for (connect_at, rpid) in msg.other_workers
        wconfig = WorkerConfig()
        wconfig.connect_at = connect_at

        let rpid=rpid, wconfig=wconfig
            if lazy
                # The constructor registers the object with a global registry.
                Worker(rpid, ()->connect_to_peer(cluster_manager, rpid, wconfig))
            else
                @async connect_to_peer(cluster_manager, rpid, wconfig)
            end
        end
    end

    send_connection_hdr(controller, false)
    send_msg_now(controller, MsgHeader(RRID(0,0), header.notify_oid), JoinCompleteMsg(Sys.CPU_THREADS, getpid()))
end

function connect_to_peer(manager::ClusterManager, rpid::Int, wconfig::WorkerConfig)
    try
        (r_s, w_s) = connect(manager, rpid, wconfig)
        w = Worker(rpid, r_s, w_s, manager; config=wconfig)
        process_messages(w.r_stream, w.w_stream, false)
        send_connection_hdr(w, true)
        send_msg_now(w, MsgHeader(), IdentifySocketMsg(myid()))
        notify(w.initialized)
    catch e
        @error "Error on $(myid()) while connecting to peer $rpid, exiting" exception=e,catch_backtrace()
        exit(1)
    end
end

function handle_msg(msg::JoinCompleteMsg, header, r_stream, w_stream, version)
    w = map_sock_wrkr[r_stream]
    environ = something(w.config.environ, Dict())
    environ[:cpu_threads] = msg.cpu_threads
    w.config.environ = environ
    w.config.ospid = msg.ospid
    w.version = version

    ntfy_channel = lookup_ref(header.notify_oid)
    put!(ntfy_channel, w.id)

    push!(default_worker_pool(), w.id)
end

max_stars_repo_name: RelationalAI-oss/julia


text:
@ms include("solvers/annealing.jl")

function loop_annealing(args)
    instances = [
        # "02",
        # "03",
        # "05",
        # "08",
        # "09",
        # "10",
        # "11",
        "12",
        # "13"
    ]
    println("Instance n° & Cost & Time")
    for instance_name in instances
        instance_path = "data/$instance_name.alp"
        inst = Instance(instance_path)
        sol = Solution(inst)

        initial_sort!(sol)

        sv = AnnealingSolver(
            inst;
            temp_init_rate=0.8,
            step_size=1,
            startsol=sol,
            temp_coef=0.985,
        )

        ms_start = ms() # nb secondes depuis démarrage avec précision à la ms
        solve(sv, swap_close_planes!, durationmax=15*60)
        ms_stop = ms()

        bestsol = sv.bestsol

        nb_calls = bestsol.solver.nb_calls
        nb_infeasable = bestsol.solver.nb_infeasable
        nb_sec = round(ms_stop - ms_start, digits=3)
        nb_call_per_sec = round(nb_calls/nb_sec, digits=3)
        println("Performance: ")
        println("  nb_calls=$nb_calls")
        println("  nb_infeasable=$nb_infeasable")
        println("  nb_sec=$nb_sec")
        println("  => nb_call_per_sec = $nb_call_per_sec call/sec")
        println("$instance_name & $(bestsol.cost) & $nb_sec")
        end
end

loop_annealing(Args.args)


text:
module TestDimensions

using Mimi
using Test

import Mimi:
    AbstractDimension, RangeDimension, Dimension, key_type,
    reset_compdefs

reset_compdefs()

dim_varargs = Dimension(:foo, :bar, :baz)   # varargs
dim_vec = Dimension([:foo, :bar, :baz]) # Vector
dim_range = Dimension(2010:2100)     # AbstractRange
rangedim = RangeDimension(2010:2100) # RangeDimension type
dim_vals = Dimension(4) # Same as 1:4

# RJP: unclear what was intended here, but you cannot instantiate an abstract type...
# dim_vals_abstract = AbstractDimension(dim_vals) # Abstract

@test key_type(dim_varargs) == Symbol
@test key_type(dim_vec) == Symbol
@test key_type(dim_range) == Int
@test key_type(dim_vals) == Int

@test length(dim_varargs) == 3
@test length(dim_vec) == 3
@test length(dim_vals) == 4
@test length(dim_range) == 91
@test length(rangedim) == 91

# Test iteration
@test [x.first for x in dim_varargs] == collect(keys(dim_varargs))
@test [x.first for x in dim_range]   == collect(keys(dim_range))

@test lastindex(dim_varargs) == :baz
@test lastindex(dim_vec) == :baz
@test lastindex(dim_range) == 2100
@test lastindex(dim_vals) == 4

@test Base.keys(rangedim) == [2010:2100...]
@test Base.values(rangedim) == [1:91...]

for i = 1:length(dim_range)
    @test getindex(dim_range, 2010 + i - 1) == i
end
@test dim_range[:] == [1:91...]
@test dim_varargs[:bar] == 2
@test dim_varargs[:] == [1,2,3]

# @test dim_vals_abstract[1:4...]== [1:4...]

# @test rangedim[2011] == 2 # TODO: this errors..

@test get(dim_varargs, :bar, 999) == 2
@test get(dim_varargs, :new, 4) == 4 #adds a key/value pair
@test get(rangedim, 2010, 1) == 1
# @test get(rangedim, 2101, 92) == 92 # TODO:  this errors ...

#test iteratable
dim_vals2 = Dimension(2:2:8)
intkeys = [2,4,6,8]
# Work around new global scoping rules
# Without the `let`, `index` is unknown in the for loop
let index = 1
    for pair in dim_vals2
        @test dim_vals2[intkeys[index]] == index
        index += 1
    end
end

rangedim2 = RangeDimension(2:2:8)
# Work around new global scoping rules
# Without the `let`, `index` is unknown in the for loop
let index = 1
    for pair in rangedim2   # uses iterate()
        @test get(rangedim2, Base.keys(rangedim2)[index]) == index
        index += 1
    end
end

@test getindex(dim_varargs, :bar) == 2
@test getindex(dim_varargs, :) == [1,2,3]

# Test resetting the time dimension

@defcomp foo2 begin
    x = Parameter(index=[time])
    y = Variable(index=[4])
end

m = Model()
set_dimension!(m, :time, 2000:2100)
@test_throws ErrorException add_comp!(m, foo2; first = 2005, last = 2105)   # Can't add a component longer than a model
add_comp!(m, foo2; first = 2005, last = 2095)

# Test that foo's time dimension is unchanged
@test_logs(
    (:warn, "Redefining dimension :time"),
    set_dimension!(m, :time, 1990:2200)
)
@test m.md.comp_defs[:foo2].first == 2005
@test m.md.comp_defs[:foo2].last == 2095

# Test parameter connections
@test_throws ErrorException set_param!(m, :foo2, :x, 1990:2200) # too long
set_param!(m, :foo2, :x, 2005:2095) # Shouldn't throw an error

# Test that foo's time dimension is updated
@test_logs(
    (:warn, "Redefining dimension :time"),
    (:warn, "Resetting foo2 component's first timestep to 2010"),
    (:warn, "Resetting foo2 component's last timestep to 2050"),
    set_dimension!(m, :time, 2010:2050)
)
@test m.md.comp_defs[:foo2].first == 2010
@test m.md.comp_defs[:foo2].last == 2050

end #module


text:
# Graph → Syntax

tocall(f, a...) = :($f($(a...)))

binding(bindings::Associative, v) =
  haskey(bindings, v) ? bindings[v] : (bindings[v] = gensym("edge"))

function syntax(head::DVertex; bindconst = !isfinal(head))
  vs = topo(head)
  ex, bs = :(;), ObjectIdDict()
  for v in vs
    x = tocall(value(v), [binding(bs, n) for n in inputs(v)]...)
    if !bindconst && isconstant(v) && nout(v) > 1
      bs[v] = value(v).value
    elseif nout(v) > 1 || (!isfinal(head) && v ≡ head)
      edge = binding(bs, v)
      push!(ex.args, :($edge = $x))
    elseif haskey(bs, v)
      if MacroTools.inexpr(ex, bs[v])
        ex = MacroTools.replace(ex, bs[v], x)
      else
        push!(ex.args, :($(bs[v]) = $x))
      end
    else
      isfinal(v) ? push!(ex.args, x) : (bs[v] = x)
    end
  end
  head ≢ vs[end] && push!(ex.args, binding(bs, head))
  return ex
end

function constructor(g)
  vertex = isa(g, DVertex) ? :dvertex : :vertex
  g = mapv(g) do v
    if isconstant(v)
      prethread!(v, typeof(v)(value(v)))
      v.value = :constant
    else
      prethread!(v, typeof(v)(Constant(value(v))))
      v.value = vertex
    end
    v
  end
  ex = syntax(g)
  decls, exs = [], []
  for x in block(ex).args
    if @capture(x, v_ = $vertex(f_, a__))
      push!(decls, :($v = $vertex($f)))
      push!(exs, :(thread!($v, $(a...))))
    else
      push!(exs, x)
    end
  end
  return :($(decls...);$(exs...))
end


text:
using Arrow, DataFrames
using PyPlot
using JuliaForHEP

include("variables.jl")
include("plot_styles.jl")
using .Vars: classes

basedir = "/work/sschaub/JuliaForHEP/final_plotting/"
output_dir = joinpath(basedir, "lola+all_1")
all_features = DataFrame(Arrow.Table(joinpath(output_dir, "all_features.arrow")))

fig1, axs1 = subplots(4, 2, figsize=(18, 30))
fig2, axs2 = subplots(3, 2, figsize=(18, 22.5))
df = groupby(all_features, :output_expected; sort=true)
for (ax, feature) in zip(permutedims([axs1; axs2]), Vars.scalar_features)
    ax.hist([i[!, feature] for i in df],
        weights=[i.weights for i in df],
        bins=30, stacked=true, label=first.(keys(df)))
    unit = any(occursin(feature), ["Evt_M", "Evt_Pt", "toplep_m"]) ? " in GeV" : ""
    ax.set_xlabel("\\verb|$feature|$unit")
    ax.set_ylabel("Events/Bin")
    ax.set_yscale(:log)

    leg = ax.legend()
    if feature in Vars.scalar_features[[end-3; end-1:end]]
        annotate_cms(ax; sel_pos=(.25, .95))
    else
        annotate_cms(ax; sel_pos=(.05, .95))
    end
end
axs2[end].axis(:off)
for (i, fig) in enumerate([fig1, fig2])
    fig.suptitle("Distributions of Scalar Features ($i)")
    fig.tight_layout()
    display(fig)
    fig.savefig(joinpath(output_dir, "scalar_features_$i.pdf"))
end

_all_features = let
    basedir = "/local/scratch/ssd/sschaub/h5files/ttH_sebastian-processed/"
    h5files = joinpath.(basedir, ["ttH", "ttbb", "ttcc", "ttlf"] .* "_lola.h5")
    names = [
        mapreduce(vcat, 0:5) do i
            ["Jet_Pt[$i]", "Jet_Eta[$i]", "Jet_Phi[$i]", "Jet_M[$i]"]
        end
        "TightLepton_" .* ["Pt", "Eta", "Phi", "M"] .* "[0]"
        "Weight_XS"
    ]
    df = DataFrame(extract_cols(h5files, names)', names)
    [df all_features]
end

fig, axs = subplots(2, 2, figsize=(18, 16))
for (ax, feature, unit) in zip(axs, [r"_Pt\[", r"_Eta\[", r"_Phi\[", r"_M\["],
                               [" in GeV", "", "", " in GeV"])
    _df = groupby(stack(_all_features, feature), :output_expected; sort=true)
    #_df = [stack(i, feature) for i in _df]
    ax.hist([i.value for i in _df],
        weights=[i.weights for i in _df],
        bins=30, stacked=true, label=first.(keys(_df)))
    ax.set_xlabel("\\verb|$(feature.pattern[2:end-2])|$unit")
    ax.set_ylabel("Events/Bin")
    ax.set_yscale(:log)

    ax.legend()
    annotate_cms(ax; sel_pos=(.05, .95))
end
fig.suptitle("Distributions of Four Momenta")
fig.tight_layout()
display(fig)
fig.savefig(joinpath(output_dir, "four_momenta.pdf"))

for i in 0:6
fig, axs = subplots(4, 2, figsize=(18, 30))
axs[2, 2].axis(:off)
for (ax, feature, (lo, hi), unit) in zip([axs[1:2, :][1:3]; axs[3:4, :][:]], Vars.names_lola_out[7i+1:end],
                                [(0, 10000); (0, 1000); (-200, 1200); fill((-5e5, 5e5), 2); fill((-5e5, 0), 2)],
                                ["GeV\$^2\$"; "GeV"; "GeV"; fill("GeV\$^2\$", 4)])
    ax.hist([i[!, feature] for i in df],
        weights=[i.weights for i in df],
        bins=range(lo, hi; length=31), stacked=true, label=first.(keys(df)))
    feature = replace(feature, r"Jet_(.*)\[([0-5])\]" => s"\\verb|\1| of Jet \2")
    feature = replace(feature, r"Lepton_(.*)" => s"\\verb|\1| of the Lepton")
    ax.set_xlabel("$feature in $unit")
    ax.set_ylabel("Events/Bin")
    ax.set_yscale(:log)

    ax.legend()
    annotate_cms(ax; sel_pos=(.05, .95))
end
fig.suptitle("LoLa outputs")
fig.tight_layout()
display(fig)
fig.savefig(joinpath(output_dir, "lola_outputs_$i.pdf"))
end


text:
module EnsembleMethods

using Statistics
import StatsBase
import IterTools: product
import MLBase

# standard included modules
using DataFrames
using Random

using ..AbsTypes
using ..Utils

import ..AbsTypes: fit!, transform!
export fit!, transform!
export VoteEnsemble, StackEnsemble, BestLearner

using ..DecisionTreeLearners

"""
    VoteEnsemble(
       Dict(
          # Output to train against
          # (:class).
          :output => :class,
          # Learners in voting committee.
          :learners => [PrunedTree(), Adaboost(), RandomForest()]
       )
    )

Set of machine learners employing majority vote to decide prediction.

Implements: `fit!`, `transform!`
"""
mutable struct VoteEnsemble <: Learner
  name::String
  model::Dict{Symbol,Any}

  function VoteEnsemble(args=Dict())
     default_args = Dict{Symbol,Any}(
      :name => "votingens",
      # Output to train against
      # (:class).
      :output => :class,
      # Learners in voting committee.
      :learners => [PrunedTree(), Adaboost(), RandomForest()]
    )
    cargs = nested_dict_merge(default_args, args)
    cargs[:name] = cargs[:name]*"_"*randstring(3)
    new(cargs[:name],cargs)
  end
end

function VoteEnsemble(learners::Vector{<:Learner}, opt::Dict)
   VoteEnsemble(Dict(:learners => learners, opt...))
end

function VoteEnsemble(learners::Vector{<:Learner}; opt...)
   VoteEnsemble(Dict(:learners => learners, Dict(pairs(opt))...))
end

function VoteEnsemble(learners::Vararg{<:Learner})
   (eltype(learners) <: Learner) || throw(ArgumentError("argument setup error"))
   v=[x for x in learners] # convert tuples to vector
   VoteEnsemble(v)
end

"""
    fit!(ve::VoteEnsemble, instances::DataFrame, labels::Vector)

Training phase of the ensemble.
"""
function fit!(ve::VoteEnsemble, instances::DataFrame, labels::Vector)
  # Train all learners
  learners = ve.model[:learners]
  for learner in learners
    fit!(learner, instances, labels)
  end
  ve.model[:learners] = learners
end

"""
    transform!(ve::VoteEnsemble, instances::DataFrame)

Prediction phase of the ensemble.
"""
function transform!(ve::VoteEnsemble, instances::DataFrame)
  # Make learners vote
  learners = ve.model[:learners]
  predictions = map(learner -> transform!(learner, instances), learners)
  # Return majority vote prediction
  return StatsBase.mode(predictions)
end

"""
    StackEnsemble(
       Dict(
          # Output to train against
          # (:class).
          :output => :class,
          # Set of learners that produce feature space for stacker.
          :learners => [PrunedTree(), Adaboost(), RandomForest()],
          # Machine learner that trains on set of learners' outputs.
          :stacker => RandomForest(),
          # Proportion of training set left to train stacker itself.
          :stacker_training_proportion => 0.3,
          # Provide original features on top of learner outputs to stacker.
          :keep_original_features => false
       )
    )

An ensemble where a 'stack' of learners is used for training and prediction.
"""
mutable struct StackEnsemble <: Learner
  name::String
  model::Dict{Symbol,Any}

  function StackEnsemble(args=Dict())
     default_args = Dict{Symbol,Any}(
      :name => "stackens",
      # Output to train against
      # (:class).
      :output => :class,
      # Set of learners that produce feature space for stacker.
      :learners => [PrunedTree(), Adaboost(), RandomForest()],
      # Machine learner that trains on set of learners' outputs.
      :stacker => RandomForest(),
      # Proportion of training set left to train stacker itself.
      :stacker_training_proportion => 0.3,
      # Provide original features on top of learner outputs to stacker.
      :keep_original_features => false
    )
    cargs = nested_dict_merge(default_args, args)
    cargs[:name] = cargs[:name]*"_"*randstring(3)
    new(cargs[:name],cargs)
  end
end

function StackEnsemble(learners::Vector{<:Learner}, opt::Dict)
   StackEnsemble(Dict(:learners => learners, opt...))
end

function StackEnsemble(learners::Vector{<:Learner}; opt...)
   StackEnsemble(Dict(:learners => learners, Dict(pairs(opt))...))
end

function StackEnsemble(learners::Vararg{<:Learner})
   (eltype(learners) <: Learner) || throw(ArgumentError("argument setup error"))
   v=[x for x in learners] # convert tuples to vector
   StackEnsemble(v)
end

"""
    fit!(se::StackEnsemble, instances::DataFrame, labels::Vector)

Training phase of the stack of learners.

- perform holdout to obtain indices for
- partition learner and stacker training sets
- partition training set for learners and stacker
- train all learners
- train stacker on learners' outputs
- build final model from the trained learners
"""
function fit!(se::StackEnsemble, instances::DataFrame, labels::Vector)
  learners = se.model[:learners]
  num_learners = size(learners, 1)
  num_instances = size(instances, 1)
  num_labels = size(labels, 1)

  # Perform holdout to obtain indices for
  # partitioning learner and stacker training sets
  shuffled_indices = randperm(num_instances)
  stack_proportion = se.model[:stacker_training_proportion]
  (learner_indices, stack_indices) = holdout(num_instances, stack_proportion)

  # Partition training set for learners and stacker
  learner_instances = instances[learner_indices, :]
  stack_instances = instances[stack_indices, :]
  learner_labels = labels[learner_indices]
  stack_labels = labels[stack_indices]

  # Train all learners
  for learner in learners
    fit!(learner, learner_instances, learner_labels)
  end

  # Train stacker on learners' outputs
  label_map = MLBase.labelmap(labels)
  stacker = se.model[:stacker]
  keep_original_features = se.model[:keep_original_features]
  stacker_instances = build_stacker_instances(learners, stack_instances,
                      label_map, keep_original_features) |> DataFrame
  fit!(stacker, stacker_instances, stack_labels)

  # Build model
  se.model[:learners] = learners
  se.model[:stacker] = stacker
  se.model[:label_map] = label_map
  se.model[:keep_original_features] = keep_original_features
end

"""
    transform!(se::StackEnsemble, instances::DataFrame)

Build stacker instances and predict
"""
function transform!(se::StackEnsemble, instances::DataFrame)
  # Build stacker instances
  learners = se.model[:learners]
  stacker = se.model[:stacker]
  label_map = se.model[:label_map]
  keep_original_features = se.model[:keep_original_features]
  stacker_instances = build_stacker_instances(learners, instances,
                      label_map, keep_original_features) |> DataFrame

  # Predict
  return transform!(stacker, stacker_instances)
end

# Build stacker instances.
function build_stacker_instances(
  learners::Vector{T}, instances::DataFrame,
  label_map, keep_original_features=false) where T<:Learner

  # Build empty stacker instance space
  num_labels = size(label_map.vs, 1)
  num_instances = size(instances, 1)
  num_learners = size(learners, 1)
  stacker_instances = zeros(num_instances, num_learners * num_labels)

  # Fill stack instances with predictions from learners
  for l_index = 1:num_learners
    predictions = transform!(learners[l_index], instances)
    for p_index in 1:size(predictions, 1)
      pred_encoding = MLBase.labelencode(label_map, predictions[p_index])
      pred_column = (l_index-1) * num_labels + pred_encoding
      stacker_instances[p_index, pred_column] =
        one(eltype(stacker_instances))
    end
  end

  # Add original features to stacker instance space if enabled
  if keep_original_features
    stacker_instances = [instances stacker_instances]
  end

  # Return stacker instances
  return stacker_instances
end

"""
    BestLearner(
       Dict(
          # Output to train against
          # (:class).
          :output => :class,
          # Function to return partitions of instance indices.
          :partition_generator => (instances, labels) -> kfold(size(instances, 1), 5),
          # Function that selects the best learner by index.
          # Arg learner_partition_scores is a (learner, partition) score matrix.
          :selection_function => (learner_partition_scores) -> findmax(mean(learner_partition_scores, dims=2))[2],
          # Score type returned by score() using respective output.
          :score_type => Real,
          # Candidate learners.
          :learners => [PrunedTree(), Adaboost(), RandomForest()],
          # Options grid for learners, to search through by BestLearner.
          # Format is [learner_1_options, learner_2_options, ...]
          # where learner_options is same as a learner's options but
          # with a list of values instead of scalar.
          :learner_options_grid => nothing
       )
    )

Selects best learner from the set by performing a
grid search on learners if grid option is indicated.
"""
mutable struct BestLearner <: Learner
  name::String
  model::Dict{Symbol,Any}

  function BestLearner(args=Dict{Symbol,Any}())
     default_args = Dict{Symbol,Any}(
      :name => "bestens",
      # Output to train against
      # (:class).
      :output => :class,
      # Function to return partitions of instance indices.
      :partition_generator => (instances, labels) -> kfold(size(instances, 1), 5),
      # Function that selects the best learner by index.
      # Arg learner_partition_scores is a (learner, partition) score matrix.
      :selection_function => (learner_partition_scores) -> findmax(mean(learner_partition_scores, dims=2))[2],
      # Score type returned by score() using respective output.
      :score_type => Real,
      # Candidate learners.
      :learners => [PrunedTree(), Adaboost(), RandomForest()],
      # Options grid for learners, to search through by BestLearner.
      # Format is [learner_1_options, learner_2_options, ...]
      # where learner_options is same as a learner's options but
      # with a list of values instead of scalar.
      :learner_options_grid => nothing
    )
    cargs = nested_dict_merge(default_args, args)
    cargs[:name] = cargs[:name]*"_"*randstring(3)
    new(cargs[:name],cargs)
  end
end

function BestLearner(learners::Vector{<:Learner}, opt::Dict)
   BestLearner(Dict(:learners => learners, opt...))
end

function BestLearner(learners::Vector{<:Learner}; opt...)
   BestLearner(Dict(:learners => learners, Dict(pairs(opt))...))
end

function BestLearner(learners::Vararg{Learner})
   (eltype(learners) <: Learner) || throw(ArgumentError("argument setup error"))
   v=[x for x in learners] # convert tuples to vector
   BestLearner(v)
end

"""
    fit!(bls::BestLearner, instances::DataFrame, labels::Vector)

Training phase:

- obtain learners as is if grid option is not present
- generate learners if grid option is present
- foreach prototype learner, generate learners with specific options found in grid
- generate partitions
- train each learner on each partition and obtain validation output
"""
function fit!(bls::BestLearner, instances::DataFrame, labels::Vector)
  # Obtain learners as is if no options grid present
  if bls.model[:learner_options_grid] == nothing
    learners = bls.model[:learners]
  # Generate learners if options grid present
  else
    # Foreach prototype learner, generate learners with specific options
    # found in grid.
    learners = Transformer[]
    for l_index in 1:length(bls.model[:learners])
      # Obtain options grid
      options_prototype = bls.model[:learner_options_grid][l_index]
      grid_list = nested_dict_to_tuples(options_prototype)
      grid_keys = map(x -> x[1], grid_list)
      grid_values = map(x -> x[2], grid_list)

      # Foreach combination of options
      # generate learner.
      for combination in product(grid_values...)
        # Assign values for each option
        learner_options = deepcopy(options_prototype)
        for g_index in 1:length(grid_list)
          keys = grid_keys[g_index]
          value = combination[g_index]
          nested_dict_set!(learner_options, keys, value)
        end

        # Generate learner
        learner_prototype = bls.model[:learners][l_index]
        learner = create_transformer(learner_prototype, learner_options)

        # Append to candidate learners
        push!(learners, learner)
      end
    end
  end

  # Generate partitions
  partition_generator = bls.model[:partition_generator]
  partitions = partition_generator(instances, labels)

  # Train each learner on each partition and obtain validation output
  num_partitions           = size(partitions, 1)
  num_learners             = size(learners, 1)
  num_instances            = size(instances, 1)
  score_type               = bls.model[:score_type]
  learner_partition_scores = Array{score_type}(undef,num_learners, num_partitions)
  for l_index = 1:num_learners, p_index = 1:num_partitions
    partition = partitions[p_index]
    rest = setdiff(1:num_instances, partition)
    learner = learners[l_index]

    training_instances   = instances[partition,:]
    training_labels      = labels[partition]
    validation_instances = instances[rest, :]
    validation_labels    = labels[rest]

    fit!(learner, training_instances, training_labels)
    predictions = transform!(learner, validation_instances)
    result = score(:accuracy, validation_labels, predictions)
    learner_partition_scores[l_index, p_index] = result
  end

  # Find best learner based on selection function
  best_learner_index =
    bls.model[:selection_function](learner_partition_scores)
  best_learner = learners[best_learner_index]

  # Retrain best learner on all training instances
  fit!(best_learner, instances, labels)

  # Create model
  bls.model[:best_learner]             = best_learner
  bls.model[:best_learner_index]       = best_learner_index
  bls.model[:learners]                 = learners
  bls.model[:learner_partition_scores] = learner_partition_scores
end

"""
    transform!(bls::BestLearner, instances::DataFrame)

Choose the best learner based on cross-validation results and use it for prediction.
"""
function transform!(bls::BestLearner, instances::DataFrame)
   transform!(bls.model[:best_learner], instances)
end

end # module


text:
module maagar

"""
    append_line(key::String, value::String)

documentation
"""
function append_line(key::String, value::String)
    len::UInt16 = ncodeunits(key) + ncodeunits(value)
    len_key::UInt8 = ncodeunits(key)
    open("./maagar", append=true) do io
        write(io, len, len_key, key, value, "\n")
    end
end # function

"""
    write_loop(c::Channel)

documentation
"""
function write_loop(c::Channel)
    for key_val in c
        append_line(key_val[1], key_val[2])
    end
end # function

"""
    find_by_key(args)

documentation
"""
function find_by_key(search_key::String)
    open("./maagar", read=true) do io
        while true
            len = read(io, 2)
            if length(len) == 0
                return nothing
            end
            len = reinterpret(UInt16, len)[1]
            len_key = reinterpret(UInt8, read(io, 1))[1]
            global key = String(read(io, len_key))
            global value = String(read(io, len - len_key))
            if key == search_key
                break
            end
        end
        return (key, value)
    end
end # function

"""
    initialize_loop()

documentation
"""
function initialize_loop()::Channel
    channel = Channel(write_loop)
    return channel
end # function

channel = initialize_loop()
put!(channel, ("1", "testing"))
put!(channel, ("2", "testing1"))
put!(channel, ("3", "testing2"))
put!(channel, ("4", "testing3"))

println(find_by_key("4"))
println(find_by_key("5"))

end # module

max_stars_repo_name: mosheduminer/maagar

max_issues_repo_name: mosheduminer/maagar

max_forks_repo_name: mosheduminer/maagar