text:
export create_fields, create_operators

"""
function create_fields(; names = (),
                         grid = nothing,
                         transform = nothing,
                         arraytype = Array,
                         floattype = ComplexF64)
# Description
Automates the constructions of fourier fields with names.

# Arguments
none

# Keyword Arguments
- `names`: tuple of strings
- `grid`: FourierGrid object
- `transfrom`: Transform object
- `arraytype`: default = Array,
- `floattype`: default = ComplexF64

# Return
- Array of fourier field objects with names
"""
function create_fields(; names = (),
                         grid = nothing,
                         arraytype = Array,
                         floattype = ComplexF64)
    transform = Transform(grid, arraytype=arraytype)
    dimensions = size(grid)
    state = []
    for name in names
        fmd = FourierMetaData(name, grid, transform)
        field_data = arraytype(zeros(floattype, dimensions...))
        push!(state, FourierField(field_data, fmd))
    end
    if length(state) == 1
        return state[1]
    else
        return state
    end
end

"""
function create_operators(g::FourierGrid;
                          names = ("x","y","z"),
                          arraytype = Array)
# Description
A convenience function to automatically put operators in global scope

# Arguments
- `g`: FourierGrid. A fourier grid object

# Keyword Arguments
- `names`: name for operators assumes at most 3 dimensional

# Return
several operators
"""
function create_operators(g::FourierGrid;
                          names = ("x","y","z"),
                          arraytype = Array)
    @assert length(names) <= 3
    ops = []
    for i in 1:length(g.grid)
        operator_name = Char(0x02202) * names[i]
        k = g.wavenumbers[i]
        op = arraytype(im .* k)
        fmd = FourierOperatorMetaData(operator_name)
        push!(ops, FourierOperator(op, fmd))
    end
    if length(ops) == 1
        return ops[1]
    else
        return ops
    end
end

## DO NOT USE THESE FUNCTIONS NORMALLY !!!!!
"""
function create_fields(mod::Module; names = (),
                         grid = nothing,
                         transform = nothing,
                         arraytype = Array,
                         floattype = ComplexF64)
# Description
Automates the constructions of fourier fields with names
# Arguments
- `mod`: Module for evaluating functions. Typically mod = @__MODULE__
# Keyword Arguments
- `names`: tuple of strings
- `grid`: FourierGrid object
- `transfrom`: Transform object
- `arraytype`: default = Array,
- `floattype`: default = ComplexF64
# Return
nothing
"""
function create_fields(mod::Module; names = (),
                         grid = nothing,
                         arraytype = Array,
                         floattype = ComplexF64)
    transform = Transform(grid, arraytype=arraytype)
    printstyled("Warning!!! ", color = :red)
    print("the name(s) ")
    printstyled("(", color = :blue)
    for name in names
        printstyled(name, ", ", color = :blue)
    end
    printstyled(")", color = :blue)
    print(" are being overwritten in the ")
    print("global scope.")
    println(" ")
    dimensions = size(grid)
    for name in names
        local fmd = FourierMetaData(name, grid, transform)
        local field_data = arraytype(zeros(floattype, dimensions...))
        local parsed_name = Meta.parse(string(name))
        Base.eval(mod, :($parsed_name = FourierField($field_data, $fmd)) )
    end
    return nothing
end

"""
function create_operators(mod::Module, g::FourierGrid;
                          names = ("x","y","z"),
                          arraytype = Array)
# Description
A convenience function to automatically put operators in global scope.
Should only be used for quick tests.

# Arguments
- `mod`: Module for evaluating functions. Typically mod = @__MODULE__
- `g`: FourierGrid. A fourier grid object
# Keyword Arguments
- `names`: name for operators assumes at most 3 dimensional
# Return
Nothing. But puts several Fourier Operators into the global scope
"""
function create_operators(mod::Module, g::FourierGrid;
                          names = ("x","y","z"),
                          arraytype = Array)
    printstyled("Warning !!!", color = :red)
    println("")
    for i in 1:length(g.grid)
        operator_name = Char(0x02202) * names[i]
        print("Introducing ")
        printstyled(operator_name, color = :blue)
        println(" into the global scope.")
        k = g.wavenumbers[i]
        parsed_name = Meta.parse(operator_name)
        op = arraytype(im .* k)
        fmd = FourierOperatorMetaData(operator_name)
        Base.eval(mod, :($parsed_name = FourierOperator($op, $fmd)))
    end
    println(" ")
    return nothing
end

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: upiterbarg/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
# function test1_well_scaled_model()
#     scale_factors = 10.0.^(0.005:0.001:0.1)
#     orig_model = MetNets.ecoli_core_model()
#     orig_min_order, orig_max_order = MetNets.nzabs_range(orig_model.S)
#     obj_ider = MetNets.ECOLI_MODEL_BIOMASS_IDER

#     prog = Progress(length(scale_factors); desc = "Testing well scaling ...")
#     for scale_factor in scale_factors

#         scl_model = MetNets.well_scaled_model(orig_model, scale_factor;
#             verbose = true)
#         scl_min_order, scl_max_order = MetNets.nzabs_range(scl_model.S)

#         @test scl_min_order >= orig_min_order
#         @test scl_max_order <= orig_max_order

#         orig_fbaout = Chemostat.LP.fba(orig_model, obj_ider)
#         scl_fbaout = Chemostat.LP.fba(scl_model, obj_ider)

#         for rxn in orig_model.rxns
#             orig_av = MetNets.av(orig_model, orig_fbaout, rxn)
#             scl_av = MetNets.av(scl_model, scl_fbaout, rxn)

#             @test isapprox(orig_av, scl_av; atol = 1e-8)
#         end

#         next!(prog; showvalues = [
#                 ("scale_factor        ", scale_factor),
#                 ("                    ", ""),
#                 (" ORIGINAL MODEL     ", ""),
#                 ("model size:         ", size(orig_model)),
#                 ("nzabs_range:        ", (orig_min_order, orig_max_order)),
#                 ("obj_val:            ", MetNets.av(orig_model, orig_fbaout, obj_ider)),
#                 ("                    ", ""),
#                 (" SCALED MODEL       ", ""),
#                 ("model size:         ", size(scl_model)),
#                 ("nzabs_range:        ", (scl_min_order, scl_max_order)),
#                 ("obj_val:            ", MetNets.av(scl_model, scl_fbaout, obj_ider))
#             ]
#         )
#     end
#     finish!(prog)
# end
# test1_well_scaled_model()

# ## ------------------------------------------------------------------
# # This test just test the same obj_fun value but test all reactions as obj_fun
# function test2_well_scaled_model()

#     scale_factors = 10.0.^(0.005:0.005:0.1)
#     orig_model = MetNets.ecoli_core_model()
#     # Open all exchanges
#     exchs = MetNets.exchanges(orig_model)
#     MetNets.lb!.([orig_model], exchs, -1000.0)

#     orig_min_order, orig_max_order = MetNets.nzabs_range(orig_model.S)

#     N = length(scale_factors) * size(orig_model, 2)
#     prog = Progress(N; desc = "Testing well scaling ...")
#     for scale_factor in scale_factors

#         scl_model = MetNets.well_scaled_model(orig_model, scale_factor;
#             verbose = true)
#         scl_min_order, scl_max_order = MetNets.nzabs_range(scl_model.S)

#         @test scl_min_order >= orig_min_order
#         @test scl_max_order <= orig_max_order

#         for obj_ider in orig_model.rxns

#             orig_fbaout = Chemostat.LP.fba(orig_model, obj_ider)
#             scl_fbaout = Chemostat.LP.fba(scl_model, obj_ider)

#             orig_av = MetNets.av(orig_model, orig_fbaout, obj_ider)
#             scl_av = MetNets.av(scl_model, scl_fbaout, obj_ider)

#             @test isapprox(orig_av, scl_av; atol = 1e-8)

#             next!(prog; showvalues = [
#                     ("scale_factor        ", scale_factor),
#                     ("obj_ider            ", obj_ider),
#                     ("                    ", ""),
#                     (" ORIGINAL MODEL     ", ""),
#                     ("model size:         ", size(orig_model)),
#                     ("nzabs_range:        ", (orig_min_order, orig_max_order)),
#                     ("obj_val:            ", MetNets.av(orig_model, orig_fbaout, obj_ider)),
#                     ("                    ", ""),
#                     (" SCALED MODEL       ", ""),
#                     ("model size:         ", size(scl_model)),
#                     ("nzabs_range:        ", (scl_min_order, scl_max_order)),
#                     ("obj_val:            ", MetNets.av(scl_model, scl_fbaout, obj_ider))
#                 ]
#             )
#         end
#     end
#     finish!(prog)
# end
# test2_well_scaled_model()

ext: jl


text:
"""
Front-End is responsible for tracking keypoints
and computing poses for the Frames.
It also decides when the system needs a new Keyframe in the map.

# Parameters

- `current_frame::Frame`: Current frame that is being processed.
    This is a shared Frame between FrontEnd, MapManager, Estimator
    and SlamManager.
- `motion_model::MotionModel`: Motion model that is used to predict
    pose for the Frame before the actual pose for it was computed.
- `map_manager::MapManager`: Map manager that is responsible for the
    creation of new Keyframes in the map.
- `params::Params`: Parameters of the system.
- `current_image::Matrix{Gray{Float64}}`: Current image that is processed.
- `previous_image::Matrix{Gray{Float64}}`: Previous processed image.
- `current_pyramid::LKPyramid`: Pre-computed pyramid
    that is used for optical flow tracking for the `current_image`.
- `previous_pyramid::LKPyramid`: Pre-computed pyramid
    that is used for optical flow tracking for the `previous_image`.
"""
mutable struct FrontEnd
    current_frame::Frame
    motion_model::MotionModel
    map_manager::MapManager
    params::Params

    current_image::Matrix{Gray{Float64}}
    previous_image::Matrix{Gray{Float64}}

    current_pyramid::LKPyramid
    previous_pyramid::LKPyramid
    geev_cache::GEEV4x4Cache
end

function FrontEnd(params::Params, frame::Frame, map_manager::MapManager)
    empty_pyr = LKPyramid(
        [Matrix{Gray{Float64}}(undef, 0, 0)],
        nothing, nothing, nothing, nothing, nothing, nothing)
    FrontEnd(
        frame, MotionModel(), map_manager, params,
        Matrix{Gray{Float64}}(undef, 0, 0), Matrix{Gray{Float64}}(undef, 0, 0),
        empty_pyr, empty_pyr, GEEV4x4Cache())
end

"""
    track!(fe::FrontEnd, image, time)

Given an image and time at which it was taken, track keypoints in it.
After tracking, decide if the system needs a new Keyframe added to the map.
If it is the first image to be tracked, then Keyframe is always needed.

# Returns

`true` if the system needs a new Keyframe, otherwise `false`.
"""
function track!(fe::FrontEnd, image, time, slam_io)
    is_kf_required = false
    lock(fe.map_manager.map_lock)
    try
        is_kf_required = track_mono!(fe, image, time, slam_io)
        is_kf_required && create_keyframe!(fe.map_manager, image)
    catch e
        showerror(stdout, e)
        display(stacktrace(catch_backtrace()))
    finally
        unlock(fe.map_manager.map_lock)
    end

    @debug "Pose $(fe.current_frame.id) WC: $(fe.current_frame.wc[1:3, 4])"
    is_kf_required
end

function track_mono!(fe::FrontEnd, image, image_time, slam_io)
    preprocess!(fe, image)
    fe.current_frame.id == 1 && return true
    # Apply motion model & update current Frame pose.
    new_pose = fe.motion_model(fe.current_frame.wc, image_time)
    set_wc!(fe.current_frame, new_pose, slam_io)

    t1 = time()
    klt_tracking!(fe)
    t2 = time()
    @debug "[FE] KLT Tracking ($(t2 - t1) sec)."

    if !fe.params.vision_initialized
        if fe.current_frame.nb_keypoints < 50
            @warn "[Front End] NB KP < 50. Reset required."
            fe.params.reset_required = true
            return false
        elseif check_ready_for_init!(fe, slam_io)
            @debug "[Front End] System ready for initialization."
            fe.params.vision_initialized = true
            return true
        else
            @debug "[Front End] System not ready for initialization."
            return false
        end
    end

    # Epipolar filtering to remove outliers.
    # In case P3P fails this pose will be used.
    t1 = time()
    pose_5pt = compute_pose_5pt!(fe; min_parallax=5.0, use_motion_model=true)
    fe.map_manager.nb_keyframes > 2 && pose_5pt ≢ nothing &&
        set_cw!(fe.current_frame, pose_5pt, slam_io)
    t2 = time()
    @debug "[FE] 5PT Pose ($(t2 - t1) sec)."

    t1 = time()
    compute_pose!(fe, slam_io)
    t2 = time()
    @debug "[FE] Compute Pose ($(t2 - t1) sec)."

    update!(fe.motion_model, fe.current_frame.wc, image_time)
    check_new_kf_required(fe)
end

"""
    compute_pose!(fe::FrontEnd)

Compute pose of a current Frame using P3P Ransac algorithm.
Pose is computed from the triangulated Keypoints (MapPoints) that are visible
in this frame.

# Returns

`true` if the pose was successfully computed and applied to current Frame,
otherwise `false`.
"""
function compute_pose!(fe::FrontEnd, slam_io)
    if fe.current_frame.nb_3d_kpts < 5
        @warn "[Front-End] Not enough 3D keypoints to compute P3P $(fe.current_frame.nb_3d_kpts)."
        return false
    end

    p3p_pixels = Vector{Point2f}(undef, fe.current_frame.nb_3d_kpts)
    p3p_pdn_positions = Vector{Point3f}(undef, fe.current_frame.nb_3d_kpts)
    p3p_3d_points = Vector{Point3f}(undef, fe.current_frame.nb_3d_kpts)
    p3p_kpids = Vector{Int64}(undef, fe.current_frame.nb_3d_kpts)
    i = 1

    for kp in values(fe.current_frame.keypoints)
        kp.is_3d || continue
        mp = get(fe.map_manager.map_points, kp.id, nothing)
        mp ≡ nothing && continue

        p3p_pdn_positions[i] = normalize(kp.position)
        # Convert pixel to `(x, y)` format, expected by P3P.
        p3p_pixels[i] = kp.undistorted_pixel[[2, 1]]
        p3p_3d_points[i] = mp.position
        p3p_kpids[i] = kp.id
        i += 1
    end

    i -= 1
    p3p_pdn_positions = @view(p3p_pdn_positions[1:i])
    p3p_pixels = @view(p3p_pixels[1:i])
    p3p_3d_points = @view(p3p_3d_points[1:i])
    p3p_kpids = @view(p3p_kpids[1:i])

    # P3P computes world → camera projection.
    res = p3p_ransac(
        p3p_3d_points, p3p_pixels, p3p_pdn_positions,
        fe.current_frame.camera.K; threshold=fe.params.max_reprojection_error,
    )
    if res ≡ nothing
        @warn "[FE] P3P Ransac returned Nothing."
        fe |> reset_frame!
        return false
    end

    n_inliers, model = res
    if n_inliers < 5 || model ≡ nothing
        @warn "[FE] P3P too few inliers - resetting!"
        fe |> reset_frame!
        return false
    end

    KP, inliers, error = model
    set_cw!(fe.current_frame, to_4x4(fe.current_frame.camera.iK * KP), slam_io)
    # Remove outliers after P3P.
    for (kpid, inlier) in zip(p3p_kpids, inliers)
        inlier || remove_obs_from_current_frame!(fe.map_manager, kpid)
    end

    # Prepare data for BA refinement.
    p3p_3d_points = @view(p3p_3d_points[inliers])
    p3p_kpids = @view(p3p_kpids[inliers])

    c = 1
    for i in 1:length(p3p_pixels)
        inliers[i] || continue
        p = p3p_pixels[i]
        p3p_pixels[c] = Point2f(p[2], p[1])
        c += 1
    end
    p3p_pixels = @view(p3p_pixels[1:n_inliers])

    # Refinement.
    new_T, init_error, res_error, outliers, n_outliers = pnp_bundle_adjustment(
        fe.current_frame.camera, fe.current_frame.cw,
        p3p_pixels, p3p_3d_points;
        iterations=10, show_trace=false,
        repr_ϵ=fe.params.max_reprojection_error)
    if length(p3p_3d_points) - n_outliers < 5 || res_error > init_error
        @warn "[FE] P3P BA too few inliers - resetting!"
        fe |> reset_frame!
        return false
    end

    for (kpid, outlier) in zip(p3p_kpids, outliers)
        outlier && remove_obs_from_current_frame!(fe.map_manager, kpid)
    end

    set_cw!(fe.current_frame, new_T, slam_io)
    true
end

"""
    compute_pose_5pt!(fe::FrontEnd; min_parallax, use_motion_model)

Copmute pose for pixel correspondences using 5-point algorithm
to recover essential matrix, then pose from it.

# Arguments

- `min_parallax`: Minimum parallax required between pixels in the current Frame
    and previous Keyframe to compute pose.
    Note, that parallax is rotation-compensated, meaning a rotation from
    current to previous frame is computed and applied to every pixel.
- `use_motion_model`: If `true`, then use constant-velocity motion model
    to predict next pose from previous frame.
    Otherwise, the computed pose will be "local".

# Returns

If successfull, 4x4 pose matrix, that transforms points
from previous Keyframe to current Frame.
Otherwise `nothing`.
"""
function compute_pose_5pt!(fe::FrontEnd; min_parallax, use_motion_model)
    if fe.current_frame.nb_keypoints < 8
        @debug "[FE] Not enough keypoints for initialization: " *
            "$(fe.current_frame.nb_keypoints)"
        return nothing
    end

    previous_keyframe = get(
        fe.map_manager.frames_map, fe.current_frame.kfid, nothing)
    previous_keyframe ≡ nothing && return nothing
    R_compensation = get_Rcw(previous_keyframe) * get_Rwc(fe.current_frame)

    n_parallax = 0
    avg_parallax = 0.0

    previous_points = Vector{Point2f}(undef, fe.current_frame.nb_keypoints)
    current_points = Vector{Point2f}(undef, fe.current_frame.nb_keypoints)
    previous_pd = Vector{Point2f}(undef, fe.current_frame.nb_keypoints)
    current_pd = Vector{Point2f}(undef, fe.current_frame.nb_keypoints)
    kp_ids = Vector{Int64}(undef, fe.current_frame.nb_keypoints)
    i = 1

    # Get all Keypoint's positions and compute rotation-compensated parallax.
    @inbounds for keypoint in values(fe.current_frame.keypoints)
        pkf_keypoint = get(previous_keyframe.keypoints, keypoint.id, nothing)
        pkf_keypoint ≡ nothing && continue

        # Convert points to `(x, y)` format as expected by five points.
        previous_points[i] = pkf_keypoint.undistorted_pixel[[2, 1]]
        current_points[i] = keypoint.undistorted_pixel[[2, 1]]
        previous_pd[i] = pkf_keypoint.position[[1, 2]]
        current_pd[i] = keypoint.position[[1, 2]]
        kp_ids[i] = keypoint.id
        i += 1

        # Compute rotation-compensated parallax.
        rot_position = R_compensation * keypoint.position
        rot_px = project(fe.current_frame.camera, rot_position)
        avg_parallax += norm(rot_px - pkf_keypoint.undistorted_pixel)
        n_parallax += 1
    end
    if n_parallax < 8
        @warn "[FE] Not enough keypoints in previous KF " *
            "to compute 5pt Essential Matrix."
        return nothing
    end

    avg_parallax /= n_parallax
    if avg_parallax < min_parallax
        @warn "[FE] Not enough parallax ($avg_parallax) " *
            "to compute 5pt Essential Matrix."
        return nothing
    end

    i -= 1
    previous_points = @view(previous_points[1:i])
    current_points = @view(current_points[1:i])
    previous_pd = @view(previous_pd[1:i])
    current_pd = @view(current_pd[1:i])
    kp_ids = @view(kp_ids[1:i])

    # `P` is `cw`: transforms from world (previous frame) to current frame.
    n_inliers, (_, P, inliers, _) = five_point_ransac(
        previous_points, current_points, previous_pd, current_pd,
        fe.current_frame.camera.K, fe.current_frame.camera.K, fe.geev_cache;
        max_repr_error=fe.params.max_reprojection_error)
    if n_inliers < 5
        @warn "[FE] Not enough inliers ($n_inliers) for the " *
            "5pt Essential Matrix."
        return nothing
    end

    if n_inliers != n_parallax
        for (i, inlier) in enumerate(inliers)
            inlier || remove_obs_from_current_frame!(fe.map_manager, kp_ids[i])
        end
    end

    if use_motion_model
        # Get motion-model translation scale from last KeyFrame.
        prev_cw = get_cw(previous_keyframe)
        current = prev_cw * get_wc(fe.current_frame)
        scale = norm(current[1:3, 4])

        R, t = P[1:3, 1:3], P[1:3, 4]
        t = scale .* normalize(t)
        return to_4x4(R, t) * prev_cw
    end
    to_4x4(P) # cw pose
end

"""
    check_ready_for_init!(fe::FrontEnd)

Check if there is enough average rotation compensated parallax
between current Frame and previous KeyFrame.

Additionally, compute Essential matrix using 5-point Ransac algorithm
to filter out outliers and check if there is enough inliers to proceed.
"""
function check_ready_for_init!(fe::FrontEnd, slam_io)
    avg_parallax = compute_parallax(
        fe, fe.current_frame.kfid;
        compensate_rotation=false, median_parallax=false)
    @debug "[FE] Initial parallax $avg_parallax vs $(fe.params.initial_parallax)."
    avg_parallax ≤ fe.params.initial_parallax && return false
    pose = compute_pose_5pt!(
        fe; min_parallax=fe.params.initial_parallax, use_motion_model=false)
    ready = pose ≢ nothing
    ready && set_cw!(fe.current_frame, pose, slam_io)
    ready
end

"""
    check_new_kf_required(fe::FrontEnd)

Check if we need to insert a new KeyFrame into the Map.
"""
function check_new_kf_required(fe::FrontEnd)
    prev_kf = get(fe.map_manager.frames_map, fe.current_frame.kfid, nothing)
    prev_kf ≡ nothing && return false

    # Id difference since last KeyFrame.
    frames_δ = fe.current_frame.id - prev_kf.id
    if fe.current_frame.nb_occupied_cells < 0.33 * fe.params.max_nb_keypoints &&
        frames_δ ≥ 5 && !fe.params.local_ba_on
        return true
    end
    if fe.current_frame.nb_3d_kpts < 20 && frames_δ ≥ 2
        return true
    end
    if fe.current_frame.nb_3d_kpts > 0.5 * fe.params.max_nb_keypoints &&
        (fe.params.local_ba_on || frames_δ < 2)
        return false
    end

    # Time difference since last KeyFrame.
    time_δ = fe.current_frame.time - prev_kf.time
    # TODO option for stereo
    median_parallax = compute_parallax(
        fe, prev_kf.kfid; compensate_rotation=true, only_2d=false,
    )
    cx = median_parallax ≥ fe.params.initial_parallax / 2.0 # TODO || stereo
    c0 = median_parallax ≥ fe.params.initial_parallax
    c1 = fe.current_frame.nb_3d_kpts < 0.75 * prev_kf.nb_3d_kpts
    c2 = fe.current_frame.nb_occupied_cells < 0.5 * fe.params.max_nb_keypoints &&
        fe.current_frame.nb_3d_kpts < 0.85 * prev_kf.nb_3d_kpts &&
        !fe.params.local_ba_on

    cx && (c0 || c1 || c2)
end

"""
    compute_parallax(
        fe::FrontEnd, current_frame_id;
        compensate_rotation = true, only_2d = true, median_parallax = true)

Compute parallax in pixels between current Frame
and the provided `current_frame_id` Frame.

# Arguments

- `compensate_rotation::Bool`:
    Compensate rotation by computing relative rotation between
    current Frame and previous Keyframe if `true`. Default is `true`.
- `only_2d::Bool`: Consider only 2d keypoints. Default is `true`.
- `median_parallax::Bool`:
    Instead of the average, compute median parallax. Default is `true`.
"""
function compute_parallax(
    fe::FrontEnd, current_frame_id;
    compensate_rotation = true, only_2d = true, median_parallax = true,
)
    frame = get(fe.map_manager.frames_map, current_frame_id, nothing)
    if frame ≡ nothing
        @warn "[Front-End] Error in `compute_parallax`! " *
            "Keyframe $current_frame_id does not exist."
        return 0.0
    end

    current_rotation = compensate_rotation ?
        get_Rcw(frame) * get_Rwc(fe.current_frame) : SMatrix{3, 3, Float64, 9}(I)

    avg_parallax = 0.0
    n_parallax = 0
    parallax_set = Float64[]

    for keypoint in values(fe.current_frame.keypoints)
        only_2d && keypoint.is_3d && continue
        keypoint.id in keys(frame.keypoints) || continue

        upx = keypoint.undistorted_pixel
        compensate_rotation && (upx =
            project(frame.camera, current_rotation * keypoint.position);)

        parallax = norm(upx - get_keypoint_unpx(frame, keypoint.id))
        avg_parallax += parallax
        n_parallax += 1

        median_parallax && push!(parallax_set, parallax)
    end
    n_parallax == 0 && return 0.0

    if median_parallax
        avg_parallax = parallax_set |> median
    else
        avg_parallax /= n_parallax
    end
    avg_parallax
end

function preprocess!(fe::FrontEnd, image)
    t1 = time()
    fe.previous_image = fe.current_image
    fe.current_image = image

    if has_gradients(fe.previous_pyramid) && has_gradients(fe.current_pyramid)
        copy!(fe.previous_pyramid, fe.current_pyramid)
        update!(fe.current_pyramid, fe.current_image)
    else
        fe.previous_pyramid = fe.current_pyramid
        fe.current_pyramid = LKPyramid(
            fe.current_image, fe.params.pyramid_levels;
            σ=fe.params.pyramid_σ, reusable=true)
    end
    t2 = time()
    @debug "[FE] Preprocess ($(t2 - t1) sec)."
end

"""
    klt_tracking!(fe::FrontEnd)

Track keypoints from previous frame to current frame.
"""
function klt_tracking!(fe::FrontEnd)
    optical_flow_matching!(
        fe.map_manager, fe.current_frame,
        fe.previous_pyramid, fe.current_pyramid, false)
end

"""
    reset_frame!(fe::FrontEnd)

Reset current Frame in Front-End.
"""
function reset_frame!(fe::FrontEnd)
    for kpid in keys(fe.current_frame.keypoints)
        remove_obs_from_current_frame!(fe.map_manager, kpid)
    end
    fe.current_frame.keypoints |> empty!
    fe.current_frame.keypoints_grid .|> empty!

    fe.current_frame.nb_2d_kpts = 0
    fe.current_frame.nb_3d_kpts = 0
    fe.current_frame.nb_keypoints = 0
    fe.current_frame.nb_occupied_cells = 0
end

"""
    reset!(fe::FrontEnd)

Reset Front-End.
"""
function reset!(fe::FrontEnd)
    empty_pyr = LKPyramid(
        [Matrix{Gray{Float64}}(undef, 0, 0)],
        nothing, nothing, nothing, nothing, nothing, nothing)
    fe.previous_pyramid = empty_pyr
    fe.current_pyramid = empty_pyr
end


text:
# 2D FWI gradient test with 1 source
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#
# Mathias Louboutin, mlouboutin3@gatech.edu
# Updated July 2020

### Model
model, model0, dm = setup_model(tti, viscoacoustic, 4)
q, srcGeometry, recGeometry, f0 = setup_geom(model)
dt = srcGeometry.dt[1]

###################################################################################################

@testset "FWI gradient test with $(nlayer) layers and tti $(tti) and viscoacoustic $(viscoacoustic) and freesurface $(fs)" begin
	# Gradient test
	h = 5f-2
	maxiter = 6
	err1 = zeros(maxiter)
	err2 = zeros(maxiter)
	h_all = zeros(maxiter)
	modelH = deepcopy(model0)

	# Observed data
	opt = Options(sum_padding=true, free_surface=fs, f0=f0)
	F = judiModeling(model, srcGeometry, recGeometry; options=opt)
	d = F*q

	# FWI gradient and function value for m0
	Jm0, grad = fwi_objective(model0, q, d; options=opt)
	# Check get same misfit as l2 misifit on forward data
	Jm01 = .5f0 * norm(F(model0)*q - d)^2
	@test Jm0 ≈ Jm01
	dJ = dot(grad, dm)

	for j=1:maxiter
		# FWI function value for m0 + h*dm
		modelH.m = model0.m + h*dm
		Jm = .5f0 * norm(F(modelH)*q - d)^2

		# Check convergence
		err1[j] = abs(Jm - Jm0)
		err2[j] = abs(Jm - (Jm0 + h*dJ))
		j == 1 ? prev = 1 : prev = j - 1
		@printf("h = %2.2e, e1 = %2.2e, rate = %2.2e", h, err1[j], err1[prev]/err1[j])
		@printf(", e2  = %2.2e, rate = %2.2e \n", err2[j], err2[prev]/err2[j])
		h_all[j] = h
		h = h * .8f0
	end

	# CHeck convergence rates
	rate_1 = sum(err1[1:end-1]./err1[2:end])/(maxiter - 1)
	rate_2 = sum(err2[1:end-1]./err2[2:end])/(maxiter - 1)

	@test isapprox(rate_1, 1.25f0; rtol=5f-2)
	@test isapprox(rate_2, 1.5625f0; rtol=5f-2)
end


text:
using LinearAlgebra
using Symbolics

# r=1 in simulations always
struct SphericalCoord
    theta::Float64 # angle from z-axis
    phi::Float64 # angle from x-axis in x-y plane
end

Base.:+(x::SphericalCoord, y::SphericalCoord) = SphericalCoord(x.theta+y.theta, x.phi+y.phi)

struct CartesianCoord
    x::Float64
    y::Float64
    z::Float64
end

LinearAlgebra.normalize( c::CartesianCoord ) = begin
    n = normalize( [c.x, c.y, c.z ] )
    CartesianCoord( n[1], n[2], n[3])
end

LinearAlgebra.norm( c::CartesianCoord ) = norm( [c.x, c.y, c.z ] )

function cartesianToSpherical( p::CartesianCoord )
    r = sqrt( p.x^2 + p.y^2 + p.z^2 )
    theta = acos(p.z/r)
    phi = atan(p.y,p.x)
    SphericalCoord( theta, phi )
end

function sphericalToCartesian( p::SphericalCoord )
    x = cos(p.phi)*sin(p.theta)
    y = sin(p.phi)*sin(p.theta)
    z = cos(p.theta)
    CartesianCoord( x, y, z )
end

# p is the position, v is the velocity
function sphericalToCartesianVelocity( p::SphericalCoord, v::SphericalCoord )
    # p.theta is dtheta, p.phi is dphi
    vx = [cos(p.phi)*cos(p.theta), - sin(p.phi)]' * [v.theta, v.phi]
    vy = [sin(p.phi)*cos(p.theta), cos(p.phi)]' * [v.theta, v.phi]
    vz = -sin(p.theta)*v.theta
    CartesianCoord( vx, vy, vz )
end

# geodesic distance
rho( p1::SphericalCoord, p2 ::SphericalCoord) = acos( cos(p1.theta)*cos(p2.theta) + sin(p1.theta)*sin(p2.theta)*cos(p1.phi - p2.phi) )
# derivative of geodesic distance \rho wrt \phi
alpha( p1::SphericalCoord, p2 ::SphericalCoord) = sin(p1.theta) * sin(p2.theta) * sin(p1.phi-p2.phi) / sin(rho(p1,p2))
# derivative of geodesic distance \rho wrt \theta
beta( p1::SphericalCoord, p2 ::SphericalCoord) = ( cos(p2.theta) * sin(p1.theta) -  cos(p1.theta) * sin(p2.theta) * cos(p1.phi-p2.phi) ) / sin(rho(p1,p2))


text:
using ModelingToolkit, LinearAlgebra, DiffEqJump, Test, Distributions
MT = ModelingToolkit

@parameters t k[1:20]
@variables A(t) B(t) C(t) D(t)
rxs = [Reaction(k[1], nothing, [A]),            # 0 -> A
       Reaction(k[2], [B], nothing),            # B -> 0
       Reaction(k[3],[A],[C]),                  # A -> C
       Reaction(k[4], [C], [A,B]),              # C -> A + B
       Reaction(k[5], [C], [A], [1], [2]),      # C -> A + A
       Reaction(k[6], [A,B], [C]),              # A + B -> C
       Reaction(k[7], [B], [A], [2], [1]),      # 2B -> A
       Reaction(k[8], [A,B], [A,C]),            # A + B -> A + C
       Reaction(k[9], [A,B], [C,D]),            # A + B -> C + D
       Reaction(k[10], [A], [C,D], [2], [1,1]), # 2A -> C + D
       Reaction(k[11], [A], [A,B], [2], [1,1]), # 2A -> A + B
       Reaction(k[12], [A,B], nothing, [3,1], nothing),           # 3A+B -> 0
       Reaction(k[13], [A,B,C], [C,D], [1,3,4], [Poisson(2), 3]),          # A+3B+4C -> k_pois * C + 3D
       Reaction(k[14], nothing, [A], nothing, [Geometric(1 // 3)]),               # 0 -> k_geom * A
       Reaction(k[15]*A/(2+A), [A], nothing; only_use_rate=true), # A -> 0 with custom rate
       Reaction(k[16], [A], [B]; only_use_rate=true),             # A -> B with custom rate.
       Reaction(k[17]*A*exp(B), [C], [D], [2], [Bernoulli(1)]),   # 2C -> D with non constant rate.
       Reaction(k[18]*B, nothing, [B], nothing, [2]),             # 0 -> 2B with non constant rate.
       Reaction(k[19]*t, [A], [B]),                                # A -> B with non constant rate.
       Reaction(k[20]*t*A, [B,C], [D],[2,1],[2])                  # 2A +B -> 2C with non constant rate.
  ]
rs = ReactionSystem(rxs,t,[A,B,C,D],k)
odesys = convert(ODESystem,rs)
sdesys = convert(SDESystem,rs)

# test show
io = IOBuffer()
show(io, rs)
str = String(take!(io))
@test count(isequal('\n'), str) < 30

# hard coded ODE rhs
function oderhs(u,k,t)
       A = u[1]; B = u[2]; C = u[3]; D = u[4];
       du = zeros(eltype(u),4)
       du[1] = k[1] - k[3]*A + k[4]*C + 2*k[5]*C - k[6]*A*B + k[7]*B^2/2 - k[9]*A*B - k[10]*A^2 - k[11]*A^2/2 - k[13]*A*B^3*C^4/144 - 3*k[12]*A^3*B/6 + 2*k[14] - k[15]*A/(2+A) - k[16] - k[19]*t*A
       du[2] = -k[2]*B + k[4]*C - k[6]*A*B - k[7]*B^2 - k[8]*A*B - k[9]*A*B + k[11]*A^2/2 - 3*k[13]*A*B^3*C^4/144 - k[12]*A^3*B/6 + k[16] + 2*k[18]*B + k[19]*t*A - 2*k[20]*t*A*B^2*C
       du[3] = k[3]*A - k[4]*C - k[5]*C  + k[6]*A*B + k[8]*A*B + k[9]*A*B + k[10]*A^2/2 - 2*k[13]*A*B^3*C^4/144 - 2*k[17]*A*exp(B)*C^2/2 - k[20]*t*A*B^2*C
       du[4] = k[9]*A*B + k[10]*A^2/2 + 3*k[13]*A*B^3*C^4/144 + k[17]*A*exp(B)*C^2/2 + 2*k[20]*t*A*B^2*C
       du
end

# sde noise coefs
function sdenoise(u,k,t)
       A = u[1]; B = u[2]; C = u[3]; D = u[4];
       G = zeros(eltype(u),length(k),length(u))
       z = zero(eltype(u))

       G = [sqrt(k[1]) z  z z;
            z -sqrt(k[2]*B) z z;
            -sqrt(k[3]*A) z sqrt(k[3]*A) z;
            sqrt(k[4]*C) sqrt(k[4]*C) -sqrt(k[4]*C) z;
            2*sqrt(k[5]*C) z -sqrt(k[5]*C) z;
            -sqrt(k[6]*A*B) -sqrt(k[6]*A*B) sqrt(k[6]*A*B) z;
            sqrt(k[7]*B^2/2) -2*sqrt(k[7]*B^2/2) z z;
            z -sqrt(k[8]*A*B) sqrt(k[8]*A*B) z;
            -sqrt(k[9]*A*B) -sqrt(k[9]*A*B) sqrt(k[9]*A*B) sqrt(k[9]*A*B);
            -2*sqrt(k[10]*A^2/2) z sqrt(k[10]*A^2/2) sqrt(k[10]*A^2/2);
            -sqrt(k[11]*A^2/2) sqrt(k[11]*A^2/2) z z;
            -3*sqrt(k[12]*A^3*B/6) -sqrt(k[12]*A^3*B/6) z z;
            -sqrt(k[13]*A*B^3*C^4/144) -3*sqrt(k[13]*A*B^3*C^4/144) -2*sqrt(k[13]*A*B^3*C^4/144) 3*sqrt(k[13]*A*B^3*C^4/144);
            2*sqrt(k[14]) z z z;
            -sqrt(k[15]*A/(2+A)) z z z;
            -sqrt(k[16]) sqrt(k[16]) z z;
            z z -2*sqrt(k[17]*A*exp(B)*C^2/2) sqrt(k[17]*A*exp(B)*C^2/2);
            z 2*sqrt(k[18]*B) z z;
            -sqrt(k[19]*t*A) sqrt(k[19]*t*A) z z;
            z -2*sqrt(k[20]*t*A*B^2*C) -sqrt(k[20]*t*A*B^2*C) +2*sqrt(k[20]*t*A*B^2*C);
            # burst variances
            z z sqrt(k[13]*A*B^3*C^4/144) * sqrt(2.) z;
            sqrt(k[14])*sqrt(6.) z z z;
            z z z z]';

       return G
end

# test by evaluating drift terms
p  = rand(length(k))
u  = rand(length(k))
t  = 0.
du = oderhs(u,p,t)
odesys = convert(ODESystem,rs)
of = ODEFunction{false}(odesys, states(rs), parameters(rs))
du2 = of.f(u,p,t)
@test norm(du-du2) < 100*eps()

# test by evaluating drift and diffusion terms
p  = rand(length(k))
u  = rand(length(k))
t  = 0.
du = oderhs(u,p,t)
G  = sdenoise(u,p,t)
sdesys = convert(SDESystem,rs)
sf = SDEFunction{false}(sdesys, states(rs), parameters(rs))
du2 = sf.f(u,p,t)
@test norm(du-du2) < 100*eps()
G2 = sf.g(u,p,t)
@test norm(G-G2) < 100*eps()

# test conversion to NonlinearSystem
ns = convert(NonlinearSystem,rs)
fnl = eval(generate_function(ns)[2])
dunl = similar(du)
fnl(dunl,u,p)
@test norm(du-dunl) < 100*eps()

# tests the noise_scaling argument.
p  = rand(length(k)+1)
u  = rand(length(k))
t  = 0.
G  = p[21]*sdenoise(u,p,t)
@variables η
sdesys_noise_scaling = convert(SDESystem,rs;noise_scaling=η)
sf = SDEFunction{false}(sdesys_noise_scaling, states(rs), parameters(sdesys_noise_scaling))
G2 = sf.g(u,p,t)
@test norm(G-G2) < 100*eps()

# tests the noise_scaling vector argument.
p  = rand(length(k)+3)
u  = rand(length(k))
t  = 0.
G  = vcat(fill(p[21],8),fill(p[22],5),fill(p[23],7),[p[22],p[23],p[23]])' .* sdenoise(u,p,t)
@variables η[1:3]
sdesys_noise_scaling = convert(SDESystem,rs;noise_scaling=vcat(fill(η[1],8),fill(η[2],5),fill(η[3],7)))
sf = SDEFunction{false}(sdesys_noise_scaling, states(rs), parameters(sdesys_noise_scaling))
G2 = sf.g(u,p,t)
@test norm(G-G2) < 100*eps()

# tests using previous parameter for noise scaling
p  = rand(length(k)+3)
p[21], p[22], p[23] = p[13], p[14], p[17]
u  = rand(length(k))
t  = 0.
G  = [p p p p]' .* sdenoise(u,p,t)
sdesys_noise_scaling = convert(SDESystem,rs;noise_scaling=vcat(k))
sf = SDEFunction{false}(sdesys_noise_scaling, states(rs), parameters(sdesys_noise_scaling))
G2 = sf.g(u,p,t)
@test norm(G-G2) < 100*eps()

# test with JumpSystem
js = convert(JumpSystem, rs)

midxs = 1:12
cidxs = 13:18
vidxs = 19:20
@test all(map(i -> typeof(equations(js)[i]) <: DiffEqJump.MassActionJump, midxs))
@test all(map(i -> typeof(equations(js)[i]) <: DiffEqJump.ConstantRateJump, cidxs))
@test all(map(i -> typeof(equations(js)[i]) <: DiffEqJump.VariableRateJump, vidxs))

pars = rand(length(k)); u0 = rand(1:10,4); ttt = rand();
jumps = Vector{Union{ConstantRateJump, MassActionJump, VariableRateJump}}(undef,length(rxs))

jumps[1] = MassActionJump(pars[1], Vector{Pair{Int,Int}}(), [1 => 1]);
jumps[2] = MassActionJump(pars[2], [2 => 1], [2 => -1]);
jumps[3] = MassActionJump(pars[3], [1 => 1], [1 => -1, 3 => 1]);
jumps[4] = MassActionJump(pars[4], [3 => 1], [1 => 1, 2 => 1, 3 => -1]);
jumps[5] = MassActionJump(pars[5], [3 => 1], [1 => 2, 3 => -1]);
jumps[6] = MassActionJump(pars[6], [1 => 1, 2 => 1], [1 => -1, 2 => -1, 3 => 1]);
jumps[7] = MassActionJump(pars[7], [2 => 2], [1 => 1, 2 => -2]);
jumps[8] = MassActionJump(pars[8], [1 => 1, 2 => 1], [2 => -1, 3 => 1]);
jumps[9] = MassActionJump(pars[9], [1 => 1, 2 => 1], [1 => -1, 2 => -1, 3 => 1, 4 => 1]);
jumps[10] = MassActionJump(pars[10], [1 => 2], [1 => -2, 3 => 1, 4 => 1]);
jumps[11] = MassActionJump(pars[11], [1 => 2], [1 => -1, 2 => 1]);
jumps[12] = MassActionJump(pars[12], [1 => 3, 2 => 1], [1 => -3, 2 => -1]);

jumps[13] = ConstantRateJump((u,p,t) -> p[13]*u[1]*binomial(u[2],3)*binomial(u[3],4),
                              integrator -> (integrator.u[1] -= 1; integrator.u[2] -= 3; integrator.u[3] += 2 - 4; integrator.u[4]+=3));
jumps[14] = ConstantRateJump((u,p,t) -> p[14], integrator -> (integrator.u[1] += 1));

jumps[15] = ConstantRateJump((u,p,t) -> p[15]*u[1]/(2+u[1]), integrator -> (integrator.u[1] -= 1))
jumps[16] = ConstantRateJump((u,p,t) -> p[16], integrator -> (integrator.u[1] -= 1; integrator.u[2] += 1;))
jumps[17] = ConstantRateJump((u,p,t) -> p[17]*u[1]*exp(u[2])*binomial(u[3],2), integrator -> (integrator.u[3] -= 2; integrator.u[4] += 1))
jumps[18] = ConstantRateJump((u,p,t) -> p[18]*u[2], integrator -> (integrator.u[2] += 2))

jumps[19] = VariableRateJump((u,p,t) -> p[19]*u[1]*t, integrator -> (integrator.u[1] -= 1; integrator.u[2] += 1))
jumps[20] = VariableRateJump((u,p,t) -> p[20]*t*u[1]*binomial(u[2],2)*u[3], integrator -> (integrator.u[2] -= 2; integrator.u[3] -= 1; integrator.u[4] += 2))

statetoid = Dict(state => i for (i,state) in enumerate(states(js)))
parammap = map((x,y)->Pair(x,y),parameters(js),pars)
for i in midxs
  maj = MT.assemble_maj(equations(js)[i], statetoid, ModelingToolkit.substituter(parammap),eltype(pars))
  @test abs(jumps[i].scaled_rates - maj.scaled_rates) < 100*eps()
  @test jumps[i].reactant_stoch == maj.reactant_stoch
  @test jumps[i].net_stoch == maj.net_stoch
end
for i in cidxs
  crj = MT.assemble_crj(js, equations(js)[i], statetoid)
  @test isapprox(crj.rate(u0,p,ttt), jumps[i].rate(u0,p,ttt))
  fake_integrator1 = (u=zeros(4),p=p,t=0); fake_integrator2 = deepcopy(fake_integrator1);
  crj.affect!(fake_integrator1);
  jumps[i].affect!(fake_integrator2);
  # Won't have equality in general due to randomness in bursts
  #@test fake_integrator1 == fake_integrator2
end
for i in vidxs
  crj = MT.assemble_vrj(js, equations(js)[i], statetoid)
  @test isapprox(crj.rate(u0,p,ttt), jumps[i].rate(u0,p,ttt))
  fake_integrator1 = (u=zeros(4),p=p,t=0.); fake_integrator2 = deepcopy(fake_integrator1);
  crj.affect!(fake_integrator1); jumps[i].affect!(fake_integrator2);
  @test fake_integrator1 == fake_integrator2
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dashbio_sequenceviewer

"""
    dashbio_sequenceviewer(;kwargs...)

A SequenceViewer component.
The sequence viewer component is used to display sequences
that represent proteins, strands of genomic information, and
more. It can apply a coverage to the sequence supplied (with
clickable coverage sections that can display specific information,
and an optional legend to describe the color codes used),
search through the sequence for specific regex, capture
mouse selection events of subparts of the sequence, display a
count of the number of nucleotides or amino acids in the
sequence,
Read more about the component here:
https://github.com/FlyBase/react-sequence-viewer
Keyword arguments:
- `id` (String; optional): The ID used to identify this component in Dash callbacks.
- `badge` (Bool; optional): The option of whether or not to display a badge showing the
amino acid count at the top of the component beside the title.
- `charsPerLine` (Real; optional): The number of amino acids that will display per line.
- `coverage` (optional): A coverage of the entire sequence; each section of the sequence
can have its own text color, background color, tooltip (on hover),
and an optional underscore. The props start and end represent the
beginning and terminating indices of the section in question.
Cannot be used at the same time as selection.. coverage has the following type: Array of lists containing elements 'start', 'end', 'color', 'bgcolor', 'tooltip', 'underscore', 'onclick'.
Those elements have the following types:
  - `start` (Real; optional)
  - `end` (Real; optional)
  - `color` (String; optional)
  - `bgcolor` (String; optional)
  - `tooltip` (String; optional)
  - `underscore` (Bool; optional)
  - `onclick` (optional)s
- `coverageClicked` (Real; optional): Contains the index of the section that was clicked last in
the coverage list supplied.
- `legend` (optional): A legend corresponding to the color codes above (optionally displayed).. legend has the following type: Array of lists containing elements 'name', 'color', 'underscore'.
Those elements have the following types:
  - `name` (String; optional)
  - `color` (String; optional)
  - `underscore` (Bool; optional)s
- `mouseSelection` (optional): Contains information about the subsequence selected
by the mouse. Start and end refer to the initial and
final indices, respectively, of the subsequence, and
"selection" contains the string that is selected.. mouseSelection has the following type: lists containing elements 'start', 'end', 'selection'.
Those elements have the following types:
  - `start` (Real; optional)
  - `end` (Real; optional)
  - `selection` (String; optional)
- `search` (Bool; optional): The option of whether or not to include a search bar in
the header. This supports regex.
- `selection` (optional): A highlighted section of the sequence; the color of the highlight
can also be defined. Takes a list of format [min, max, color] where
min is a number that represents the starting index of the selection,
max is a number that represents the stopping index of the selection,
and color is a string that defines the highlight color.
Cannot be used at the same time as coverage.
- `sequence` (String; optional): The amino acid sequence that will be displayed.
- `sequenceMaxHeight` (String; optional): The maximum height of the sequence.
- `showLineNumbers` (Bool; optional): The option of whether or not to display line numbers.
- `subpartSelected` (optional): A list of the subparts selected using the
"search" function or the "selection" property.. subpartSelected has the following type: Array of lists containing elements 'start', 'end', 'sequence'.
Those elements have the following types:
  - `start` (Real; optional)
  - `end` (Real; optional)
  - `sequence` (String; optional)s
- `title` (String; optional): A string that displays at the top of the component.
- `toolbar` (Bool; optional): The option of whether or not to display a toolbar at the top
that allows the user to choose the number of letters per line.
- `wrapAminoAcids` (Bool; optional): The option of whether or not to display the list of amino acids
as broken up into separate lines of a fixed length set by
charsPerLine.
"""
function dashbio_sequenceviewer(; kwargs...)
        available_props = Symbol[:id, :badge, :charsPerLine, :coverage, :coverageClicked, :legend, :mouseSelection, :search, :selection, :sequence, :sequenceMaxHeight, :showLineNumbers, :subpartSelected, :title, :toolbar, :wrapAminoAcids]
        wild_props = Symbol[]
        return Component("dashbio_sequenceviewer", "SequenceViewer", "dash_bio", available_props, wild_props; kwargs...)
end


text:
export SyncFileObservable, file

import Base: ==
import Base: show

"""
    SyncFileObservable(path::String)

File observable, which synchronously emits content of the file line by line as a `String` objects on subscription.

See also: [`file`](@ref), [`Subscribable`](@ref)
"""
struct SyncFileObservable <: Subscribable{String}
    path :: String
end

function on_subscribe!(observable::SyncFileObservable, actor)
    f = open(observable.path, "r")
    for line in eachline(f)
        next!(actor, line)
    end
    complete!(actor)
    close(f)
    return voidTeardown
end

"""
    file(path::String)

Creation operator for the `SyncFileObservable` with a given path.

See also: [`SyncFileObservable`](@ref)
"""
file(path::String) = SyncFileObservable(path)

Base.:(==)(f1::SyncFileObservable, f2::SyncFileObservable) = f1.path == f2.path

Base.show(io::IO, observable::SyncFileObservable) = print(io, "SyncFileObservable()")


text:
using Compat.Test

include("testexamples.jl")
include("mathprogtest.jl")


text:
const TeamViewPokemonBoxRects = [
    Rect(
        (29 + 280 * (y - 1), 121 + 880 * (x - 1)),
        (264, 800),
        resolution = DesignResolution
    )
    for x in 1:2 for y in 1:3
]
const TeamViewLineHeight = 50
const TeamViewNameRect = Rect((108, 1), (TeamViewLineHeight, 410), resolution = DesignResolution)
const TeamViewAbilityRect = Rect((158, 1), (TeamViewLineHeight, 410), resolution = DesignResolution)
const TeamViewItemRect = Rect((208, 1), (TeamViewLineHeight, 410), resolution = DesignResolution)
const TeamViewMoveRects = [
    Rect((13 + 64 * (i - 1), 476), (TeamViewLineHeight, 320), resolution = DesignResolution)
    for i in 1:4
]
const TeamViewAuthorRect = Rect((877, 512), (64, 334), resolution = DesignResolution)
const TeamViewTitleRect = Rect((877, 875), (64, 535), resolution = DesignResolution)
const TeamViewRentalCodeRect = Rect((941, 820), (64, 560), resolution = DesignResolution)
const TeamViewMoveTypeRects = [
    Rect((5 + 64 * (i - 1), 414), (64, 64), resolution = DesignResolution)
    for i in 1:4
]
const TypeIconSheetRect = Rect((8, 8), (49, 49))

struct TeamView <: AbstractPokemonScene
    team::Team
end

feature_image_name(::Type{TeamView}, ::PokemonBattle) = "team-view"

should_include(::Type{TeamView}, ctx::PokemonContext) =
    !isnothing(ctx.data.context)  # should not be included when in a battle

function _get_resized_gender_icons(ctx, height)
    data = ctx.data
    icons, icon_masks = if !isnothing(data.resized_gender_icons)
        data.resized_gender_icons, data.resized_gender_icon_masks
    else
        sheet = gender_icons()
        icons = [block(sheet, i) for i in 2:3]
        icon_masks = [block(sheet.mask, i, sheet.table_size, sheet.block_size) for i in 2:3]
        icons, icon_masks
    end
    h, w = size(icons[1])  # (50, 50)
    if h != height
        new_size = height, w * height ÷ h
        icons = [imresize(icon, new_size) for icon in icons]
        icon_masks = [imresize(mask, new_size) for mask in icon_masks]
        data.resized_gender_icons = icons
        data.resized_gender_icon_masks = icon_masks
    end
    icons, icon_masks
end

function parse_gender(line, ctx::PokemonContext)
    line_h, line_w = size(line)
    icons, icon_masks = _get_resized_gender_icons(ctx, line_h)
    ch = cld(line_h, 2)
    indices = (ch:ch, 1:line_w)
    (female_pos, female_dist), (male_pos, male_dist) = (
        template_match(
            icon, line;
            indices = indices,
            σ = 0,
            mask = mask
        )
        for (icon, mask) in zip(icons, icon_masks)
    )
    r = female_dist / male_dist  # Use the ratio to determine if one is dominant.
    w = size(icons[1], 2)
    if r > 10
        GENDER_MALE, male_pos ± w ÷ 2
    elseif r < 0.1
        GENDER_FEMALE, female_pos ± w ÷ 2
    else
        GENDER_NULL, 0:0
    end
end

function parse_name_level(line, gender_pos, ctx::PokemonContext)
    plist = poke_list(ctx)
    pid, level = if length(gender_pos) > 0  # has gender
        len = length(gender_pos)
        gender_start, gender_end = first(gender_pos) + len ÷ 10, last(gender_pos) - len ÷ 10
        name_box = prepare_text_for_ocr(@view line[:, 1:gender_start - 1])
        level_box = prepare_text_for_ocr(@view line[:, gender_end + 1:end])
        # TODO: Check if these are needed
        # pids = data_search_n(plist, name; n = 5)
        # dist1, id1 = pids[1]
        # dist1 < 2 && return id1, 50
        name = ocr(name_box, ctx)
        level = ocr(Int, level_box, ctx; language = "eng")
        data_search(plist, name), level
    else
        name_box = prepare_text_for_ocr(line)
        name = ocr(name_box, ctx)
        m = match(r"^(.+)Lv(.+)$", name)
        name, level_string = if isnothing(m)
            name, name
        else
            m.captures[1], m.captures[2]
        end
        data_search(plist, name), parse_int(level_string)
    end
    PokemonID(pid), level
end

function parse_name_line(box, ctx::PokemonContext)
    # id, gender, level
    line = box[TeamViewNameRect]
    gender, gender_pos = parse_gender(line, ctx)
    id, level = parse_name_level(line, gender_pos, ctx)
    id, gender, level
end

function _parse_simple_line(::Type{T}, line, dex, ctx) where T
    line = prepare_text_for_ocr(line)
    text = ocr(line, ctx)
    data_search(dex, text) |> T
end

function parse_ability_line(box, ctx::PokemonContext; poke_id = missing)
    alist = ability_list(ctx)
    if !ismissing(poke_id)
        dex = poke_dex()
        dex_poke = dex[poke_id]
        abilities = dex_poke.abilities
        length(abilities) == 1 && return abilities[1]
        alist = Dict{AbilityID, String}(
            ability => alist[ability]
            for ability in abilities
        )
    end
    _parse_simple_line(AbilityID, box[TeamViewAbilityRect], alist, ctx)
end

parse_item_line(box, ctx::PokemonContext) =
    _parse_simple_line(ItemID, box[TeamViewItemRect], item_list(ctx), ctx)

function get_move_with_type(moves, box, i, ctx::PokemonContext)
    moves = MoveID.(moves)
    icons = type_icons()
    icon = box[TeamViewMoveTypeRects[i]]
    t = table_search(
        icon, icons,
        rect = TypeIconSheetRect
    )
    t = enum_from_string(t, PokemonType)
    mdex = move_dex()
    for move in moves
        mdex[move].type == t && return move
    end
    moves[1]
end

function parse_move_line(box, i, ctx::PokemonContext)
    line = box[TeamViewMoveRects[i]]
    line = Gray.(line)
    # TODO: check if this threshold is OK
    Float32(stdmult(⊙, line)) < 0.05f0 && return missing
    mlist = move_list(ctx)
    text = ocr(line, ctx)
    dists = data_search_n(mlist, text, n = 5)
    dists = [x for x in dists if x[1] == dists[1][1]]
    move = if length(dists) == 1
        MoveID(dists[1][2])
    else
        moves = [x[2] for x in dists]
        get_move_with_type(moves, box, i, ctx::PokemonContext)
    end
end

function parse_moves(box, ctx::PokemonContext)
    moves = MoveID[]
    for i in 1:4
        move = parse_move_line(box, i, ctx)
        ismissing(move) && break
        push!(moves, move)
    end
    moves
end

function parse_pokemon_box(box, ctx::PokemonContext)
    poke = Pokemon()
    poke.id, poke.gender, poke.level = parse_name_line(box, ctx)
    poke.ability = parse_ability_line(box, ctx; poke_id = poke.id)
    poke.item = parse_item_line(box, ctx)
    poke.moves = parse_moves(box, ctx)
    poke
end

parse_pokemon_box(img, i, ctx::PokemonContext) = parse_pokemon_box(img[TeamViewPokemonBoxRects[i]], ctx)

parse_author(img, ctx::PokemonContext) = ocr_multiple_lang(prepare_text_for_ocr(img[TeamViewAuthorRect]), ctx)[1]
parse_title(img, ctx::PokemonContext) = ocr_multiple_lang(prepare_text_for_ocr(img[TeamViewTitleRect]), ctx)[1]
function parse_rental_code(img, ctx::PokemonContext)
    text = ocr(img[TeamViewRentalCodeRect], ctx; language = "eng") |> remove_spaces
    if length(text) == 14
        text = "$(text[1:4]) $(text[5:8]) $(text[9:12]) $(text[13:14])"
    end
    text
end

function _parse_scene(::Type{TeamView}, frame::VsFrame, ctx::PokemonContext)
    img = image(frame)
    team = Team()
    for i in 1:6
        push!(team.pokemons, parse_pokemon_box(img, i, ctx))
    end
    team.author = parse_author(img, ctx)
    team.title = parse_title(img, ctx)
    team.rental_code = parse_rental_code(img, ctx)
    TeamView(team)
end


text:
using RecyclePermanent
using BenchmarkTools
using JLD2
include("../models/gkModel.jl")

@load "data/gkData.jld" ys θ
ys = ys[1:100]
model = makeGKModel(θ, ys)

N = 80 * model.n

rpio = RPIO{model.Tx}(N, model.n, Threads.nthreads())
@btime simpleEstimate($model, $rpio)
@btime biasedEstimate($model, $rpio)
@btime rpEstimate($model, $rpio)

rpio = RPIO{model.Tx}(N, model.n, 1)
@btime simpleEstimate($model, $rpio)
@btime biasedEstimate($model, $rpio)
@btime rpEstimate($model, $rpio)


text:
splitobs(lm::Dict{T,Vector{I}}; at = 0.7) where {T,I<:Integer} =
    splitobs(lm, at)

function splitobs(lm::Dict{T,Vector{I}}, at::AbstractFloat) where {T,I<:Integer}
    0 < at < 1 || throw(ArgumentError("the parameter \"at\" must be in interval (0, 1)"))
    n = mapreduce(length, +, values(lm))
    k = nlabel(lm)
    # preallocate the indices vectors
    idx1 = Vector{I}()
    idx2 = Vector{I}()
    # sizehint will save us a few heavy memory allocations
    # we specify "+ k" to deal with trailing observations when
    # the number of observations from a class isn't divideable
    # by "at" or "1-at"
    sizehint!(idx1, ceil(Int, n * at     + k))
    sizehint!(idx2, ceil(Int, n * (1-at) + k))
    # loop through all label indices
    for indices in values(lm)
        i1, i2 = splitobs(indices, at)
        append!(idx1, i1)
        append!(idx2, i2)
    end
    idx1, idx2
end

@generated function splitobs(lm::Dict{T,Vector{I}}, at::NTuple{N,AbstractFloat}) where {T,I<:Integer,N}
    quote
        (all(map(_ispos, at)) && sum(at) < 1) || throw(ArgumentError("all elements in \"at\" must be positive and their sum must be smaller than 1"))
        n = mapreduce(length, +, values(lm))
        k = nlabel(lm)
        # preallocate the indices vectors
        @nexprs $(N+1) i -> idx_i = Vector{I}()
        # sizehint will save us a few heavy memory allocations
        # we specify "+ k" to deal with trailing observations when
        # the number of observations from a class isn't divideable
        # by "at" or "1-at"
        @nexprs $(N) i -> sizehint!(idx_i, ceil(Int, n*at[i] + k))
        sizehint!($(Symbol(:idx_, Symbol(N+1))), ceil(Int, n*(1-sum(at)) + k))
        # loop through all label indices
        for indices in values(lm)
            tup = splitobs(indices, at)
            @nexprs $(N+1) i -> append!(idx_i, tup[i])
        end
        # return a tuple of all indices vectors
        @ntuple $(N+1) idx
    end
end

"""
    stratifiedobs([f], data, [p = 0.7], [shuffle = true], [obsdim]) -> Tuple

Partition the `data` into multiple disjoint subsets proportional
to the value(s) of `p`. The observations are assignmed to a data
subset using stratified sampling without replacement. These
subsets are then returned as a `Tuple` of subsets, where the
first element contains the fraction of observations of `data`
that is specified by the first float in `p`.

For example, if `p` is a `Float64` itself, then the return-value
will be a tuple with two elements (i.e. subsets), in which the
first element contains the fraction of observations specified by
`p` and the second element contains the rest. In the following
code the first subset `train` will contain around 70% of the
observations and the second subset `test` the rest. The key
difference to [`splitobs`](@ref) is that the class distribution
in `y` will actively be preserved in `train` and `test`.

```julia
train, test = stratifiedobs(y, p = 0.7)
```

If `p` is a tuple of `Float64` then additional subsets will be
created. In this example `train` will contain about 50% of the
observations, `val` will contain around 30%, and `test` the
remaining 20%.

```julia
train, val, test = stratifiedobs(y, p = (0.5, 0.3))
```

It is also possible to call `stratifiedobs` with multiple data
arguments as tuple, which all must have the same number of total
observations. Note that if `data` is a tuple, then it will be
assumed that the last element of the tuple contains the targets.

```julia
train, test = stratifiedobs((X, y), p = 0.7)
(X_train,y_train), (X_test,y_test) = stratifiedobs((X, y), p = 0.7)
```

The optional parameter `shuffle` determines if the resulting data
subsets should be shuffled. If `false`, then the observations in
the subsets will be grouped together according to their labels.

```julia
julia> y = ["a", "b", "b", "b", "b", "a"] # 2 imbalanced classes
6-element Array{String,1}:
 "a"
 "b"
 "b"
 "b"
 "b"
 "a"

julia> train, test = stratifiedobs(y, p = 0.5, shuffle = false)
(String["b","b","a"],String["b","b","a"])
```

The optional parameter `obsdim` can be used to specify which
dimension denotes the observations, if that concept makes sense
for the type of `data`. See `?ObsDim` for more information.

```julia
# 2 imbalanced classes in one-of-k encoding
julia> X = [1 0; 1 0; 1 0; 1 0; 0 1; 0 1]
6×2 Array{Int64,2}:
 1  0
 1  0
 1  0
 1  0
 0  1
 0  1

julia> train, test = stratifiedobs(argmax, X, p = 0.5, obsdim = 1)
([1 0; 1 0; 0 1], [0 1; 1 0; 1 0])
```

For this function to work, the type of `data` must implement
[`nobs`](@ref) and [`getobs`](@ref). For example, the following
code allows `stratifiedobs` to work on a `DataTable`.

```julia
# Make DataTables.jl work
LearnBase.getobs(data::DataTable, i) = data[i,:]
LearnBase.nobs(data::DataTable) = nrow(data)
```

You can use the parameter `f` to specify how to extract or
retrieve the targets from each observation of the given `data`.

```julia
julia> data = DataTable(Any[rand(6), rand(6), [:a,:b,:b,:b,:b,:a]], [:X1,:X2,:Y])
6×3 DataTables.DataTable
│ Row │ X1        │ X2          │ Y │
├─────┼───────────┼─────────────┼───┤
│ 1   │ 0.226582  │ 0.0443222   │ a │
│ 2   │ 0.504629  │ 0.722906    │ b │
│ 3   │ 0.933372  │ 0.812814    │ b │
│ 4   │ 0.522172  │ 0.245457    │ b │
│ 5   │ 0.505208  │ 0.11202     │ b │
│ 6   │ 0.0997825 │ 0.000341996 │ a │

julia> train, test = stratifiedobs(row->row[:Y], data, 0.5);

julia> getobs(train)
3×3 DataTables.DataTable
│ Row │ X1        │ X2          │ Y │
├─────┼───────────┼─────────────┼───┤
│ 1   │ 0.933372  │ 0.812814    │ b │
│ 2   │ 0.522172  │ 0.245457    │ b │
│ 3   │ 0.0997825 │ 0.000341996 │ a │

julia> getobs(test)
3×3 DataTables.DataTable
│ Row │ X1       │ X2        │ Y │
├─────┼──────────┼───────────┼───┤
│ 1   │ 0.504629 │ 0.722906  │ b │
│ 2   │ 0.226582 │ 0.0443222 │ a │
│ 3   │ 0.505208 │ 0.11202   │ b │
```

see [`DataSubset`](@ref) for more information on data subsets.

see also [`undersample`](@ref), [`oversample`](@ref), [`splitobs`](@ref).
"""
function stratifiedobs(data; p = 0.7, shuffle = true, obsdim = default_obsdim(data))
    stratifiedobs(identity, data, p, shuffle, convert(ObsDimension, obsdim))
end

function stratifiedobs(f, data; p = 0.7, shuffle = true, obsdim = default_obsdim(data))
    stratifiedobs(f, data, p, shuffle, convert(ObsDimension, obsdim))
end

function stratifiedobs(data, p::AbstractFloat, args...)
    stratifiedobs(identity, data, p, args...)
end

function stratifiedobs(data, p::NTuple{N,AbstractFloat}, args...) where N
    stratifiedobs(identity, data, p, args...)
end

function stratifiedobs(f, data, p::Union{NTuple,AbstractFloat}, shuffle::Bool = true, obsdim = default_obsdim(data))
    # The given data is always shuffled to qualify as performing
    # stratified sampling without replacement.
    data_shuf = shuffleobs(data, obsdim)
    # FIXME: if i put the following line before `data_shuf` is
    #        defined, it breaks type inference for `data_shuf`.
    allowcontainer(stratifiedobs, data) || throw(MethodError(stratifiedobs, (f,data,shuffle,obsdim)))
    idx_tup = splitobs(labelmap(eachtarget(f, data_shuf, obsdim)), p)
    # Setting the parameter "shuffle = false" specifies that the
    # classes are ordered in the resulting subsets respectively.
    shuffle && foreach(x->isempty(x) || shuffle!(x), idx_tup)
    map(idx -> datasubset(data_shuf, idx, obsdim), idx_tup)
end


text:
"""
$(TYPEDEF)
$(TYPEDFIELDS)

"""
struct Uniform <: AbstractDistDelayed
    low::NestedFloat
    high::NestedFloat
end

function uniform(low::Float64, high::Float64)::Float64
    return first(rand(Distributions.Uniform(low, high), 1))
end


text:
@testset "FASTA" begin
    @testset "Record" begin
        record = FASTA.Record()
        @test !isfilled(record)
        @test_throws ArgumentError FASTA.identifier(record)

        record = FASTA.Record(">foo\nACGT\n")
        @test isfilled(record)
        @test hasseqname(record)
        @test FASTA.hasidentifier(record)
        @test seqname(record) == FASTA.identifier(record) == "foo"
        @test !FASTA.hasdescription(record)
        @test_throws BioCore.Exceptions.MissingFieldException FASTA.description(record)
        @test hassequence(record)
        @test FASTA.hassequence(record)
        @test FASTA.sequence(record) == dna"ACGT"
        @test FASTA.sequence(record, 2:3) == dna"CG"
        @test FASTA.sequence(String, record) == "ACGT"
        @test FASTA.sequence(String, record, 2:3) == "CG"
        @test record == FASTA.Record(">foo\nACGT\n")

        record = FASTA.Record("""
        >CYS1_DICDI fragment
        SCWSFSTTGNVEGQHFISQNKL
        VSLSEQNLVDCDHECMEYEGE
        """)
        @test isfilled(record)
        @test FASTA.identifier(record) == "CYS1_DICDI"
        @test FASTA.description(record) == "fragment"
        @test FASTA.sequence(record) == aa"SCWSFSTTGNVEGQHFISQNKLVSLSEQNLVDCDHECMEYEGE"
        @test FASTA.sequence(record, 10:15) == aa"NVEGQH"
    end

    output = IOBuffer()
    writer = FASTA.Writer(output, 5)
    write(writer, FASTA.Record("seq1", dna"TTA"))
    write(writer, FASTA.Record("seq2", "some description", dna"ACGTNN"))
    flush(writer)
    @test String(take!(output)) == """
    >seq1
    TTA
    >seq2 some description
    ACGTN
    N
    """

    reader = FASTA.Reader(IOBuffer(
    """
    >seqA some description
    QIKDLLVSSSTDLDTTLKMK
    ILELPFASGDLSM
    >seqB
    VLMALGMTDLFIPSANLTG*
    """))
    record = FASTA.Record()
    @test read!(reader, record) === record
    @test FASTA.identifier(record) == "seqA"
    @test FASTA.description(record) == "some description"
    @test FASTA.sequence(record) == aa"QIKDLLVSSSTDLDTTLKMKILELPFASGDLSM"
    @test read!(reader, record) === record
    @test FASTA.identifier(record) == "seqB"
    @test !FASTA.hasdescription(record)
    @test FASTA.sequence(record) == aa"VLMALGMTDLFIPSANLTG*"

    function test_fasta_parse(filename, valid)
        # Reading from a stream
        stream = open(FASTA.Reader, filename)
        @test eltype(stream) == FASTA.Record
        if valid
            for seqrec in stream end
            @test true  # no error
            @test close(stream) === nothing
        else
            @test_throws Exception begin
                for seqrec in stream end
            end
            return
        end

        # in-place parsing
        stream = open(FASTA.Reader, filename)
        entry = eltype(stream)()
        while !eof(stream)
            read!(stream, entry)
        end

        # Check round trip
        output = IOBuffer()
        writer = FASTA.Writer(output, width = 60)
        outputB = IOBuffer()
        writerB = FASTA.Writer(outputB, width = -1)
        expected_entries = Any[]
        for seqrec in open(FASTA.Reader, filename)
            write(writer, seqrec)
            write(writerB, seqrec)
            push!(expected_entries, seqrec)
        end
        flush(writer)
        flush(writerB)

        seekstart(output)
        seekstart(outputB)
        read_entries = FASTA.Record[]
        read_entriesB = FASTA.Record[]
        for seqrec in FASTA.Reader(output)
            push!(read_entries, seqrec)
        end
        for seqrec in FASTA.Reader(outputB)
            push!(read_entriesB, seqrec)
        end
        @test expected_entries == read_entries
        @test expected_entries == read_entriesB
    end

    get_bio_fmt_specimens()
    path = joinpath(dirname(@__FILE__), "..", "BioFmtSpecimens", "FASTA")
    for specimen in YAML.load_file(joinpath(path, "index.yml"))
        tags = specimen["tags"]
        valid = get(specimen, "valid", true)
        if occursin("comments", tags)
            # currently comments are not supported
            continue
        end
        test_fasta_parse(joinpath(path, specimen["filename"]), valid)
    end

    @testset "Faidx" begin
        fastastr = """
        >chr1
        CCACACCACACCCACACACC
        >chr2
        ATGCATGCATGCAT
        GCATGCATGCATGC
        >chr3
        AAATAGCCCTCATGTACGTCTCCTCCAAGCCCTGTTGTCTCTTACCCGGA
        TGTTCAACCAAAAGCTACTTACTACCTTTATTTTATGTTTACTTTTTATA
        >chr4
        TACTT
        """
        # generated with `samtools faidx`
        faistr = """
        chr1	20	6	20	21
        chr2	28	33	14	15
        chr3	100	69	50	51
        chr4	5	177	5	6
        """
        mktempdir() do dir
            filepath = joinpath(dir, "test.fa")
            write(filepath, fastastr)
            write(filepath * ".fai", faistr)
            open(FASTA.Reader, filepath, index=filepath * ".fai") do reader
                chr3 = reader["chr3"]
                @test FASTA.identifier(chr3) == "chr3"
                @test FASTA.sequence(chr3) == dna"""
                AAATAGCCCTCATGTACGTCTCCTCCAAGCCCTGTTGTCTCTTACCCGGA
                TGTTCAACCAAAAGCTACTTACTACCTTTATTTTATGTTTACTTTTTATA
                """

                chr2 = reader["chr2"]
                @test FASTA.identifier(chr2) == "chr2"
                @test FASTA.sequence(chr2) == dna"""
                ATGCATGCATGCAT
                GCATGCATGCATGC
                """

                chr4 = reader["chr4"]
                @test FASTA.identifier(chr4) == "chr4"
                @test FASTA.sequence(chr4) == dna"""
                TACTT
                """

                chr1 = reader["chr1"]
                @test FASTA.identifier(chr1) == "chr1"
                @test FASTA.sequence(chr1) == dna"""
                CCACACCACACCCACACACC
                """

                @test_throws ArgumentError reader["chr5"]
            end
        end

        # invalid index
        @test_throws ArgumentError FASTA.Reader(IOBuffer(fastastr), index=π)
    end

    @testset "append" begin
        intempdir() do
            filepath = "test.fa"
            writer = open(FASTA.Writer, filepath)
            write(writer, FASTA.Record("seq1", dna"AAA"))
            close(writer)
            writer = open(FASTA.Writer, filepath, append=true)
            write(writer, FASTA.Record("seq2", dna"CCC"))
            close(writer)
            seqs = open(collect, FASTA.Reader, filepath)
            @test length(seqs) == 2
        end
    end
end


text:
using Gurobi
function printIIS(anyM::anyModel)

    # computes iis
    compute_conflict!(anyM.optModel)

    if MOI.get(anyM.optModel, MOI.ConflictStatus()) != MOI.ConflictStatusCode(3) return end
    # loops over constraint tables to find constraints within iis
    allCns_pair = vcat(collect.(vcat(anyM.parts.bal.cns, anyM.parts.trd.cns, anyM.parts.lim.cns, map(x -> x.cns,values(anyM.parts.tech))...))...)

    for cns in allCns_pair
        if cns[1] == :objEqn continue end

        allConstr_arr = findall(map(x -> MOI.ConflictParticipationStatusCode(0) != MOI.get(anyM.optModel.moi_backend, MOI.ConstraintConflictStatus(), x.index),cns[2][!,:cns]))
        # prints constraints within iis
        if !isempty(allConstr_arr)
            println("$(length(allConstr_arr)) of IIS in $(cns[1]) constraints.")
            colSet_dic = Dict(x => Symbol(split(string(x),"_")[1]) for x in intCol(cns[2]))
            for iisConstr in allConstr_arr
                row = cns[2][iisConstr,:]
                dimStr_arr = map(x -> row[x] == 0 ?  "" : string(x,": ",join(getUniName(row[x], anyM.sets[colSet_dic[x]])," < ")),collect(keys(colSet_dic)))
                println("$(join(filter(x -> x != "",dimStr_arr),", ")), constraint: $(row[:cns])")
            end
        end
    end
end
export printIIS

max_stars_repo_name: leonardgoeke/anyMOD

max_issues_repo_name: leonardgoeke/anyMOD

max_forks_repo_name: leonardgoeke/anyMOD


text:
###################################################################################################
# STRUCTS
###################################################################################################

# Alignment
struct AlignTemp
    strand::String                  # Methylation call strand
    R1::BAM.Record                  # First record from left to right
    R2::BAM.Record                  # Second record from left to right
end
mutable struct AllAlignTemps
    paired::Bool                    # Boolean indicating if record has pair
    templates::Array{AlignTemp,1}   # All alignment templates mapping to a region
end

# CpelTdm Configuration
mutable struct CpeltdmConfig
    pe::Bool                                # Paired/Single end
    min_cov::Union{Int64,Float64}           # Minimum coverage accepted
    matched::Bool                           # Matched vs unmatched group comparison
    bound_check::Bool                       # Parameter space boundary check
    trim::NTuple{4,Union{Int64,Float64}}    # Trimming of reads
    max_size_subreg::Union{Int64,Float64}   # Maximum size subregion (determines K)
    max_size_anal_reg::Union{Int64,Float64} # Maximum size of analysis region
    # Initializing method
    CpeltdmConfig(pe,min_cov,matched,bound_check,trim,max_size_subreg,max_size_anal_reg) =
        new(pe,min_cov,matched,bound_check,trim,max_size_subreg,max_size_anal_reg)
end

# ROI data
mutable struct RoiData
    # Chromosome properties
    chr::String
    chrst::Int64
    chrend::Int64
    n_vec::Vector{Int64}
    cpg_pos::Vector{Int64}
    # Analysis flags
    diff_anal_done::Bool
    analyzed1::Vector{Bool}
    analyzed2::Vector{Bool}
    # Data
    mmls1::Vector{Float64}
    mmls2::Vector{Float64}
    nmes1::Vector{Float64}
    nmes2::Vector{Float64}
    θ1s::Vector{Vector{Float64}}
    θ2s::Vector{Vector{Float64}}
    mml_test::NTuple{2,Float64}
    nme_test::NTuple{2,Float64}
    pdm_test::NTuple{2,Float64}
    # Initializing method
    RoiData(s1,s2) = new("",0,0,[],[],false,fill(false,s1),fill(false,s2),fill(NaN,s1),fill(NaN,s2),
        fill(NaN,s1),fill(NaN,s2),fill([],s1),fill([],s2),(NaN,NaN),(NaN,NaN),(NaN,NaN))
end


text:
module Functors

using MacroTools

export @functor, fmap, fmapstructure, fcollect

include("functor.jl")

end # module


text:
# Dispatcher.jl - Adaptive hash-graph persistency mechanism for computational
#                 task graphs written at 0x0α Research by Corneliu Cofaru, 2019

"""
DispatcherCache.jl is a `hash-chain` optimizer for Dispatcher delayed
execution graphs. It employes a hashing mechanism to check wether the
state associated to a node is the `DispatchGraph`  that is to
be executed has already been `hashed` (and hence, an output is available)
or, it is new or changed. Depending on the current state (by 'state'
one understands the called function source code, input arguments and
other input node dependencies), the current task becomes a load-from-disk
or an execute-and-store-to-disk operation. This is done is such a manner that
the minimimum number of load/execute operations are performed, minimizing
both persistency and computational demands.
"""
module DispatcherCache

    using Serialization
    using Dispatcher
    using IterTools
    using JSON
    using TranscodingStreams
    using CodecBzip2
    using CodecZlib

    import Dispatcher: run!
    export add_hash_cache!

    include("constants.jl")
    include("utils.jl")
    include("wrappers.jl")
    include("hash.jl")
    include("core.jl")

end # module


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

eltype(::Type{AbstractSet{T}}) where {T} = T

mutable struct Set{T} <: AbstractSet{T}
    dict::Dict{T,Nothing}

    Set{T}() where {T} = new(Dict{T,Nothing}())
    Set{T}(s::Set{T}) where {T} = new(Dict{T,Nothing}(s.dict))
end

Set{T}(itr) where {T} = union!(Set{T}(), itr)
Set() = Set{Any}()

"""
    Set([itr])

Construct a [`Set`](@ref) of the values generated by the given iterable object, or an
empty set. Should be used instead of [`BitSet`](@ref) for sparse integer sets, or
for sets of arbitrary objects.
"""
Set(itr) = Set{eltype(itr)}(itr)
function Set(g::Generator)
    T = @default_eltype(g)
    (_isleaftype(T) || T === Union{}) || return grow_to!(Set{T}(), g)
    return Set{T}(g)
end

empty(s::Set{T}, ::Type{U}=T) where {T,U} = Set{U}()

# return an empty set with eltype T, which is mutable (can be grown)
# by default, a Set is returned
emptymutable(s::AbstractSet{T}, ::Type{U}=T) where {T,U} = Set{U}()

_similar_for(c::AbstractSet, T, itr, isz) = empty(c, T)

function show(io::IO, s::Set)
    print(io, "Set(")
    show_vector(io, s)
    print(io, ')')
end

isempty(s::Set) = isempty(s.dict)
length(s::Set)  = length(s.dict)
in(x, s::Set) = haskey(s.dict, x)
push!(s::Set, x) = (s.dict[x] = nothing; s)
pop!(s::Set, x) = (pop!(s.dict, x); x)
pop!(s::Set, x, deflt) = x in s ? pop!(s, x) : deflt

function pop!(s::Set)
    isempty(s) && throw(ArgumentError("set must be non-empty"))
    idx = start(s.dict)
    val = s.dict.keys[idx]
    _delete!(s.dict, idx)
    val
end

delete!(s::Set, x) = (delete!(s.dict, x); s)

copy(s::Set) = copymutable(s)

# Set is the default mutable fall-back
copymutable(s::AbstractSet{T}) where {T} = Set{T}(s)

sizehint!(s::Set, newsz) = (sizehint!(s.dict, newsz); s)
empty!(s::Set) = (empty!(s.dict); s)
rehash!(s::Set) = (rehash!(s.dict); s)

start(s::Set)       = start(s.dict)
done(s::Set, state) = done(s.dict, state)
# NOTE: manually optimized to take advantage of Dict representation
next(s::Set, i)     = (s.dict.keys[i], skip_deleted(s.dict, i+1))

"""
    union(s, itrs...)
    ∪(s, itrs...)

Construct the union of sets. Maintain order with arrays.

# Examples
```jldoctest
julia> union([1, 2], [3, 4])
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> union([1, 2], [2, 4])
3-element Array{Int64,1}:
 1
 2
 4

julia> union([4, 2], 1:2)
3-element Array{Int64,1}:
 4
 2
 1

julia> union(Set([1, 2]), 2:3)
Set([2, 3, 1])
```
"""
function union end

_in(itr) = x -> x in itr

union(s, sets...) = union!(emptymutable(s, promote_eltype(s, sets...)), s, sets...)
union(s::AbstractSet) = copy(s)

const ∪ = union

"""
    union!(s::Union{AbstractSet,AbstractVector}, itrs...)

Construct the union of passed in sets and overwrite `s` with the result.
Maintain order with arrays.

# Examples
```jldoctest
julia> a = Set([1, 3, 4, 5]);

julia> union!(a, 1:2:8);

julia> a
Set([7, 4, 3, 5, 1])
```
"""
union!(s::AbstractSet, sets...) = foldl(union!, s, sets)

# default generic 2-args implementation with push!
union!(s::AbstractSet, itr) = foldl(push!, s, itr)

function union!(s::Set{T}, itr) where T
    haslength(itr) && sizehint!(s, length(itr))
    for x=itr
        push!(s, x)
        length(s) == max_values(T) && break
    end
    s
end

"""
    intersect(s, itrs...)
    ∩(s, itrs...)

Construct the intersection of sets.
Maintain order with arrays.

# Examples
```jldoctest
julia> intersect([1, 2, 3], [3, 4, 5])
1-element Array{Int64,1}:
 3

julia> intersect([1, 4, 4, 5, 6], [4, 6, 6, 7, 8])
2-element Array{Int64,1}:
 4
 6

julia> intersect(Set([1, 2]), BitSet([2, 3]))
Set([2])
```
"""
intersect(s::AbstractSet, itr, itrs...) = intersect!(intersect(s, itr), itrs...)
intersect(s) = union(s)
intersect(s::AbstractSet, itr) = mapfilter(_in(s), push!, itr, emptymutable(s))

const ∩ = intersect

"""
    intersect!(s::Union{AbstractSet,AbstractVector}, itrs...)

Intersect all passed in sets and overwrite `s` with the result.
Maintain order with arrays.
"""
intersect!(s::AbstractSet, itrs...) = foldl(intersect!, s, itrs)
intersect!(s::AbstractSet, s2::AbstractSet) = filter!(_in(s2), s)
intersect!(s::AbstractSet, itr) = intersect!(s, union!(emptymutable(s), itr))

"""
    setdiff(s, itrs...)

Construct the set of elements in `s` but not in any of the iterables in `itrs`.
Maintain order with arrays.

# Examples
```jldoctest
julia> setdiff([1,2,3], [3,4,5])
2-element Array{Int64,1}:
 1
 2
```
"""
setdiff(s::AbstractSet, itrs...) = setdiff!(copymutable(s), itrs...)
setdiff(s) = union(s)

"""
    setdiff!(s, itrs...)

Remove from set `s` (in-place) each element of each iterable from `itrs`.
Maintain order with arrays.

# Examples
```jldoctest
julia> a = Set([1, 3, 4, 5]);

julia> setdiff!(a, 1:2:6);

julia> a
Set([4])
```
"""
setdiff!(s::AbstractSet, itrs...) = foldl(setdiff!, s, itrs)
setdiff!(s::AbstractSet, itr) = foldl(delete!, s, itr)

"""
    symdiff(s, itrs...)

Construct the symmetric difference of elements in the passed in sets.
When `s` is not an `AbstractSet`, the order is maintained.
Note that in this case the multiplicity of elements matters.

# Examples
```jldoctest
julia> symdiff([1,2,3], [3,4,5], [4,5,6])
3-element Array{Int64,1}:
 1
 2
 6

julia> symdiff([1,2,1], [2, 1, 2])
2-element Array{Int64,1}:
 1
 2

julia> symdiff(unique([1,2,1]), unique([2, 1, 2]))
0-element Array{Int64,1}
```
"""
symdiff(s, sets...) = symdiff!(emptymutable(s, promote_eltype(s, sets...)), s, sets...)
symdiff(s) = symdiff!(copy(s))

"""
    symdiff!(s::Union{AbstractSet,AbstractVector}, itrs...)

Construct the symmetric difference of the passed in sets, and overwrite `s` with the result.
When `s` is an array, the order is maintained.
Note that in this case the multiplicity of elements matters.
"""
symdiff!(s::AbstractSet, itrs...) = foldl(symdiff!, s, itrs)

function symdiff!(s::AbstractSet, itr)
    for x in itr
        x in s ? delete!(s, x) : push!(s, x)
    end
    s
end

==(l::Set, r::Set) = (length(l) == length(r)) && (l <= r)
<( l::Set, r::Set) = (length(l) < length(r)) && (l <= r)
<=(l::Set, r::Set) = issubset(l, r)

"""
    issubset(a, b)
    ⊆(a,b) -> Bool
    ⊈(a,b) -> Bool
    ⊊(a,b) -> Bool

Determine whether every element of `a` is also in `b`, using [`in`](@ref).

# Examples
```jldoctest
julia> issubset([1, 2], [1, 2, 3])
true

julia> issubset([1, 2, 3], [1, 2])
false
```
"""
function issubset(l, r)
    for elt in l
        if !in(elt, r)
            return false
        end
    end
    return true
end

# use the implementation below when it becoms as efficient
# issubset(l, r) = all(_in(r), l)

const ⊆ = issubset
⊊(l::Set, r::Set) = <(l, r)
⊈(l::Set, r::Set) = !⊆(l, r)
⊇(l, r) = issubset(r, l)
⊉(l::Set, r::Set) = !⊇(l, r)
⊋(l::Set, r::Set) = <(r, l)

⊊(l::T, r::T) where {T<:AbstractSet} = <(l, r)
⊈(l::T, r::T) where {T<:AbstractSet} = !⊆(l, r)
⊉(l::T, r::T) where {T<:AbstractSet} = !⊇(l, r)
⊋(l::T, r::T) where {T<:AbstractSet} = <(r, l)

"""
    unique(itr)

Return an array containing only the unique elements of collection `itr`,
as determined by [`isequal`](@ref), in the order that the first of each
set of equivalent elements originally appears. The element type of the
input is preserved.

# Examples
```jldoctest
julia> unique([1, 2, 6, 2])
3-element Array{Int64,1}:
 1
 2
 6

julia> unique(Real[1, 1.0, 2])
2-element Array{Real,1}:
 1
 2
```
"""
function unique(itr)
    T = @default_eltype(itr)
    out = Vector{T}()
    seen = Set{T}()
    i = start(itr)
    if done(itr, i)
        return out
    end
    x, i = next(itr, i)
    if !_isleaftype(T) && iteratoreltype(itr) == EltypeUnknown()
        S = typeof(x)
        return _unique_from(itr, S[x], Set{S}((x,)), i)
    end
    push!(seen, x)
    push!(out, x)
    return unique_from(itr, out, seen, i)
end

_unique_from(itr, out, seen, i) = unique_from(itr, out, seen, i)
@inline function unique_from(itr, out::Vector{T}, seen, i) where T
    while !done(itr, i)
        x, i = next(itr, i)
        S = typeof(x)
        if !(S === T || S <: T)
            R = typejoin(S, T)
            seenR = convert(Set{R}, seen)
            outR = convert(Vector{R}, out)
            if !in(x, seenR)
                push!(seenR, x)
                push!(outR, x)
            end
            return _unique_from(itr, outR, seenR, i)
        end
        if !in(x, seen)
            push!(seen, x)
            push!(out, x)
        end
    end
    return out
end

"""
    unique(f, itr)

Returns an array containing one value from `itr` for each unique value produced by `f`
applied to elements of `itr`.

# Examples
```jldoctest
julia> unique(x -> x^2, [1, -1, 3, -3, 4])
3-element Array{Int64,1}:
 1
 3
 4
```
"""
function unique(f::Callable, C)
    out = Vector{eltype(C)}()
    seen = Set()
    for x in C
        y = f(x)
        if !in(y, seen)
            push!(seen, y)
            push!(out, x)
        end
    end
    out
end

# If A is not grouped, then we will need to keep track of all of the elements that we have
# seen so far.
function _unique!(A::AbstractVector)
    seen = Set{eltype(A)}()
    idxs = eachindex(A)
    i = state = start(idxs)
    for x in A
        if x ∉ seen
            push!(seen, x)
            i, state = next(idxs, state)
            A[i] = x
        end
    end
    resize!(A, i - first(idxs) + 1)
end

# If A is grouped, so that each unique element is in a contiguous group, then we only
# need to keep track of one element at a time. We replace the elements of A with the
# unique elements that we see in the order that we see them. Once we have iterated
# through A, we resize A based on the number of unique elements that we see.
function _groupedunique!(A::AbstractVector)
    isempty(A) && return A
    idxs = eachindex(A)
    y = first(A)
    state = start(idxs)
    i, state = next(idxs, state)
    for x in A
        if !isequal(x, y)
            i, state = next(idxs, state)
            y = A[i] = x
        end
    end
    resize!(A, i - first(idxs) + 1)
end

"""
    unique!(A::AbstractVector)

Remove duplicate items as determined by [`isequal`](@ref), then return the modified `A`.
`unique!` will return the elements of `A` in the order that they occur. If you do not care
about the order of the returned data, then calling `(sort!(A); unique!(A))` will be much
more efficient as long as the elements of `A` can be sorted.

# Examples
```jldoctest
julia> unique!([1, 1, 1])
1-element Array{Int64,1}:
 1

julia> A = [7, 3, 2, 3, 7, 5];

julia> unique!(A)
4-element Array{Int64,1}:
 7
 3
 2
 5

julia> B = [7, 6, 42, 6, 7, 42];

julia> sort!(B);  # unique! is able to process sorted data much more efficiently.

julia> unique!(B)
3-element Array{Int64,1}:
  6
  7
 42
```
"""
function unique!(A::Union{AbstractVector{<:Real}, AbstractVector{<:AbstractString},
                          AbstractVector{<:Symbol}})
    if isempty(A)
        return A
    elseif issorted(A) || issorted(A, rev=true)
        return _groupedunique!(A)
    else
        return _unique!(A)
    end
end
# issorted fails for some element types, so the method above has to be restricted to
# elements with isless/< defined.
function unique!(A)
    if isempty(A)
        return A
    else
        return _unique!(A)
    end
end

"""
    allunique(itr) -> Bool

Return `true` if all values from `itr` are distinct when compared with [`isequal`](@ref).

# Examples
```jldoctest
julia> a = [1; 2; 3]
3-element Array{Int64,1}:
 1
 2
 3

julia> allunique([a, a])
false
```
"""
function allunique(C)
    seen = Set{eltype(C)}()
    for x in C
        if in(x, seen)
            return false
        else
            push!(seen, x)
        end
    end
    true
end

allunique(::Set) = true

allunique(r::AbstractRange{T}) where {T} = (step(r) != zero(T)) || (length(r) <= 1)

filter(pred, s::AbstractSet) = mapfilter(pred, push!, s, emptymutable(s))
filter!(f, s::Set) = unsafe_filter!(f, s)

# it must be safe to delete the current element while iterating over s:
unsafe_filter!(pred, s::AbstractSet) = mapfilter(!pred, delete!, s, s)

# TODO: delete mapfilter in favor of comprehensions/foldl/filter when competitive
function mapfilter(pred, f, itr, res)
    for x in itr
        pred(x) && f(res, x)
    end
    res
end

const hashs_seed = UInt === UInt64 ? 0x852ada37cfe8e0ce : 0xcfe8e0ce
function hash(s::Set, h::UInt)
    hv = hashs_seed
    for x in s
        hv ⊻= hash(x)
    end
    hash(hv, h)
end

convert(::Type{T}, s::T) where {T<:AbstractSet} = s
convert(::Type{T}, s::AbstractSet) where {T<:AbstractSet} = T(s)

## replace/replace! ##

"""
    replace!(A, old_new::Pair...; [count::Integer])

For each pair `old=>new` in `old_new`, replace all occurrences
of `old` in collection `A` by `new`.
If `count` is specified, then replace at most `count` occurrences in total.
See also [`replace`](@ref replace(A, old_new::Pair...)).

# Examples
```jldoctest
julia> replace!([1, 2, 1, 3], 1=>0, 2=>4, count=2)
4-element Array{Int64,1}:
 0
 4
 1
 3

julia> replace!(Set([1, 2, 3]), 1=>0)
Set([0, 2, 3])
```
"""
replace!(A, old_new::Pair...; count::Integer=typemax(Int)) = _replace!(A, eltype(A), count, old_new)

# we use this wrapper because using directly eltype(A) as the type
# parameter below for Some degrades performance
function _replace!(A, ::Type{K}, count::Integer, old_new::Tuple{Vararg{Pair}}) where K
    @inline function prednew(x)
        for o_n in old_new
            first(o_n) == x && return Some{K}(last(o_n))
        end
    end
    replace!(prednew, A, count=count)
end

"""
    replace!(pred::Function, A, new; [count::Integer])

Replace all occurrences `x` in collection `A` for which `pred(x)` is true
by `new`.

# Examples
```jldoctest
julia> A = [1, 2, 3, 1];

julia> replace!(isodd, A, 0, count=2)
4-element Array{Int64,1}:
 0
 2
 0
 1
```
"""
replace!(pred::Callable, A, new; count::Integer=typemax(Int)) =
    replace!(x -> if pred(x) Some(new) end, A, count=count)

"""
    replace!(prednew::Function, A; [count::Integer])

For each value `x` in `A`, `prednew(x)` is called and must
return either `nothing`, in which case no replacement occurs,
or a value, possibly wrapped as a [`Some`](@ref) object, which
will be used as a replacement for `x`.

# Examples
```jldoctest
julia> replace!(x -> isodd(x) ? 2x : nothing, [1, 2, 3, 4])
4-element Array{Int64,1}:
 2
 2
 6
 4

julia> replace!(Union{Int,Nothing}[0, 1, 2, nothing, 4], count=2) do x
           x !== nothing && iseven(x) ? Some(nothing) : nothing
       end
5-element Array{Union{Nothing,Int64},1}:
  nothing
 1
  nothing
  nothing
 4

julia> replace!(Dict(1=>2, 3=>4)) do kv
           if first(kv) < 3; first(kv)=>3 end
       end
Dict{Int64,Int64} with 2 entries:
  3 => 4
  1 => 3

julia> replace!(x->2x, Set([3, 6]))
Set([6, 12])
```
"""
replace!(prednew::Callable, A; count::Integer=typemax(Int)) =
    replace!(prednew, A, count=clamp(count, typemin(Int), typemax(Int)) % Int)

"""
    replace(A, old_new::Pair...; [count::Integer])

Return a copy of collection `A` where, for each pair `old=>new` in `old_new`,
all occurrences of `old` are replaced by `new`.
If `count` is specified, then replace at most `count` occurrences in total.
See also [`replace!`](@ref).

# Examples
```jldoctest
julia> replace([1, 2, 1, 3], 1=>0, 2=>4, count=2)
4-element Array{Int64,1}:
 0
 4
 1
 3
```
"""
replace(A, old_new::Pair...; count::Integer=typemax(Int)) =
    _replace!(copy(A), eltype(A), count, old_new)

"""
    replace(pred::Function, A, new; [count::Integer])

Return a copy of collection `A` where all occurrences `x` for which
`pred(x)` is true are replaced by `new`.

# Examples
```jldoctest
julia> replace(isodd, [1, 2, 3, 1], 0, count=2)
4-element Array{Int64,1}:
 0
 2
 0
 1
```
"""
replace(pred::Callable, A, new; count::Integer=typemax(Int)) =
    replace!(x -> if pred(x) Some(new) end, copy(A), count=count)

"""
    replace(prednew::Function, A; [count::Integer])

Return a copy of `A` where for each value `x` in `A`, `prednew(x)` is called
and must return  either `nothing`, in which case no replacement occurs,
or a value, possibly wrapped as a [`Some`](@ref) object, which
will be used as a replacement for `x`.

# Examples
```jldoctest
julia> replace(x -> isodd(x) ? 2x : nothing, [1, 2, 3, 4])
4-element Array{Int64,1}:
 2
 2
 6
 4

julia> replace(Union{Int,Nothing}[0, 1, 2, nothing, 4], count=2) do x
           x !== nothing && iseven(x) ? Some(nothing) : nothing
       end
5-element Array{Union{Nothing,Int64},1}:
  nothing
 1
  nothing
  nothing
 4

julia> replace(Dict(1=>2, 3=>4)) do kv
           if first(kv) < 3; first(kv)=>3 end
       end
Dict{Int64,Int64} with 2 entries:
  3 => 4
  1 => 3
```
"""
replace(prednew::Callable, A; count::Integer=typemax(Int)) = replace!(prednew, copy(A), count=count)

# Handle ambiguities
replace!(a::Callable, b::Pair; count::Integer=-1) = throw(MethodError(replace!, (a, b)))
replace!(a::Callable, b::Pair, c::Pair; count::Integer=-1) = throw(MethodError(replace!, (a, b, c)))
replace(a::Callable, b::Pair; count::Integer=-1) = throw(MethodError(replace, (a, b)))
replace(a::Callable, b::Pair, c::Pair; count::Integer=-1) = throw(MethodError(replace, (a, b, c)))
replace(a::AbstractString, b::Pair, c::Pair) = throw(MethodError(replace, (a, b, c)))

### replace! for AbstractDict/AbstractSet

askey(k, ::AbstractDict) = k.first
askey(k, ::AbstractSet) = k

function _replace_update_dict!(repl::Vector{<:Pair}, x, y::Some)
    push!(repl, x => y.value)
    true
end

_replace_update_dict!(repl::Vector{<:Pair}, x, ::Nothing) = false
_replace_update_dict!(repl::Vector{<:Pair}, x, y) = _replace_update_dict!(repl, x, Some(y))

function replace!(prednew::Callable, A::Union{AbstractDict,AbstractSet}; count::Int=typemax(Int))
    count < 0 && throw(DomainError(count, "`count` must not be negative"))
    count == 0 && return A
    repl = Pair{eltype(A),eltype(A)}[]
    c = 0
    for x in A
        c += _replace_update_dict!(repl, x, prednew(x))
        c == count && break
    end
    for oldnew in repl
        pop!(A, askey(first(oldnew), A))
    end
    for oldnew in repl
        push!(A, last(oldnew))
    end
    A
end

### AbstractArray

function _replace_update!(A::AbstractArray, i::Integer, y::Some)
    @inbounds A[i] = y.value
    true
end

_replace_update!(A::AbstractArray, i::Integer, ::Nothing) = false
_replace_update!(A::AbstractArray, i::Integer, y) = _replace_update!(A, i, Some(y))

function replace!(prednew::Callable, A::AbstractArray; count::Int=typemax(Int))
    count < 0 && throw(DomainError(count, "`count` must not be negative"))
    count == 0 && return A
    c = 0
    for i in eachindex(A)
        c += _replace_update!(A, i, prednew(A[i]))
        c == count && break
    end
    A
end

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
using SafeTestsets

@safetestset "Utils" begin include("utils_test.jl") end

@safetestset "RkMatrix" begin include("rkmatrix_test.jl") end

@safetestset "Compressors" begin include("compressor_test.jl") end

@safetestset "HMatrix" begin include("hmatrix_test.jl") end

@safetestset "Addition" begin include("addition_test.jl") end

@safetestset "Multiplication" begin include("multiplication_test.jl") end

@safetestset "Triangular" begin include("triangular_test.jl") end

@safetestset "LU" begin include("lu_test.jl") end

# @safetestset "DHMatrix" begin include("dhmatrix_test.jl") end

max_stars_repo_name: WaveProp/HMatrices

max_forks_repo_name: WaveProp/HMatrices


text:
"""
    IndicesArrayStyle{S}

This is a `BroadcastStyle` for IndicesArray's It preserves the dimension
names. `S` should be the `BroadcastStyle` of the wrapped type.
"""
struct IndicesArrayStyle{S <: BroadcastStyle} <: AbstractArrayStyle{Any} end
IndicesArrayStyle(::S) where {S} = IndicesArrayStyle{S}()
IndicesArrayStyle(::S, ::Val{N}) where {S,N} = IndicesArrayStyle(S(Val(N)))
IndicesArrayStyle(::Val{N}) where N = IndicesArrayStyle{DefaultArrayStyle{N}}()
function IndicesArrayStyle(a::BroadcastStyle, b::BroadcastStyle)
    inner_style = BroadcastStyle(a, b)

    # if the inner_style is Unknown then so is the outer-style
    if inner_style isa Unknown
        return Unknown()
    else
        return IndicesArrayStyle(inner_style)
    end
end

function Base.BroadcastStyle(::Type{T}) where {T<:IndicesArray}
    return IndicesArrayStyle{typeof(BroadcastStyle(parent_type(T)))}()
end

Base.BroadcastStyle(::IndicesArrayStyle{A}, ::IndicesArrayStyle{B}) where {A, B} = IndicesArrayStyle(A(), B())
Base.BroadcastStyle(::IndicesArrayStyle{A}, b::B) where {A, B} = IndicesArrayStyle(A(), b)
Base.BroadcastStyle(a::A, ::IndicesArrayStyle{B}) where {A, B} = IndicesArrayStyle(a, B())
Base.BroadcastStyle(::IndicesArrayStyle{A}, b::DefaultArrayStyle) where {A} = IndicesArrayStyle(A(), b)
Base.BroadcastStyle(a::AbstractArrayStyle{M}, ::IndicesArrayStyle{B}) where {B,M} = IndicesArrayStyle(a, B())

"""
    unwrap_broadcasted

Recursively unwraps `IndicesArray`s and `IndicesArrayStyle`s.
replacing the `IndicesArray`s with the wrapped array,
and `IndicesArrayStyle` with the wrapped `BroadcastStyle`.
"""
function unwrap_broadcasted(bc::Broadcasted{IndicesArrayStyle{S}}) where S
    return Broadcasted{S}(bc.f, map(unwrap_broadcasted, bc.args))
end
#unwrap_broadcasted(bc::Broadcasted{BroadcastIndexStyle{S}}) where S = Broadcasted{S}(bc.f, map(unwrap_broadcasted, bc.args))
unwrap_broadcasted(a::IndicesArray) = parent(a)
#unwrap_broadcasted(a::AbstractIndex) = parent(a)
unwrap_broadcasted(x) = x

# We need to implement copy because if the wrapper array type does not support setindex
# then the `similar` based default method will not work
# TODO: need to make this flexible to other IndicesArrays types
function Broadcast.copy(bc::Broadcasted{IndicesArrayStyle{S}}) where S
    data = unwrap_broadcasted(bc)
    return rebuild(IArray, copy(data), combine_axes(bc.args...))
end
# TODO: copyto! for broadcasting

# TODO
@inline function Broadcast.combine_axes(A::IndicesArray, B::IndicesArray, C...)
    return combine_axes(combine_axes(axes(A), axes(B)), C...)
end

@inline function Broadcast.combine_axes(A::AbstractArray, B::IndicesArray, C...)
    return combine_axes(combine_axes(axes(A), axes(B)), C...)
end

@inline function Broadcast.combine_axes(A::IndicesArray, B::AbstractArray, C...)
    return combine_axes(combine_axes(axes(A), axes(B)), C...)
end

Broadcast.combine_axes(A::IndicesArray, B::IndicesArray) = _combine_axes(indices(A), indices(B))
Broadcast.combine_axes(A::IndicesArray, B::AbstractArray) = _combine_axes(indices(A), indices(B))
Broadcast.combine_axes(A::AbstractArray, B::IndicesArray) = _combine_axes(indices(A), indices(B))
Broadcast.combine_axes(A::IndicesArray) = indices(A)

_combine_axes(a::Tuple{Any,Vararg{Any}}, b::Tuple{Any,Vararg{Any}}) = combine_indices(a, b)
#    (combine_indices(first(a), first(b))..., _combine_axes(tail(a), tail(b))...)
#end
_combine_axes(a::Tuple{Any,Vararg{Any}}, b::Tuple{}) = a
_combine_axes(a::Tuple{}, b::Tuple{Any,Vararg{Any}}) = b
_combine_axes(a::Tuple{}, b::Tuple{}) = ()

function Broadcast.combine_axes(
    A::Tuple{<:AbstractIndex, Vararg{Any}},
    B::Tuple{<:AbstractIndex, Vararg{Any}}
   )
    return (combine_indices(first(A), first(B))..., combine_axes(tail(A), tail(B))...)
end

function Broadcast.combine_axes(
    A::Tuple{<:AbstractIndex, Vararg{Any}},
    B::Tuple{Any, Vararg{Any}}
   )
    return (combine_indices(first(A), first(B))..., combine_axes(tail(A), tail(B))...)
end

function Broadcast.combine_axes(
    A::Tuple{Any, Vararg{Any}},
    B::Tuple{<:AbstractIndex, Vararg{Any}}
   )
    return (combine_indices(first(A), first(B))..., combine_axes(tail(A), tail(B))...)
end

function Broadcast.combine_axes(A::Tuple{}, B::Tuple{<:AbstractIndex, Vararg{Any}})
    return (combine_indices(first(B))..., combine_axes(A, tail(B))...)
end

function Broadcast.combine_axes(A::Tuple{<:AbstractIndex, Vararg{Any}}, B::Tuple{})
    return (combine_indices(first(A))..., combine_axes(tail(A), B)...)
end


text:
#!/usr/bin/env julia

# Resources:
#  https://learnxinyminutes.com/docs/julia/
#  https://en.wikibooks.org/wiki/Introducing_Julia/Dictionaries_and_sets
#  https://en.wikibooks.org/wiki/Introducing_Julia/Arrays_and_tuples#Tuples
#  https://docs.julialang.org/en/v1/manual/types/#Type-Aliases-1

module Wires

const Coordinates = Tuple{Int64, Int64}

const Wire = Dict{Coordinates, Int64}

function parse_direction(direction :: AbstractString)
  len = parse(Int64, direction[2:end])
  letter = direction[1]

  if letter == 'U'
    (0, len)
  elseif letter == 'R'
    (len, 0)
  elseif letter == 'D'
    (0, -len)
  elseif letter == 'L'
    (-len, 0)
  else
    error("Invalid direction '$(direction)'")
  end
end

function gen_range(start :: Int64, offset :: Int64)
  step = offset < 0 ? -1 : 1

  start:step:start+offset
end

function gen_segment(origin :: Coordinates, direction :: Coordinates, index :: Int64)
  segment = Wire()

  ox, oy = origin
  px, py = direction

  for x in gen_range(ox, px)
    for y in gen_range(oy, py)
      if x == ox && y == oy
        continue
      end

      segment[x, y] = index
      index += 1
    end
  end

  segment
end

function parse_wire(line :: String)
  wire = Wire()
  position :: Coordinates = (0, 0)
  index = 1

  for direction_string = split(line, ",")
    direction = parse_direction(direction_string)
    segment = gen_segment(position, direction, index)

    for (coords, index) in segment
      if !haskey(wire, coords)
        wire[coords] = index
      end
    end

    index += length(segment)
    position = (position[1]+direction[1], position[2]+direction[2])
  end

  wire
end

function manhattan_distance_to_origin(cell :: Coordinates)
  x, y = cell
  abs(x) + abs(y)
end

function wire_coordinates_set(wire :: Wire)
  Set([k for (k, _) in wire])
end

function wires_crossings_set(wire1 :: Wire, wire2 :: Wire)
  intersect(wire_coordinates_set(wire1), wire_coordinates_set(wire2))
end

function wires_crossing_distance(wire1 :: Wire, wire2 :: Wire)
  crossings = wires_crossings_set(wire1, wire2)
  minimum(manhattan_distance_to_origin, crossings)
end

function wires_crossing_steps(wire1 :: Wire, wire2 :: Wire)
  crossings = wires_crossings_set(wire1, wire2)
  best_steps = length(wire1) + length(wire2)

  for coords in crossings
    steps = wire1[coords] + wire2[coords]

    if steps < best_steps
      best_steps = steps
    end
  end

  best_steps
end

function problem1(line1 :: String, line2 :: String)
  wire1 = parse_wire(line1)
  wire2 = parse_wire(line2)

  wires_crossing_distance(wire1, wire2)
end

function problem2(line1 :: String, line2 :: String)
  wire1 = parse_wire(line1)
  wire2 = parse_wire(line2)

  wires_crossing_steps(wire1, wire2)
end

export parse_wire, parse_direction, gen_segment, wires_crossing_distance,
  problem1, problem2

end # end module

if PROGRAM_FILE != "" && realpath(@__FILE__) == realpath(PROGRAM_FILE)
  using Printf

  if isempty(ARGS)
    println("Usage:\n  $(PROGRAM_FILE) 1|2 input.txt")
  else
    using .Wires

    open(ARGS[2]) do f
      line1, line2 = eachline(f)

      if ARGS[1] == "1"
        println(problem1(line1, line2))
      elseif ARGS[1] == "2"
        println(problem2(line1, line2))
      else
        error("Invalid problem number '$(ARGS[1])")
      end
    end
  end
end

max_stars_repo_name: bfontaine/advent-of-code

max_issues_repo_name: bfontaine/advent-of-code

max_forks_repo_name: bfontaine/advent-of-code