text:
struct Observation
    pixel::Point2f
    point::Point3f
    pose::NTuple{6, Float64}

    point_order::Int64
    pose_order::Int64

    constant::Bool
    in_covmap::Bool

    kfid::Int64
    mpid::Int64
end

struct LocalBACache
    observations::Vector{Observation}
    outliers::Vector{Bool} # same length as observations
    bad_keypoints::Set{Int64}

    θ::Vector{Float64}
    θconst::Vector{Bool}
    pixels::Matrix{Float64}

    poses_ids::Vector{Int64}
    points_ids::Vector{Int64}

    poses_remap::Vector{Int64} # order id → kfid
    points_remap::Vector{Int64} # order id → mpid
end

function LocalBACache(
    observations, bad_keypoints, θ, θconst, pixels,
    poses_ids, points_ids, poses_remap, points_remap,
)
    outliers = Vector{Bool}(undef, length(observations))
    LocalBACache(
        observations, outliers, bad_keypoints, θ, θconst, pixels,
        poses_ids, points_ids, poses_remap, points_remap)
end

"""
Estimator is responsible for performing Bundle-Adjustment procedure
and map filtering that removes low-information Keyframes from the map.

# Arguments

- `map_manager::MapManager`: Map manager with all the map information.
- `params::Params`: Parameters of the system.
- `frame_queue::Vector{Frame}`: Queue of the Frames to process.
- `new_kf_available::Bool`: Whether or not new Keyframe in the queue
    is available.
- `exit_required::Bool`: Whether the estimator should end its work and exit.
- `queue_lock:ReentrantLock`: Lock for the queue.
"""
mutable struct Estimator{V}
    map_manager::MapManager
    params::Params
    slam_io::V

    frame_queue::Vector{Frame}
    new_kf_available::Bool
    exit_required::Bool

    queue_lock::ReentrantLock
end

function Estimator(map_manager::MapManager, params::Params, slam_io)
    Estimator(
        map_manager, params, slam_io, Frame[],
        false, false, ReentrantLock())
end

"""
    run!(estimator::Estimator)

Main routine that starts Estimator.
"""
function run!(estimator::Estimator)
    while !estimator.exit_required
        new_kf = get_new_kf!(estimator)
        if new_kf ≡ nothing
            sleep(1e-2)
            continue
        end

        if estimator.params.do_local_bundle_adjustment && new_kf.kfid ≥ 2
            lock(estimator.map_manager.optimization_lock)
            try
                t1 = time()
                local_bundle_adjustment!(estimator, new_kf)
                t2 = time()
                @debug "[ES] BA Time ($(t2 - t1) seconds)."
            catch e
                showerror(stdout, e); println()
                display(stacktrace(catch_backtrace())); println()
            finally
                unlock(estimator.map_manager.optimization_lock)
            end
        end

        if estimator.params.map_filtering
            t1 = time()
            map_filtering!(estimator, new_kf)
            t2 = time()
            @debug "[ES] Map filtering ($(t2 - t1) seconds)."
        end
    end
    @debug "[ES] Exit required."
end

"""
    add_new_kf!(estimator::Estimator, frame::Frame)

Add `frame` to the Estimator queue to be processed.
"""
function add_new_kf!(estimator::Estimator, frame::Frame)
    lock(estimator.queue_lock) do
        push!(estimator.frame_queue, frame)
        estimator.new_kf_available = true
    end
end

"""
    get_new_kf!(estimator::Estimator)

Get `frame` from the queue if it is available.
"""
function get_new_kf!(estimator::Estimator)
    lock(estimator.queue_lock) do
        if isempty(estimator.frame_queue)
            estimator.new_kf_available = false
            return nothing
        end

        # TODO if more than 1 frame in queue, add them to ba anyway.
        estimator.new_kf_available = false
        kf = popfirst!(estimator.frame_queue)
        kf
    end
end

function _get_ba_parameters(
    map_manager::MapManager, frame::Frame,
    covisibility_map::OrderedDict{Int64, Int64}, min_cov_score,
)
    # poses: kfid → (order id, θ).
    poses = Dict{Int64, Tuple{Int64, NTuple{6, Float64}}}()
    constant_poses = Set{Int64}()
    # map_points: mpid → (order id, point).
    map_points = Dict{Int64, Tuple{Int64, Point3f}}()

    processed_keypoints_ids = Set{Int64}()
    bad_keypoints = Set{Int64}()

    observations = Vector{Observation}(undef, 0)
    poses_remap = Vector{Int64}(undef, 0)
    points_remap = Vector{Int64}(undef, 0)
    sizehint!(observations, 1000)
    sizehint!(poses_remap, 10)
    sizehint!(points_remap, 1000)

    for (co_kfid, score) in covisibility_map
        co_frame = get_keyframe(map_manager, co_kfid)
        co_frame ≡ nothing && (remove_covisible_kf!(frame, co_kfid); continue)
        (co_kfid > frame.kfid || get_3d_keypoints_nb(co_frame) == 0 ||
            score == 0) && continue

        if !(co_kfid in keys(poses))
            if !(co_kfid in constant_poses)
                is_constant = score < min_cov_score || co_kfid == 0
                is_constant && (push!(constant_poses, co_kfid); continue)
            end
        end

        for kpid in get_3d_keypoints_ids(co_frame)
            kpid in processed_keypoints_ids && continue
            push!(processed_keypoints_ids, kpid)

            mp = get_mappoint(map_manager, kpid)
            mp ≡ nothing && continue
            is_bad!(mp) && (push!(bad_keypoints, kpid); continue)

            mp_order_id = length(map_points) + 1
            mp_position = get_position(mp)
            map_points[kpid] = (mp_order_id, mp_position)
            push!(points_remap, kpid)

            # For each observer, add observation: px ← (mp, pose).
            for ob_kfid in get_observers(mp)
                ob_kfid > frame.kfid && continue

                ob_frame = get_keyframe(map_manager, ob_kfid)
                if ob_frame ≡ nothing
                    remove_mappoint_obs!(map_manager, kpid, ob_kfid)
                    continue
                end
                ob_pixel = get_keypoint_unpx(ob_frame, kpid)
                if ob_pixel ≡ nothing
                    remove_mappoint_obs!(map_manager, kpid, ob_kfid)
                    continue
                end

                in_processed = ob_kfid in keys(poses)
                in_constants = ob_kfid in constant_poses
                in_covmap = ob_kfid in keys(covisibility_map)

                is_constant = ob_kfid == 0 || in_constants || !in_covmap
                !is_constant && in_covmap && (is_constant =
                    covisibility_map[ob_kfid] < min_cov_score;)

                if in_processed
                    pose_order_id, ob_pose = poses[ob_kfid]
                else
                    ob_pose = get_cw_ba(ob_frame)
                    pose_order_id = length(poses) + 1
                    poses[ob_kfid] = (pose_order_id, ob_pose)

                    push!(poses_remap, ob_kfid)
                    is_constant && push!(constant_poses, ob_kfid)
                end

                push!(observations, Observation(
                    ob_pixel, mp_position, ob_pose,
                    mp_order_id, pose_order_id,
                    is_constant, in_covmap, ob_kfid, kpid))
            end
        end
    end

    n_observations = length(observations)
    n_poses, n_points = length(poses), length(map_points)
    point_shift = n_poses * 6

    θ = Vector{Float64}(undef, point_shift + n_points * 3)
    θconst = Vector{Bool}(undef, n_poses)
    poses_ids = Vector{Int64}(undef, n_observations)
    points_ids = Vector{Int64}(undef, n_observations)
    pixels = Matrix{Float64}(undef, 2, n_observations)

    processed_poses = fill(false, n_poses)
    processed_points = fill(false, n_points)

    for (oi, observation) in enumerate(observations)
        pixels[:, oi] .= observation.pixel
        poses_ids[oi] = observation.pose_order
        points_ids[oi] = observation.point_order

        if !processed_poses[observation.pose_order]
            processed_poses[observation.pose_order] = true
            p = (observation.pose_order - 1) * 6

            θ[(p + 1):(p + 6)] .= observation.pose
            θconst[observation.pose_order] = observation.constant
        end
        if !processed_points[observation.point_order]
            processed_points[observation.point_order] = true
            p = point_shift + (observation.point_order - 1) * 3
            θ[(p + 1):(p + 3)] .= observation.point
        end
    end

    LocalBACache(
        observations, bad_keypoints, θ, θconst, pixels,
        poses_ids, points_ids, poses_remap, points_remap)
end

function _update_ba_parameters!(
    map_manager::MapManager, cache::LocalBACache, current_kfid, slam_io,
)
    points_shift = length(cache.poses_remap) * 6

    for (i, kfid) in enumerate(cache.poses_remap)
        p = (i - 1) * 6
        kf = get_keyframe(map_manager, kfid)
        set_cw_ba!(kf, @view(cache.θ[(p + 1):(p + 6)]), slam_io)
    end

    for i in 1:length(cache.observations)
        cache.outliers[i] || continue

        obs = cache.observations[i]
        obs.in_covmap && remove_mappoint_obs!(map_manager, obs.mpid, obs.kfid)
        obs.kfid == current_kfid &&
            remove_obs_from_current_frame!(map_manager, obs.mpid)

        push!(cache.bad_keypoints, obs.mpid)
    end

    for (i, mpid) in enumerate(cache.points_remap)
        mp = get_mappoint(map_manager, mpid)
        if is_bad!(mp)
            remove_mappoint!(map_manager, mpid)
            mpid in cache.bad_keypoints && pop!(cache.bad_keypoints, mpid)
        else
            p = points_shift + (i - 1) * 3
            set_position!(mp, @view(cache.θ[(p + 1):(p + 3)]))
        end
    end

    for bad_kpid in cache.bad_keypoints
        mp = get_mappoint(map_manager, bad_kpid)
        mp ≡ nothing && continue
        is_bad!(mp) && remove_mappoint!(map_manager, mp)
    end
end

"""
    local_bundle_adjustment!(estimator::Estimator, new_frame::Frame)

Perform Bundle-Adjustment on the new frame and its covisibility graph.

Minimize error function over all KeyFrame's extrinsic parameters
in the covisibility graph and their corresponding MapPoint's positions.
Afterwards, update these parameters.
"""
function local_bundle_adjustment!(estimator::Estimator, new_frame::Frame)
    if new_frame.nb_3d_kpts < estimator.params.min_cov_score
        @warn "[ES] Not enough 3D keypoints for BA: $(new_frame.nb_3d_kpts)."
        return
    end

    estimator.params.local_ba_on = true
    covisibility_map = get_covisible_map(new_frame)
    covisibility_map[new_frame.kfid] = new_frame.nb_3d_kpts

    # Get up to 5 latest KeyFrames.
    co_kfids = sort!(collect(keys(covisibility_map)); rev=true)
    co_kfids = co_kfids[1:min(5, length(co_kfids))]
    covisibility_map = OrderedDict{Int64, Int64}(
        kfid => covisibility_map[kfid] for kfid in co_kfids)

    cache = _get_ba_parameters(
        estimator.map_manager, new_frame, covisibility_map,
        estimator.params.min_cov_score)
    bundle_adjustment!(cache, new_frame.camera; show_trace=false)

    lock(estimator.map_manager.map_lock)
    try
        _update_ba_parameters!(
            estimator.map_manager, cache, new_frame.kfid, estimator.slam_io)
    catch e
        showerror(stdout, e); println()
        display(stacktrace(catch_backtrace())); println()
    finally
        unlock(estimator.map_manager.map_lock)
    end

    estimator.params.local_ba_on = false
end

"""
    map_filtering!(estimator::Estimator, new_keyframe::Frame)

Filter out KeyFrames that share too many MapPoints with other KeyFrames
in the covisibility graph. Since they are not informative.
"""
function map_filtering!(estimator::Estimator, new_keyframe::Frame)
    estimator.params.filtering_ratio ≥ 1 && return
    new_keyframe.kfid < 20 && return

    n_removed = 0
    for kfid in keys(get_covisible_map(new_keyframe))
        estimator.new_kf_available && break
        kfid == 0 && break
        kfid ≥ new_keyframe.kfid && continue

        if !has_keyframe(estimator.map_manager, kfid)
            remove_covisible_kf!(new_keyframe, kfid)
            continue
        end

        kf = get_keyframe(estimator.map_manager, kfid)
        if kf.nb_3d_kpts < estimator.params.min_cov_score ÷ 2
            lock(estimator.map_manager.map_lock) do
                remove_keyframe!(estimator.map_manager, kfid)
            end
            n_removed += 1
            continue
        end

        n_good, n_total = 0, 0
        for kp in get_3d_keypoints(kf)
            if !(kp.id in keys(estimator.map_manager.map_points))
                remove_mappoint_obs!(estimator.map_manager, kp.id, kfid)
                continue
            end
            mp = get_mappoint(estimator.map_manager, kp.id)
            mp ≡ nothing && continue

            get_observers_number(mp) > 4 && (n_good += 1;)
            n_total += 1;

            estimator.new_kf_available && break
        end

        ratio = n_good / n_total
        if ratio > estimator.params.filtering_ratio
            lock(estimator.map_manager.map_lock) do
                remove_keyframe!(estimator.map_manager, kfid)
            end
            n_removed += 1
        end
    end
    @debug "[ES] Removed $n_removed KeyFrames."
end

"""
    reset!(estimator::Estimator)

Reset Estimator.
"""
function reset!(estimator::Estimator)
    lock(estimator.queue_lock) do
        estimator.new_kf_available = false
        empty!(estimator.frame_queue)
    end
end


text:
import RowEchelon

#=
Homework 3
Section 1.7
Problem 5
=#

# Create the matrix
A = [0 -3 9 ;
	 2 1 -7 ;
	 -1 4 -5 ;
	 1 -4 -2]

# Row Reduce the matrix
A_rref = RowEchelon.rref(A)

# Check the diagonal of the matrix for all 1s.
index = 1
missingOne = false
while index <= size(A_rref)[2]
	global index
	global missingOne
	if A_rref[index,index] != 1
		missingOne = true
	end
	index += 1
end

if missingOne
	# Then the columns are not
	println("The columns of the matrix are dependent")
else
	# Then the columns are independent
	println("The columns of the matrix are independent")
end

ext: jl


text:
export DcSpfMcmcSampler

struct DcSpfMcmcSampler
    ν::Int
    desired_samples::Int
    burn_in::Int
    lag::Int
    par_info::Dict
    curr::Vector{Int} # Storage for current val
    prop::Vector{Int} # Storage for proposed val
    ind_del::Vector{Int} # Storage for indexing of deletions
    ind_add::Vector{Int} # Storage for indexing of additions
    vals::Vector{Int} # Storgae for new values to insert
    function DcSpfMcmcSampler(
        ;ν=4, desired_samples=1000, burn_in=0, lag=1
    )
    # req_samples = burn_in + 1 + (desired_samples - 1) * lag
    curr = Int[]
    prop = Int[]
    ind_del = zeros(Int, ν)
    ind_add = zeros(Int, ν)
    vals = zeros(Int, ν)
    par_info = Dict()
    par_info[:ν] = "(maximum number of edit operations)"
    new(
        ν,
        desired_samples, burn_in, lag, par_info,
        curr, prop, ind_del, ind_add, vals
        )
    end
end

function rand_flip(x::Int, V::UnitRange)
    tmp = rand(1:(length(V)-1))
    if tmp >= x
        return tmp+1
    else
        return tmp
    end
end

function rand_flip!(x::Path{Int}, i::Int, V::UnitRange)
    tmp = rand(1:(length(V)-1))
    if tmp >= x[i]
        x[i] = tmp+1
    else
        x[i] = tmp
    end
end

function flip_accept_reject!(
    I_curr::Path{Int},
    I_prop::Path{Int},
    model::DcSPF,
    mcmc::DcSpfMcmcSampler
    )

    V = model.V
    n = length(I_prop)
    n_flip = rand(1:min(mcmc.ν, n))
    k = n_flip
    i = 0
    while k > 0
        u = rand()
        q = (n - k) / n
        while q > u  # skip
            i += 1
            n -= 1
            q *= (n - k) / n
        end
        i+=1
        # i is now index to flip

        I_prop[i] = rand_flip!(I_prop, i, V)
        n -= 1
        k -= 1
    end

    d, I_mode, γ = (model.d, model.mode, model.γ)
    log_α = γ * (d(I_curr, I_mode) - d(I_prop, I_mode))
    if exp(rand()) < log_α
        # Accept
        copy!(I_curr, I_prop)
        return 1
    else
        # Reject
        copy!(I_prop, I_curr)
        return 0
    end
end

function rand_resize!(
    x::Path, m::Int, V::UnitRange
    )
    n = length(x)
    if n == m
        return x
    elseif m < n
        rand_delete!(x, n-m)
    else
        rand_insert!(x, m-n, V)
    end
end

function draw_sample!(
    sample_out::Vector{Path{Int}},
    mcmc::DcSpfMcmcSampler,
    model::DcSPF;
    lag::Int=mcmc.lag,
    init::Path=model.mode
    )

    p = model.p # length distributuion
    V = model.V
    I_curr = copy(init)
    I_prop = copy(init)

    acc_count = 0
    count = 0
    for i in eachindex(sample_out)
        n = rand(p) # Sample length

        # Resize the paths
        rand_resize!(I_curr, n, V)
        copy!(I_prop, I_curr)

        # Now do lag number of mcmc samples and store last
        for j in 1:lag
            acc_count += flip_accept_reject!(
                I_curr, I_prop, model, mcmc
            )
            count += 1
        end
        copy!(sample_out[i], I_curr)

    end
    return acc_count/count
end

function draw_sample(
    mcmc::DcSpfMcmcSampler,
    model::DcSPF;
    lag::Int=mcmc.lag,
    desired_samples::Int=mcmc.desired_samples,
    init::Path=model.mode
    )
    sample_out = [Int[] for i in 1:desired_samples]
    a = draw_sample!(sample_out, mcmc, model, lag=lag)
    println("Acceptance prob: $a")
    return sample_out
end


text:
# ---
# title: 139. Word Break
# id: problem139
# author: zhwang
# date: 2022-02-18
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/word-break/description/>
# hidden: true
# ---
#
# Given a **non-empty** string _s_ and a dictionary _wordDict_ containing a list
# of **non-empty** words, determine if _s_ can be segmented into a space-
# separated sequence of one or more dictionary words.
#
# **Note:**
#
#   * The same word in the dictionary may be reused multiple times in the segmentation.
#   * You may assume the dictionary does not contain duplicate words.
#
# **Example 1:**
#
#
#
#     Input: s = "leetcode", wordDict = ["leet", "code"]
#     Output: true
#     Explanation: Return true because "leetcode" can be segmented as "leet code".
#
#
# **Example 2:**
#
#
#
#     Input: s = "applepenapple", wordDict = ["apple", "pen"]
#     Output: true
#     Explanation: Return true because "applepenapple" can be segmented as "apple pen apple".
#                  Note that you are allowed to reuse a dictionary word.
#
#
# **Example 3:**
#
#
#
#     Input: s = "catsandog", wordDict = ["cats", "dog", "sand", "and", "cat"]
#     Output: false
#
#
#
## @lc code=start
using LeetCode

## BFS
function word_break_bfs(s::String, word_dict::Vector{String})::Bool
    n, word_dict = length(s), sort!(word_dict; by=x -> length(x))
    valids, valid_pos = fill(false, n), [0]
    while !isempty(valid_pos)
        pos = popfirst!(valid_pos)
        for word in word_dict
            (new_pos = pos + length(word)) > n && break
            new_pos == n && s[(pos + 1):new_pos] == word && return true
            if !valids[new_pos] && s[(pos + 1):new_pos] == word
                valids[new_pos] = true
                push!(valid_pos, new_pos)
            end
        end
    end
    return false
end

## Dynamic Programming
function word_break(s::String, word_dict::Vector{String})::Bool
    lens = sort!(unique!(length.(word_dict)))
    dp = append!([true], fill(false, length(s)))
    for i in eachindex(dp)
        for len in lens
            i > len || break
            dp[i - len] && s[(i - len):(i - 1)] ∈ word_dict && (dp[i] = true)
        end
    end
    return last(dp)
end
## @lc code=end
## @lc test=start
@testset "139.word-break.jl" begin
    tcase1 = ("leetcode", ["leet", "code"])
    tcase2 = ("applepenapple", ["apple", "pen"])
    tcase3 = ("bb", ["a", "b", "bbb", "bbbb"])
    fcase1 = ("catsandog", ["cats", "dog", "sand", "and", "cat"])
    fcase2 = ("a", ["b"])
    fcase3 = ("a"^150 * "b", ["a"^i for i in 1:10])
    @test word_break(tcase1...) && word_break_bfs(tcase1...)
    @test word_break(tcase2...) && word_break_bfs(tcase2...)
    @test word_break(tcase3...) && word_break_bfs(tcase3...)
    @test !word_break(fcase1...) && !word_break_bfs(fcase1...)
    @test !word_break(fcase2...) && !word_break_bfs(fcase2...)
    @test !word_break(fcase3...) && !word_break_bfs(fcase3...)
end
## @lc test=end

## @lc info=start
# link: [solution to 139](https://leetcode-cn.com/problems/word-break/solution/139-dan-ci-chai-fen-bfs-sou-suo-by-rex-r-jgjc/)
## @lc info=end


text:
macro gdtrace()
    quote
        if tracing
            dt = Dict()
            if o.extended_trace
                dt["x"] = copy(x)
                dt["g(x)"] = copy(g)
            end
            g_norm = vecnorm(g, Inf)
            update!(tr,
                    iteration,
                    f_x,
                    g_norm,
                    dt,
                    o.store_trace,
                    o.show_trace,
                    o.show_every,
                    o.callback)
        end
    end
end

immutable GradientDescent{T} <: Optimizer
    linesearch!::Function
    P::T
    precondprep!::Function
end

GradientDescent(; linesearch!::Function = hz_linesearch!,
                P = nothing, precondprep! = (P, x) -> nothing) =
                    GradientDescent(linesearch!, P, precondprep!)

function optimize{T}(d::DifferentiableFunction,
                     initial_x::Array{T},
                     mo::GradientDescent,
                     o::OptimizationOptions)
    # Print header if show_trace is set
    print_header(o)

    # Maintain current state in x and previous state in x_previous
    x, x_previous = copy(initial_x), copy(initial_x)

    # Count the total number of iterations
    iteration = 0

    # Track calls to function and gradient
    f_calls, g_calls = 0, 0

    # Count number of parameters
    n = length(x)

    # Maintain current gradient in g
    g = similar(x)

    # The current search direction
    s = similar(x)

    # Buffers for use in line search
    x_ls = similar(x)
    g_ls = similar(x)

    # Store f(x) in f_x
    f_x_previous, f_x = NaN, d.fg!(x, g)
    f_calls, g_calls = f_calls + 1, g_calls + 1

    # Keep track of step-sizes
    alpha = alphainit(one(T), x, g, f_x)

    # TODO: How should this flag be set?
    mayterminate = false

    # Maintain a cache for line search results
    lsr = LineSearchResults(T)

    # Trace the history of states visited
    tr = OptimizationTrace{typeof(mo)}()
    tracing = o.store_trace || o.show_trace || o.extended_trace || o.callback != nothing
    @gdtrace

    # Assess multiple types of convergence
    x_converged, f_converged = false, false
    g_converged = vecnorm(g, Inf) < o.g_tol

    # Iterate until convergence
    converged = g_converged
    while !converged && iteration < o.iterations
        # Increment the number of steps we've had to perform
        iteration += 1

        # Search direction is always the negative preconditioned gradient
        mo.precondprep!(mo.P, x)
        A_ldiv_B!(s, mo.P, g)
        @simd for i in 1:n
            @inbounds s[i] = -s[i]
        end

        # Refresh the line search cache
        dphi0 = vecdot(g, s)
        clear!(lsr)
        push!(lsr, zero(T), f_x, dphi0)

        # Determine the distance of movement along the search line
        alpha, f_update, g_update =
          mo.linesearch!(d, x, s, x_ls, g_ls, lsr, alpha, mayterminate)
        f_calls, g_calls = f_calls + f_update, g_calls + g_update

        # Maintain a record of previous position
        copy!(x_previous, x)

        # Update current position # x = x + alpha * s
        LinAlg.axpy!(alpha, s, x)

        # Update the function value and gradient
        f_x_previous, f_x = f_x, d.fg!(x, g)
        f_calls, g_calls = f_calls + 1, g_calls + 1

        x_converged,
        f_converged,
        g_converged,
        converged = assess_convergence(x,
                                       x_previous,
                                       f_x,
                                       f_x_previous,
                                       g,
                                       o.x_tol,
                                       o.f_tol,
                                       o.g_tol)

        @gdtrace
    end

    return MultivariateOptimizationResults("Gradient Descent",
                                           initial_x,
                                           x,
                                           Float64(f_x),
                                           iteration,
                                           iteration == o.iterations,
                                           x_converged,
                                           o.x_tol,
                                           f_converged,
                                           o.f_tol,
                                           g_converged,
                                           o.g_tol,
                                           tr,
                                           f_calls,
                                           g_calls)
end


text:
mev(Σ::Matrix{T}, c, ls::Vector{Bool}) where T <: AbstractFloat = det(reduceband(Σ, ls))

function mormbased(Σ::Matrix{T}, c::Array{T, N}, ls::Vector{Bool}) where {T <: AbstractFloat, N}
  norm(reduceband(c, ls))/norm(reduceband(Σ, ls))^(N/2)
end

"""
  detoverdetfitfunction(a::Array{N}, b::Array{N})

computes the maximizing function det(C_n)/det(C_2)^(n/2). It assumes, that product
of singular values from HOSVD of tensor is a good approximation of hyperdeterminant
of the tensor (whatever that means).
Returns the value of the maximizin function
"""
function hosvdapprox(Σ::Matrix{T}, c::Array{T,N}, fibres::Vector{Bool} = [fill(true, size(Σ, 1))...]) where {T <: AbstractFloat, N}
  c = reduceband(c, fibres)
  Σ = reduceband(Σ, fibres)
  cunf = unfoldsym(c)
  eigc = abs.(eigvals(cunf*cunf'))
  eigΣ = abs.(eigvals(Σ*Σ'))
  sum(log.(eigc)-N/2*log.(eigΣ))/2
end

"""
  reduceband(ar::Array{N}, k::Vector{Bool})

Returns n-array without values at indices in ind
```jldoctest
julia>  reshape(collect(1.:27.),(3,3,3))
3×3×3 Array{Float64,3}:
[:, :, 1] =
 1.0  4.0  7.0
 2.0  5.0  8.0
 3.0  6.0  9.0

[:, :, 2] =
 10.0  13.0  16.0
 11.0  14.0  17.0
 12.0  15.0  18.0

[:, :, 3] =
 19.0  22.0  25.0
 20.0  23.0  26.0
 21.0  24.0  27.0

julia> reduceband(reshape(collect(1.:27.),(3,3,3)), [true, false, false])
1×1×1 Array{Float64,3}:
[:, :, 1] =
 1.0
```
TODO reimplement in blocks
"""
reduceband(ar::Array{T,N}, fibres::Vector{Bool}) where {T <: AbstractFloat, N} =
  ar[fill(fibres, N)...]

"""
  function unfoldsym{T <: Real, N}(ar::Array{T,N})

Returns a matrix of size (i, k^(N-1)) that is an unfold of symmetric array ar
"""
function unfoldsym(ar::Array{T,N}) where {T <: AbstractFloat, N}
  i = size(ar, 1)
  return reshape(ar, i, i^(N-1))
end

"""
TODO reimplement in blocks
"""
function unfoldsym(t::SymmetricTensor{T, N}) where {T <: AbstractFloat, N}
   t = unfoldsym(Array(t))
   t*t'
end

#greedy algorithm

"""
  greedestep(c::Vector{Array{Float}}, maxfunction::Function, ls::Vector{Bool})

Returns vector of bools that determines bands that maximise a function. True means include
a band, false exclude a band. It changes one true to false in input ls

```jldoctest
julia> a = reshape(collect(1.:9.), 3,3);

julia> b = reshape(collect(1.: 27.), 3,3,3);

julia> testf(ar,bool)= det(ar[1][bool,bool])

julia> greedestep(ar, testf, [true, true, true])
3-element Array{Bool,1}:
true
true
false
```
"""
function greedestep(Σ::Matrix{T}, c::Array{T, N}, maxfunction::Function,
                    ls::Vector{Bool}) where {T <: AbstractFloat, N}
  inds = findall(ls)
  bestval = SharedArray{T}(length(ls))
  bestval .= -Inf
  bestls = copy(ls)
  @sync @distributed for i in inds
    templs = copy(ls)
    templs[i] = false
    bestval[i] = maxfunction(Σ, c, templs)
  end
  v, i = findmax(bestval)
  bestls[i] = false
  return bestls, v, i
end

"""
  greedesearchdata(Σ::SymmetricTensor{T,2}, c::SymmetricTensor{T, N}, maxfunction::Function, k::Int)

returns array of bools that are non-outliers features
"""
function greedesearchdata(Σ::SymmetricTensor{T,2}, c::SymmetricTensor{T, N}, maxfunction::Function, k::Int) where {T <: AbstractFloat, N}
  ls =  [true for i=1:Σ.dats]
  Σ = Array(Σ)
  c = Array(c)
  result = []
  for i = 1:k
    ls, value, j = greedestep(Σ, c, maxfunction, ls)
    push!(result, (ls,value,j))
    value != -Inf || throw(AssertionError(" for k = $(k) optimisation does not work"))
  end
  result
end

"""
  function cumfsel(Σ::SymmetricTensor{T,2}, c::SymmetricTensor{T, N}, f::String, k::Int = Σ.dats) where {T <: AbstractFloat, N}

Returns an Array of tuples (ind::Array{Bool}, fval::Float64, i::Int). Given
k-th Array ind are marginals removed after k -steps as those with low N'th order
dependency, fval, the value of the target function at step k and i, a feature removed
at step k.

Uses Σ - the covariance matrix and c - the N'th cumulant tensor to measure the
N'th order dependencies between marginals.
Function f is the optimization function, ["hosvd", "norm", "mev"] are supported.

```jldoctest

julia> Random.seed!(42);

julia> x = rand(12,10);

julia> c = cumulants(x, 4);

julia> cumfsel(c[2], c[4], "hosvd")
10-element Array{Any,1}:
 (Bool[true, true, true, false, true, true, true, true, true, true], 27.2519, 4)
 (Bool[true, true, false, false, true, true, true, true, true, true], 22.6659, 3)
 (Bool[true, true, false, false, false, true, true, true, true, true], 18.1387, 5)
 (Bool[false, true, false, false, false, true, true, true, true, true], 14.4492, 1)
 (Bool[false, true, false, false, false, true, true, false, true, true], 11.2086, 8)
 (Bool[false, true, false, false, false, true, true, false, true, false], 7.84083, 10)
 (Bool[false, false, false, false, false, true, true, false, true, false], 5.15192, 2)
 (Bool[false, false, false, false, false, false, true, false, true, false], 2.56748, 6)
 (Bool[false, false, false, false, false, false, true, false, false, false], 0.30936, 9)
 (Bool[false, false, false, false, false, false, false, false, false, false], 0.0, 7)
```
"""
function cumfsel(Σ::SymmetricTensor{T,2}, c::SymmetricTensor{T, N}, f::String, k::Int = Σ.dats) where {T <: AbstractFloat, N}
  if f == "hosvd"
    return greedesearchdata(Σ, c, hosvdapprox, k)
  elseif f == "norm"
    return greedesearchdata(Σ, c, mormbased, k)
  elseif f == "mev"
    return greedesearchdata(Σ, c ,mev, k)
  end
  throw(AssertionError("$(f) not supported use hosvd, norm or mev"))
end

"""
  cumfsel(Σ::Matrix{T}, k::Int = size(Σ, 1))

cumfsel that uses as default the mev method
"""
cumfsel(Σ::SymmetricTensor{T,2}, k::Int = Σ.dats) where T <: AbstractFloat = cumfsel(Σ, SymmetricTensor(convert(Array{T}, ones(2,2,2))), "mev", k)


text:
"""
    search(graph, queries, n_neighbors; max_candidates) -> indices, distances

Search the kNN `graph` for the nearest neighbors of the points in `queries`.
`max_candidates` controls how large the candidate queue should be (min `n_neighbors`);
larger values increase accuracy at the cost of speed.
"""
function search(graph::G,
                queries::AbstractVector,
                n_neighbors::Integer;
                max_candidates=max(n_neighbors, 20),
                ) where {V, U, G <: ApproximateKNNGraph{V, U}}

    length(queries) ≥ 1 || error("queries must have at least 1 element")
    n_neighbors ≥ 1 || error("n_neighbors must be at least 1")
    max_candidates ≥ 5 || error("max_candidates must be at least 5")

    data = graph.data
    metric = graph.metric
    # lists of candidates, sorted by distance
    candidates = [BinaryMaxHeap{Tuple{U, V, Bool}}() for _ in 1:length(queries)]
    # a set of seen candidates per thread
    seen_sets = [BitVector(undef, length(data)) for _ in 1:Threads.nthreads()]
    Threads.@threads for i in eachindex(queries)
        # zero out seen
        seen = seen_sets[Threads.threadid()]
        seen .= false
        # initialize with random
        init_candidates!(candidates[i], seen, graph, queries[i], max_candidates)
        while true
            next_candidate = get_next_candidate!(candidates[i])
            if isnothing(next_candidate)
                break
            end

            for v in outneighbors(graph, next_candidate[2])
                if !seen[v]
                    dist = evaluate(metric, queries[i], data[v])
                    if dist ≤ top(candidates[i])[1]
                        pop!(candidates[i]) # pop maximum
                        push!(candidates[i], (dist, v, false))
                    end
                    seen[v] = true
                end
            end
        end
    end

    return deheap_knns(candidates, n_neighbors)
end

function search(graph::ApproximateKNNGraph,
                queries::AbstractMatrix,
                n_neighbors::Integer;
                max_candidates=max(n_neighbors, 20),
                )
    query_cols = collect(eachcol(queries))
    return search(graph, query_cols, n_neighbors; max_candidates=max_candidates)
end

function get_next_candidate!(candidates::BinaryMaxHeap{Tuple{U, V, Bool}}) where {U <: Real, V}
    min_idx = -1
    min_dist = typemax(U)
    for (i, t) in enumerate(candidates.valtree)
        if t[1] < min_dist && !t[3]
            min_idx = i
            min_dist = t[1]
        end
    end
    if min_idx != -1 # found an unvisited candidate
        dist, node, _ = candidates.valtree[min_idx]
        cand = (dist, node, true)
        candidates.valtree[min_idx] = cand # mark visited
        return cand
    end
    return nothing
end

function init_candidates!(candidates, seen, graph, query, max_candidates)
    for v in KNNGraphs.sample_neighbors(nv(graph), max_candidates)
        dist = evaluate(graph.metric, query, graph.data[v])
        push!(candidates, (dist, v, false))
        seen[v] = true
    end
    return candidates
end

"""
Remove the `k` nearest neighbors from each heap in `knn_heaps`.
Return two k x length(knn_heaps) arrays for the indices and
distances to each point's kNN.
"""
function deheap_knns(heaps::Vector{BinaryMaxHeap{Tuple{U, V, Bool}}}, k) where {U, V}

    ids = Array{V}(undef, (k, length(heaps)))
    dists = Array{U}(undef, (k, length(heaps)))

    for i in 1:length(heaps)
        len = length(heaps[i])
        for j in 1:len
            # NOTE: these are max heaps, so we only want the last k
            node_dist, node_idx, _ = pop!(heaps[i])
            neighbor_idx = 1 + len - j
            if neighbor_idx <= k
                ids[neighbor_idx, i] = node_idx
                dists[neighbor_idx, i] = node_dist
            end
        end
    end
    return ids, dists
end


text:
push!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))

using Test
using StaticArrays, IntervalSets, LinearAlgebra, UnPack

import ClimaCore: ClimaCore, slab, Spaces, Domains, Meshes, Geometry, Topologies, Spaces, Fields, Operators
using ClimaCore.Geometry

using Logging: global_logger
using TerminalLoggers: TerminalLogger
global_logger(TerminalLogger())

# set up function space
function hvspace_2D(xlim = (-π, π), zlim = (0, 4π), helem = 20, velem = 20, npoly = 1)
    FT = Float64
    vertdomain = Domains.IntervalDomain(
        Geometry.ZPoint{FT}(zlim[1]),
        Geometry.ZPoint{FT}(zlim[2]);
        boundary_tags = (:bottom, :top),
    )
    vertmesh = Meshes.IntervalMesh(vertdomain, nelems = velem)
    vert_center_space = Spaces.CenterFiniteDifferenceSpace(vertmesh)

    horzdomain = Domains.IntervalDomain(Geometry.XPoint{FT}(xlim[1]) .. Geometry.XPoint{FT}(xlim[2]), periodic = true)
    horzmesh = Meshes.IntervalMesh(horzdomain; nelems = helem)
    horztopology = Topologies.IntervalTopology(horzmesh)

    quad = Spaces.Quadratures.GLL{npoly + 1}()
    horzspace = Spaces.SpectralElementSpace1D(horztopology, quad)

    hv_center_space = Spaces.ExtrudedFiniteDifferenceSpace(horzspace, vert_center_space)
    hv_face_space = Spaces.FaceExtrudedFiniteDifferenceSpace(hv_center_space)
    return (hv_center_space, hv_face_space)
end

function pressure(ρθ)
    if ρθ >= 0
        return MSLP * (R_d * ρθ / MSLP)^γ
    else
        return NaN
    end
end

Φ(z) = grav * z

abstract type BCtag end
struct ZeroFlux <: BCtag end

bc_divF2C_bottom!(::ZeroFlux, dY, Y, p, t) = Operators.SetValue(Geometry.WVector(0.0))
bc_divF2C_top!(::ZeroFlux, dY, Y, p, t) = Operators.SetValue(Geometry.WVector(0.0))

function init_bubble_2d(x, z)
    θ₀ = atm_T_ini
    cp_d = C_p
    cv_d = C_v
    p₀ = MSLP
    g = grav
    γ = cp_d / cv_d
    x_c = 0.0
    z_c = 350.0
    r_c = 250.0
    θ_b = atm_T_ini
    θ_c = 0.5

    # auxiliary quantities
    r = sqrt((x - x_c)^2 + (z - z_c)^2)
    θ_p = r < r_c ? 0.5 * θ_c * (1.0 + cospi(r / r_c)) : 0.0 # potential temperature perturbation

    θ = θ_b + θ_p # potential temperature
    π_exn = 1.0 - g * z / cp_d / θ # exner function
    T = π_exn * θ # temperature
    p = p₀ * π_exn^(cp_d / R_d) # pressure
    ρ = p / R_d / T # density
    ρθ = ρ * θ # potential temperature density

    return (ρ = ρ, ρθ = ρθ, ρuₕ = ρ * Geometry.UVector(0.0))
end

function init_sea_breeze_2d(x, z)
    θ₀ = atm_T_ini
    cp_d = C_p
    cv_d = C_v
    p₀ = MSLP
    g = grav
    γ = cp_d / cv_d
    z_c = 100.0
    θ_b = atm_T_ini
    θ_p = z < z_c ? rand() - 0.5 : 0.0 # potential temperature perturbation
    θ = θ_b + θ_p # potential temperature
    π_exn = 1.0 - g * z / cp_d / θ # exner function
    T = π_exn * θ # temperature
    p = p₀ * π_exn^(cp_d / R_d) # pressure
    ρ = p / R_d / T # density
    ρθ = ρ * θ # potential temperature density
    return (ρ = ρ, ρθ = ρθ, ρuₕ = ρ * Geometry.UVector(0.0))
end

function atm_rhs!(dY, Y, params, t)
    ρw = Y.ρw
    Yc = Y.Yc
    dYc = dY.Yc
    dρw = dY.ρw

    center_coords = Fields.coordinate_field(params.domain.hv_center_space)
    face_coords = Fields.coordinate_field(params.domain.hv_face_space)

    # spectral horizontal operators
    hdiv = Operators.Divergence()
    hgrad = Operators.Gradient()
    hwdiv = Operators.WeakDivergence()
    hwgrad = Operators.WeakGradient()

    # vertical FD operators with BC's
    vdivf2c = Operators.DivergenceF2C(
        bottom = Operators.SetValue(Geometry.WVector(0.0)),
        top = Operators.SetValue(Geometry.WVector(0.0)),
    )
    vvdivc2f = Operators.DivergenceC2F(
        bottom = Operators.SetDivergence(Geometry.WVector(0.0)),
        top = Operators.SetDivergence(Geometry.WVector(0.0)),
    )
    uvdivf2c = Operators.DivergenceF2C(
        bottom = Operators.SetValue(Geometry.WVector(0.0) ⊗ Geometry.UVector(0.0)),
        top = Operators.SetValue(Geometry.WVector(0.0) ⊗ Geometry.UVector(0.0)),
    )
    If = Operators.InterpolateC2F(bottom = Operators.Extrapolate(), top = Operators.Extrapolate())
    Ic = Operators.InterpolateF2C()
    ∂ = Operators.DivergenceF2C(
        bottom = Operators.SetValue(Geometry.WVector(0.0)),
        top = Operators.SetValue(Geometry.WVector(0.0)),
    )
    ∂f = Operators.GradientC2F()
    ∂c = Operators.GradientF2C()
    B = Operators.SetBoundaryOperator(
        bottom = Operators.SetValue(Geometry.WVector(0.0)),
        top = Operators.SetValue(Geometry.WVector(0.0)),
    )

    fcc = Operators.FluxCorrectionC2C(bottom = Operators.Extrapolate(), top = Operators.Extrapolate())
    fcf = Operators.FluxCorrectionF2F(bottom = Operators.Extrapolate(), top = Operators.Extrapolate())
    ∇_z_ρθ = Operators.DivergenceF2C(
        bottom = bc_divF2C_bottom!(params.bc.ρθ.bottom, dY, Y, params, t),
        top = bc_divF2C_top!(params.bc.ρθ.top, dY, Y, params, t),
    )

    uₕ = @. Yc.ρuₕ / Yc.ρ
    w = @. ρw / If(Yc.ρ)
    wc = @. Ic(ρw) / Yc.ρ
    p = @. pressure(Yc.ρθ)
    θ = @. Yc.ρθ / Yc.ρ
    Yfρ = @. If(Yc.ρ)

    ### HYPERVISCOSITY
    # 1) compute hyperviscosity coefficients
    @. dYc.ρθ = hwdiv(hgrad(θ))
    @. dYc.ρuₕ = hwdiv(hgrad(uₕ))
    @. dρw = hwdiv(hgrad(w))
    Spaces.weighted_dss!(dYc)
    Spaces.weighted_dss!(dρw)

    κ₄ = 0.0 # m^4/s
    @. dYc.ρθ = -κ₄ * hwdiv(Yc.ρ * hgrad(dYc.ρθ))
    @. dYc.ρuₕ = -κ₄ * hwdiv(Yc.ρ * hgrad(dYc.ρuₕ))
    @. dρw = -κ₄ * hwdiv(Yfρ * hgrad(dρw))

    # density
    @. dYc.ρ = -∂(ρw)
    @. dYc.ρ -= hdiv(Yc.ρuₕ)

    # potential temperature
    @. dYc.ρθ += -(∇_z_ρθ(ρw * If(Yc.ρθ / Yc.ρ)))
    @. dYc.ρθ -= hdiv(uₕ * Yc.ρθ)

    # horizontal momentum
    Ih = Ref(Geometry.Axis2Tensor((Geometry.UAxis(), Geometry.UAxis()), @SMatrix [1.0]))
    @. dYc.ρuₕ += -uvdivf2c(ρw ⊗ If(uₕ))
    @. dYc.ρuₕ -= hdiv(Yc.ρuₕ ⊗ uₕ + p * Ih)

    # vertical momentum
    @. dρw +=
        B(Geometry.transform(Geometry.WAxis(), -(∂f(p)) - If(Yc.ρ) * ∂f(Φ(center_coords.z))) - vvdivc2f(Ic(ρw ⊗ w)))
    uₕf = @. If(Yc.ρuₕ / Yc.ρ) # requires boundary conditions
    @. dρw -= hdiv(uₕf ⊗ ρw)

    ### UPWIND FLUX CORRECTION
    upwind_correction = false
    if upwind_correction
        @. dYc.ρ += fcc(w, Yc.ρ)
        @. dYc.ρθ += fcc(w, Yc.ρθ)
        @. dYc.ρuₕ += fcc(w, Yc.ρuₕ)
        @. dρw += fcf(wc, ρw)
    end

    ### DIFFUSION
    κ₂ = 5.0 # m^2/s
    #  1a) horizontal div of horizontal grad of horiz momentun
    @. dYc.ρuₕ += hwdiv(κ₂ * (Yc.ρ * hgrad(Yc.ρuₕ / Yc.ρ)))
    #  1b) vertical div of vertical grad of horiz momentun
    @. dYc.ρuₕ += uvdivf2c(κ₂ * (Yfρ * ∂f(Yc.ρuₕ / Yc.ρ)))

    #  1c) horizontal div of horizontal grad of vert momentum
    @. dρw += hwdiv(κ₂ * (Yfρ * hgrad(ρw / Yfρ)))
    #  1d) vertical div of vertical grad of vert momentun
    @. dρw += vvdivc2f(κ₂ * (Yc.ρ * ∂c(ρw / Yfρ)))

    #  2a) horizontal div of horizontal grad of potential temperature
    @. dYc.ρθ += hwdiv(κ₂ * (Yc.ρ * hgrad(Yc.ρθ / Yc.ρ)))
    #  2b) vertical div of vertial grad of potential temperature
    @. dYc.ρθ += ∇_z_ρθ(κ₂ * (Yfρ * ∂f(Yc.ρθ / Yc.ρ)))

    Spaces.weighted_dss!(dYc)
    Spaces.weighted_dss!(dρw)
    return dY
end

# init simulation
function atm_init(; xmin = -500, xmax = 500, zmin = 0, zmax = 1000, npoly = 3, helem = 20, velem = 20, bc = nothing)

    # construct domain spaces
    hv_center_space, hv_face_space = hvspace_2D((xmin, xmax), (zmin, zmax), helem, velem, npoly) # [m]
    center_coords = Fields.coordinate_field(hv_center_space)
    face_coords = Fields.coordinate_field(hv_face_space)
    domain = (hv_center_space = hv_center_space, hv_face_space = hv_face_space)

    # initialize prognostic variables
    Yc = map(center_coords) do coord
        sea_breeze = init_sea_breeze_2d(coord.x, coord.z)
        sea_breeze
    end

    ρw = map(face_coords) do coord
        Geometry.WVector(0.0)
    end

    Y = Fields.FieldVector(Yc = Yc, ρw = ρw)

    # select boundary conditions
    if bc === nothing
        bc = (
            ρθ = (bottom = CoupledFlux(), top = ZeroFlux()),
            ρu = nothing, # for now BCs are hard coded, except for ρθ
        )
    end

    return Y, bc, domain
end

using OrdinaryDiffEq
function atm_run!(Y, bc, domain)
    dYdt = similar(Y)
    params = (aux_params = 0.0, T_sfc = 1.0, bc = bc, domain = domain)
    atm_rhs!(dYdt, Y, params, 0.0)
    prob = ODEProblem(atm_rhs!, Y, (0.0, 250.0), params)
    Δt = 0.025
    sol = solve(prob, SSPRK33(), dt = Δt, saveat = 1.0, progress = true, progress_message = (dt, u, params, t) -> t)
end


text:
# ---
# title: 448. Find All Numbers Disappeared in an Array
# id: problem448
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/find-all-numbers-disappeared-in-an-array/description/>
# hidden: true
# ---
#
# Given an array of integers where 1 ≤ a[i] ≤ _n_ ( _n_ = size of array), some
# elements appear twice and others appear once.
#
# Find all the elements of [1, _n_ ] inclusive that do not appear in this array.
#
# Could you do it without extra space and in O( _n_ ) runtime? You may assume
# the returned list does not count as extra space.
#
# **Example:**
#
#
#
#     Input:
#     [4,3,2,7,8,2,3,1]
#
#     Output:
#     [5,6]
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#=
    https://github.com/JuliaArrays/EndpointRanges.jl
    by Tim Holy
    [this has minor edits for v0.7-DEV]
=#

module EndpointRanges

#using Compat

import Base: +, -, *, /, ÷, %
using Base: ViewIndex, tail, indices1

export ibegin, iend

abstract type Endpoint end
struct IBegin <: Endpoint end
struct IEnd   <: Endpoint end
const ibegin = IBegin()
const iend   = IEnd()

(::IBegin)(b::Integer, e::Integer) = b
(::IEnd  )(b::Integer, e::Integer) = e
(::IBegin)(r::AbstractRange) = first(r)
(::IEnd  )(r::AbstractRange) = last(r)

struct IndexFunction{F<:Function} <: Endpoint
    index::F
end
(f::IndexFunction)(r::AbstractRange) = f.index(r)

for op in (:+, :-)
    @eval $op(x::Endpoint) = IndexFunction(r->x(r))
end
for op in (:+, :-, :*, :/, :÷, :%)
    @eval $op(x::Endpoint, y::Endpoint) = IndexFunction(r->$op(x(r), y(r)))
    @eval $op(x::Endpoint, y::Number) = IndexFunction(r->$op(x(r), y))
    @eval $op(x::Number, y::Endpoint) = IndexFunction(r->$op(x, y(r)))
end

# deliberately not <: AbstractUnitRange{Int}
abstract type EndpointRange{T} end
struct EndpointUnitRange{F<:Union{Int,Endpoint},L<:Union{Int,Endpoint}} <: EndpointRange{Int}
    start::F
    stop::L
end
struct EndpointStepRange{F<:Union{Int,Endpoint},L<:Union{Int,Endpoint}} <: EndpointRange{Int}
    start::F
    step::Int
    stop::L
end

(r::EndpointUnitRange)(s::AbstractRange) = r.start(s):r.stop(s)
(r::EndpointUnitRange{Int,E})(s::AbstractRange) where {E<:Endpoint} = r.start:r.stop(s)
(r::EndpointUnitRange{E,Int})(s::AbstractRange) where {E<:Endpoint} = r.start(s):r.stop

(r::EndpointStepRange)(s::AbstractRange) = r.start(s):r.step:r.stop(s)
(r::EndpointStepRange{Int,E})(s::AbstractRange) where {E<:Endpoint} = r.start:r.step:r.stop(s)
(r::EndpointStepRange{E,Int})(s::AbstractRange) where {E<:Endpoint} = r.start(s):r.step:r.stop

Base.colon(start::Endpoint, stop::Endpoint) = EndpointUnitRange(start, stop)
Base.colon(start::Endpoint, stop::Int) = EndpointUnitRange(start, stop)
Base.colon(start::Int, stop::Endpoint) = EndpointUnitRange(start, stop)

Base.colon(start::Endpoint, step::Int, stop::Endpoint) = EndpointStepRange(start, step, stop)
Base.colon(start::Endpoint, step::Int, stop::Int) = EndpointStepRange(start, step, stop)
Base.colon(start::Int, step::Int, stop::Endpoint) = EndpointStepRange(start, step, stop)

function Base.getindex(r::UnitRange, s::EndpointRange)
    getindex(r, newindex(indices1(r), s))
end

function Base.getindex(r::AbstractUnitRange, s::EndpointRange)
    getindex(r, newindex(indices1(r), s))
end

function Base.getindex(r::StepRange, s::EndpointRange)
    getindex(r, newindex(indices1(r), s))
end

function Base.getindex(r::StepRangeLen, s::EndpointRange)
    getindex(r, newindex(indices1(r), s))
end

function Base.getindex(r::LinRange, s::EndpointRange)
    getindex(r, newindex(indices1(r), s))
end

@inline function Base.to_indices(A, inds, I::Tuple{Union{Endpoint, EndpointRange}, Vararg{Any}})
    (newindex(inds[1], I[1]), to_indices(A, Base._maybetail(inds), Base.tail(I))...)
end

@inline newindices(indsA, inds) = (newindex(indsA[1], inds[1]), newindices(tail(indsA), tail(inds))...)
newindices(::Tuple{}, ::Tuple{}) = ()

newindex(indA, i::Union{Real, AbstractArray, Colon}) = i
newindex(indA, i::EndpointRange) = i(indA)
newindex(indA, i::IBegin) = first(indA)
newindex(indA, i::IEnd)   = last(indA)

end # module


text:
using Distributions
using LinearAlgebra
using Random

# Define the DGP
function simulate_sample(N0, N1)
    dy1 = Beta(2, 3);
    dy2 = Beta(2, 4);
    X0 = [ones(N0) rand([0, 1], N0, 5)]
    X1 = kron([ones(2, 5) [1, 0]], ones(N1))
    y0 = [X0[i, end] == 1 ? rand(dy2) : rand(dy1) for i in 1:N0]
    y1 = LinRange(0, 1, N1) |> collect |> x -> repeat(x, 2)
    return dy1, dy2, y0, X0, y1, X1
end

# # Not run
# Random.seed!(1);
# N0, N1 = 500, 0;
# dy1, dy2, y0, X0, y1, X1 = simulate_sample(N0, N1);

# # Not run
# begin
#     using Gadfly
#     Random.seed!(1)
#     N0, N1 = 500, 50
#     dy1, dy2, y0, X0, y1, X1 = simulate_sample(N0, N1)
#     plot(
#         layer(x = y1, y = pdf.(dy1, y1), Geom.line()),
#         layer(x = y1, y = pdf.(dy2, y1), Geom.line())
#     )
# end


text:
include("make_literate.jl")

using Pkg
current_env = Base.load_path()[1]
Pkg.activate(@__DIR__)

using Morbit
using Documenter

DocMeta.setdocmeta!(Morbit, :DocTestSetup, :(using Morbit); recursive=true)

makedocs(;
    modules=[Morbit],
    authors="Manuel Berkemeier",
    repo="https://github.com/manuelbb-upb/Morbit.jl/blob/{commit}{path}#{line}",
    sitename="Morbit.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://manuelbb-upb.github.io/Morbit.jl",
        #assets=["/custom_assets/iframeResizer.min.js",],
    ),
    pages=[
        "Home" => "index.md",
        "Examples" => [
            "Two Parabolas" => "example_two_parabolas.md",
            "ZDT3" => "example_zdt.md"
        ],
        "Models" => [
            "ExactModels" => "ExactModel.md",
            "RbfModels" => "RbfModel.md",
            "TaylorModels" => "TaylorModel.md",
            "LagrangeModels" => "LagrangeModel.md"
        ],
        "Random Notebooks" => [
            "Finite Differences" => "notebook_finite_differences.md",
            "Lagrange Interpolation" => "notebook_polynomial_interpolation.md"
        ],
        "Pretty Printing" => "custom_logging.md",
        "Internals" => [
            "DocStrings" => "dev_man.md",
            "Interfaces" => "Interfaces.md",
            "`AbstractResult` Interface" => "AbstractResultInterface.md",
            "`AbstractIteraData` Interface" => "AbstractIterDataInterface.md",
            "`AbstractDB` Interface" => "AbstractDBInterface.md",
        ]
    ],
)

deploydocs(;
    repo="github.com/manuelbb-upb/Morbit.jl",
    #devbranch="stronger_typing" # TODO Make this "master" after merging
)

Pkg.activate(current_env)


text:
module D
export d
d = 4
end


text:
@static if VERSION < v"1.2"
    function hasfield(T::Type, name::Symbol)
        Base.@_pure_meta
        Base.fieldindex(T, name, false) > 0
    end
end

max_stars_repo_name: anowacki/QuakeML


text:
function heatbathtest_4D(NX,NY,NZ,NT,β,NC)
    Dim = 4
    Nwing = 0

    #=
    u1 = IdentityGauges(NC,Nwing,NX,NY,NZ,NT)
    U = Array{typeof(u1),1}(undef,Dim)
    U[1] = u1
    for μ=2:Dim
        U[μ] = IdentityGauges(NC,Nwing,NX,NY,NZ,NT)
    end
    =#

    U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "cold")

    h = Heatbath(U,β)

    gauge_action = GaugeAction(U)
    plaqloop = make_loops_fromname("plaquette",Dim=Dim)
    append!(plaqloop,plaqloop')
    βinp = β/2
    push!(gauge_action,βinp,plaqloop)
    rectloop = make_loops_fromname("rectangular",Dim=Dim)
    append!(rectloop,rectloop')
    βinp = β/2
    push!(gauge_action,βinp,rectloop)
    hnew = Heatbath_update(U,gauge_action)

    temp1 = similar(U[1])
    temp2 = similar(U[1])
    temp3 = similar(U[1])

    comb = 6
    factor = 1/(comb*U[1].NV*U[1].NC)
    @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
    println("plaq_t = $plaq_t")
    poly = calculate_Polyakov_loop(U,temp1,temp2)
    println("polyakov loop = $(real(poly)) $(imag(poly))")
    #hnew = Heatbath_update(U,gauge_action)

    numhb = 200
    numOR = 3
    plaq_ave = 0.0
    for itrj = 1:numhb
        heatbath!(U,hnew)
        for ior=1:numOR
            overrelaxation!(U,hnew)
        end
        #heatbath!(U,h)
        #=
        if NC == 2
            heatbath_SU2!(U,NC,[temp1,temp2,temp3],β)
        elseif NC == 3
            heatbath_SU3!(U,NC,[temp1,temp2,temp3],β)
        else
            heatbath_SUN!(U,NC,[temp1,temp2,temp3],β)
        end
        =#
        plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
        plaq_ave += plaq_t

        if itrj % 40 == 0
            @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
            println("$itrj plaq_t = $plaq_t average: $(plaq_ave/itrj)")
            #println("$itrj plaq_t = $plaq_t")
            poly = calculate_Polyakov_loop(U,temp1,temp2)
            println("$itrj polyakov loop = $(real(poly)) $(imag(poly))")
        end
    end

    return  plaq_ave/numhb #plaq_t

end

function heatbathtest_2D(NX,NT,β,NC)
    Dim = 2
    Nwing = 1

    #=
    u1 = RandomGauges(NC,Nwing,NX,NT)
    U = Array{typeof(u1),1}(undef,Dim)
    U[1] = u1
    for μ=2:Dim
        U[μ] = RandomGauges(NC,Nwing,NX,NT)
    end
    =#

    U = Initialize_Gaugefields(NC,Nwing,NX,NT,condition = "hot",randomnumber="Reproducible")

    temp1 = similar(U[1])
    temp2 = similar(U[1])
    temp3 = similar(U[1])

    #comb = 6
    if Dim == 4
        comb = 6 #4*3/2
    elseif Dim == 3
        comb = 3
    elseif Dim == 2
        comb = 1
    else
        error("dimension $Dim is not supported")
    end

    factor = 1/(comb*U[1].NV*U[1].NC)
    @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
    println("plaq_t = $plaq_t")
    poly = calculate_Polyakov_loop(U,temp1,temp2)
    println("polyakov loop = $(real(poly)) $(imag(poly))")

    gauge_action = GaugeAction(U)
    plaqloop = make_loops_fromname("plaquette",Dim=Dim)
    append!(plaqloop,plaqloop')
    βinp = β/2
    push!(gauge_action,βinp,plaqloop)
    rectloop = make_loops_fromname("rectangular",Dim=Dim)
    append!(rectloop,rectloop')
    βinp = β/2
    push!(gauge_action,βinp,rectloop)
    hnew = Heatbath_update(U,gauge_action)

    numhb = 1000
    numOR = 3
    plaq_ave = 0.0
    for itrj = 1:numhb
        heatbath!(U,hnew)
        for ior=1:numOR
            overrelaxation!(U,hnew)
        end
        #heatbath!(U,[temp1,temp2,temp3],β)
        #=
        if NC == 2
            heatbath_SU2!(U,NC,[temp1,temp2,temp3],β,Dim)
        elseif NC == 3
            heatbath_SU3!(U,NC,[temp1,temp2,temp3],β,Dim)
        else
            heatbath_SUN!(U,NC,[temp1,temp2,temp3],β,Dim)
        end
        =#
        plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
        plaq_ave += plaq_t

        if itrj % 200 == 0
            #@time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
            println("$itrj plaq_t = $plaq_t average: $(plaq_ave/itrj)")
            poly = calculate_Polyakov_loop(U,temp1,temp2)
            println("$itrj polyakov loop = $(real(poly)) $(imag(poly))")
        end
    end

    return plaq_ave/numhb

end

println("2D system")
@testset "2D" begin
    NX = 4
    #NY = 4
    #NZ = 4
    NT = 4
    Nwing = 1

    @testset "NC=2" begin
        β = 2.3
        NC = 2
        println("NC = $NC")
        #val =  0.8275345965584693
        #val =0.6414596466929057
        val = 0.8237678929985711
        @time plaq_ave = heatbathtest_2D(NX,NT,β,NC)
        @test abs(plaq_ave-val)/abs(val) < eps
    end

    @testset "NC=3" begin
        β = 5.7
        NC = 3
        println("NC = $NC")
        val = 0.8070960744356934
        @time plaq_ave = heatbathtest_2D(NX,NT,β,NC)
        @test abs(plaq_ave-val)/abs(val) < eps
    end

    @testset "NC=4" begin
        β = 5.7
        NC = 4
        println("NC = $NC")
        val = 0.38564608114572807
        @time plaq_ave = heatbathtest_2D(NX,NT,β,NC)
        @test abs(plaq_ave-val)/abs(val) < eps
    end

end

println("4D system")
@testset "4D" begin
    NX = 4
    NY = 4
    NZ = 4
    NT = 4
    Nwing = 1

    @testset "NC=2" begin
        β = 2.3
        NC = 2
        println("NC = $NC")
        #val =0.9042749888688881
        val = 0.9040744318206326
        @time plaq_ave = heatbathtest_4D(NX,NY,NZ,NT,β,NC)
        @test abs(plaq_ave-val)/abs(val) < eps
    end

    @testset "NC=3" begin

        NC = 3
        β = 2.3*NC #5.7
        println("NC = $NC")
        #val = 0.9148612416401057
        val = 0.9148494229319207
        #@time plaq_t = heatbathtest_4D(NX,NY,NZ,NT,β,NC)
        @time plaq_ave = heatbathtest_4D(NX,NY,NZ,NT,β,NC)
        @test abs(plaq_ave-val)/abs(val) < eps
        #@test abs(plaq_t-val)/abs(val) < eps
    end

    @testset "NC=4" begin
        β = 5.7
        NC = 4
        println("NC = $NC")
        #val  =0.791346082492281
        val = 0.792025880016283
        #@time plaq_t = heatbathtest_4D(NX,NY,NZ,NT,β,NC)
        @time plaq_ave = heatbathtest_4D(NX,NY,NZ,NT,β,NC)
        @test abs(plaq_ave-val)/abs(val) < eps
        #@test abs(plaq_t-val)/abs(val) < eps
    end

end


text:
module DocumenterCitations

using Documenter
using Documenter.Anchors
using Documenter.Builder
using Documenter.Documents
using Documenter.Selectors
using Documenter.Utilities
using Documenter.Expanders

using Markdown
using Bibliography
using Bibliography: xnames, xyear, xlink, xtitle, xin

export CitationBibliography
struct CitationBibliography <: Documenter.Plugin
    bib::Dict
end
function CitationBibliography(filename::AbstractString="")
    filename == "" && return CitationBibliography(Dict())
    bf = import_bibtex(filename)
    return CitationBibliography(bf)
end

include("citations.jl")
include("bibliography.jl")

end


text:
@testset "1637.widest-vertical-area-between-two-points-containing-no-points.jl" begin
    @test max_width_of_vertical_area([[8,7],[9,9],[7,4],[9,7]]) == 1
    @test max_width_of_vertical_area([[3,1],[9,0],[1,0],[1,4],[5,3],[8,8]]) == 3
end


text:
using SIAMFANLEquations
using SIAMFANLEquations.TestProblems
using SIAMFANLEquations.Examples
using Test
using LinearAlgebra

include("Chapter1/nsolsc_solution_test.jl")
include("Chapter1/ptcsolsc_test.jl")
include("Chapter2/basic2d_test.jl")
include("Chapter2/heq_lu_test.jl")
include("Chapter2/bvp_test.jl")
include("Chapter2/beam_test.jl")
include("Chapter2/pde_lin_test.jl")
include("Chapter2/nsolpde_test.jl")
include("Chapter2/knowsdt_test.jl")
include("Chapter3/gmres_test.jl")
include("Chapter3/mgs_test.jl")
include("Chapter3/bicgstab_test.jl")
include("Chapter3/Krylov_pde_test.jl")
include("Chapter3/ptcKrylovTest.jl")
include("Chapter3/ptcKrylovTestB.jl")
include("Chapter3/nk_test.jl")
include("Chapter3/nk_pde.jl")
include("Chapter3/nk_heq.jl")
include("Chapter4/reldiff.jl")
include("Chapter4/alex_test.jl")
include("Chapter4/ci_pde_aa.jl")
include("Chapter4/heq_aa.jl")
include("Chapter4/linear_aa.jl")

@testset "Scalar Equations: Chapter 1" begin
    @test nsolsc_solution_test()
    @test ptcsolsc_test()
end
@testset "nsol and ptcsol: Chapter 2" begin
    @test basic2d_test()
    @test bvp_test(201)
    @test beam_test()
    @test heq_lu_test()
    @test pde_lin_test(31)
    @test nsolpde_test(31)
    @test knowsdt_test()
end
@testset "Newton-Krylov solvers: Chapter 3" begin
    @test nk_test()
    @test nk_pde()
    @test nk_heq()
    @test ptcKrylovTest()
    @test ptcKrylovTestB()
end
@testset "Krylov solvers: Chapter 3" begin
    @test gmres_test()
    @test mgs_test()
    @test bicgstab_test()
    @test gmres_test_pde(31)
    @test bicgstab_test_pde(31)
end
@testset "Anderson Acceleration: Chapter 4" begin
    @test ci_pde_aa()
    @test heq_aa()
    @test linear_aa()
    @test alex_test()
end


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export linearfst
"""
`linearfst(ilabels,olabels,w,[isym,[osym=isym]]; W=TropicalWeight{Float32})`

Creates a linear WFST. `ilabels`, `olabels` and `w` are the input labels, output labels and weights respectively which must be vectors of the same size.
If input and output tables `isym` and `osym` are not provided, the symbol tables are derived from `ilabels` and `olabels`.

```julia
julia> A = linearfst([1,1,3,4],[:a,:b,:c,:d],ones(TropicalWeight{Float64},4))
WFST #states: 5, #arcs: 4, #isym: 3, #osym: 4
|1/0.0|
1:a/0.0 → (2)
(2)
1:b/0.0 → (3)
(3)
3:c/0.0 → (4)
(4)
4:d/0.0 → (5)
((5/0.0))

```
"""
function linearfst(ilabels::AbstractVector,
                   olabels::AbstractVector,
                   w::AbstractVector,
                   isym::AbstractDict,
                   osym::AbstractDict=isym; W = TropicalWeight{Float32}
                  )
  N = length(w)
  if !(length(ilabels) == length(olabels) == N)
    throw(ErrorException("`ilabels`, `olabels`, and `w` must have all the same lengths"))
  end
  fst = WFST(isym,osym; W = W)
  for i = 1:N
    add_arc!(fst, i, i+1, ilabels[i], olabels[i], W(w[i]))
  end
  initial!(fst,1)
  final!(fst,length(w)+1)
  return fst
end

linearfst(ilabels,olabels,w) =
linearfst(ilabels,olabels,w,
          Dict(s=>k for (k,s) in enumerate(unique(ilabels) )),
          Dict(s=>k for (k,s) in enumerate(unique(olabels) ))
              )

export matrix2wfst
"""
`matrix2wfst(isym,[osym,] X; W = TropicalWeight{Float32})`

Creates a WFST with weight type `W` and input output tables `isym` and `osym` using the matrix `X`.
If `X` is a matrix of dimensions `Ns`x`Nt`, the resulting WFST will have `Nt+1` states.
Arcs form the `t`-th state of `fst` will go only to the `t+1`-th state and will correspond to the non-zero element of the `t`-th column of `X`.
State `1` and `Nt+1` are labelled as initial and final state, respectively.

```julia
julia> matrix2wfst(["a","b","c"],[1 2 3; 1 2 3; 1 2 3])
WFST #states: 4, #arcs: 9, #isym: 3, #osym: 3
|1/0.0f0|
a:a/1.0f0 → (2)
b:b/1.0f0 → (2)
c:c/1.0f0 → (2)
(2)
a:a/2.0f0 → (3)
b:b/2.0f0 → (3)
c:c/2.0f0 → (3)
(3)
a:a/3.0f0 → (4)
b:b/3.0f0 → (4)
c:c/3.0f0 → (4)
((4/0.0f0))

```
"""
function matrix2wfst(isym::AbstractDict{I,D},
                     osym::AbstractDict{O,D},
                     X::AbstractMatrix; W = TropicalWeight{Float32}) where {D,I,O}
  fst = WFST(isym, osym; W = W)
  Ns, Nt = size(X)
  add_states!(fst,Nt+1)
  for t = 1:Nt, j = 1:Ns
    if W(X[j,t]) != zero(W)
      arc = Arc{W,D}(j,j,W(X[j,t]),t+1)
      push_arc!(fst,t,arc)
    end
  end
  initial!(fst,1)
  final!(fst,Nt+1)
  return fst
end

matrix2wfst(isym,X;kwargs...) = matrix2wfst(isym,isym,X;kwargs...)
matrix2wfst(isym::Vector,osym::Vector,X;kwargs...) =
matrix2wfst(Dict(s=>k for (k,s) in enumerate(unique(isym) )),
            Dict(s=>k for (k,s) in enumerate(unique(osym) )),X; kwargs...)


text:
using MathToolkit
using Base.Test

@test recvec([1,2,3,4,5,6,7,8,9,10]) == [1,-2,1]
@test recvec([0,1,1,2,3,5,8,13], output=Float64) == [1.0,-1.0,-1.0]
@test recvec([2^k for k=0:7], compute=Int128) == [1,-2]
@test recvec([1//2^k for k=0:7], compute=Int128, output=Float64) == [2.0,-1.0]
@test recvec([1]) == []
@test recvec([1,2,3,4,5,6,7,8,9,10,42]) == []
@test recvec([BigInt(2)^k + k%16 for k=1:64]) ==
        [1,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,2]

@test abs(transpose([0.75,BigFloat(1)])*pslq([0.75, 1],1e-10))[1] < 1e-10
@test abs(transpose([0.5,BigFloat(1)])*pslq([0.5, 1],1e-10))[1] < 1e-10


text:
"""
R-Tree variants.
"""
@enum RTreeVariant::Int RTreeLinear = 0 RTreeQuadratic = 1 RTreeStar = 2

"""
Pool of the `T1` R-tree nodes (`Leaf` or `Branch`) with `T2` children
(`Branch`, `Leaf` or `Elem`).
It allows reusing the deleted nodes and reduce the stress on GC.
"""
mutable struct NodePool{T1<:Node, T2} <: AbstractPool{T1}
    objs::Vector{T1}

    capacity::Int       # pool capacity
    node_capacity::Int  # node capacity (max number of children)

    NodePool{T1,T2}(pool_capacity::Integer, node_capacity::Integer) where {T1<:Node, T2} =
        new{T1,T2}(sizehint!(Vector{T1}(), pool_capacity), pool_capacity, node_capacity)
end

const BranchPool{T,N,V} = NodePool{Branch{T,N,V}, Branch{T,N,V}}
const TwigPool{T,N,V} = NodePool{Branch{T,N,V}, Leaf{T,N,V}}
const LeafPool{T,N,V} = NodePool{Leaf{T,N,V}, V}

node_capacity(pool::NodePool) = pool.node_capacity
newelem(pool::NodePool{T1,T2}) where {T1, T2} = T1(nothing, T2)

"""
R-Tree: `N`-dimensional spatial data index [guttman84].

R-tree groups data elements (`V`) into leaves (`Leaf`) and leaves into branches (`Branch`).
It uses various heuristics to ensure that the minimal bounding rectangles (MBRs)
of the nodes (`Rect{T,N}` rectangles that encompass the data elements attached to these nodes)
stay compact and that the MBRs of the nodes that are on the same level of
R-tree hierarchy have minimal overlap with each other.
This property makes R-trees efficient for spatial queries.

To facilitate spatial indexing, the `V` data elements need to support `HasMBR`
trait (i.e. define `mbrtype(V)` and `mbr(v::V)` methods) and, optionally,
`HasID` trait (via `idtype(V)` and `id(v::V)` methods). `mbr(v::V)` should
return minimal bounding rectangle (MBR) of type `Rect{T,N}` that contains `v`.
`SpatialElem{T,N,D}` type provides generic implementation of spatial data element
that explicitly stores `id`, `mbr` and data object of type `D` and implements
`HasMBR` and `HasID` traits.

# Parameters

The behaviour of `RTree` is defined by the parameters supplied at its creation:
  * `T`: the numeric type for the spatial coordinate
  * `N`: the number of spatial dimensions
  * `variant`: one of `RTreeLinear`, `RTreeQuadratic`, or `RTreeStar` (default)
  * `tight_mbrs`: recalculate node MBR when the child is removed (default is `true`)
  * `branch_capacity`: capacity of branch nodes (default is `100`)
  * `leaf_capacity`: capacity of leaf nodes (default is `100`)
  * `leafpool_capacity`: How many detached 1st level nodes (leaves) should be kept for reuse (default is `100`)
  * `twigpool_capacity`: How many detached 2nd level nodes should be kept for reuse (default is `100`)
  * `branchpool_capacity`: How many other (level > 2) detached branch nodes should be kept for reuse (default is `100`)
  * `nearmin_overlap`: How many candidates to consider when identifying the node with minimal overlap (default is `32`)
  * `fill_factor`: How much should the node be filled (fraction of its capacity) after splitting (default is `0.7`)
  * `split_factor`: How much can the sizes of the two nodes differ after splitting (default is `0.4`)
  * `reinsert_factor`: How much should the node be underfilled (fraction of its capacity)
    to consider removing it and redistributing its children to other nodes (default is `0.3`)

# Performance

The nodes in R-tree have limited capacity (maximual number of children)
specified at `RTree` creation (`leaf_capacity` and `branch_capacity`).
Larger capacities results in shorter trees, but they time required to locate
the specific spatial region grows linearly with the capacity.

# References
[guttman84] “R-Trees: A Dynamic Index Structure for Spatial Searching”
    A. Guttman, Proc. 1984 ACM-SIGMOD Conference on Management of
    Data (1985), 47-57.
[beckmann90] "The R*-tree: an efficient and robust access method for points and rectangles"
    N. Beckmann, H.P. Kriegel, R. Schneider, B. Seeger, Proc. 1990 ACM SIGMOD
    international conference on Management of data (1990), p.322
"""
mutable struct RTree{T,N,V} <: SpatialIndex{T,N,V}
    # R-tree insert/update policy
    variant::RTreeVariant
    tight_mbrs::Bool        # whether MBR are always updated upon child change/removal

    nearmin_overlap::Int
    fill_factor::Float64
    split_factor::Float64
    reinsert_factor::Float64

    root::Node{T,N}     # root node

    nelems::Int         # number of data elements in a tree (`length()`)
    nnodes_perlevel::Vector{Int}    # number of nodes (branches and leaves) at each tree level
    nelem_insertions::Int           # number of elements insertions
    nelem_deletions::Int            # number of elements deletions
    nnode_splits::Int               # number of node splits
    nnode_reinsertions::Int         # number of node reinsertions

    branchpool::BranchPool{T,N,V}
    twigpool::TwigPool{T,N,V}
    leafpool::LeafPool{T,N,V}

    function RTree{T,N,V}(;
        variant::RTreeVariant = RTreeStar,
        tight_mbrs::Bool = true,
        branch_capacity::Integer = 100,
        leaf_capacity::Integer = 100,
        leafpool_capacity::Integer = 100,
        branchpool_capacity::Integer = 100,
        nearmin_overlap::Integer = floor(Int, 0.32 * leaf_capacity),
        fill_factor::Real = 0.7,
        split_factor::Real = 0.4,
        reinsert_factor::Real = 0.3
    ) where {T<:Number,N,V}
        check_eltype_rtree(V, Leaf{T,N,V})
        leafpool = LeafPool{T,N,V}(leafpool_capacity, leaf_capacity)
        new{T,N,V}(variant, tight_mbrs,
                   nearmin_overlap, fill_factor,
                   split_factor, reinsert_factor,
                   acquire!(leafpool),
                   0, [1], # init with 1 leaf == parent
                   0, 0, 0, 0,
                   BranchPool{T,N,V}(branchpool_capacity, branch_capacity),
                   TwigPool{T,N,V}(branchpool_capacity, branch_capacity),
                   leafpool)
    end
end

"""
    RTree{T,N}([K], V; kwargs...) where {T,N,V}

Construct `N`-dimensional [`RTree`](@ref) object that stores
[`SpatialElem{T,N,K,V}`](@ref `SpatialElem`) elements.
If `K` (the type of spatial element identifier) if not specified,
the stored elements will not have ids.
"""
RTree{T,N}(::Type{V}; kwargs...) where {T,N,V} =
    RTree{T,N,SpatialElem{T,N,Nothing,V}}(; kwargs...)

RTree{T,N}(::Type{K}, ::Type{V}; kwargs...) where {T,N,K,V} =
    RTree{T,N,SpatialElem{T,N,K,V}}(; kwargs...)

# create an empty R-tree of the same type and with the same settings as `tree`
Base.similar(tree::RTree) =
    typeof(tree)(variant=tree.variant,
                 tight_mbrs=tree.tight_mbrs,
                 branch_capacity=capacity(Branch, tree),
                 leaf_capacity=capacity(Leaf, tree),
                 leafpool_capacity=capacity(tree.leafpool),
                 branchpool_capacity=capacity(tree.branchpool),
                 nearmin_overlap=tree.nearmin_overlap,
                 fill_factor=tree.fill_factor,
                 split_factor=tree.split_factor,
                 reinsert_factor=tree.reinsert_factor)

mbrtype(::Type{<:RTree{T,N}}) where {T,N} = Rect{T,N}
mbrtype(tree::RTree) = mbrtype(typeof(tree))
regiontype(R::Type{<:RTree}) = mbrtype(R)
leaftype(::Type{RTree{T,N,V}}) where {T,N,V} = Leaf{T,N,V}
leaftype(tree::RTree) = leaftype(typeof(tree))
branchtype(::Type{RTree{T,N,V}}) where {T,N,V} = Branch{T,N,V}
branchtype(tree::RTree) = branchtype(typeof(tree))
nodetype(::Type{RTree{T,N,V}}) where {T,N,V} = Node{T,N,V} # not a concrete type
nodetype(tree::RTree) = nodetype(typeof(tree))

check_eltype_rtree(el::Any, tree::RTree) =
    check_eltype_rtree(typeof(el), leaftype(tree))

height(tree::RTree) = level(tree.root)
variant(tree::RTree) = tree.variant
mbr(tree::RTree) = mbr(tree.root)

capacity(::T, tree::RTree) where T<:Node = capacity(T, tree)
capacity(::Type{<:Leaf}, tree::RTree) = node_capacity(tree.leafpool)
capacity(::Type{<:Branch}, tree::RTree) = node_capacity(tree.branchpool)

# acquire Leaf/Branch from the corresponding node pool
function acquire(tree::RTree, ::Type{<:Leaf}, lev::Int = 1, br::Rect = empty(mbrtype(tree)))
    @assert lev == 1
    leaf = acquire!(tree.leafpool)
    @assert isempty(leaf)
    leaf.mbr = br
    return leaf
end

function acquire(tree::RTree, ::Type{<:Branch}, lev::Int, br::Rect = empty(mbrtype(tree)))
    @assert lev > 1
    node = acquire!(lev == 2 ? tree.twigpool : tree.branchpool)
    @assert isempty(node)
    node.mbr = br
    node.level = lev
    return node
end

# detach the Leaf/Branch from the tree and release it to the corresponding node pool
# it's up to the calling context to take care of the R-tree integrity (parent and children nodes)
function release(tree::RTree, node::Node)
    node.parent = nothing # don't refer to the parent (so it could be GCed)
    empty!(node.children) # don't refer to the children (so they could be GCed)
    if node isa Leaf
        release!(tree.leafpool, node)
    else
        release!(eltype(node.children) <: Leaf ? tree.twigpool : tree.branchpool, node)
    end
end

# low-level insertion of a child into the back of the list of the node children
# `child` parent is set to `node`, MBR of `node` is extended with child MBR
# no subtree balancing or anything like that
function _attach!(node::Node, child::Any, tree::RTree; force::Bool = false)
    force || (length(node) < capacity(node, tree)) ||
        throw(SpatialIndexException("Node is full ($(capacity(node, tree)) items)"))
    if child isa Node
        (level(node) == level(child) + 1) ||
            throw(SpatialIndexException("Parent ($(level(node)) and child ($(level(child))) levels don't match"))
        child.parent = node
    end
    push!(node.children, child)
    node.mbr = length(node) > 1 ? combine(node.mbr, mbr(child)) : mbr(child)
    return node
end

# low-level removal of a child node from the specified position, no subtree balancing
# return true if node MBR was updated
function _detach!(node::Node, pos::Integer, tree::RTree; updatembr = tree.tight_mbrs)
    if length(node) > 1
        delmbr = mbr(node[pos])   # cache it, since its required for "touches" later
        # copy the last child into i-th position and pop it
        # this should be more efficient that deleting the i-th node
        if pos < length(node)
            node.children[pos] = pop!(node.children)
        else
            pop!(node.children)
        end
        if updatembr && touches(mbr(node), delmbr)
            syncmbr!(node)
            return true
        end
    else # node becomes empty
        @assert pos == 1
        pop!(node.children)
        node.mbr = empty(mbrtype(node)) # reset br
        return true
    end
    return false
end

# The context for R-tree node insertion/balancing.
struct SubtreeContext{T,N,V,O}
    tree::RTree{T,N,V}
    level_overflow::O       # nothing or an indicator of the R-tree level overflow

    function SubtreeContext(tree::RTree{T,N,V}) where {T,N,V}
        O = variant(tree) == RTreeStar ? BitVector : Nothing
        new{T,N,V,O}(tree, O == BitVector ? falses(height(tree)) : nothing)
    end
end

isoverflow(con::SubtreeContext, lev::Int) = con.level_overflow[lev]
setoverflow!(con::SubtreeContext, lev::Int) = con.level_overflow[lev] = true


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# No aggregation
# ------------------------------------------------------------
"""
    NoAggregation

Empty functor object for running an L-shaped algorithm without aggregation (multi-cut L-shaped).

"""
struct NoAggregation <: AbstractAggregation end

function aggregate_cut!(lshaped::AbstractLShaped, ::NoAggregation, cut::HyperPlane)
    return add_cut!(lshaped, cut)
end

function aggregate_cut!(cutqueue::CutQueue, ::NoAggregation, ::MetaDataChannel, t::Integer, cut::HyperPlane, x::AbstractArray)
    put!(cutqueue, (t, cut))
    return nothing
end

function num_thetas(num_subproblems::Integer, ::NoAggregation)
    return num_subproblems
end

function num_thetas(num_subproblems::Integer, ::NoAggregation, ::AbstractScenarioProblems)
    return num_subproblems
end

function flush!(::AbstractLShaped, ::NoAggregation)
    return false
end

function flush!(::CutQueue, ::NoAggregation, ::MetaDataChannel, ::Integer, ::AbstractArray)
    return false
end

# API
# ------------------------------------------------------------
"""
    DontAggregate

Factory object for [`NoAggregation`](@ref). Passed by default to `aggregate` in `LShaped.Optimizer`.

"""
struct DontAggregate <: AbstractAggregator end

function (::DontAggregate)(::Integer, ::Type{<:AbstractFloat})
    return NoAggregation()
end

function remote_aggregator(::NoAggregation, ::AbstractScenarioProblems, ::Integer)
    return DontAggregate()
end

function str(::DontAggregate)
    return "disaggregate cuts"
end


text:
struct _enum_TTransportExceptionTypes
    UNKNOWN::Int32
    NOT_OPEN::Int32
    ALREADY_OPEN::Int32
    TIMED_OUT::Int32
    END_OF_FILE::Int32
end

const TransportExceptionTypes = _enum_TTransportExceptionTypes(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4))

struct TTransportException <: Exception
    typ::Int32
    message::AbstractString

    TTransportException(typ=TransportExceptionTypes.UNKNOWN, message::AbstractString="") = new(typ, message)
end

# TODO: Thrift SASL server transport
# Thrift SASL client transport
mutable struct TSASLClientTransport <: TTransport
    tp::TTransport
    mech::String
    callback::Function
    rbuff::IOBuffer
    wbuff::IOBuffer

    function TSASLClientTransport(tp::TTransport, mech::String=SASL_MECH_PLAIN, callback::Function=sasl_callback_default)
        validate_sasl_mech(mech)
        new(TFramedTransport(tp), mech, callback, PipeBuffer(), PipeBuffer())
    end
end

rawio(t::TSASLClientTransport)  = rawio(t.tp)
close(t::TSASLClientTransport)  = close(t.tp)
isopen(t::TSASLClientTransport) = isopen(t.tp)
flush(t::TSASLClientTransport)  = flush(t.tp)

read!(t::TSASLClientTransport, buff::Vector{UInt8}) = read!(t.tp, buff)
read(t::TSASLClientTransport, UInt8) = read(t.tp, UInt8)
function write(t::TSASLClientTransport, buff::Vector{UInt8})
    @debug("TSASLClientTransport buffering bytes", len=length(buff))
    write(t.tp, buff)
end
function write(t::TSASLClientTransport, b::UInt8)
    @debug("TSASLClientTransport buffering 1 byte")
    write(t.tp, b)
end

function open(t::TSASLClientTransport)
    open(t.tp)
    sasl_negotiate(rawio(t), t.mech, t.callback)
end

# Thrift Framed Transport
mutable struct TFramedTransport <: TTransport
    tp::TTransport
    rbuff::IOBuffer
    wbuff::IOBuffer
    TFramedTransport(tp::TTransport) = new(tp, PipeBuffer(), PipeBuffer())
end
rawio(t::TFramedTransport)  = rawio(t.tp)
open(t::TFramedTransport)   = open(t.tp)
close(t::TFramedTransport)  = close(t.tp)
isopen(t::TFramedTransport) = isopen(t.tp)

readframesz(t::TFramedTransport) = _read_fixed(t.tp, UInt32(0), 4, true)
function readframe(t::TFramedTransport)
    @debug("TFramedTransport reading frame")
    sz = readframesz(t)
    @debug("TFramedTransport reading frame", sz)
    write(t.rbuff, read!(t.tp, Vector{UInt8}(undef,sz)))
    @debug("TFramedTransport read frame", sz)
    nothing
end

function read!(t::TFramedTransport, buff::Vector{UInt8})
    ntotal = length(buff)
    nread = 0

    while nread < ntotal
        navlb = bytesavailable(t.rbuff)
        nremain = ntotal - nread
        if navlb < nremain
            @debug("reading new frame", navlb, nremain)
            readframe(t)
            navlb = bytesavailable(t.rbuff)
        end
        nbuff = min(navlb, nremain)
        Base.read_sub(t.rbuff, buff, nread+1, nbuff)
        nread += nbuff
    end
    buff
end
function read(t::TFramedTransport, UInt8)
    navlb = bytesavailable(t.rbuff)
    if navlb == 0
        readframe(t)
    end
    return read(t.rbuff, UInt8)
end

function write(t::TFramedTransport, buff::Vector{UInt8})
    @debug("TFramedTransport buffering bytes", len=length(buff))
    write(t.wbuff, buff)
end
function write(t::TFramedTransport, b::UInt8)
    @debug("TFramedTransport buffering 1 byte")
    write(t.wbuff, b)
end
function flush(t::TFramedTransport)
    szbuff = IOBuffer()
    navlb = bytesavailable(t.wbuff)
    @debug("sending data", navlb)
    _write_fixed(szbuff, UInt32(navlb), true)
    nbyt = write(t.tp, take!(szbuff))
    nbyt += write(t.tp, take!(t.wbuff))
    @debug("wrote frame", nbyt)
    flush(t.tp)
end

# Thrift Socket Transport
mutable struct TSocket <: TTransport
    host::AbstractString
    port::Integer

    io::TCPSocket

    TSocket(host::AbstractString, port::Integer) = new(host, port)
    TSocket(port::Integer) = TSocket("127.0.0.1", port)
end

mutable struct TServerSocket <: TServerTransport
    host::AbstractString
    port::Integer

    io::TCPServer

    TServerSocket(host::AbstractString, port::Integer) = new(host, port)
    TServerSocket(port::Integer) = TServerSocket("", port)
end

const TSocketBase = Union{TSocket, TServerSocket}

open(tsock::TServerSocket) = nothing
open(tsock::TSocket) = (!isopen(tsock) && (tsock.io = connect(tsock.host, tsock.port)); nothing)

listen(tsock::TServerSocket) = (tsock.io = isempty(tsock.host) ? listen(tsock.port) : listen(parseip(tsock.host), tsock.port); nothing)
function accept(tsock::TServerSocket)
    accsock = TSocket(tsock.host, tsock.port)
    accsock.io = accept(tsock.io)
    accsock
end

close(tsock::TSocketBase) = (isopen(tsock.io) && close(tsock.io); nothing)
rawio(tsock::TSocketBase) = tsock.io
read!(tsock::TSocketBase, buff::Vector{UInt8}) = read!(tsock.io, buff)
read(tsock::TSocketBase, UInt8) = read(tsock.io, UInt8)
write(tsock::TSocketBase, buff::Vector{UInt8}) = write(tsock.io, buff)
write(tsock::TSocketBase, b::UInt8) = write(tsock.io, b)
flush(tsock::TSocketBase)   = flush(tsock.io)
isopen(tsock::TSocketBase)  = (isdefined(tsock, :io) && isreadable(tsock.io) && iswritable(tsock.io))

# Thrift Memory Transport
mutable struct TMemoryTransport <: TTransport
    buff::IOBuffer

    TMemoryTransport() = new(PipeBuffer())
    TMemoryTransport(buff::Array{UInt8}) = new(PipeBuffer(buff))
end

rawio(t::TMemoryTransport)  = t.buff
open(t::TMemoryTransport)   = nothing
close(t::TMemoryTransport)  = nothing
isopen(t::TMemoryTransport) = true
flush(t::TMemoryTransport)  = nothing
read!(t::TMemoryTransport, buff::Vector{UInt8}) = read!(t.buff, buff)
read(t::TMemoryTransport, UInt8) = read(t.buff, UInt8)
write(t::TMemoryTransport, buff::Vector{UInt8}) = write(t.buff, buff)
write(t::TMemoryTransport, b::UInt8) = write(t.buff, b)

# Thrift File IO Transport
mutable struct TFileTransport <: TTransport
    handle::IO
end

rawio(t::TFileTransport)  = t.handle
open(t::TFileTransport)   = nothing
close(t::TFileTransport)  = nothing
isopen(t::TFileTransport) = true
flush(t::TFileTransport)  = flush(t.handle)
read!(t::TFileTransport, buff::Vector{UInt8}) = read!(t.handle, buff)
read(t::TFileTransport, UInt8) = read(t.handle, UInt8)
write(t::TFileTransport, buff::Vector{UInt8}) = write(t.handle, buff)
write(t::TFileTransport, b::UInt8) = write(t.handle, b)


text:
module mmmmPoiss_06122017
using FinEtools
using FinEtoolsHeatDiff
using Test
import LinearAlgebra: cholesky, Symmetric
function test()

  # println("""
  #
  # Heat conduction example described by Amuthan A. Ramabathiran
  # http://www.codeproject.com/Articles/579983/Finite-Element-programming-in-Julia:
  # Unit square, with known temperature distribution along the boundary,
  # and uniform heat generation rate inside.  Mesh of regular linear TRIANGLES,
  # in a grid of 1000 x 1000 edges (2M triangles, 1M degrees of freedom).
  # Version: 05/29/2017
  # """
  # )
  t0 = time()

  A = 1.0 # dimension of the domain (length of the side of the square)
  thermal_conductivity = [i==j ? one(FFlt) : zero(FFlt) for i=1:2, j=1:2]; # conductivity matrix
  Q = -6.0; # internal heat generation rate
  function getsource!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    forceout[1] = Q; #heat source
  end
  tempf(x) = (1.0 .+ x[:,1].^2 .+ 2*x[:,2].^2);#the exact distribution of temperature
  N = 10;# number of subdivisions along the sides of the square domain

  # println("Mesh generation")
  fens,fes =T3block(A, A, N, N)

  geom = NodalField(fens.xyz)
  Temp = NodalField(zeros(size(fens.xyz,1),1))

  # println("Searching nodes  for BC")
  l1 = selectnode(fens; box=[0. 0. 0. A], inflate = 1.0/N/100.0)
  l2 = selectnode(fens; box=[A A 0. A], inflate = 1.0/N/100.0)
  l3 = selectnode(fens; box=[0. A 0. 0.], inflate = 1.0/N/100.0)
  l4 = selectnode(fens; box=[0. A A A], inflate = 1.0/N/100.0)
  List = vcat(l1, l2, l3, l4)
  setebc!(Temp, List, true, 1, tempf(geom.values[List,:])[:])
  applyebc!(Temp)
  numberdofs!(Temp)

  t1 = time()

  material = MatHeatDiff(thermal_conductivity)

  femm = FEMMHeatDiff(IntegDomain(fes, TriRule(1), 100.), material)

  # println("Conductivity")
  K = conductivity(femm, geom, Temp)
  # println("Nonzero EBC")
  @show F2 = nzebcloadsconductivity(femm, geom, Temp);
  # println("Internal heat generation")
  # fi = ForceIntensity(FFlt, getsource!);# alternative  specification
  fi = ForceIntensity(FFlt[Q]);
    F1 = distribloads(femm, geom, Temp, fi, 3);

  # println("Factorization")
  K = cholesky(Symmetric(K))
  # println("Solution of the factorized system")
  U = K\(F1+F2)
  scattersysvec!(Temp,U[:])

  # println("Total time elapsed = $(time() - t0) [s]")
  # println("Solution time elapsed = $(time() - t1) [s]")

  Error= 0.0
  for k=1:size(fens.xyz,1)
    Error = Error.+abs.(Temp.values[k,1].-tempf(reshape(fens.xyz[k,:], (1,2))))
  end
  # println("Error =$Error")

  # File =  "a.vtk"
  # MeshExportModule.vtkexportmesh (File, fes.conn, [geom.values Temp.values], MeshExportModule.T3; scalars=Temp.values, scalars_name ="Temperature")

  @test Error[1]<1.e-5

  true
end
end
using .mmmmPoiss_06122017
mmmmPoiss_06122017.test()

max_stars_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff

max_issues_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff

max_forks_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff


text:
mutable struct OBOPDataset
    file_name::String
    total_voters::Int64
    total_itens::Int64
    item_names::Array{String,1}
    C::Array{Int64,2}
    D::Array{Int64,2}
    interval_init::Array{Float64,1}
end

mutable struct OBOPSolution
    id::Int64
    objective::Int64
    bucket::Array{Int64,1}
    start_time::Float64
    total_time::Float64
    generation::Int64
    local_search::Bool
    function OBOPSolution(n,start_time)
        new(0,typemin(Int64),zeros(Int64,n),start_time,0.0,0,false)
    end
end

mutable struct OBOPStatistics
    elite_mean::Array{Float64,1}
    elite_std::Array{Float64,1}
    pop_mean::Array{Float64,1}
    pop_std::Array{Float64,1}
    best_generations::Array{Int64,1}
    time_generations::Array{Float64,1}
    time_total::Float64
    time_best::Float64
    total_local_search::Int64
    total_local_search_effective::Int64
end

max_stars_repo_name: LuizHNLorena/AbrkgaObop

max_issues_repo_name: LuizHNLorena/AbrkgaObop

max_forks_repo_name: LuizHNLorena/AbrkgaObop


text:
using BEAST
using CompScienceMeshes
using LinearAlgebra

# define geometry & material parameters
f₀= 3.0e6 # 3 GHz
c = 3.0e8
ϵ₀ = 8.85418782e-12
μ₀ = 1.25663706e-6
ϵᵣ=1; μᵣ=1
l = w = 1.0 #Length and width of capacitor plates
d = 0.1     #seperation of plates
h = 1/12    #size of meshes
ω = 2π*f₀
vₑ = c/sqrt(ϵᵣ)
λₑ = vₑ/f₀
η = sqrt(μ₀*μᵣ/ϵ₀*ϵᵣ)
κ = ω/vₑ

println("wavenumber = ", κ)
println("freq = ", f₀, " Hz")
println("λ = ", λₑ, "m; l = ", l ,"m")

# Build and mesh the structure
Γ₁ = meshrectangle(l,w,h);
CompScienceMeshes.translate!(Γ₁, point(0.0,0.0,d))
Γ₀ = meshrectangle(l,w,h)

γ₁ = meshsegment(l, h, 3)
CompScienceMeshes.translate!(γ₁, point(0.0,0.0,d))
γ₀ = meshsegment(l, h, 3)

Γ = weld(Γ₁, Γ₀)

# define the excitation
V₀ = 1.0
f = ScalarTrace(p -> V₀)

#define basis function
RT = rt_ports(Γ,[γ₁,γ₀])

# define the equations
@hilbertspace j
@hilbertspace k
T = MWSingleLayer3D(im*κ)
trc = X->ntrace(X,γ₁)
EFIE = @varform η*T[k,j] == f[trc(k)]

# discretise & solve the equation
efie = @discretise EFIE j∈RT k∈RT
u = solve(efie)

#Calculate & plot face currents
fcr,geo = facecurrents(u, RT); println("Face currents calculated")
# patch_mat(geo, fcr); println("Facecurrents plotted")

# Compute the Scalar Potential across the ports
z = range(-0.5, stop=0.5, length=100)
pts1 = point.(0.5,0.5,z)
# pts1 = [point(0.5,0.5,a) for a in range(-0.5,stop=0.5,length=100)]
volt = η*potential(MWSingleLayerPotential3D(κ), pts1, u, RT, type=ComplexF64)

# Plot the Scalar potential
using Plots
Φ = real.(getindex.(volt,1))
display(plot(z,Φ,xlabel="height",ylabel="scalar potential",label=""))

# Compare numerical capacitance with parallel plate approximation
idx = getindex_rtg(RT) #get index of global RWG
Q = norm(-u[idx]/(im*ω))
C = (ϵ₀*ϵᵣ*(l*w))/d
println("Relative different w.r.t. large plate capacitance: ", 2*abs(Q-C)/abs(Q+C))


text:
# Utils for MathOptInterface
#
##################################################
# Import legacy from LinQuadOptInterface to ease the integration
# of KNITRO quadratic and linear facilities.
##################################################
# URL: https://github.com/JuliaOpt/LinQuadOptInterface.jl
#
# LICENSE:
# MIT License
# Copyright (c) 2017 Oscar Dowson, Joaquim Dias Garcia and contributors
##################################################

function reduce_duplicates!(
    rows::Vector{T},
    cols::Vector{T},
    vals::Vector{S},
) where {T} where {S}
    @assert length(rows) == length(cols) == length(vals)
    for i in 1:length(rows)
        if rows[i] > cols[i]
            tmp = rows[i]
            rows[i] = cols[i]
            cols[i] = tmp
        end
    end
    return findnz(sparse(rows, cols, vals))
end

"""
    canonical_quadratic_reduction(func::ScalarQuadraticFunction)

Reduce a ScalarQuadraticFunction into three arrays, returned in the following
order:
 1. a vector of quadratic row indices
 2. a vector of quadratic column indices
 3. a vector of quadratic coefficients

Warning: we assume in this function that all variables are correctly
ordered, that is no deletion or swap has occured.
"""
function canonical_quadratic_reduction(func::MOI.ScalarQuadraticFunction)
    quad_columns_1, quad_columns_2, quad_coefficients = (
        Cint[term.variable_1.value for term in func.quadratic_terms],
        Cint[term.variable_2.value for term in func.quadratic_terms],
        Cdouble[term.coefficient for term in func.quadratic_terms],
    )
    # Take care of difference between MOI standards and KNITRO ones.
    for i in 1:length(quad_coefficients)
        @inbounds if quad_columns_1[i] == quad_columns_2[i]
            quad_coefficients[i] *= 0.5
        end
    end
    return reduce_duplicates!(quad_columns_1, quad_columns_2, quad_coefficients)
end

"""
    canonical_linear_reduction(func::Quad)

Reduce a ScalarQuadraticFunction into two arrays, returned in the following
order:
 1. a vector of linear column indices
 2. a vector of linear coefficients

Warning: we assume in this function that all variables are correctly
ordered, that is no deletion or swap has occured.
"""
function canonical_linear_reduction(func::MOI.ScalarQuadraticFunction)
    affine_columns = Cint[term.variable.value - 1 for term in func.affine_terms]
    affine_coefficients = Cdouble[term.coefficient for term in func.affine_terms]
    return affine_columns, affine_coefficients
end
function canonical_linear_reduction(func::MOI.ScalarAffineFunction)
    affine_columns = Cint[term.variable.value - 1 for term in func.terms]
    affine_coefficients = Cdouble[term.coefficient for term in func.terms]
    return affine_columns, affine_coefficients
end

function canonical_vector_affine_reduction(func::MOI.VectorAffineFunction)
    index_cols = Cint[]
    index_vars = Cint[]
    coefs = Cdouble[]

    for t in func.terms
        push!(index_cols, t.output_index - 1)
        push!(index_vars, t.scalar_term.variable.value - 1)
        push!(coefs, t.scalar_term.coefficient)
    end
    return index_cols, index_vars, coefs
end

# Convert Julia'Inf to KNITRO's Inf.
function check_value(val::Float64)
    if val > KN_INFINITY
        return KN_INFINITY
    elseif val < -KN_INFINITY
        return -KN_INFINITY
    end
    return val
end


text:
"
    SimpleEnv does nothing fancy to ensure the correctness of the environments it returns.
    Supports setleftenv! and setrightenv!
    Only used internally (in idmrg); no public constructor is provided
"
struct SimpleEnv{H<:Operator,V} <:Cache
    opp :: H
    lw :: PeriodicArray{V,2}
    rw :: PeriodicArray{V,2}
end

function SimpleEnv(state,pars::AbstractInfEnv)
    lw = similar(pars.lw)
    rw = similar(pars.rw)
    for i = 1:length(state)
        lw[i,:] = leftenv(pars,i,state)
        rw[i,:] = rightenv(pars,i,state)
    end

    return SimpleEnv(pars.opp,lw,rw)
end

leftenv(pars::SimpleEnv,pos::Int,state) = pars.lw[pos,:];
function setleftenv!(pars::SimpleEnv,pos,mps,lw)
    for i in 1:length(lw)
        pars.lw[pos,i] = lw[i]
    end
end

rightenv(pars::SimpleEnv,pos::Int,state) = pars.rw[pos,:];
function setrightenv!(pars::SimpleEnv,pos,mps,rw)
    for i in 1:length(rw)
        pars.rw[pos,i] = rw[i]
    end
end


text:
module ConfigCheck

    using Formatting

    export ConfigEntry, validateConfigEntries, validateByConfigGroup

    mutable struct ConfigEntry
        name      :: String
        required  :: Symbol
        valid_vts :: AbstractArray
        default   :: Any
        desc      :: String
        function ConfigEntry(
            name      :: String,
            required  :: Symbol,
            valid_vts :: AbstractArray,
            default   :: Any = nothing;
            desc      :: String = "",
        )

            return new(
                name,
                required,
                valid_vts,
                default,
                desc,
            )
        end

    end
#=
    function validateByConfigGroup(
        cfg                :: Dict,
        cfg_desc           :: Dict{Symbol, AbstractArray{ConfigEntry, 1}},
        target_group_names :: Union{AbstractArray{Symbol}, Nothing} = nothing;
        verbose :: Bool = false,
    )

        # cfg_desc is a Dict whose key is the name of a config group and
        # value is an Array of ConfigEntry

        all_group_names = keys(cfg_desc.groups)
        if target_group_names == nothing
            target_group_names = all_group_names
        end

        new_cfg = Dict{Symbol,Any}()
        for group_name in all_group_names
            if group_name in target_group_names
                new_cfg[group_name] = validateConfig(cfg[group_name], cfg_desc[group_name])
            else
                # simply transfer it
                new_cfg[group_name] = cfg[group_name]
            end
        end

        return new_cfg
    end

=#
    function validateConfigEntries(
        cfg         :: Dict,
        cfg_entries :: AbstractArray{ConfigEntry, 1};
        verbose     :: Bool = false,
    )

        new_cfg = Dict{String, Any}()

        for entry in cfg_entries

            name      = entry.name
            required  = entry.required
            valid_vts = entry.valid_vts
            default   = entry.default

            if ! ( required in [:optional, :required] )
                throw(ErrorException("The `required` only takes :optional or :required"))
            end

            if haskey(cfg, name)
                pass = false
                for valid_vt in valid_vts

                    if typeof(valid_vt) <: Union{DataType, UnionAll}
                        if typeof(cfg[name]) <: valid_vt
                            pass = true
                        end
                    else
                        if cfg[name] == valid_vt
                            pass = true
                        end
                    end

                    if pass
                        break
                    end

                end

                if pass

                    new_cfg[name] = cfg[name]
                    verbose && println(format("[Validation] Config `{:s}` : {:s}", string(name), string(new_cfg[name])))
                else
                    throw(ErrorException(format(
                        "[Error] Invalid value of key `{:s}`: {:s}. Valid values/types: `{:s}`.",
                        string(name),
                        string(cfg[name]),
                        join(string.(valid_vts), "` ,`")
                    )))
                end

            else

                msg = format(
                    "Missing config: `{:s}`. Valid values/types: `{:s}`.",
                    string(name),
                    join(string.(valid_vts), "` ,`")
                )

                if required == :required
                    throw(ErrorException(format("[Required] {:s}", msg)))
                else
                    new_cfg[name] = default
                    verbose && println(format("[Optional] {:s} is set to default: {:s}", string(name), string(new_cfg[name])))

                end
            end
        end

        #dropped_names = filter(x -> !( haskey(new_cfg, x) ), keys(cfg))
        #for dropped_name in dropped_names
        #    println(format("The config `{:s}` is not used.", string(dropped_name)))
        #end

        return new_cfg

    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using Patents
using PatentsLens
using Graphs
using Test

apps = PatentsLens.read("data/portfolio_carbios.jsonl")
fams = aggregate_families(apps)
g = citationgraph(fams)

@testset "application" begin
    @test apps isa Vector{Application}
    @test fams isa Vector{Family}

end

@testset "mainpath" begin
    # start = rand(eachindex(fams), 10)
    # mp = mainpath(g, SPCEdge(:log), ForwardBackwardLocal(start))
    # mc = Patents.component(mp, 1)
    # @test all(f in fams[start] for f in fams[mc.vertices[mc.start]])
end


text:
@doc """
### ODLCommonTools.getAllTypeParams

Default fallback method for getAllTypeParams.

  Input:
    mesh
    eqn
    opts

  Output:
    nothing
"""->
function getAllTypeParams(mesh::AbstractMesh, eqn::AbstractSolutionData, opts)

  error("getAllTypeParams default fallback reached. A physics-specific definition for this is needed.")

  return nothing

end