text:
### Common preprocessing

struct OnColRow{T}
    row::Int
    cols::T
    h::Vector{UInt}

    OnColRow(row::Union{Signed,Unsigned},
             cols::NTuple{<:Any, AbstractVector}, h::Vector{UInt}) =
        new{typeof(cols)}(Int(row), cols, h)
end

struct OnCol{T,N} <: AbstractVector{OnColRow{T}}
    len::Int
    cols::T
    h::Vector{UInt}

    function OnCol(cs::AbstractVector...)
        @assert length(cs) > 1
        len = length(cs[1])
        @assert all(x -> firstindex(x) == 1, cs)
        @assert all(x -> lastindex(x) == len, cs)
        new{typeof(cs), length(cs)}(len, cs, UInt[])
    end
end

Base.IndexStyle(::Type{<:OnCol}) = Base.IndexLinear()

@inline Base.size(oc::OnCol) = (oc.len,)

@inline function Base.getindex(oc::OnCol, i::Int)
    @boundscheck checkbounds(oc, i)
    return OnColRow(i, oc.cols, oc.h)
end

Base.hash(ocr1::OnColRow, h::UInt) = throw(MethodError(hash, (ocr1, h)))
@inline Base.hash(ocr1::OnColRow) = @inbounds ocr1.h[ocr1.row]

# Hashing one column at a time is faster since it can use SIMD
function _prehash(oc::OnCol)
    h = oc.h
    resize!(h, oc.len)
    fill!(h, Base.tuplehash_seed)
    for col in reverse(oc.cols)
        h .= hash.(col, h)
    end
end

# TODO: rewrite isequal and isless to use @generated
# or some other approach that would keep them efficient and avoid code duplication

Base.:(==)(x::OnColRow, y::OnColRow) = throw(MethodError(==, (x, y)))

@inline function Base.isequal(ocr1::OnColRow{<:NTuple{2, AbstractVector}},
                              ocr2::OnColRow{<:NTuple{2, AbstractVector}})
    r1 = ocr1.row
    c11, c12 = ocr1.cols
    r2 = ocr2.row
    c21, c22 = ocr2.cols

    return @inbounds isequal(c11[r1], c21[r2]) && isequal(c12[r1], c22[r2])
end

Base.isequal(ocr1::OnColRow{<:NTuple{N,AbstractVector}},
             ocr2::OnColRow{<:NTuple{N,AbstractVector}}) where {N} =
    isequal(ntuple(i -> @inbounds(ocr1.cols[i][ocr1.row]), N),
            ntuple(i -> @inbounds(ocr2.cols[i][ocr2.row]), N))

@inline function Base.isless(ocr1::OnColRow{<:NTuple{2, AbstractVector}},
                             ocr2::OnColRow{<:NTuple{2, AbstractVector}})
    r1 = ocr1.row
    c11, c12 = ocr1.cols
    r2 = ocr2.row
    c21, c22 = ocr2.cols

    c11r = @inbounds c11[r1]
    c12r = @inbounds c12[r1]
    c21r = @inbounds c21[r2]
    c22r = @inbounds c22[r2]

    isless(c11r, c21r) || (isequal(c11r, c21r) && isless(c12r, c22r))
end

@inline Base.isless(ocr1::OnColRow{<:NTuple{N,AbstractVector}},
                    ocr2::OnColRow{<:NTuple{N,AbstractVector}}) where {N} =
    isless(ntuple(i -> @inbounds(ocr1.cols[i][ocr1.row]), N),
           ntuple(i -> @inbounds(ocr2.cols[i][ocr2.row]), N))

prepare_on_col() = throw(ArgumentError("at least one on column required when joining"))
prepare_on_col(c::AbstractVector) = c
prepare_on_col(cs::AbstractVector...) = OnCol(cs...)

# Return if it is allowed to use refpool instead of the original array for joining.
# There are multiple conditions that must be met to allow for this.
# If it is allowed we are sure that missing can be used as a sentinel
check_mapping_allowed(short::AbstractVector, refarray_long::AbstractVector,
                      refpool_long, invrefpool_long) =
    !isempty(short) && !isnothing(refpool_long) && !isnothing(invrefpool_long) &&
        eltype(refarray_long) <: Union{Signed, Unsigned}

@noinline map_refarray(mapping::AbstractVector, refarray::AbstractVector,
                       ::Val{fi})  where {fi} =
    [@inbounds mapping[r - fi + 1] for r in refarray]

function map2refs(x::AbstractVector, invrefpool)
    x_refpool = DataAPI.refpool(x)
    if x_refpool isa AbstractVector{<:Integer} && 0 <= firstindex(x_refpool) <= 1
        # here we know that x_refpool is AbstractVector that allows integer indexing
        # and its firstindex must be an integer
        # if firstindex is not 0 or 1 then we fallback to slow path for safety reasons
        # all refpool we currently know have firstindex 0 or 1
        # if there is some very strange firstindex we might run into overflow issues
        # below use function barrier as mapping is not type stable
        mapping = [get(invrefpool, v, missing) for v in x_refpool]
        return map_refarray(mapping, DataAPI.refarray(x), Val(Int(firstindex(x_refpool))))
    else
        return [get(invrefpool, v, missing) for v in x]
    end
end

function preprocess_columns(joiner::DataFrameJoiner)
    right_len = length(joiner.dfr_on[!, 1])
    left_len = length(joiner.dfl_on[!, 1])
    right_shorter = right_len < left_len

    left_cols = collect(eachcol(joiner.dfl_on))
    right_cols = collect(eachcol(joiner.dfr_on))

    # if column of the longer table supports DataAPI.refpool and DataAPI.invrefpool
    # remap matching left and right columns to use refs
    if right_shorter
        for i in eachindex(left_cols, right_cols)
            rc = right_cols[i]
            lc = left_cols[i]

            lc_refs = DataAPI.refarray(lc)
            lc_refpool = DataAPI.refpool(lc)
            lc_invrefpool = DataAPI.invrefpool(lc)
            if check_mapping_allowed(rc, lc_refs, lc_refpool, lc_invrefpool)
                right_cols[i] = map2refs(rc, lc_invrefpool)
                left_cols[i] = lc_refs
            end
        end
    else
        for i in eachindex(left_cols, right_cols)
            rc = right_cols[i]
            lc = left_cols[i]

            rc_refs = DataAPI.refarray(rc)
            rc_refpool = DataAPI.refpool(rc)
            rc_invrefpool = DataAPI.invrefpool(rc)
            if check_mapping_allowed(lc, rc_refs, rc_refpool, rc_invrefpool)
                right_cols[i] = rc_refs
                left_cols[i] = map2refs(lc, rc_invrefpool)
            end
        end
    end

    disallow_sorted = false

    for (lc, rc) in zip(left_cols, right_cols)
        @assert length(lc) == left_len
        @assert length(rc) == right_len
        lc_et = nonmissingtype(eltype(lc))
        rc_et = nonmissingtype(eltype(rc))

        # special case common safe scenarios when eltype between left and right
        # column can be different or non-concrete
        lc_et <: Real && rc_et <: Real && continue
        lc_et <: AbstractString && rc_et <: AbstractString && continue

        # otherwise we require non-missing eltype of both sides to be the same and concrete
        lc_et === rc_et && isconcretetype(lc_et) && continue

        # we disallow using sorted branch for some columns that theoretically
        # could be safely sorted (e.g. having Any eltype but holding strings)
        # for safety reasons assuming that such cases will be rare in practice
        disallow_sorted = true
    end

    # TODO:
    # If DataAPI.invrefpool vectors are found in the "on" columns
    # then potentially the following optimizations can be done:
    # 1. identify rows in shorter table that should be dropped
    # 2. develop custom _innerjoin_sorted and _innerjoin_unsorted that drop rows
    #    from shorter table that do not match rows from longer table based on
    #    PooledArray refpool check
    # This optimization would significantly complicate the code (especially
    # sorted path). It should be added if in practice we find that the use case
    # is often enough and that the benefits are significant. The two cases when
    # the benefits should be expected are:
    # 1. Shorter table is sorted when we drop rows not matching longer table rows
    # 2. Shorter table does not have duplicates when we drop rows not matching
    #    longer table rows

    left_col = prepare_on_col(left_cols...)
    right_col = prepare_on_col(right_cols...)

    return left_col, right_col, right_shorter, disallow_sorted
end

### innerjoin logic

@inline function find_next_range(x::AbstractArray, start::Int, start_value)
    stop_value = start_value
    n = length(x)
    stop = start + 1
    while stop <= n
        @inbounds stop_value = x[stop]
        isequal(start_value, stop_value) || break
        stop += 1
    end
    return stop, stop_value
end

function _innerjoin_sorted(left::AbstractArray, right::AbstractArray)
    left_n = length(left)
    right_n = length(right)

    left_ixs = Int[]
    right_ixs = Int[]

    (left_n == 0 || right_n == 0) && return left_ixs, right_ixs

    # lower bound assuming we get matches
    sizehint!(left_ixs, min(left_n, right_n))
    sizehint!(right_ixs, min(left_n, right_n))

    left_cur = 1
    left_val = left[left_cur]
    left_new, left_tmp = find_next_range(left, left_cur, left_val)

    right_cur = 1
    right_val = right[right_cur]
    right_new, right_tmp = find_next_range(right, right_cur, right_val)

    while left_cur <= left_n && right_cur <= right_n
        if isequal(left_val, right_val)
            if left_new - left_cur == right_new - right_cur == 1
                push!(left_ixs, left_cur)
                push!(right_ixs, right_cur)
            else
                idx = length(left_ixs)
                left_range = left_cur:left_new - 1
                right_range = right_cur:right_new - 1
                to_grow = Base.checked_add(idx, Base.checked_mul(length(left_range),
                                                                 length(right_range)))
                resize!(left_ixs, to_grow)
                resize!(right_ixs, to_grow)
                @inbounds for right_i in right_range, left_i in left_range
                    idx += 1
                    left_ixs[idx] = left_i
                    right_ixs[idx] = right_i
                end
            end
            left_cur, left_val = left_new, left_tmp
            left_new, left_tmp = find_next_range(left, left_cur, left_val)
            right_cur, right_val = right_new, right_tmp
            right_new, right_tmp = find_next_range(right, right_cur, right_val)
        elseif isless(left_val, right_val)
            left_cur, left_val = left_new, left_tmp
            left_new, left_tmp = find_next_range(left, left_cur, left_val)
        else
            right_cur, right_val = right_new, right_tmp
            right_new, right_tmp = find_next_range(right, right_cur, right_val)
        end
    end

    return left_ixs, right_ixs
end

# optimistically assume that shorter table does not have duplicates in on column
# if this is not the case we call _innerjoin_dup
# which efficiently uses the work already done and continues with the more
# memory expensive algorithm that allows for duplicates
function _innerjoin_unsorted(left::AbstractArray, right::AbstractArray{T}) where {T}
    dict = Dict{T, Int}()

    right_len = length(right)
    # we make sure that:
    # * we do not preallocate dict of size larger than half of size of Int
    #   (this is relevant in 32 bit architectures)
    # * dict has at least 2x more slots than the number of values we
    #   might store in it to avoid reallocations of internal structures when
    #   we populate it later and to minimize the number of hash collisions;
    #   typically Dict allows for 16 probes;
    #   the value of multiplier is heuristic was determined by empirical tests
    sizehint!(dict, 2 * min(right_len, typemax(Int) >> 2))

    right isa OnCol && _prehash(right)
    left isa OnCol && _prehash(left)

    for (idx_r, val_r) in enumerate(right)
        haskey(dict, val_r) && return _innerjoin_dup(left, right, dict, idx_r)
        dict[val_r] = idx_r
    end

    left_ixs = Int[]
    right_ixs = Int[]

    # lower bound assuming we get matches
    sizehint!(left_ixs, right_len)
    sizehint!(right_ixs, right_len)

    for (idx_l, val_l) in enumerate(left)
        # we know that dict contains only positive values
        idx_r = get(dict, val_l, -1)
        if idx_r != -1
            push!(left_ixs, idx_l)
            push!(right_ixs, idx_r)
        end
    end
    return left_ixs, right_ixs
end

extrema_missing(x::AbstractVector{Missing}) = (1, 0)

function extrema_missing(x::AbstractVector{T}) where {T<:Union{Integer, Missing}}
    try
        return extrema(skipmissing(x))
    catch
        S = nonmissingtype(T)
        return S(1), S(0)
    end
end

function _innerjoin_unsorted_int(left::AbstractVector{<:Union{Integer, Missing}},
                                 right::AbstractVector{<:Union{Integer, Missing}})
    minv, maxv = extrema_missing(right)

    val_range = big(maxv) - big(minv)
    if val_range > typemax(Int) - 3 || val_range ÷ 2 > max(64, length(right)) ||
       minv < typemin(Int) + 2 || maxv > typemax(Int) - 3
       return _innerjoin_unsorted(left, right)
    end

    offset = 1 - Int(minv) # we are now sure it does not overflow
    len = Int(maxv) - Int(minv) + 2
    group_map = zeros(Int, len)

    @inbounds for (idx_r, val_r) in enumerate(right)
        i = val_r === missing ? length(group_map) : Int(val_r) + offset
        if group_map[i] > 0
            return _innerjoin_dup_int(left, right, group_map, idx_r, offset,
                                      Int(minv), Int(maxv))
        end
        group_map[i] = idx_r
    end

    left_ixs = Int[]
    right_ixs = Int[]

    right_len = length(right)
    sizehint!(left_ixs, right_len)
    sizehint!(right_ixs, right_len)

    @inbounds for (idx_l, val_l) in enumerate(left)
        if val_l === missing
            idx_r = group_map[end]
            if idx_r > 0
                push!(left_ixs, idx_l)
                push!(right_ixs, idx_r)
            end
        elseif minv <= val_l <= maxv
            idx_r = group_map[Int(val_l) + offset]
            if idx_r > 0
                push!(left_ixs, idx_l)
                push!(right_ixs, idx_r)
            end
        end
    end
    return left_ixs, right_ixs
end

# we fall back to general case if we have duplicates
# normally it should happen fast as we reuse work already done
function _innerjoin_dup(left::AbstractArray, right::AbstractArray{T},
                        dict::Dict{T, Int}, idx_r_start::Int) where {T}
    ngroups = idx_r_start - 1
    right_len = length(right)
    groups = Vector{Int}(undef, right_len)
    groups[1:ngroups] = 1:ngroups

    @inbounds for idx_r in idx_r_start:right_len
        val_r = right[idx_r]
        # we know that group ids are positive
        group_id = get(dict, val_r, -1)
        if group_id == -1
            ngroups += 1
            groups[idx_r] = ngroups
            dict[val_r] = ngroups
        else
            groups[idx_r] = group_id
        end
    end

    @assert ngroups > 0 # we should not get here with 0-length right
    return _innerjoin_postprocess(left, dict, groups, ngroups, right_len)
end

function _innerjoin_dup_int(left::AbstractVector{<:Union{Integer, Missing}},
                            right::AbstractVector{<:Union{Integer, Missing}},
                            group_map::Vector{Int}, idx_r_start::Int, offset::Int,
                            minv::Int, maxv::Int)
    ngroups = idx_r_start - 1
    right_len = length(right)
    groups = Vector{Int}(undef, right_len)
    groups[1:ngroups] = 1:ngroups

    @inbounds for idx_r in idx_r_start:right_len
        val_r = right[idx_r]
        i = val_r === missing ? length(group_map) : Int(val_r) + offset
        group_map_val = group_map[i]
        if group_map_val > 0
            groups[idx_r] = group_map_val
        else
            ngroups += 1
            groups[idx_r] = ngroups
            group_map[i] = ngroups
        end
    end

    @assert ngroups > 0 # we should not get here with 0-length right
    return _innerjoin_postprocess_int(left, group_map, groups, ngroups, right_len,
                                      offset, minv, maxv)
end

function compute_join_indices!(groups::Vector{Int},
                               starts::Vector, rperm::Vector)
    @inbounds for gix in groups
        starts[gix] += 1
    end

    cumsum!(starts, starts)

    @inbounds for (i, gix) in enumerate(groups)
        rperm[starts[gix]] = i
        starts[gix] -= 1
    end
    push!(starts, length(groups))
    return nothing
end

function _innerjoin_postprocess(left::AbstractArray, dict::Dict{T, Int},
                                groups::Vector{Int}, ngroups::Int,
                                right_len::Int) where {T}
    starts = zeros(Int, ngroups)
    rperm = Vector{Int}(undef, right_len)

    left_ixs = Int[]
    right_ixs = Int[]

    # lower bound assuming we get matches
    sizehint!(left_ixs, right_len)
    sizehint!(right_ixs, right_len)

    compute_join_indices!(groups, starts, rperm)

    n = 0
    @inbounds for (idx_l, val_l) in enumerate(left)
        group_id = get(dict, val_l, -1)
        if group_id != -1
            ref_stop = starts[group_id + 1]
            l = ref_stop - starts[group_id]
            newn = n + l
            resize!(left_ixs, newn)
            for i in n+1:n+l
                left_ixs[i] = idx_l
            end
            resize!(right_ixs, newn)
            for i in 1:l
                right_ixs[n + i] = rperm[ref_stop - i + 1]
            end
            n = newn
        end
    end

    return left_ixs, right_ixs
end

function _innerjoin_postprocess_int(left::AbstractVector{<:Union{Integer, Missing}},
                                    group_map::Vector{Int},
                                    groups::Vector{Int}, ngroups::Int, right_len::Int,
                                    offset::Int, minv::Int, maxv::Int)
    starts = zeros(Int, ngroups)
    rperm = Vector{Int}(undef, right_len)

    left_ixs = Int[]
    right_ixs = Int[]

    sizehint!(left_ixs, right_len)
    sizehint!(right_ixs, right_len)

    compute_join_indices!(groups, starts, rperm)

    n = 0
    @inbounds for (idx_l, val_l) in enumerate(left)
        if val_l === missing
            group_id = group_map[end]
        elseif minv <= val_l <= maxv
            group_id = group_map[Int(val_l) + offset]
        else
            group_id = 0
        end

        if group_id > 0
            ref_stop = starts[group_id + 1]
            l = ref_stop - starts[group_id]
            newn = n + l
            resize!(left_ixs, newn)
            for i in n+1:n+l
                left_ixs[i] = idx_l
            end
            resize!(right_ixs, newn)
            for i in 1:l
                right_ixs[n + i] = rperm[ref_stop - i + 1]
            end
            n = newn
        end
    end

    return left_ixs, right_ixs
end

function find_inner_rows(joiner::DataFrameJoiner)

    left_col, right_col, right_shorter, disallow_sorted = preprocess_columns(joiner)

    # we treat this case separately so we know we have at least one element later
    (isempty(left_col) || isempty(right_col)) && return Int[], Int[]

    # if sorting is not disallowed try using a fast algorithm that works
    # on sorted columns; if it is not run or errors fall back to the unsorted case
    # the try-catch is used to handle the case when columns on which we join
    # contain values that are not comparable
    if !disallow_sorted
        try
            if issorted(left_col) && issorted(right_col)
                return _innerjoin_sorted(left_col, right_col)
            end
        catch
            # nothing to do - one of the columns is not sortable
        end
    end

    if right_shorter
        if left_col isa AbstractVector{<:Union{Integer, Missing}} &&
           right_col isa AbstractVector{<:Union{Integer, Missing}}
            return _innerjoin_unsorted_int(left_col, right_col)
        else
            return _innerjoin_unsorted(left_col, right_col)
        end
    else
        if left_col isa AbstractVector{<:Union{Integer, Missing}} &&
           right_col isa AbstractVector{<:Union{Integer, Missing}}
            return reverse(_innerjoin_unsorted_int(right_col, left_col))
        else
            return reverse(_innerjoin_unsorted(right_col, left_col))
        end
    end

    error("unreachable reached")
end

### semijoin logic

function _semijoin_sorted(left::AbstractArray, right::AbstractArray,
                          seen_rows::AbstractVector{Bool})
    left_n = length(left)
    right_n = length(right)

    @assert left_n > 0 && right_n > 0

    left_cur = 1
    left_val = left[left_cur]
    left_new, left_tmp = find_next_range(left, left_cur, left_val)

    right_cur = 1
    right_val = right[right_cur]
    right_new, right_tmp = find_next_range(right, right_cur, right_val)

    while left_cur <= left_n && right_cur <= right_n
        if isequal(left_val, right_val)
            seen_rows[left_cur:left_new - 1] .= true
            left_cur, left_val = left_new, left_tmp
            left_new, left_tmp = find_next_range(left, left_cur, left_val)
            right_cur, right_val = right_new, right_tmp
            right_new, right_tmp = find_next_range(right, right_cur, right_val)
        elseif isless(left_val, right_val)
            left_cur, left_val = left_new, left_tmp
            left_new, left_tmp = find_next_range(left, left_cur, left_val)
        else
            right_cur, right_val = right_new, right_tmp
            right_new, right_tmp = find_next_range(right, right_cur, right_val)
        end
    end

    return seen_rows
end

# optimistically assume that shorter table does not have duplicates in on column
# if this is not the case we call _semijoin_dup
# which efficiently uses the work already done and continues with the more
# memory expensive algorithm that allows for duplicates
# note that in semijoin and antijoin we do not have to do it if right table is
# shorter as then we process left table row by row anyway
function _semijoin_unsorted(left::AbstractArray, right::AbstractArray{T},
                            seen_rows::AbstractVector{Bool},
                            right_shorter::Bool) where {T}
    right_len = length(right)
    right isa OnCol && _prehash(right)
    left isa OnCol && _prehash(left)

    if right_shorter
        @assert length(left) == length(seen_rows)
        set = Set{T}()
        sizehint!(set, 2 * min(right_len, typemax(Int) >> 2))
        for val_r in right
            push!(set, val_r)
        end
        @inbounds for (idx_l, val_l) in enumerate(left)
            seen_rows[idx_l] = val_l in set
        end
    else
        @assert length(right) == length(seen_rows)
        dict = Dict{T, Int}()
        sizehint!(dict, 2 * min(right_len, typemax(Int) >> 2))
        for (idx_r, val_r) in enumerate(right)
            haskey(dict, val_r) && return _semijoin_dup(left, right, dict, idx_r,
                                                        seen_rows)
            dict[val_r] = idx_r
        end
        @inbounds for (idx_l, val_l) in enumerate(left)
            # we know that dict contains only positive values
            idx_r = get(dict, val_l, -1)
            if idx_r != -1
                seen_rows[idx_r] = true
            end
        end
    end

    return seen_rows
end

function _semijoin_unsorted_int(left::AbstractVector{<:Union{Integer, Missing}},
                                right::AbstractVector{<:Union{Integer, Missing}},
                                seen_rows::AbstractVector{Bool},
                                right_shorter::Bool)
    minv, maxv = extrema_missing(right)

    val_range = big(maxv) - big(minv)
    if val_range > typemax(Int) - 3 || val_range ÷ 2 > max(64, length(right)) ||
       minv < typemin(Int) + 2 || maxv > typemax(Int) - 3
       return _semijoin_unsorted(left, right, seen_rows, right_shorter)
    end

    offset = 1 - Int(minv) # we are now sure it does not overflow
    len = Int(maxv) - Int(minv) + 2
    group_map = zeros(Int, len)

    if right_shorter
        @inbounds for (idx_r, val_r) in enumerate(right)
            i = val_r === missing ? length(group_map) : Int(val_r) + offset
            group_map[i] = idx_r
        end
        @inbounds for (idx_l, val_l) in enumerate(left)
            if val_l === missing
                idx_r = group_map[end]
                seen_rows[idx_l] = idx_r > 0
            elseif minv <= val_l <= maxv
                idx_r = group_map[Int(val_l) + offset]
                seen_rows[idx_l] = idx_r > 0
            end
        end
    else
        @inbounds for (idx_r, val_r) in enumerate(right)
            i = val_r === missing ? length(group_map) : Int(val_r) + offset
            if group_map[i] > 0
                return _semijoin_dup_int(left, right, group_map, idx_r, offset,
                                         Int(minv), Int(maxv), seen_rows)
            end
            group_map[i] = idx_r
        end
        @inbounds for (idx_l, val_l) in enumerate(left)
            if val_l === missing
                idx_r = group_map[end]
                if idx_r > 0
                    seen_rows[idx_r] = true
                end
            elseif minv <= val_l <= maxv
                idx_r = group_map[Int(val_l) + offset]
                if idx_r > 0
                    seen_rows[idx_r] = true
                end
            end
        end
    end

    return seen_rows
end

# we fall back to general case if we have duplicates
# normally it should happen fast as we reuse work already done
function _semijoin_dup(left::AbstractArray, right::AbstractArray{T},
                       dict::Dict{T, Int}, idx_r_start::Int,
                       seen_rows::AbstractVector{Bool}) where {T}
    ngroups = idx_r_start - 1
    right_len = length(right)
    groups = Vector{Int}(undef, right_len)
    groups[1:ngroups] = 1:ngroups

    @inbounds for idx_r in idx_r_start:right_len
        val_r = right[idx_r]
        # we know that group ids are positive
        group_id = get(dict, val_r, -1)
        if group_id == -1
            ngroups += 1
            groups[idx_r] = ngroups
            dict[val_r] = ngroups
        else
            groups[idx_r] = group_id
        end
    end

    @assert ngroups > 0 # we should not get here with 0-length right
    @assert length(right) == length(seen_rows)
    return _semijoin_postprocess(left, dict, groups, ngroups, right_len,
                                 seen_rows)
end

function _semijoin_dup_int(left::AbstractVector{<:Union{Integer, Missing}},
                           right::AbstractVector{<:Union{Integer, Missing}},
                           group_map::Vector{Int}, idx_r_start::Int, offset::Int,
                           minv::Int, maxv::Int, seen_rows::AbstractVector{Bool})
    ngroups = idx_r_start - 1
    right_len = length(right)
    groups = Vector{Int}(undef, right_len)
    groups[1:ngroups] = 1:ngroups

    @inbounds for idx_r in idx_r_start:right_len
        val_r = right[idx_r]
        i = val_r === missing ? length(group_map) : Int(val_r) + offset
        group_map_val = group_map[i]
        if group_map_val > 0
            groups[idx_r] = group_map_val
        else
            ngroups += 1
            groups[idx_r] = ngroups
            group_map[i] = ngroups
        end
    end

    @assert ngroups > 0 # we should not get here with 0-length right
    @assert length(right) == length(seen_rows)
    return _semijoin_postprocess_int(left, group_map, groups, ngroups, right_len,
                                     offset, minv, maxv, seen_rows)
end

function _semijoin_postprocess(left::AbstractArray, dict::Dict{T, Int},
                               groups::Vector{Int}, ngroups::Int, right_len::Int,
                               seen_rows::AbstractVector{Bool}) where {T}
    starts = zeros(Int, ngroups)
    rperm = Vector{Int}(undef, right_len)

    compute_join_indices!(groups, starts, rperm)

    @inbounds for (idx_l, val_l) in enumerate(left)
        group_id = get(dict, val_l, -1)
        if group_id != -1
            ref_stop = starts[group_id + 1]
            l = ref_stop - starts[group_id]
            for i in 1:l
                seen_rows[rperm[ref_stop - i + 1]] = true
            end
        end
    end

    return seen_rows
end

function _semijoin_postprocess_int(left::AbstractVector{<:Union{Integer, Missing}},
                                   group_map::Vector{Int},
                                   groups::Vector{Int}, ngroups::Int, right_len::Int,
                                   offset::Int, minv::Int, maxv::Int,
                                   seen_rows::AbstractVector{Bool})
    starts = zeros(Int, ngroups)
    rperm = Vector{Int}(undef, right_len)

    compute_join_indices!(groups, starts, rperm)

    @inbounds for (idx_l, val_l) in enumerate(left)
        if val_l === missing
            group_id = group_map[end]
        elseif minv <= val_l <= maxv
            group_id = group_map[Int(val_l) + offset]
        else
            group_id = 0
        end

        if group_id > 0
            ref_stop = starts[group_id + 1]
            l = ref_stop - starts[group_id]
            for i in 1:l
                seen_rows[rperm[ref_stop - i + 1]] = true
            end
        end
    end

    return seen_rows
end

function find_semi_rows(joiner::DataFrameJoiner)
    left_col, right_col, right_shorter, disallow_sorted = preprocess_columns(joiner)

    seen_rows = falses(length(left_col))

    # we treat this case separately so we know we have at least one element later
    (isempty(left_col) || isempty(right_col)) && return falses(length(left_col))

    # if sorting is not disallowed try using a fast algorithm that works
    # on sorted columns; if it is not run or errors fall back to the unsorted case
    # the try-catch is used to handle the case when columns on which we join
    # contain values that are not comparable
    if !disallow_sorted
        try
            if issorted(left_col) && issorted(right_col)
                return _semijoin_sorted(left_col, right_col, seen_rows)
            end
        catch
            # nothing to do - one of the columns is not sortable
        end
    end

    if right_shorter
        if left_col isa AbstractVector{<:Union{Integer, Missing}} &&
           right_col isa AbstractVector{<:Union{Integer, Missing}}
            return _semijoin_unsorted_int(left_col, right_col, seen_rows, right_shorter)
        else
            return _semijoin_unsorted(left_col, right_col, seen_rows, right_shorter)
        end
    else
        if left_col isa AbstractVector{<:Union{Integer, Missing}} &&
           right_col isa AbstractVector{<:Union{Integer, Missing}}
            return _semijoin_unsorted_int(right_col, left_col, seen_rows, right_shorter)
        else
            return _semijoin_unsorted(right_col, left_col, seen_rows, right_shorter)
        end
    end
    error("unreachable reached")
end


text:
module Simulate_Pendulum

using  Modia3D
import Modia3D.ModiaMath

# Solids
groundMaterial = Modia3D.Material(color="DarkGreen", transparency=0.5)
# massProperties = Modia3D.MassProperties()
vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)
vmat2 = Modia3D.Material(color="Red")

rCM = zeros(3)
#rCM[2] = 5.0
#I2 = fill(100.0,(9))
I = fill(0.0,(3,3))
#I[3] = 50.0
#I[7] = 80.0
#println("I = ",I)
massProp = Modia3D.MassProperties(5.0, rCM, I)

@assembly Pendulum(;Lx = 1.0, Ly=0.2*Lx, Lz=0.2*Lx, m=1.0) begin
   world  = Modia3D.Object3D(Modia3D.CoordinateSystem(0.5*Lx))

   # Pendulum
   body   = Modia3D.Object3D(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), massProp, vmat1))
   frame1 = Modia3D.Object3D(body; r=[-Lx/2, 0.0, 0.0])
   cyl    = Modia3D.Object3D(frame1,Modia3D.Cylinder(Ly/2,1.2*Ly; material=vmat2))

   # Connect pendulum to world
   rev    = Modia3D.Revolute(world, frame1)
end

const Lx = 1.6
const m  = 0.5
# Modia3D.visualizeAssembly!( Pendulum(Lx=Lx, m=m) )

pendulum = Pendulum(Lx=Lx, m=m, sceneOptions=Modia3D.SceneOptions(visualizeFrames=true, defaultFrameLength=0.3))
model    = Modia3D.SimulationModel( pendulum )
result   = ModiaMath.simulate!(model, stopTime=10.0, tolerance=1e-6,log=false)

ModiaMath.plot(result, ["rev.phi", "rev.w", "rev.a"])

println("... success of Simulate_Pendulum.jl!")
end

ext: jl


text:
using DifferentialEquations

function time_stepping_example(x_min, x_max, M, T, N, rho, sigma_bar)
    M_bar = M + 2
    x = linspace(x_min, x_max, M)
    Delta = x[2]-x[1]
    mu_tilde_T = -0.1+T+0.1*x
    sigma_tilde_T = sigma_bar*x
    c_tilde = exp.(x)

    L_1_plus = Tridiagonal(zeros(M_bar-1), -1 * ones(M_bar), ones(M_bar-1))[2: end-1, :]/Delta
    L_1_plus[:, 2] = L_1_plus[:, 2] + L_1_plus[:, 1]
    L_1_plus[:, end-1] = L_1_plus[:, end-1] + L_1_plus[:, end]
    L_1_plus = L_1_plus[:, 2: end-1]

    L_2 = Tridiagonal(ones(M_bar-1), -2 * ones(M_bar), ones(M_bar-1))[2: end-1, :]/(Delta^2)
    L_2[:, 2] = L_2[:, 2] + L_2[:, 1]
    L_2[:, end-1] = L_2[:, end-1] + L_2[:, end]
    L_2 = L_2[:, 2: end-1]

    L_tilde = Diagonal(rho*ones(M)) - (Diagonal(mu_tilde_T)*L_1_plus + Diagonal(sigma_tilde_T.^2/2)*L_2)

    u_T = L_tilde\c_tilde

    function f(du,u,p,t)
        mu_tilde = -0.1+0.01*t+0.1*x#mu_tilde_T#-0.1+t+0.1*x
        L = Diagonal(rho*ones(M)) - (Diagonal(mu_tilde)*L_1_plus + Diagonal(sigma_tilde_T.^2/2)*L_2)
        A_mul_B!(du,L,u)
        du .-= c_tilde
    end

    tstops=(0.0, T)
    prob = ODEProblem(f, u_T, tstops)
    solve(prob)
end

x_min = 0.01
x_max = 1
M = 40
T = 0.3
N = 10 #This is not being used right now.
rho = 0.05
sigma_bar = 0.1

sol = time_stepping_example(x_min, x_max, M, T, N, rho, sigma_bar)

using Plots
plot(sol)


text:
include("./genericLogs.jl")
include("./sopIO.jl")


text:
module JuliaOnLambda

using AWS
using Mocking

@service ECR

const DEFAULT_PROJECT_PATH = joinpath(@__DIR__, "..")
const DEFAULT_DOCKERFILE = joinpath(DEFAULT_PROJECT_PATH, "Dockerfile")
const REPO_NOT_FOUND_EXCEPTION = "RepositoryNotFoundException"

include("utilities.jl")

"""
Get the ECR Repo called `repository_name`. If it does not exist, create a new one.
"""
function _get_create_ecr_repo(repository_name::AbstractString)
    try
        response = @mock ECR.describe_repositories(
            Dict("repositoryNames" => [repository_name])
        )
        return response["repositories"][1]["repositoryUri"]
    catch e
        if e.code == REPO_NOT_FOUND_EXCEPTION
            @info "Creating ECR repository, $repository_name"
            response = @mock ECR.create_repository(repository_name)
            return response["repository"]["repositoryUri"]
        else
            rethrow(e)
        end
    end
end

"""
Build and tag the Docker image.
"""
function _create_docker_image(
    image::AbstractString,
    tag::AbstractString;
    dockerfile_path::AbstractString=DEFAULT_DOCKERFILE,
)
    image_uri = _image_uri(image, tag)

    _build_docker_image(image_uri; dockerfile_path=dockerfile_path)

    return nothing
end

function _tag_docker_image(
    image::AbstractString, tag::AbstractString, repository_uri::AbstractString
)
    return run(`docker tag $(image):$(tag) $repository_uri`)
end

"""
Upload Docker image
"""
function _upload_docker_image(repository_name::AbstractString)
    return run(`docker push $repository_name`)
end

end


text:
function get_margin(obj::FieldOfView)
    return jcall(obj, "getMargin", jdouble, ())
end


text:
@enum DirectionType begin
    Forward
    Backward
end

"""
    Track{BCFwd,BOut,DFwd,DBwd,T<:Real}

Represents a neutron trajectory across the domain, with certain azimuthal angle `ϕ`, entry
and exit points `p` and `q`, length `ℓ` and formed by `segments` coming from its
segmentation.
"""
mutable struct Track{BCFwd,BCBwd,DFwd,DBwd,T<:Real}
    uid::Int
    azim_idx::Int
    track_idx::Int

    p::Point2D{T}
    q::Point2D{T}

    ϕ::T
    ℓ::T

    ABC::SVector{3,T}
    segments::Vector{Segment{T}}

    next_track_fwd::Track
    next_track_bwd::Track

    function Track{BCFwd,BCBwd,DFwd,DBwd}(
        uid::Int,
        azim_idx::Int,
        track_idx::Int,
        p::Point2D{T},
        q::Point2D{T},
        ϕ::T,
        ℓ::T,
        ABC::SVector{3,T},
        segments::Vector{Segment{T}}
    ) where {BCFwd,BCBwd,DFwd,DBwd,T}
        track = new{BCFwd,BCBwd,DFwd,DBwd,T}(
            uid, azim_idx, track_idx, p, q, ϕ, ℓ, ABC, segments
        )
        track.next_track_fwd = track
        track.next_track_bwd = track
        return track
    end
end

universal_id(track::Track) = track.uid
azim_idx(track::Track) = track.azim_idx
track_idx(track::Track) = track.track_idx
bc_fwd(::Track{BCFwd}) where {BCFwd} = BCFwd
bc_bwd(::Track{BCFwd,BCBwd}) where {BCFwd,BCBwd} = BCBwd
dir_next_track_fwd(::Track{BCFwd,BCBwd,DFwd}) where {BCFwd,BCBwd,DFwd} = DFwd
dir_next_track_bwd(::Track{BCFwd,BCBwd,DFwd,DBwd}) where {BCFwd,BCBwd,DFwd,DBwd} = DBwd

function show(io::IO, track::Track)
    @unpack ϕ, p, q, ℓ, segments = track
    # println(io, typeof(t))
    println(io, "  Azimuthal angle: ", round(rad2deg(ϕ), digits=2))
    println(io, "  Entry point: ", p)
    println(io, "  Exit point: ", q)
    println(io, "  Length: ", ℓ)
    println(io, "  # of segments: ", length(segments)) # if it is zero, run segmentize!
    println(io, "  Boundary fwd: ", bc_fwd(track))
    print(io,   "  Boundary bwd: ", bc_bwd(track))
    # avoid printing circular references (since it is a cyclic ray tracing...)
    # println(io, "  Next track fwd: ", next_track_fwd)
    # print(io,   "  Next track fwd: ", next_track_bwd)
end

# since we are in active development of the package and there might be unconsidered cases in
# the segmentation of tracks, let's have this maximum number of iterations.
const MAX_ITER = 10_000

function _segmentize_track!(t, track::Track, k::Int, rtol::Real) # t::TrackGenerator
    @unpack mesh, tiny_step = t
    @unpack ϕ, segments = track

    # since we are going to compute them, empty!
    empty!(segments)

    # move a tiny step in ϕ direction to get inside the mesh
    xp = advance_step(track.p, tiny_step, ϕ)

    i = 0
    element = -1
    prev_element = -1
    while i < MAX_ITER

        # find the element or cell where `xp` lies
        element = find_element(mesh, xp)

        # we might be at the boundary of the mesh
        if inboundary(mesh, xp, tiny_step)
            if isempty(segments)
                # if we just started to segmentize, move a tiny step forward
                xp = advance_step(xp, tiny_step, ϕ)
                continue
            else
                # or we just finished, so leave
                break
            end
        end

        # if we are inside the domain and `find_element` did not return an element, we might
        # be dealing with a deformed mesh, so we might need to increase the knn search.
        if isequal(element, -1)
            element = find_element(mesh, xp, k)
            if isequal(element, -1)
                error("Try increasing `k`. If the problem persists, raise an issue, this " *
                      "might be a case that hasn't been presented before.")
            end
        end

        # move a step if the new element is the same as the previous one
        if isequal(prev_element, element)
            xp = advance_step(xp, tiny_step, ϕ)
            continue
        end

        # compute intersections between track and the element
        p, q = intersections(mesh, element, track)

        # we might be in a vertex, so we just continue
        if isapprox(p, q)
            xp = advance_step(xp, tiny_step, ϕ)
            continue
        end

        segment = Segment(p, q, element)
        push!(segments, segment)

        # update new starting point and previous element
        xp = advance_step(q, tiny_step, ϕ)
        prev_element = element

        i += 1 # just for safety
    end

    if !isapprox(track.ℓ, sum(ℓ.(track.segments)); rtol=rtol)
        error("Track with `uid` $(string(track.uid)) has a length that do not match the " *
              "sum of its segments lengths with the provided tolerance `rtol`. Check " *
              "whether this is an actual error or increase `rtol`.")
    end

    return nothing
end


text:
# cuSPARSE functions for managing the library

function cusparseCreate()
    handle = Ref{cusparseHandle_t}()
    res = @retry_reclaim err->isequal(err, CUSPARSE_STATUS_ALLOC_FAILED) ||
                              isequal(err, CUSPARSE_STATUS_NOT_INITIALIZED) begin
        unsafe_cusparseCreate(handle)
    end
    if res != CUSPARSE_STATUS_SUCCESS
         throw_api_error(res)
    end
    handle[]
end

function cusparseGetProperty(property::libraryPropertyType)
    value_ref = Ref{Cint}()
    cusparseGetProperty(property, value_ref)
    value_ref[]
end

version() = VersionNumber(cusparseGetProperty(CUDA.MAJOR_VERSION),
                          cusparseGetProperty(CUDA.MINOR_VERSION),
                          cusparseGetProperty(CUDA.PATCH_LEVEL))


text:
"""
Get global settings.
"""
function get_options()

    opt = Dict()

    opt["base_folder"] = "/data02/Ican/vic_sim/past_1km"
    opt["target_folder"] = "/data02/Ican/vic_sim/jan_past_1031"
    opt["eval_folder"] = "/data02/Ican/vic_sim/jan_eval_1031"

    opt["stat_sel"] = ["191.2", "122.11", "2.32", "2.279", "224.1", "12.70", "62.5", "22.4"]

    opt["resolutions"] = ["1km", "5km", "10km", "25km", "50km"]

    opt["startyear"] = 1982
    opt["endyear"] = 2012
    opt["timestep"] = 3
    opt["output_force"] = "FALSE"
    opt["full_energy"] = "FALSE"
    opt["output_binary"] = "TRUE"

    return opt

end

"""
Get wsh_info file.
"""
function get_wsh_info()

    @load Pkg.dir("IcanProj", "data", "wsh_info.jld2") wsh_info

    return wsh_info

end


text:
##########################################################################
# Geography
immutable Position
    x::Float64
    y::Float64
end

+(A::Position, B::Position) = Position(A.x + B.x, A.y + B.y)
-(A::Position, B::Position) = Position(A.x - B.x, A.y - B.y)
rotate(p::Position, θ) = Position(cos(θ)*p.x - sin(θ)*p.y, sin(θ)*p.x + cos(θ)*p.y)

get_distance(A::Position, B::Position) = sqrt((A.x - B.x)^2 + (A.y - B.y)^2)
get_angle(A::Position, B::Position) = atan2(A.y - B.y, A.x - B.x)

immutable Velocity
    x::Float64
    y::Float64
end

##########################################################################
# Systems
abstract System

typealias AuxPrecodingParams Dict{ASCIIString, Any} # belongs in precoding.jl
typealias AuxAssignmentParams Dict{ASCIIString, Any} # belongs in assignment.jl

type SinglecarrierSystem <: System
    carrier_frequency::Float64
    bandwidth::Float64

    aux_precoding_params::AuxPrecodingParams
    aux_assignment_params::AuxAssignmentParams
end
SinglecarrierSystem() =
    SinglecarrierSystem(0, 0, AuxPrecodingParams(), AuxAssignmentParams())
SinglecarrierSystem(carrier_frequency, bandwidth) =
    SinglecarrierSystem(carrier_frequency, bandwidth, AuxPrecodingParams(), AuxAssignmentParams())

type MulticarrierSystem <: System
    carrier_frequency::Float64
    bandwidth::Float64

    num_subcarriers::Int

    aux_precoding_params::AuxPrecodingParams
    aux_assignment_params::AuxAssignmentParams
end
MulticarrierSystem() =
    MulticarrierSystem(0, 0, 1, AuxPrecodingParams(), AuxAssignmentParams())
MulticarrierSystem(carrier_frequency, bandwidth, num_subcarriers) =
    MulticarrierSystem(carrier_frequency, bandwidth, num_subcarriers, AuxPrecodingParams(), AuxAssignmentParams())

##########################################################################
# Antenna params
abstract AntennaParams

##########################################################################
# Propagation environments
abstract PropagationEnvironment
abstract PropagationEnvironmentState

immutable SimpleLargescaleFadingEnvironment <: PropagationEnvironment
    pathloss_alpha::Float64
    pathloss_beta::Float64
    penetration_loss_dB::Float64
    shadow_sigma_dB::Float64
end

immutable SimpleLargescaleFadingEnvironmentState <: PropagationEnvironmentState
    shadow_realizations_dB::Vector{Float64}
    LoSs::BitArray
end

##########################################################################
# Nodes
abstract Node

abstract BS <: Node

type CanonicalBS <: BS
    num_antennas::Int

    transmit_power::Float64
end

type PhysicalBS{AntennaParams_t <: AntennaParams} <: BS
    num_antennas::Int
    position::Position

    transmit_power::Float64

    antenna_params::AntennaParams_t
end

abstract MS <: Node

type CanonicalMS <: MS
    num_antennas::Int

    user_priority::Float64
    num_streams::Int

    receiver_noise_power::Float64
end

type PhysicalMS{PropagationEnvironmentState_t <: PropagationEnvironmentState} <: MS
    num_antennas::Int
    position::Position
    velocity::Velocity

    user_priority::Float64
    num_streams::Int

    antenna_gain_dB::Float64
    noise_figure::Float64
    propagation_environment_state::PropagationEnvironmentState_t
end

get_num_antennas(node) = node.num_antennas

get_user_priority(MS) = MS.user_priority
set_user_priority!(MS, α) = (MS.user_priority = α)

get_num_streams(MS::MS) = MS.num_streams
set_num_streams!(MS::MS, d) = (MS.num_streams = d)

get_receiver_noise_power(MS::CanonicalMS, network) = MS.receiver_noise_power
set_receiver_noise_power!(MS::CanonicalMS, sigma2, network) = (MS.receiver_noise_power = sigma2)
get_receiver_noise_power_dBm(MS::CanonicalMS, network) = 10*log10(get_receiver_noise_power(MS, network))
set_receiver_noise_power_dBm!(MS::CanonicalMS, sigma2_dBm, network) = set_receiver_noise_power(MS, 10^(sigma2_dBm/10), network)

get_receiver_noise_power(MS::PhysicalMS, network) = 10^(get_receiver_noise_power_dBm(MS, network)/10)
set_receiver_noise_power!(MS::PhysicalMS, sigma2, network) = (MS.noise_figure = 174. + 10*log10(sigma2) - 10*log10(network.system.bandwidth))
get_receiver_noise_power_dBm(MS::PhysicalMS, network) = (-174. + 10*log10(network.system.bandwidth) + MS.noise_figure)
set_receiver_noise_power_dBm!(MS::PhysicalMS, sigma2_dBm, network) = (MS.noise_figure = 174. + sigma2_dBm - 10*log10(network.system.bandwidth))

get_transmit_power(BS) = BS.transmit_power
set_transmit_power!(BS, P) = (BS.transmit_power = P)
get_transmit_power_dBm(BS) = 10*log10(get_transmit_power(BS))
set_transmit_power_dBm!(BS, P_dBm) = set_transmit_power!(BS, 10^(P_dBm/10))

get_distance(MS, BS) = get_distance(MS.position, BS.position)

##########################################################################
# Networks
abstract Network
abstract CanonicalNetwork <: Network
abstract PhysicalNetwork <: Network

typealias AuxNetworkParams Dict{ASCIIString, Any}

get_num_MSs(network) = length(network.MSs)
get_num_BSs(network) = length(network.BSs)

get_assignment(network) = network.assignment
set_assignment!(network, assignment) = (network.assignment = assignment)

has_aux_network_param(network, k) = haskey(network.aux_network_params, k)
get_aux_network_param(network, k) = (network.aux_network_params[k])
set_aux_network_param!(network, v, k) = (network.aux_network_params[k] = v)
get_aux_network_params(network) = network.aux_network_params
set_aux_network_params!(network, additional) = merge!(network.aux_network_params, additional)

has_aux_precoding_param(network, k) = haskey(network.aux_precoding_params, k)
get_aux_precoding_param(network, k) = (network.system.aux_precoding_params[k])
set_aux_precoding_param!(network, v, k) = (network.system.aux_precoding_params[k] = v)
get_aux_precoding_params(network) = network.system.aux_precoding_params
set_aux_precoding_params!(network, additional) = merge!(network.system.aux_precoding_params, additional)

has_aux_assignment_param(network, k) = haskey(network.aux_assignment_params, k)
get_aux_assignment_param(network, k) = (network.system.aux_assignment_params[k])
set_aux_assignment_param!(network, v, k) = (network.system.aux_assignment_params[k] = v)
get_aux_assignment_params(network) = network.system.aux_assignment_params
set_aux_assignment_params!(network, additional) = merge!(network.system.aux_assignment_params, additional)

get_num_MS_antennas(network) = [ get_num_antennas(network.MSs[k]) for k = 1:get_num_MSs(network) ]
get_num_BS_antennas(network) = [ get_num_antennas(network.BSs[i]) for i = 1:get_num_BSs(network) ]

get_transmit_powers(network) = [ get_transmit_power(network.BSs[i]) for i = 1:get_num_BSs(network) ]
set_transmit_powers!(network, P::Real) = (for BS in network.BSs; set_transmit_power!(BS, P); end)
set_transmit_powers!(network, Ps::Vector) = (for i = 1:length(Ps); set_transmit_power!(network.BSs[i], Ps[i]); end)
get_transmit_powers_dBm(network) = [ get_transmit_power_dBm(network.BSs[i]) for i = 1:get_num_BSs(network) ]
set_transmit_powers_dBm!(network, P_dBm::Real) = (for BS in network.BSs; set_transmit_power_dBm!(BS, P_dBm); end)
set_transmit_powers_dBm!(network, Ps_dBm::Vector) = (for i = 1:length(Ps_dBm); set_transmit_power_dBm!(network.BSs[i], Ps_dBm[i]); end)

get_receiver_noise_powers(network) = [ get_receiver_noise_power(network.MSs[k], network) for k = 1:get_num_MSs(network) ]
set_receiver_noise_powers!(network, sigma2::Real) = (for MS in network.MSs; set_receiver_noise_power!(MS, sigma2, network); end)
set_receiver_noise_powers!(network, sigma2s::Vector) = (for k = 1:length(sigma2s); set_receiver_noise_power!(network.MSs[k], sigma2s[k], network); end)
get_receiver_noise_powers_dBm(network) = [ get_receiver_noise_power_dBm(network.MSs[k], network) for k = 1:get_num_MSs(network) ]
set_receiver_noise_powers_dBm!(network, sigma2_dBm::Real) = (for MS in network.MSs; set_receiver_noise_power_dBm!(MS, sigma2_dBm, network); end)
set_receiver_noise_powers_dBm!(network, sigma2s_dBm::Vector) = (for k = 1:length(sigma2s_dBm); set_receiver_noise_power_dBm!(network.MSs[k], sigma2s_dBm[k], network); end)

get_user_priorities(network) = [ get_user_priority(network.MSs[k]) for k = 1:get_num_MSs(network) ]
set_user_priorities!(network, α::Real) = (for MS in network.MSs; set_user_priority!(MS, α); end)
set_user_priorities!(network, α::Vector) = (for k = 1:length(α); set_user_priority!(network.MSs[k], α[k]); end)

get_num_streams(network::Network) = [ get_num_streams(network.MSs[k]) for k = 1:get_num_MSs(network) ]
set_num_streams!(network::Network, d::Int) = (for MS in network.MSs; set_num_streams(MS, d); end)
set_num_streams!(network::Network, ds::Vector) = (for k = 1:length(ds); set_num_streams(network.MSs[k], ds[k]); end)

require_equal_num_MS_antennas(network) = (Ns = get_num_MS_antennas(network); all(Ns .== Ns[1]) || Lumberjack.error("MSs must all have the same number of antennas."))
require_equal_num_BS_antennas(network) = (Ms = get_num_BS_antennas(network); all(Ms .== Ms[1]) || Lumberjack.error("BSs must all have the same number of antennas."))
require_single_antenna_MSs(network) = (Ns = get_num_MS_antennas(network); all(Ns .== 1) || Lumberjack.error("MSs must not have multiple antennas."))
require_single_antenna_BSs(network) = (Ms = get_num_BS_antennas(network); all(Ms .== 1) || Lumberjack.error("BSs must not have multiple antennas."))
require_equal_num_streams(network) = (ds = get_num_streams(network); all(ds .== ds[1]) || Lumberjack.error("MSs must all have the same number of streams."))
require_single_stream(network) = (ds = get_num_streams(network); all(ds .== 1) || Lumberjack.error("MSs must not have multiple streams."))

get_distances(network::PhysicalNetwork) = [ get_distance(network.MSs[k], network.BSs[i]) for k = 1:get_num_MSs(network), i = 1:get_num_BSs(network) ]
get_angles(network::PhysicalNetwork) = [ get_angle(network.MSs[k], network.BSs[i]) for k = 1:get_num_MSs(network), i = 1:get_num_BSs(network) ]

##########################################################################
# OmnidirectionalAntennaParams (not angle dependent)
immutable OmnidirectionalAntennaParams <: AntennaParams
    antenna_gain_dB::Float64
end

get_antenna_gain(antenna_params::OmnidirectionalAntennaParams) =
    10^(antenna_params.antenna_gain_dB/10)
get_antenna_gain(antenna_params::OmnidirectionalAntennaParams, angle) =
    get_antenna_gain(antenna_params)

get_angle{AntennaParams_t <: OmnidirectionalAntennaParams}(MS, BS::PhysicalBS{AntennaParams_t}) =
    get_angle(MS.position, BS.position)


text:
# %%
using Plots
# %%
function test()
    plot()
	plot!(rand(10);l1...)
    plot!(rand(10);l2...)
    plot!(rand(10);l3...)
    plot!(rand(10);l4...)
	plot!(rand(10);l5...)
end
# %%
module MyPlot
export l1b, l1, l2, l3, l4, l5
using Plots
using Measurements
	cp = palette(:default)
	l1b = Dict(:color => :black, :linewidth => 4.0)
	l1 = Dict(:color => :black, :linewidth => 2.0)
	l2 = Dict(:color => cp[1])
	l3 = Dict(:color => cp[2], :linestyle => :dash)
	l4 = Dict(:color => cp[3], :linestyle => :dashdot)
	l5 = Dict(:color => cp[4], :linestyle => :dashdotdot)
	function set()
		default(size=(600, 400), foreground_color_legend=nothing, margin=3Plots.mm, linewidth=1.5, legendtitlefontsize=10.0, minorgrid = true)
		default(:bglegend,plot_color(default(:bg), 0.2))
		default(:fglegend,plot_color(default(:bg), 0.2))
	end
end
# %%


text:
module DataUtils

export sampleimages
function sampleimages(images::Array{Float64,3},patchwidth::Int,patchheight::Int,numsamples::Int; scalevariance=true)
    width, height = size(images[:,:,1])
    array::Array{Float64,2} = zeros(patchwidth*patchheight,numsamples)
    for index=1:numsamples
        image_index = rand(1:size(images,3))
        x = rand(1:width-patchwidth+1)
        y = rand(1:height-patchheight+1)
        sample = images[x:x+patchwidth-1,y:y+patchheight-1,image_index]
        array[:,index] = reshape(sample,patchwidth*patchheight)
        array[:,index] -= mean(array[:,index]) #subtract mean
    end

    if scalevariance
        # rescale images to fit in range 0.1 to 0.9
        stddev = std(array)
        array = max(min(array,3*stddev),-3*stddev) / (3*stddev)
        array = (array + 1.0) * 0.4 + 0.1
    end
    return array
end

import Plotly
export displaynetwork_plotly
function displaynetwork_plotly(A,filename,username,userkey)
        m,n = size(A)
        sz = int(sqrt(m))
        A -= mean(A)
        layout = [
            "autosize" => false,
            "width" => 500,
            "height"=> 500
        ]

        gridsize = int(ceil(sqrt(n)))
        buffer = 1
        griddata = ones(gridsize*(sz+1)+1,gridsize*(sz+1)+1)
        index = 1
        for i = 1:gridsize
                for j = 1:gridsize
                        if index > n
                                continue
                        end
                        columnlimit = maximum(abs(A[:,index]))
                        griddata[buffer+(i-1)*(sz+buffer)+(1:sz),buffer+(j-1)*(sz+buffer)+(1:sz)] = reshape(A[:,index],sz,sz)/columnlimit
                        index += 1
                end
        end

        Plotly.signin(username, userkey)
        data = [
          [
            "z" => griddata,
            "colorscale" => "Greys",
            "type" => "heatmap"
          ]
        ]
        response = Plotly.plot(data, ["layout" => layout, "filename" => filename, "fileopt" => "overwrite"])
        plot_url = response["url"]
end

end


text:
using Test
import ArchGDAL; const AG = ArchGDAL
using Tables

@testset "Tables Support" begin
    dataset = AG.read(joinpath(@__DIR__, "data/point.geojson"))
    dataset1 = AG.read(joinpath(@__DIR__, "data/multi_geom.csv"), options = ["GEOM_POSSIBLE_NAMES=point,linestring", "KEEP_GEOM_COLUMNS=NO"])
    dataset2 = AG.read(joinpath(@__DIR__, "data/missing_testcase.csv"), options = ["GEOM_POSSIBLE_NAMES=point,linestring", "KEEP_GEOM_COLUMNS=NO"])
    @test dataset isa ArchGDAL.IDataset
    @test dataset1 isa ArchGDAL.IDataset
    @test dataset2 isa ArchGDAL.IDataset
    layer = AG.getlayer(dataset, 0)
    layer1 = AG.getlayer(dataset1, 0)
    layer2 = AG.getlayer(dataset2, 0)
    gt = AG.Table(layer)
    gt1 = AG.Table(layer1)
    gt2 = AG.Table(layer2)

    @testset "read layer to table" begin
        @test AG.getlayer(gt) === layer
        @test AG.getlayer(gt1) === layer1
        @test AG.getlayer(gt2) === layer2
        @test sprint(print, gt) == "Table with 4 features\n"
        @test sprint(print, gt1) == "Table with 2 features\n"
        @test sprint(print, gt2) == "Table with 9 features\n"
    end

    @testset "Tables methods" begin
        @test Tables.schema(layer).names == propertynames(gt)
        @test Tables.schema(layer1).names == propertynames(gt1)
        @test Tables.schema(layer2).names == propertynames(gt2)
        @test Tables.istable(AG.Table) == true
        @test Tables.rows(gt) == AG.Table(layer)
        @test Tables.rows(gt1) == AG.Table(layer1)
        @test Tables.rows(gt2) == AG.Table(layer2)
    end

    @testset "Misc. methods" begin
        @test Base.size(gt) == 4
        @test Base.size(gt1) == 2
        @test Base.size(gt2) == 9
        @test Base.length(gt) == 4
        @test Base.length(gt1) == 2
        @test Base.length(gt2) == 9
        @test Base.IteratorSize(typeof(gt)) == Base.HasLength()
        @test Base.IteratorEltype(typeof(gt1)) == Base.HasEltype()
        @test propertynames(gt) == (:FID, :pointname, Symbol(""))
        @test propertynames(gt1) == (:id, :zoom, :location, :point, :linestring)
        @test propertynames(gt2) == (:id, :zoom, :location, :point, :linestring)
        @test getproperty(gt, :FID) == [iterate(gt, i)[1].FID for i in 0:size(gt)-1]
        @test getproperty(gt1, :zoom) == [iterate(gt1, i)[1].zoom for i in 0:size(gt1)-1]
        @test sprint(print, gt2[5].linestring) == sprint(print, gt2[3].point)
        @test sprint(print, gt2[9].linestring) == sprint(print, gt2[7].point)
        @test collect(findall(x->x=="missing", getproperty(gt2, i)) for i in [:id, :zoom, :location]) == [[6], [4, 8], [3, 7, 8]]
        @test iterate(gt, 5) === nothing
        @test iterate(gt1, 3) === nothing
        @test typeof([getindex(gt, i) for i in 1:size(gt)]) == typeof([iterate(gt, i)[1] for i in 0:size(gt)-1])
        @test typeof([getindex(gt1, i) for i in 1:size(gt1)]) == typeof([iterate(gt1, i)[1] for i in 0:size(gt1)-1])

        AG.resetreading!(layer)
        AG.resetreading!(layer1)

        @test AG.nextnamedtuple(layer) isa NamedTuple{(:FID, :pointname, Symbol("")),Tuple{Float64,String,ArchGDAL.IGeometry}}
        @test AG.nextnamedtuple(layer1) isa NamedTuple{(:id, :zoom, :location, :point, :linestring),Tuple{String,String,String,ArchGDAL.IGeometry,ArchGDAL.IGeometry}}
        for i in 1:4
            @test AG.schema_names(layer)[i] isa Base.Generator || AG.schema_names(layer)[i] isa ArchGDAL.IFeatureDefnView
            @test AG.schema_names(layer1)[i] isa Base.Generator || AG.schema_names(layer1)[i] isa ArchGDAL.IFeatureDefnView
        end
    end
end


text:
module NiancatGames

using Niancat.Gameface
using Niancat.Languages
using Niancat.Users
using Niancat.Scores
import Niancat.Gameface: gamecommand
using UUIDs

#
# Commands
#

struct SetPuzzle <: GameCommand
    puzzle::String
end

struct Guess <: GameCommand
    word::Word

    Guess(s::String) = new(Word(s))
end

struct GetPuzzle <: GameCommand end

#
# Command responses
#

struct PuzzleIs <: Response
    puzzle::String
end

struct NoPuzzleSet <: Response end

struct Incorrect <: Response
    word::Word

    Incorrect(s::String) = new(Word(s))
end

struct Correct <: Response
    word::Word

    Correct(s::String) = new(Word(s))
end

struct Rejected <: Response end

struct NewPuzzle <: Response
    puzzle::String
end

#
# Notifications
#

struct CorrectNotification <: GameNotification
    user::User
end

#
# Game logic
#

mutable struct NiancatGame <: Game
    gameservice::GameService
    round::UUID
    dictionary::SwedishDictionary
    puzzle::Union{Nothing, String}

    NiancatGame(d::SwedishDictionary, gameservice::GameService) = new(gameservice, uuid4(), d, nothing)
    function NiancatGame(state::String, gameservice::GameService)
        # The state looks like
        #     PUZZLE#dictionary#round
        stateparts = split(state, '#')

        # The first part is the puzzle
        puzzle = stateparts[1]
        if puzzle == ""
            puzzle = nothing
        end

        # The second part is the dictionary id. We can use the id to fetch a dictionary
        # from the GameService object.
        dictionaryid = stateparts[2]
        dictionary = finddictionary(gameservice, dictionaryid)

        # The third part is the round UUID.
        round = tryparse(UUID, stateparts[3])
        if round === nothing
            round = uuid4()
        end

        new(gameservice, round, dictionary, puzzle)
    end
end

Gameface.gameround(game::NiancatGame) :: String = string(game.round)

function Gameface.gamecommand(game::NiancatGame, ::User, setpuzzle::SetPuzzle) :: Response
    isfirstround = game.puzzle === nothing
    if !isfirstround && isanagram(game.puzzle, setpuzzle.puzzle)
        return Rejected()
    end
    solutions = findanagrams(game.dictionary, setpuzzle.puzzle)
    if solutions == []
        return Rejected()
    end

    # Show a solution board to the teams, showing which user solves which word.
    # Show the scoreboard as well.
    # Only send it if there was a previous round.
    if !isfirstround
        notify!(game.gameservice, SolutionboardNotification(gameinstanceid(game.gameservice), gameround(game)))
        notify!(game.gameservice, ScoreboardNotification(gameinstanceid(game.gameservice), gameround(game)))
    end

    game.puzzle = setpuzzle.puzzle
    game.round = uuid4()

    for solution in solutions
        event!(game.gameservice, GameEvent(Gameface.GameEvent_Solution, string(game.round), solution))
    end

    NewPuzzle(game.puzzle)
end

function Gameface.gamecommand(game::NiancatGame, ::User, getpuzzle::GetPuzzle) :: Response
    if game.puzzle === nothing
        NoPuzzleSet()
    else
        PuzzleIs(game.puzzle)
    end
end

function Gameface.gamecommand(game::NiancatGame, user::User, guess::Guess) :: Response
    if guess.word.w in game.dictionary && isanagram(guess.word.w, game.puzzle)
        solution = normalizedword(game.dictionary, guess.word.w)
        score!(game.gameservice, user, Score(1, string(game.round), solution, 1))
        event!(game.gameservice, GameUserEvent(Gameface.UserEvent_Solution, user, string(game.round), solution))
        notify!(game.gameservice, CorrectNotification(user))
        Correct(guess.word.w)
    else
        Incorrect(guess.word.w)
    end
end

include("NiancatFormatting.jl")

export NiancatGame

end


text:
# based on https://github.com/chkwon/TrafficAssignment.jl

include(joinpath("Julia_files/load_network_uni-class.jl"))

using Graphs

function create_graph(start_node, end_node)
    @assert Base.length(start_node)==Base.length(end_node)

    no_node = max(maximum(start_node), maximum(end_node))
    no_arc = Base.length(start_node)

    graph = simple_inclist(no_node)
    for i=1:no_arc
        add_edge!(graph, start_node[i], end_node[i])
    end
    return graph
end

function get_vector(state, origin, destination, link_dic)
    current = destination
    parent = -1
    x = zeros(Int, maximum(link_dic))

    while parent != origin
        parent = state.parents[current]

        link_idx = link_dic[parent,current]

        if link_idx != 0
            x[link_idx] = 1
        end

        current = parent
    end

    return x
end

#root_dir = "G:/My Drive/Github/PoA/Price_of_Anarchy_for_Transportation_Networks"
out_dir = "../results/"
files_ID ="_cdc_all_comp_apr_2012"
month_w = "Apr"
instances_ = ["AM", "MD", "PM", "NT"]

for instance in instances_

    ta_data = load_ta_network_(out_dir, files_ID, month_w, instance)

    # unpacking data from ta_data
    network_name = ta_data.network_name

    number_of_zones = ta_data.number_of_zones
    number_of_nodes = ta_data.number_of_nodes
    first_thru_node = ta_data.first_thru_node
    number_of_links = ta_data.number_of_links

    start_node = ta_data.start_node
    end_node = ta_data.end_node
    capacity = ta_data.capacity
    link_length = ta_data.link_length

    free_flow_time = ta_data.free_flow_time
    B = ta_data.B
    power = ta_data.power
    speed_limit = ta_data.speed_limit
    toll = ta_data.toll
    link_type = ta_data.link_type
    number_of_zones = ta_data.number_of_zones
    total_od_flow = ta_data.total_od_flow
    travel_demand = ta_data.travel_demand
    od_pairs = ta_data.od_pairs

    toll_factor = ta_data.toll_factor
    distance_factor = ta_data.distance_factor

    best_objective = ta_data.best_objective

    # preparing a graph
    graph = create_graph(start_node, end_node)
    link_dic = sparse(start_node, end_node, 1:number_of_links);

    function MSA(coeffs)
        polyEval(coeffs, pt) = sum([coeffs[i] * pt^(i-1) for i = 1:length(coeffs)])

        function BPR(x)
            bpr = similar(x)
            for i=1:length(bpr)
                bpr[i] = free_flow_time[i] * polyEval( coeffs, (x[i]/capacity[i]) )
            end
            return bpr
        end

        function all_or_nothing(travel_time)
            state = []
            path = []
            x = zeros(size(start_node))

            for r=1:size(travel_demand)[1]
                # for each origin node r, find shortest paths to all destination nodes
                state = dijkstra_shortest_paths(graph, travel_time, r)

                for s=1:size(travel_demand)[2]
                    # for each destination node s, find the shortest-path vector
                    # load travel demand
                    x = x + travel_demand[r,s] * get_vector(state, r, s, link_dic)
                end
            end

            return x
        end

        # Finding a starting feasible solution
        travel_time = BPR(zeros(number_of_links))
        xl = all_or_nothing(travel_time)

        max_iter_no = 1e3
        l = 1
        #average_excess_cost = 1
        tol = 1e-5

        while l < max_iter_no
            l += 1

            xl_old = xl

            # Finding yl
            travel_time = BPR(xl)

            yl = all_or_nothing(travel_time)

            xl = xl + (yl - xl)/l

            xl_new = xl

            relative_gap = norm(xl_new - xl_old, 1) / norm(xl_new, 1)

            if relative_gap < tol
                break
            end
        end

        return xl
    end

    # getting the coefficients of the costs
    coeffs_dict_Apr_AM = readall(out_dir * files_ID * "/coeffs_dict_" * month_w * "_" * instance * ".json")
    coeffs_dict_Apr_AM = JSON.parse(out_dir * files_ID * "/coeffs_dict_" * month_w * "_" * instance * ".json")

    xl = Dict{(Int64,Float64,Float64,Int64),Array{Float64,1}}()

    deg_grid = [4:8]
    c_grid = .5:.5:3.
    lamb_grid = [10. .^(-3:4)]

    for deg in deg_grid
        for c in c_grid
            for lam in lamb_grid
                coeffs_1 = coeffs_dict_Apr_AM["($(deg),$(c),$(lam),1)"]
                coeffs_2 = coeffs_dict_Apr_AM["($(deg),$(c),$(lam),2)"]
                coeffs_3 = coeffs_dict_Apr_AM["($(deg),$(c),$(lam),3)"]
                xl[(deg, c, lam, 1)] = MSA(coeffs_1)
                xl[(deg, c, lam, 2)] = MSA(coeffs_2)
                xl[(deg, c, lam, 3)] = MSA(coeffs_3)
            end
        end
    end

    outfile = open(out_dir * files_ID * "uni-class_traffic_assignment_MSA_flows_" * month_w * "_" * instance * ".json", "w")

    JSON.print(outfile, xl)

    close(outfile)
end


text:
# VARIOUS TYPES FOR CALCULATIONS

type LapStatus
    currentLap::Int64       # current lap number
    currentIt::Int64        # current iteration in current lap
    switchLap::Bool
    nextLap::Bool
    s_lapTrigger::Float64
end

# Structure of coeffConst:
# 1st dimension is the polynomial coefficient
# 2nd dimension specifies the state (1 = eY, 2 = ePsi, 3 = v or 1 = xDot, 2 = yDot, 3 = psiDot, 4 = ePsi, 5 = eY)
# 3th dimension specifies one of the two lap numbers between which are iterated

type MpcCoeff           # coefficients for trajectory approximation
    coeffCost::Array{Float64}
    coeffConst::Array{Float64}
    order::Int64
    pLength::Int64      # small values here may lead to numerical problems since the functions are only approximated in a short horizon
                        # "small" values are about 2*N, good values about 4*N
                        # numerical problems occur at the edges (s=0, when v is almost 0 and s does not change fast and at s=s_target)
    c_Vx::Array{Float64,1}
    c_Vy::Array{Float64,1}
    c_Psi::Array{Float64,1}
    MpcCoeff(coeffCost=Float64[], coeffConst=Float64[], order=4, pLength=0,c_Vx=Float64[],c_Vy=Float64[],c_Psi=Float64[]) = new(coeffCost, coeffConst, order, pLength, c_Vx, c_Vy, c_Psi)
end

type OldTrajectory      # information about previous trajectories
    oldTraj::Array{Float64}             # contains all states over all laps
    oldInput::Array{Float64}            # contains all inputs
    oldTimes::Array{Float64}            # contains times related to states and inputs
    oldCost::Array{Int64}               # contains costs of laps
    count::Array{Int64}                 # contains the counter for each lap
    prebuf::Int64
    postbuf::Int64
    idx_start::Array{Int64}             # index of the first measurement with s > 0
    idx_end::Array{Int64}               # index of the last measurement with s < s_target
    OldTrajectory(oldTraj=Float64[],oldInput=Float64[],oldTimes=Float64[],oldCost=Float64[],count=Int64[],prebuf=50,postbuf=50,idx_start=Int64[],idx_end=Int64[]) = new(oldTraj,oldInput,oldTimes,oldCost,count,prebuf,postbuf,idx_start,idx_end)
end

type SafeSetData
    oldSS::Array{Float64}           # contains data from previous laps usefull to build the safe set
    cost2target::Array{Float64}     # cost to arrive at the target, i.e. how many iterations from the start to the end of the lap
    oldCost::Array{Int64}               # contains costs of laps
    count::Array{Int64}                 # contains the counter for each lap
    prebuff::Int64
    postbuff::Int64
    idx_start::Array{Int64}             # index of the first measurement with s > 0
    idx_end::Array{Int64}               # index of the last measurement with s < s_target
    oldSS_xy::Array{Float64}

    SafeSetData(oldSS=Float64[],cost2target=Float64[],oldCost=Int64[],count=Int64[],prebuf=50,postbuf=50,idx_start=Int64[],idx_end=Int64[],oldSS_xy=Float64[]) =
    new(oldSS,cost2target,oldCost,count,prebuf,postbuf,idx_start,idx_end,oldSS_xy)
end

type SelectedStates                 # Values needed for the convex hull formulation
    selStates::Array{Float64}       # selected states from previous laps ...
    statesCost::Array{Float64}      # ... and their related costs
    Np::Int64                       # number of states to select from each previous lap
    Nl::Int64                       # number of previous laps to include in the convex hull
    version::Bool                   # false if you want to use convex hull
    simulator::Bool                 # true to use tuning made for simulator,  false to use tuning made for BARC
    shift::Int64
    SelectedStates(selStates=Float64[],statesCost=Float64[],Np=6,Nl=2,version=false, simulator=true,shift=3) = new(selStates,statesCost,Np,Nl,version,simulator,shift)
end

type MpcParams          # parameters for MPC solver
    N::Int64
    nz::Int64
    OrderCostCons::Int64
    Q::Array{Float64,1}
    Q_term::Array{Float64,1}
    R::Array{Float64,1}
    vPathFollowing::Float64
    QderivZ::Array{Float64,1}
    QderivU::Array{Float64,1}
    Q_term_cost::Float64
    delay_df::Int64
    delay_a::Int64
    Q_lane::Float64                # weight on the soft constraint for the lane
    Q_vel::Float64                 # weight on the soft constraint for the maximum velocity
    Q_slack::Array{Float64,1}               # weight on the slack variables for the terminal constraint
    Q_obs::Array{Float64}          # weight used to esclude some of the old trajectories from the optimization problem

    MpcParams(N=0,nz=0,OrderCostCons=0,Q=Float64[],Q_term=Float64[],R=Float64[],vPathFollowing=1.0,QderivZ=Float64[],QderivU=Float64[],Q_term_cost=1.0,delay_df=0,delay_a=0,Q_lane=1.0,Q_vel=1.0,Q_slack=Float64[],Q_obs=Float64[]) = new(N,nz,OrderCostCons,Q,Q_term,R,vPathFollowing,QderivZ,QderivU,Q_term_cost,delay_df,delay_a,Q_lane,Q_vel,Q_slack,Q_obs)
end

type PosInfo            # current position information
    s::Float64
    s_target::Float64
    PosInfo(s=0,s_target=0) = new(s,s_target)
end

type MpcSol             # MPC solution output
    a_x::Float64
    d_f::Float64
    solverStatus::Symbol
    u::Array{Float64}
    z::Array{Float64}
    cost::Array{Float64}
    eps_alpha::Array{Float64}
    costSlack::Array{Float64}
    MpcSol(a_x=0.0,d_f=0.0,solverStatus=Symbol(),u=Float64[],z=Float64[],cost=Float64[],eps_alpha=Float64[],costSlack=Float64[]) = new(a_x,d_f,solverStatus,u,z,cost,eps_alpha,costSlack)
end

type TrackCoeff         # coefficients of track
    coeffAngle::Array{Float64,1}
    coeffCurvature::Array{Float64,1}
    nPolyCurvature::Int64      # order of the interpolation polynom
    width::Float64               # lane width -> is used in cost function as soft constraints (to stay on track)
    TrackCoeff(coeffAngle=Float64[],coeffCurvature=Float64[],nPolyCurvature=4,width=1.0) = new(coeffAngle,coeffCurvature,nPolyCurvature,width)
end

type ModelParams
    l_A::Float64
    l_B::Float64
    m::Float64
    I_z::Float64
    dt::Float64
    u_lb::Array{Float64}        # lower bounds for u
    u_ub::Array{Float64}        # upper bounds
    z_lb::Array{Float64}
    z_ub::Array{Float64}
    c0::Array{Float64}
    c_f::Float64
    ModelParams(l_A=0.25,l_B=0.25,m=1.98,I_z=0.24,dt=0.1,u_lb=Float64[],u_ub=Float64[],z_lb=Float64[],z_ub=Float64[],c0=Float64[],c_f=0.0) = new(l_A,l_B,m,I_z,dt,u_lb,u_ub,z_lb,z_ub,c0,c_f)
end

type Obstacle
    obstacle_active::Bool       # true if we have to consider the obstacles in the optimization problem
    lap_active::Int64           # number of the first lap in which the obstacles are used
    obs_detect::Float64         # maximum distance at which we can detect obstacles (in terms of s!!)
    n_obs::Int64                # number of obstacles in the track
    s_obs_init::Array{Float64}  # initial s coordinate of each obstacle
    ey_obs_init::Array{Float64} # initial ey coordinate of each obstacle
    v_obs_init::Array{Float64}  # initial velocity of each obstacles
    r_s::Float64                # radius on the s coordinate of the ellipse describing the obstacles
    r_ey::Float64               # radius on the ey coordinate of the ellipse describing the obstacle
    inv_step::Int64             # number of step of invariance required for the safe set

    Obstacle(obstacle_active=false,lap_active=10,obs_detect=1.0,n_obs=1,s_obs_init=Float64[],ey_obs_init=Float64[],v_obs_init=Float64[],r_s=0.5,r_ey=0.3,inv_step=1) = new(obstacle_active,lap_active,obs_detect,n_obs,s_obs_init,ey_obs_init,v_obs_init,r_s,r_ey,inv_step)
end

max_stars_repo_name: RoyaFiroozi/barc

max_issues_repo_name: RoyaFiroozi/barc

max_forks_repo_name: RoyaFiroozi/barc


text:
#start boundary conditions
println("calculating right-hand-side...\n")

F = zeros(dofs,1)

(nbc,ni,ix) = IDs(bc) #nbc has number of conditions
                   #ix has edges nodes ids for each condition
global bc,edge,elem,nodes,i,F
i=1

while true
    global i,ix,elem
    if i>nbc
        break
    end
    rows = [j for k=1:length(ni[i]) for j=1:length(elem[:,1]) if issubset(ix[i][k,1],elem[j,:]) &
                                                                 issubset(ix[i][k,2],elem[j,:]) &
                                                                 issubset(ix[i][k,3],elem[j,:])] #contains superset id element
    xedge = nodes[elem[rows,:],:][:,:,1]; yedge = nodes[elem[rows,:],:][:,:,2] #contains x,y of edges

    superelem = elem[rows,:]
    mapx = [sum(dims=2,superelem[k,:] .== ix[i][k,:]') for k=1:length(superelem[:,1])]#zero the non edge indexes of element
    for k=1:length(mapx)
        g(x) = f(x)[i]
        a = BC(g,xedge[k,:]',yedge[k,:]',mapx[k])
        id = [ix[i][k,:]' (ix[i][k,:]' .+ Int64(0.5*dofs))]
        id = sort(id,dims=2)

        idlocal = [mapx[k]' mapx[k]'] #has the global degrees of freedoms
        idlocal = [idlocal[i] .==1 for i=1:12]
        F[id] .= F[id] .+ a[idlocal]'
    end
    i = i+1
end

max_stars_repo_name: billyon/first

max_issues_repo_name: billyon/first

max_forks_repo_name: billyon/first


text:
# Load in, combine, and plot chains:
using DelimitedFiles
using JLD2
using Statistics
using Printf
using MCMCDiagnostics
using PyPlot
include("../../../src/nlog_prior.jl")
if !@isdefined(CGS)
  include("../../../src/CGS.jl")
  using Main.CGS
end

n = 8
nplanet = n-1
nparam = (n-1)*5+2
nstep = 2000

#fnames = readdlm("jld2_list.txt")
#fnames = readdlm("jld2_list_eps0.1.txt")
fnames = readdlm("../../../data/jld2_list_nstep2000_eps0.1.txt")

nfile = size(fnames)[1]
#ncurr = 250
#ncurr = 350
ncurr = 2000
state_total = zeros(nparam*2+2,ncurr*nfile)

nacc_all = zeros(Int64,nfile)
neff_tot = zeros(nfile,37)
for i=1:nfile
  i1 = (i-1)*ncurr+1
  i2 = i*ncurr
  @load string("../../../data/output_files/",fnames[i]) state nacc
  nacc_all[i] = nacc
  state_total[:,i1:i2] = state[:,1:ncurr]
  for j=1:37
    neff_tot[i,j] += effective_sample_size(state[j,1:ncurr])
  end
end
for j=1:37
  println(j," ",sum(neff_tot[:,j]))
end

clf()
plot(nacc_all)
#read(stdin,Char)

@load "../../../data/output_files/T1_run_hmc_student_ecc_lgndof_V1exp2nuinv_nstep2000_eps0.1_nleap20_201.jld2"

x_opt = [4.598791344049799e-5, 1.5108213441190745, 7257.550484415569, -0.00495015358026971, 0.004671261296998595,
         4.40588192304326e-5, 2.4219505838480635, 7258.587230337616, -0.001188890950465861, 0.0014051713506332184,
         1.305145694403463e-5, 4.049193371144702, 7257.06806273418, -0.0064708229088599216, 0.003529480254615964,
         2.3551037452505444e-5, 6.101020579519802, 7257.827698633612, 0.003217590595460454, -0.003935834596847923,
         3.4925310671129774e-5, 9.207561405004297, 7257.073658382135, -0.009581969253521241, -0.00016501113996023978,
         4.417497904782645e-5, 12.352405680378052, 7257.715102776943, 0.00288626733667452, 0.0016434273040851228,
         1.0931051156428014e-5, 18.772977044007867, 7249.605686691679, -0.004439910317349172, 0.0002702876058394804,
         1.2119892906238703, 0.8222963482043892]

planet = ["b","c","d","e","f","g","h"]
param = ["mu","P","t0","ecos","esin"]
prior = zeros(nfile*ncurr)
for j=1:nfile*ncurr
  prior[j] = nlog_prior(state_total[1:37,j])
end

# compute f0:
#f0total = zeros(nfile*ncurr)
#for j=1:nfile*ncurr
#  f0total[j] = f0(state_total[1:37,j])
#end

for iparam=1:nparam
  clf()
  cparam = state_total[iparam,:]
#  plot(cparam[cparam .== circshift(cparam,100)],".")
  plot(cparam,".")
  plot([0,ncurr*nfile],x_opt[iparam] .+ [0,0])
  plot([0,ncurr*nfile],x_opt[iparam] + sqrt(cov_save[iparam,iparam]) .+ [0,0],linestyle="--")
  plot([0,ncurr*nfile],x_opt[iparam] - sqrt(cov_save[iparam,iparam]) .+ [0,0],linestyle="--")
#  axis([0,700,minimum(cparam),maximum(cparam)])
  if iparam < nparam-1
    pname = string(planet[div(iparam-1,5)+1]," ",param[mod(iparam-1,5)+1])
  elseif iparam == nparam-1
    pname = "log(nu)"
  else
    pname = "V1 e^{1/(2nu)}"
  end
  println(iparam," ",pname," chains: ",@sprintf("%6.3g",mean(state_total[iparam,:])),"+-",@sprintf("%6.3g",std(state_total[iparam,:]))," Fisher: ",@sprintf("%6.3g",x_opt[iparam]),"+-",@sprintf("%6.3g",sqrt(cov_save[iparam,iparam]))," ratio: ",@sprintf("%6.3g",std(state_total[iparam,:])/sqrt(cov_save[iparam,iparam])))
#  read(stdin,Char)
end

clf()
# Probability distribution of e*sin(omega) vs. e*cos(omega)
plot(state_total[4,:],state_total[5,:],".",alpha=0.02)
plot(state_total[9,:],state_total[10,:],".",alpha=0.02)
plot([0,0],[0,0],".")
xlabel(L"$e_b \cos{\omega_b}$")
ylabel(L"$e_b \sin{\omega_b}$")
#read(stdin,Char)

include("../../../src/histogram_code.jl")

# plot histogram of each

esinb_bin, esinb_hist, esinb_bin_square,esinb_hist_square = histogram(state_total[4,:],50)
ecosb_bin, ecosb_hist, ecosb_bin_square,ecosb_hist_square = histogram(state_total[5,:],50)
esinc_bin, esinc_hist, esinc_bin_square,esinc_hist_square = histogram(state_total[9,:],50)
ecosc_bin, ecosc_hist, ecosc_bin_square,ecosc_hist_square = histogram(state_total[10,:],50)

clf();
plot(esinb_bin_square,esinb_hist_square ./ maximum(esinb_hist_square),label=L"$e_b \sin{\omega_b}$")
plot(ecosb_bin_square,ecosb_hist_square ./ maximum(ecosb_hist_square),label=L"$e_b \cos{\omega_b}$")
plot(esinc_bin_square,esinc_hist_square ./ maximum(esinc_hist_square),label=L"$e_c \sin{\omega_c}$")
plot(ecosc_bin_square,ecosc_hist_square ./ maximum(ecosc_hist_square),label=L"$e_c \cos{\omega_c}$")
axis([-0.01,0.01,0,1.05])
legend()
#read(stdin,Char)

# Now plot histogram of longitudes of inner two planets:

#delta_omega = atan.(state_total[10,:],state_total[9,:]) .- atan.(state_total[5,:],state_total[4,:]) .-pi
ecosc = state_total[9,:]
esinc = state_total[10,:]
ecosb = state_total[4,:]
esinb = state_total[5,:]
delta_omega = atan.(ecosc .*esinb .-ecosb .*esinc, -ecosb .*ecosc .-esinb .*esinc)

delta_omega[delta_omega .< -pi] .+= 2\pi
delta_omega[delta_omega .< -pi] .+= 2\pi
delta_omega[delta_omega .>  pi] .-= 2\pi
delta_omega[delta_omega .>  pi] .-= 2\pi

domega_bin, domega_hist, domega_bin_square, domega_hist_square = histogram(delta_omega,50)

#domega_hist_square[domega_hist_square .== 0.0] .+= 200.0
clf(); plot(domega_bin_square*180/pi,domega_hist_square./maximum(domega_hist_square),linewidth=2)

xlabel(L"$\omega_c - \omega_b - 180^\circ$",fontsize=20)
ylabel("Probability",fontsize=20)
axis([-180,180,0,1.05])
tight_layout()
savefig("../delta_omega_bc.pdf",bbox_inches="tight")
# Now print mass ratios:
xfit = zeros(2,nplanet)
fac = 0.09*MSUN/MEARTH
for i=1:nplanet
  xfit[:,i] = [mean(state_total[(i-1)*5+1,:])*fac,std(state_total[(i-1)*5+1,:])*fac]
  println("planet: ",i," mass: ",@sprintf("%6.4f",xfit[1,i]),"+-",@sprintf("%6.4f",xfit[2,i]))
end

#@save "T1_hmc_total_02212020.jld2" state_total


text:
using MatplotWrap

x = -pi:0.01:pi
y = sin.(x)

plot(x, y, "-r") # matplot::plot(x, y)
matplot_show() # matplot::show()


text:
##########################
# Geomagnetic Index Data #
##########################

export GeomagneticIndexData
"""
GeomagneticIndexData stores geomagnetic index data stores the geomagnetic
activity Kp index and the dervied Ap indices.

# Attributes
- `data::Dict`
"""
struct GeomagneticIndexData
    data::Dict{Int, Tuple{Array{Float64, 1}, Float64, Array{Float64, 1}, Float64}}
end

function get_kp_decimal(number, digit)
    if isempty(strip(string(number)))
        number = 0.0
    else
        number = parse(Int, number)
    end
    decimal = 0.0
    if parse(Int, digit) == 3
        decimal = 1.0/3.0
    elseif parse(Int, digit) == 7
        decimal = 2.0/3.0
    end

    return number + decimal
end

function GeomagneticIndexData(filepath::String)
    data = Dict{Int, Tuple{Array{Float64, 1}, Float64, Array{Float64, 1}, Float64}}()

    # Load file and read in data
    open(filepath, "r") do fp
        for line in readlines(fp)
            if !isempty(line)
                # Day of observation
                year  = 2000 + parse(Int, line[1:2])
                month = parse(Int, line[3:4])
                day   = parse(Int, line[5:6])

                # Data
                kp_data = Float64[]

                for i in 13:2:27
                    push!(kp_data, get_kp_decimal(line[i], line[i+1]))
                end

                kp_daily = get_kp_decimal(line[29:30], line[31])

                # Data
                ap_data = Float64[]

                for i in 32:3:55
                    push!(ap_data, parse(Float64, line[i:(i+2)]))
                end

                ap_daily = parse(Float64, line[56:58])

                mjd = round(Int, caldate_to_mjd(year, month, day))
                data[mjd] = (kp_data, kp_daily, ap_data, ap_daily)
            end
        end
    end

    GeomagneticIndexData(data)
end

export GEOMAGNETIC_DATA
"""
Module-wide global GeomagneticIndexData. This data object is used as the
default source of geomagnetic data by dynamics models if no explicit
GeomagneticIndexData file is provided to those transformations.

This value can be overridden in your own code as follows:

```julia
SatelliteDynamics.SpaceWeather.GEOMAGNETIC_DATA = GeomagneticIndexData(index_file)
```

This global variable defaults used are pull from GFZ Potsdam FTP server:
ftp://ftp.gfz-potsdam.de/pub/home/obs/kp-ap/wdc
"""
global GEOMAGNETIC_DATA = GeomagneticIndexData(abspath(@__DIR__, "../../data/kpall.wdc"))

export KpIndices
"""
Retrieve Kp Indices for day in question

Arguments:
- `mjd::Real` Modified Julian date of desired data. Time system of input is UT1
- `epc::Epoch` Epoch of desired input Time system of input is UT1
"""
function KpIndices(geoindices::GeomagneticIndexData, mjd::Real)
    return geoindices.data[floor(Int, mjd)][1]
end

KpIndices(mjd::Real)  = KpIndices(GEOMAGNETIC_DATA, mjd)
KpIndices(epc::Epoch) = KpIndices(GEOMAGNETIC_DATA, mjd(epc, tsys="UT1"))

export KpIndex
"""
Kp index for the given epoch.

Arguments:
- `mjd::Real` Modified Julian date of desired data. Time system of input is UT1
- `epc::Epoch` Epoch of desired input Time system of input is UT1
"""
function KpIndex(geoindices::GeomagneticIndexData, mjd::Real)
    # Base UT1 day for index
    mjd_ut1  = floor(Int, mjd)

    # Get Index into vector based on 3-hour window
    hour_idx = floor(Int, (mjd - mjd_ut1)*8)
    return geoindices.data[mjd_ut1][1][hour_idx+1]
end

KpIndex(mjd::Real)  = KpIndex(GEOMAGNETIC_DATA, mjd)
KpIndex(epc::Epoch) = KpIndex(GEOMAGNETIC_DATA, mjd(epc, tsys="UT1"))

export KpDailyIndex
"""
Retrieve daily Kp index for the day in question.

Arguments:
- `mjd::Real` Modified Julian date of desired data. Time system of input is UT1
- `epc::Epoch` Epoch of desired input Time system of input is UT1
"""
function KpDailyIndex(geoindices::GeomagneticIndexData, mjd::Real)
    return geoindices.data[floor(Int, mjd)][2]
end

KpDailyIndex(mjd::Real)  = KpDailyIndex(GEOMAGNETIC_DATA, mjd)
KpDailyIndex(epc::Epoch) = KpDailyIndex(GEOMAGNETIC_DATA, mjd(epc, tsys="UT1"))

export ApIndices
"""
Retrieve Ap Indices for day in question

Arguments:
- `mjd::Real` Modified Julian date of desired data. Time system of input is UT1
- `epc::Epoch` Epoch of desired input Time system of input is UT1
"""
function ApIndices(geoindices::GeomagneticIndexData, mjd::Real)
    return geoindices.data[floor(Int, mjd)][3]
end

ApIndices(mjd::Real)  = ApIndices(GEOMAGNETIC_DATA, mjd)
ApIndices(epc::Epoch) = ApIndices(GEOMAGNETIC_DATA, mjd(epc, tsys="UT1"))

export ApIndex
"""
Ap index for the given epoch.

Arguments:
- `mjd::Real` Modified Julian date of desired data. Time system of input is UT1
- `epc::Epoch` Epoch of desired input Time system of input is UT1
"""
function ApIndex(geoindices::GeomagneticIndexData, mjd::Real)
    # Base UT1 day for index
    mjd_ut1  = floor(Int, mjd)

    # Get Index into vector based on 3-hour window
    hour_idx = floor(Int, (mjd - mjd_ut1)*8)
    return geoindices.data[mjd_ut1][3][hour_idx+1]
end

ApIndex(mjd::Real)  = ApIndex(GEOMAGNETIC_DATA, mjd)
ApIndex(epc::Epoch) = ApIndex(GEOMAGNETIC_DATA, mjd(epc, tsys="UT1"))

export ApDailyIndex
"""
Retrieve daily Ap index for the day in question.

Arguments:
- `mjd::Real` Modified Julian date of desired data. Time system of input is UT1
- `epc::Epoch` Epoch of desired input Time system of input is UT1
"""
function ApDailyIndex(geoindices::GeomagneticIndexData, mjd::Real)
    return geoindices.data[floor(Int, mjd)][4]
end

ApDailyIndex(mjd::Real)  = ApDailyIndex(GEOMAGNETIC_DATA, mjd)
ApDailyIndex(epc::Epoch) = ApDailyIndex(GEOMAGNETIC_DATA, mjd(epc, tsys="UT1"))

###################
# Solar Flux Data #
###################

export SolarFluxData
"""
SolarFluxData stores solar flux index data.

The data is a dictionary keyed to the modified julian date in UT of the day of
observation. All observations are taken at 20:00 UT at the of question. The values
are a tuple storing the solar radio flux the F10.7cm solar radio flux. The first
entry is the flux observed directly at the Ottowa radio observatory station,
while the second entry is the flux value adjusted to 1 AU. The final two entries
is the 81-day average, cenerted on the day entry, of the observed and adjusted
solar radio flux.

# Attributes
- `data::Dict{Int, Tuple{Float64, Float64, Float64}` Stores observed F10.7cm
solar radio flux, the adjusted solar radio flux, and the 81-day average centered
on the day in question.
"""
struct SolarFluxData
    data::Dict{Int, Tuple{Float64, Float64, Float64, Float64}}
end

function SolarFluxData(filepath::String)
    data = Dict{Int, Tuple{Float64, Float64, Float64, Float64}}()

    # Used to calculate the average solar radio flux for the days in question
    avg  = []

    # Load file and read in data
    open(filepath, "r") do fp
        for i in 1:2
            # Skip first two lines
            readline(fp)
        end

        for line in readlines(fp)
            split_line = split(line)

            # Day of observation
            year  = parse(Int, split_line[1][1:4])
            month = parse(Int, split_line[1][5:6])
            day   = parse(Int, split_line[1][7:8])

            # Data
            f107_obs = parse(Float64, split_line[5])
            f107_adj = parse(Float64, split_line[6])

            # Push data to array for computing 81 day average
            push!(avg, [year month day f107_obs f107_adj])
        end
    end

    # Format data in matrix to make it possible to compute 81-day average
    avg = vcat(avg...)

    day_avg = 81            # Number of days in average
    day_os  = round(Int, (day_avg-1)/2) # One sided average direction
    n_data  = size(avg)[1]
    for i in 1:n_data
        min_index = max(1, i-day_os)
        max_index = min(n_data, i+day_os)

        n_avg   = 0 # Number of data points in average
        obs_avg = 0 # Observed flux average
        adj_avg = 0 # Adjusted flux average
        for idx in min_index:max_index
            n_avg += 1

            obs_avg += avg[idx, 4]
            adj_avg += avg[idx, 5]
        end

        # Calculate data average
        obs_avg = obs_avg/n_avg
        adj_avg = adj_avg/n_avg

        # Insert entry into dictionary
        year      = round(Int, avg[i, 1])
        month     = round(Int, avg[i, 2])
        day       = round(Int, avg[i, 3])
        mjd       = round(Int, caldate_to_mjd(year, month, day))
        data[mjd] = (avg[i, 4], avg[i, 5], obs_avg, adj_avg)
    end

    SolarFluxData(data)
end

export SOLAR_FLUX_DATA
"""
Module-wide global SolarFluxData object. This data object is used as the
default source of Solar Flux Data by dynamics models if no explicit
SolarFluxData file is provided to those transformations.

This value can be overridden in your own code as follows:

```julia
SatelliteDynamics.SpaceWeather.SOLAR_FLUX_DATA = SolarFluxData(flux_file)
```

This global variable defaults uses flux data provide by National Resources Canada:
ftp://ftp.seismo.nrcan.gc.ca/spaceweather/solar_flux/daily_flux_values/fluxtable.txt
"""
global SOLAR_FLUX_DATA = SolarFluxData(abspath(@__DIR__, "../../data/fluxtable.txt"))

export f107Data
"""
F10.7 cm solar flux data on the day in question.

Arguments:
- `mjd::Real` Modified Julian date of desired data. Time system of input is UT1
- `epc::Epoch` Epoch of desired input Time system of input is UT1

Returns:
- `data::Tuple{Float64, Float64, Float64, Float64}` Flux data on the day in question
Elements are observed, adjusted, 81-observed average, 81-adjusted average.
"""
function f107Data(solarflux::SolarFluxData, mjd::Real)
    return solarflux.data[floor(Int, mjd)]
end

f107Data(mjd::Real)  = f107Data(SOLAR_FLUX_DATA, mjd)
f107Data(epc::Epoch) = f107Data(SOLAR_FLUX_DATA, mjd(epc, tsys="UT1"))

export f107Observed
"""
Retrieve the F10.7 cm solar flux data observed at the day in question.

Arguments:
- `mjd::Real` Modified Julian date of desired data. Time system of input is UT1
- `epc::Epoch` Epoch of desired input Time system of input is UT1
"""
function f107Observed(solarflux::SolarFluxData, mjd::Real)
    return solarflux.data[floor(Int, mjd)][1]
end

f107Observed(mjd::Real)  = f107Observed(SOLAR_FLUX_DATA, mjd)
f107Observed(epc::Epoch) = f107Observed(SOLAR_FLUX_DATA, mjd(epc, tsys="UT1"))

export f107Adjusted
"""
Retrieve the F10.7 cm solar flux data adjusted to 1 AU distance on the day in
question.

Arguments:
- `mjd::Real` Modified Julian date of desired data. Time system of input is UT1
- `epc::Epoch` Epoch of desired input Time system of input is UT1
"""
function f107Adjusted(solarflux::SolarFluxData, mjd::Real)
    return solarflux.data[floor(Int, mjd)][2]
end

f107Adjusted(mjd::Real)  = f107Adjusted(SOLAR_FLUX_DATA, mjd)
f107Adjusted(epc::Epoch) = f107Adjusted(SOLAR_FLUX_DATA, mjd(epc, tsys="UT1"))

export f107ObservedAvg
"""
Retrieve the 81-day average of the F10.7 cm solar flux data observed at the
observatory centered on the day in question.

Arguments:
- `mjd::Real` Modified Julian date of desired data. Time system of input is UT1
- `epc::Epoch` Epoch of desired input Time system of input is UT1
"""
function f107ObservedAvg(solarflux::SolarFluxData, mjd::Real)
    return solarflux.data[floor(Int, mjd)][3]
end

f107ObservedAvg(mjd::Real)  = f107ObservedAvg(SOLAR_FLUX_DATA, mjd)
f107ObservedAvg(epc::Epoch) = f107ObservedAvg(SOLAR_FLUX_DATA, mjd(epc, tsys="UT1"))

export f107AdjustedAvg
"""
Retrieve the 81-day average of the F10.7 cm solar flux data adjusted to 1 AU
centered on the day in question.

Arguments:
- `mjd::Real` Modified Julian date of desired data. Time system of input is UT1
- `epc::Epoch` Epoch of desired input Time system of input is UT1
"""
function f107AdjustedAvg(solarflux::SolarFluxData, mjd::Real)
    return solarflux.data[floor(Int, mjd)][4]
end

f107AdjustedAvg(mjd::Real)  = f107AdjustedAvg(SOLAR_FLUX_DATA, mjd)
f107AdjustedAvg(epc::Epoch) = f107AdjustedAvg(SOLAR_FLUX_DATA, mjd(epc, tsys="UT1"))


text:
"""
    loc_max(arr::Array{Float64,1})

Shorcut to Julia's findmax.
"""
function loc_max(arr::Array{Float64,1})
    return findmax(arr)[2]
end

"""
    function r5(arg_::float64)

Shorcut to Julia's round with 5 digits.
"""
function r5(arg_::Float64)
    return round(arg_;digits=5)
end

"""
    Complex_NFT(x_arr::Array{Float64,1},y_arr::Array{Float64,1},t_::Float64)

Fastest implementation of type III NUDFT that I could think of right now.

Calculates the square of the magnitude of the fourier transform of data secuence (x_arr,y_arr) at PERIOD t_

I have disabled all possible security checks or assumptions on arrays, namely:

- function is inlined
- loops can be performed in arbitrary order, so every cicle is independent of the others.
- never waits on a previous iteration to make forward progress
- Does not check for memory dependences on iterations
- only unit increments on loop index
- non nested loops

And have the following numeric assumptions
- no nan values
- no inf values
- no signed zeros
- no division, using reciprocal multiplication  (a/b -> a*1/b)
- fused operations allowed
- common math functions are approximated when convenient (sin, exp, log, etc)
- re-association of associative operations is allowed (!)

Benchmark results with same data (300 random numbers each) on 1000 random periods, 100 evals/period:

    memory estimate:  0 bytes
    allocs estimate:  0
    --------------
    minimum time:     9.149  μs (0.00% GC)
    median time:      9.389  μs (0.00% GC)
    mean time:        10.179 μs (0.00% GC)
    maximum time:     36.571 μs (0.00% GC)

and mapping over an array of 1000 random periods, 5 runs, 100 evals/run:

    memory estimate:  7.95 KiB
    allocs estimate:  2
    --------------
    minimum time:     10.768 ms (0.00% GC)
    median time:      11.237 ms (0.00% GC)
    mean time:        11.173 ms (0.00% GC)
    maximum time:     11.529 ms (0.00% GC)

The benchmraks does count random number generation time.

"""
@inline @fastmath function Complex_NFT(x_arr::Array{Float64,1},y_arr::Array{Float64,1},t_::Float64)
    s::ComplexF64=0.0+0.0im
    @simd ivdep for i=1:length(x_arr)
        @inbounds s+=y_arr[i]*exp(-6.283185307179586im*x_arr[i]/t_)
    end
    return abs2(s)
end

"""
    find_period_ephemeris(X_::Array{Float64,1},Y_::Array{Float64,1},T_::{Float64:Float64:Float64})

Attempt to find a period and ephemeris of the signal data (X_,Y_) on the Julia range T_ via maximum of fourier transform modulus.

The ephemeris is found with the formula: period*(arg(F)/2pi+0.25), where the 0.25 counts for displaying the phase as a sin wave.

Returns a tuple (period, ephemeris, magnitude of max fourier transform)
"""
@inline function find_period_ephemeris(X_::Array{Float64,1},Y_::Array{Float64,1},T_::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}})
    F::Array{ComplexF64}=map(t->Complex_NFT(X_,Y_,t),T_)
    i=loc_max(abs2.(F))
    return (T_[i], T_[i]*(angle(F[i])/2pi+0.25), abs2(F[i]))
end

"""
    normalize101!(arr_::Array{Float64})

Normalize the array in-place to mean=0 and std=1

Is used to feed find_period_ephemeris function; the normalized magnitude eliminates period aliasing due to cuasi-constant sample on the time series due to phase cancellation on the fourier transform sum.
"""
function normalize101!(arr_::Array{Float64})
    arr_.= (arr_.-mean(arr_))./std(arr_)
end

max_stars_repo_name: santiagohenao/NUDFT-III

max_issues_repo_name: santiagohenao/NUDFT-III

max_forks_repo_name: santiagohenao/NUDFT-III


text:
# Functions to create standard digraph matrices

# Adjaceny Matrix
function adjacency(G::SimpleDigraph)
    n = NV(G)
    A = zeros(Int, (n, n))

    # create a table from V to 1:n
    d = vertex2idx(G)
    E = elist(G)

    for e in E
        i = d[e[1]]
        j = d[e[2]]
        A[i, j] = 1
    end

    return A
end

# incidence matrix
function incidence(G::SimpleDigraph)
    n = NV(G)
    m = NE(G)
    M = spzeros(Int, n, m)
    d = vertex2idx(G)
    E = elist(G)

    idx = 0
    for e in E
        i = d[e[1]]
        j = d[e[2]]
        idx += 1
        M[i, idx] = 1
        M[j, idx] = -1
    end

    return M
end