text:
function searchsortedfirst_join(_f, v::AbstractVector, x, lo::T, hi::T, o::Ordering, ::Val{T2})::keytype(v) where T<:Integer where T2
    u = T(1)
    lo = lo - u
    hi = hi + u
    @inbounds while lo < hi - u
        m = midpoint(lo, hi)
        fvm = _f(v[m])::T2
        if lt(o, fvm , x)
            lo = m
        else
            hi = m
        end
    end
    return hi
end

# index of the last value of vector a that is less than or equal to x;
# returns 0 if x is less than all values of v.
function searchsortedlast_join(_f, v::AbstractVector, x, lo::T, hi::T, o::Ordering, ::Val{T2})::keytype(v) where T<:Integer where T2
    u = T(1)
    lo = lo - u
    hi = hi + u
    @inbounds while lo < hi - u
        m = midpoint(lo, hi)
        fvm = _f(v[m])::T2
        if lt(o, x, fvm)
            hi = m
        else
            lo = m
        end
    end
    return lo
end

# returns the range of indices of v equal to x
# if v does not contain x, returns a 0-length range
# indicating the insertion point of x
function searchsorted_join(_f, v::AbstractVector, x, ilo::T, ihi::T, o::Ordering, ::Val{T2})::UnitRange{keytype(v)} where T<:Integer where T2
    u = T(1)
    lo = ilo - u
    hi = ihi + u
    @inbounds while lo < hi - u
        m = midpoint(lo, hi)
        fvm = _f(v[m])::T2
        if lt(o, fvm, x)
            lo = m
        elseif lt(o, x, fvm)
            hi = m
        else
            a = searchsortedfirst_join(_f, v, x, max(lo,ilo), m, o, Val(T2))
            b = searchsortedlast_join(_f, v, x, m, min(hi,ihi), o, Val(T2))
            return a : b
        end
    end
    return (lo + 1) : (hi - 1)
end

function _fill_range_for_accelerated_join!(ranges, starts, loc, x, f, sz, chunk; threads = true)
    loc_cumsum = cumsum(loc)
    @_threadsfor threads for i in 1:length(x)
        _index_ = hash(f(x[i])) % chunk + 1
        if loc[_index_]
            actual_start = loc_cumsum[_index_]
            st = starts[actual_start]
            actual_start == length(starts) ? en = sz : en = starts[actual_start+1]-1
            ranges[i] = st:en
        else
            ranges[i] = 1:0
        end
    end
end
# TODO how the hashing behave for Categorical Arrays?
function _find_permute_and_fill_range_for_join!(ranges, dsr, dsl, oncols_right, oncols_left, stable, alg, mapformats, accelerate, chunk = 2^10; nsfpaj = true, threads = true)
    if isempty(dsr)
        idx = []
        fill!(ranges, 1:nrow(dsr))
        last_valid_range = -1
    else
        if accelerate
            if mapformats[2]
                _fr = getformat(dsr, oncols_right[1])
            else
                _fr = identity
            end
            grng = _divide_for_fast_join(_columns(dsr)[oncols_right[1]], _fr, chunk; threads = threads)
            if mapformats[1]
                _fl = getformat(dsl, oncols_left[1])
            else
                _fl = identity
            end
            _fill_range_for_accelerated_join!(ranges, grng.starts, grng.starts_loc, _columns(dsl)[oncols_left[1]], _fl, nrow(dsr), chunk; threads = threads)
            if dsr isa SubDataset
                starts, idx, last_valid_range =  _sortperm_v(dsr, oncols_right, stable = stable, a = alg, mapformats = mapformats[2], notsortpaforjoin = nsfpaj, givenrange = grng, threads = threads)

            else
                starts, idx, last_valid_range =  _sortperm(dsr, oncols_right, stable = stable, a = alg, mapformats = mapformats[2], notsortpaforjoin = nsfpaj, givenrange = grng, threads = threads)
            end
        else
            if dsr isa SubDataset
                starts, idx, last_valid_range =  _sortperm_v(dsr, oncols_right, stable = stable, a = alg, mapformats = mapformats[2], notsortpaforjoin = nsfpaj, threads = threads)
            else
                starts, idx, last_valid_range =  _sortperm(dsr, oncols_right, stable = stable, a = alg, mapformats = mapformats[2], notsortpaforjoin = nsfpaj, threads = threads)
            end
            fill!(ranges, 1:nrow(dsr))
        end
    end
    idx, last_valid_range == length(idx)
end

function _sort_for_join_after_hash(dsr, oncols_right, stable, alg, mapformats, nsfpaj, grng; threads = true)
    if dsr isa SubDataset
        starts, idx, last_valid_range =  _sortperm_v(dsr, oncols_right, stable = stable, a = alg, mapformats = mapformats[2], notsortpaforjoin = nsfpaj, givenrange = grng, threads = threads)
    else
        starts, idx, last_valid_range =  _sortperm(dsr, oncols_right, stable = stable, a = alg, mapformats = mapformats[2], notsortpaforjoin = nsfpaj, givenrange = grng, threads = threads)
    end
end

# find count of each groups id
function _find_counts_for_join(groups, ngroups)
   res = zeros(Int, ngroups)
   for i in 1:length(groups)
       res[groups[i]] += 1
   end
   res
end

function _find_range_for_join!(ranges, ldata, gslots, reps, where, minval, sz; threads = true)
    @_threadsfor threads for i in 1:length(ldata)
        ranges[i] = _query_dictionary_for_join_int(identity, ldata[i], gslots, reps, where, minval, sz)
    end
end

Base.@propagate_inbounds function _fill_val_join!(x, r, val)
    @simd for i in r
        x[i] = val
    end
end

function _fill_val_join!(x, r2, val, inbits, r)
    cnt = 1
    lo = r2.start
    for i in r
        if inbits[i]
            x[cnt+lo-1] = val
            cnt += 1
        end
    end
end
# F1 and F2 are here for type stability when threads = false
function _find_ranges_for_join!(ranges, x, y, _fl::F1, _fr::F2, ::Val{T1}, ::Val{T2}; type = :both, threads = true) where {T1, T2, F1, F2}
    if type == :both
        @_threadsfor threads for i in 1:length(x)
            ranges[i] = searchsorted_join(_fr, y, DataAPI.unwrap(_fl(x[i]))::T1, ranges[i].start, ranges[i].stop, Base.Order.Forward, Val(T2))
        end
    # TODO having another elseif branch is better for performance.
    elseif type == :left || type == :rightstrict
        @_threadsfor threads for i in 1:length(x)
             _flx = _fl(DataAPI.unwrap(x[i]))::T1
            hi = searchsortedlast_join(_fr, y, _flx, ranges[i].start, ranges[i].stop, Base.Order.Forward, Val(T2))
            lo = ranges[i].start
            if type === :rightstrict
                ranges[i] = hi+1:ranges[i].stop
            else
                ranges[i] = lo:hi
            end
        end
    elseif type == :right || type == :leftstrict
        @_threadsfor threads for i in 1:length(x)
            _flx = _fl(DataAPI.unwrap(x[i]))::T1
            lo = searchsortedfirst_join(_fr, y, _flx, ranges[i].start, ranges[i].stop, Base.Order.Forward, Val(T2))
            hi = ranges[i].stop
            if type === :leftstrict
                    ranges[i] = ranges[i].start:lo-1
            else
                ranges[i] = lo:hi
            end
        end
    end
end

function _find_ranges_for_join_pa!(ranges, x, invpool, y, _fl::F1, _fr::F2, ::Val{T1}, ::Val{T2}; type = :both, threads = true) where {T1, T2, F1, F2}
    if type == :both
        @_threadsfor threads for i in 1:length(x)
            revmap_paval_ref = get(invpool, DataAPI.unwrap(_fl(x[i]))::T1, missing)
            if ismissing(revmap_paval_ref)
                ranges[i] = 1:0
            else
                #_fr is identity
                ranges[i] = searchsorted_join(identity, y, revmap_paval_ref, ranges[i].start, ranges[i].stop, Base.Order.Forward, Val(T2))
            end
        end
    elseif type == :left || type == :rightstrict
        @_threadsfor threads for i in 1:length(x)
            revmap_paval_ref = get(invpool, _fl(DataAPI.unwrap(x[i]))::T1, missing)
            if ismissing(revmap_paval_ref)
                ranges[i] = 1:0
            else
                #_fr is identity
                hi = searchsoredlast_join(identity, y, revmap_paval_ref, ranges[i].start, ranges[i].stop, Base.Order.Forward, Val(T2))
                lo = ranges[i].start
                if type === :rightstrict
                    ranges[i] = hi+1:ranges[i].stop
                else
                    ranges[i] = lo:hi
                end
            end
        end
    elseif type == :right || type == :leftstrict
        @_threadsfor threads for i in 1:length(x)
            revmap_paval_ref = get(invpool, _fl(DataAPI.unwrap(x[i]))::T1, missing)
            if ismissing(revmap_paval_ref)
                ranges[i] = 1:0
            else
                #_fr is identity
                lo = searchsortedfirst_join(identity, y, revmap_paval_ref, ranges[i].start, ranges[i].stop, Base.Order.Forward, Val(T2))
                hi = ranges[i].stop
                if type === :leftstrict
                    ranges[i] = ranges[i].start:lo-1
                else
                    range[i] = lo:hi
                end
            end
        end
    end

end

function _fill_oncols_left_table_left!(_res, x, ranges, en, total, fill_val; threads = true)
    @_threadsfor threads for i in 1:length(x)
        i == 1 ? lo = 1 : lo = en[i - 1] + 1
        hi = en[i]
        _fill_val_join!(_res, lo:hi, x[i])
    end
    @_threadsfor threads for i in en[length(x)]+1:total
        _res[i] = fill_val
    end
end

function _fill_oncols_left_table_inner!(_res, x, ranges, en, total; inbits = nothing, en2 = nothing, threads = true)
    if inbits === nothing
        @_threadsfor threads for i in 1:length(x)
            length(ranges[i]) == 0 && continue
            i == 1 ? lo = 1 : lo = en[i - 1] + 1
            hi = en[i]
            _fill_val_join!(_res, lo:hi, x[i])
        end
    else
        @_threadsfor threads for i in 1:length(x)
            length(ranges[i]) == 0 && continue
            if i == 1
                lo = 1
                lo2 = 1
            else
                lo = en[i - 1] + 1
                lo2 = en2[i-1] + 1
            end
            hi = en[i]
            # @show sum(view(inbits, lo:hi))
            # sum(view(inbits, lo:hi)) == 0 && continue
            hi2 = en2[i]
            _fill_val_join!(_res, lo2:hi2, x[i], inbits, lo:hi)
        end
    end
end

function _fill_oncols_left_table_anti!(_res, x, ranges, en, total; threads = true)
    @_threadsfor threads for i in 1:length(x)
        length(ranges[i]) != 0 && continue
        i == 1 ? lo = 1 : lo = en[i - 1] + 1
        hi = en[i]
        _fill_val_join!(_res, lo:hi, x[i])
    end
end

function _fill_right_cols_table_left!(_res, x, ranges, en, total, fill_val; threads = true)
    @_threadsfor threads for i in 1:length(ranges)
        i == 1 ? lo = 1 : lo = en[i - 1] + 1
        hi = en[i]
        length(ranges[i]) == 0 ? _fill_val_join!(_res, lo:hi, fill_val) : copyto!(_res, lo, x, ranges[i].start, length(ranges[i]))
    end
end

function _fill_right_col_range!(_res, r2, x, ranges, inbits, r)
    cnt = 1
    cnt_r = 1
    lo = r2.start
    for i in r
        if inbits[i]
            _res[lo+cnt-1] = x[ranges[cnt_r]]
            cnt += 1
        end
        cnt_r += 1
    end
end

function _fill_right_cols_table_inner!(_res, x, ranges, en, total; inbits = nothing, en2 = nothing, threads = true)
    if inbits === nothing
        @_threadsfor threads for i in 1:length(ranges)
            length(ranges[i]) == 0 && continue
            i == 1 ? lo = 1 : lo = en[i - 1] + 1
            hi = en[i]
            copyto!(_res, lo, x, ranges[i].start, length(ranges[i]))
        end
    else
        @_threadsfor threads for i in 1:length(ranges)
            length(ranges[i]) == 0 && continue
            if i == 1
                lo = 1
                lo2 = 1
            else
                lo = en[i - 1] + 1
                lo2 = en2[i-1] + 1
            end
            hi = en[i]

            hi2 = en2[i]
            _fill_right_col_range!(_res, lo2:hi2, x, ranges[i], inbits, lo:hi)
        end
    end
end

function _create_multiple_match_col_left(ranges, total_length)
    res = allocatecol(Bool, total_length)
    cnt = 0
    for i in 1:length(ranges)
        if length(ranges[i]) == 0
            cnt += 1
            res[cnt] = false
        else
            if length(ranges[i]) == 1
                cnt += 1
                res[cnt] = false
            else
                for j in ranges[i]
                    cnt += 1
                    res[cnt] = true
                end
            end
        end
    end
    res
end
function _create_multiple_match_col_inner(ranges, en, total_length)
    res = allocatecol(Bool, total_length)
    cnt = 0
    if en === nothing
        for i in 1:length(ranges)
            if length(ranges[i]) == 0
                nothing
            else
                if length(ranges[i]) == 1
                    cnt += 1
                    res[cnt] = false
                else
                    for j in ranges[i]
                        cnt += 1
                        res[cnt] = true
                    end
                end
            end
        end
    else
        for i in 1:length(ranges)
            if i == 1
                lo = 1
            else
                lo = en[i - 1] + 1
            end
            hi = en[i]
            if length(lo:hi) == 0
                nothing
            elseif length(lo:hi) == 1
                cnt += 1
                res[cnt] = false
            else
                for j in lo:hi
                    cnt += 1
                    res[cnt] = true
                end
            end
        end
    end

    res
end

ISLE(x, ::Missing) = true
ISLE(x, y) = (x <= y)
ISLE(::Missing, y) = false
ISLE(::Missing, ::Missing) = false

function _mark_lt_part!(inbits, x_l, x_r, _fl::F1, _fr::F2, ranges, r_perms, en, ::Val{T}; strict = false, threads = true) where {T, F1, F2}
    revised_ends = zeros(T, length(en))
    @_threadsfor threads for i in 1:length(ranges)
        if length(ranges[i]) == 0
            if i !== 1
                revised_ends[i] = 0
            end
            continue
        end
        i == 1 ? lo = 1 : lo = en[i - 1] + 1
        hi = en[i]
        total = 0
        cnt = 1
        for j in ranges[i]
            if strict
                inbits[lo + cnt - 1] = isless(_fl(x_l[i]), _fr(x_r[r_perms[j]]))
            else
                inbits[lo + cnt - 1] = ISLE(_fl(x_l[i]), _fr(x_r[r_perms[j]]))
            end
            total += inbits[lo + cnt - 1]
            cnt += 1
        end
        revised_ends[i] = total
        if total == 0
            ranges[i] = 1:0
        end
    end
    cumsum!(revised_ends, revised_ends)
end

function _change_refpool_find_range_for_join!(ranges, dsl, dsr, r_perms, oncols_left, oncols_right, lmf, rmf, j; type = :both, nsfpaj = true, threads = true)
    var_l = _columns(dsl)[oncols_left[j]]
    var_r = _columns(dsr)[oncols_right[j]]
    l_idx = oncols_left[j]
    r_idx = oncols_right[j]
    if lmf
        _fl = getformat(dsl, l_idx)
    else
        _fl = identity
    end
    if rmf
        _fr = getformat(dsr, r_idx)
    else
        _fr = identity
    end

    T1 = Core.Compiler.return_type(DataAPI.unwrap∘_fl, Tuple{eltype(var_l)})

    if DataAPI.refpool(var_r) !== nothing && nsfpaj
        # sort taken care for refs ordering of modified values, but we still need to change refs
        if _fr == identity
            var_r_cpy = var_r
        else
            var_r_cpy = map(_fr, var_r)
        end
        T2 = eltype(DataAPI.refarray(var_r_cpy))
        # now _fr must be identity
        _fr = identity
        # we should use invpool of right column
        _find_ranges_for_join_pa!(ranges, var_l, DataAPI.invrefpool(var_r_cpy), view(DataAPI.refarray(var_r_cpy), r_perms), _fl, _fr, Val(T1), Val(T2); type = type, threads = threads)
    else
        T2 = Core.Compiler.return_type(_fr, Tuple{eltype(var_r)})
        _find_ranges_for_join!(ranges, var_l, view(var_r, r_perms), _fl, _fr, Val(T1), Val(T2); type = type, threads = threads)
    end
end

function _join_left(dsl, dsr, ::Val{T}; onleft, onright, makeunique = false, mapformats = [true, true], stable = false, alg = HeapSort, check = true, accelerate = false, method = :sort, threads = true, multiple_match::Bool = false, multiple_match_name = :multiple, obs_id = [false, false], obs_id_name = :obs_id) where T
    isempty(dsl) && return copy(dsl)
    if method == :hash
        ranges, a, idx, minval, reps, sz, right_cols = _find_ranges_for_join_using_hash(dsl, dsr, onleft, onright, mapformats, makeunique, Val(T); threads = threads)
    elseif method == :sort
        oncols_left = onleft
        oncols_right = onright
        right_cols = setdiff(1:length(index(dsr)), oncols_right)
        if !makeunique && !isempty(intersect(_names(dsl), _names(dsr)[right_cols]))
            throw(ArgumentError("duplicate column names, pass `makeunique = true` to make them unique using a suffix automatically." ))
        end

        ranges = Vector{UnitRange{T}}(undef, nrow(dsl))
        if length(oncols_left) == 1 && nrow(dsr)>1
            success, result = _join_left_dict(dsl, dsr, ranges, oncols_left, oncols_right, right_cols, Val(T); makeunique = makeunique, mapformats = mapformats, check = check, threads = threads, multiple_match = multiple_match, multiple_match_name = multiple_match_name, obs_id = obs_id, obs_id_name = obs_id_name)
            if success
                return result
            end
        end
        idx, uniquemode = _find_permute_and_fill_range_for_join!(ranges, dsr, dsl, oncols_right, oncols_left, stable, alg, mapformats, accelerate; threads = threads)

        for j in 1:length(oncols_left)
            _change_refpool_find_range_for_join!(ranges, dsl, dsr, idx, oncols_left, oncols_right, mapformats[1], mapformats[2], j; threads = threads)
        end
    end
    new_ends = map(x -> max(1, length(x)), ranges)
    cumsum!(new_ends, new_ends)
    total_length = new_ends[end]

    if check
        @assert total_length < 10*nrow(dsl) "the output data set will be very large ($(total_length)×$(ncol(dsl)+length(right_cols))) compared to the left data set size ($(nrow(dsl))×$(ncol(dsl))), make sure that the `on` keyword is selected properly, alternatively, pass `check = false` to ignore this error."
    end

    if multiple_match
        multiple_match_col = _create_multiple_match_col_left(ranges, total_length)
    end

    res = []
    for j in 1:length(index(dsl))
        addmissing = false
        _res = allocatecol(_columns(dsl)[j], total_length, addmissing = false)
        if DataAPI.refpool(_res) !== nothing
            # fill_val = DataAPI.invrefpool(_res)[missing]
            _fill_oncols_left_table_left!(_res.refs, DataAPI.refarray(_columns(dsl)[j]), ranges, new_ends, total_length, missing;  threads = threads)
        else
            _fill_oncols_left_table_left!(_res, _columns(dsl)[j], ranges, new_ends, total_length, missing; threads = threads)
        end
        push!(res, _res)

    end
    if dsl isa SubDataset
        newds = Dataset(res, copy(index(dsl)), copycols = false)
    else
        newds = Dataset(res, Index(copy(index(dsl).lookup), copy(index(dsl).names), copy(index(dsl).format)), copycols = false)
    end

    for j in 1:length(right_cols)
        _res = allocatecol(_columns(dsr)[right_cols[j]], total_length)
        if DataAPI.refpool(_res) !== nothing
            fill_val = DataAPI.invrefpool(_res)[missing]
            _fill_right_cols_table_left!(_res.refs, view(DataAPI.refarray(_columns(dsr)[right_cols[j]]), idx), ranges, new_ends, total_length, fill_val; threads = threads)
        else
            _fill_right_cols_table_left!(_res, view(_columns(dsr)[right_cols[j]], idx), ranges, new_ends, total_length, missing; threads = threads)
        end
        push!(_columns(newds), _res)

        new_var_name = make_unique([_names(dsl); _names(dsr)[right_cols[j]]], makeunique = makeunique)[end]
        push!(index(newds), new_var_name)
        setformat!(newds, index(newds)[new_var_name], getformat(dsr, _names(dsr)[right_cols[j]]))
    end
    if multiple_match
        insertcols!(newds, ncol(newds)+1, multiple_match_name => multiple_match_col, unsupported_copy_cols = false, makeunique = makeunique)
    end
    if obs_id[1]
        obs_id_name1 = Symbol(obs_id_name, "_left")
        obs_id_left = allocatecol(nrow(dsl) < typemax(Int32) ? Int32 : Int64, total_length)
        _fill_oncols_left_table_left!(obs_id_left, 1:nrow(dsl), ranges, new_ends, total_length, missing; threads = threads)
        insertcols!(newds, ncol(newds)+1, obs_id_name1 => obs_id_left, unsupported_copy_cols = false, makeunique = makeunique)
    end
    if obs_id[2]
        obs_id_name2 = Symbol(obs_id_name, "_right")
        obs_id_right = allocatecol(T, total_length)
        _fill_right_cols_table_left!(obs_id_right, idx, ranges, new_ends, total_length, missing, threads = threads)
        insertcols!(newds, ncol(newds)+1, obs_id_name2 => obs_id_right, unsupported_copy_cols = false, makeunique = makeunique)
    end
    newds

end

function _join_left!(dsl::Dataset, dsr::AbstractDataset, ::Val{T}; onleft, onright, makeunique = false, mapformats = [true, true], stable = false, alg = HeapSort, check = true, accelerate = false, method = :sort, threads = true, multiple_match = false, multiple_match_name = :multiple, obs_id = [false, false], obs_id_name = :obs_id) where T
    isempty(dsl) && return dsl
    if method == :hash
        ranges, a, idx, minval, reps, sz, right_cols = _find_ranges_for_join_using_hash(dsl, dsr, onleft, onright, mapformats, makeunique, Val(T); threads = threads)
    elseif method == :sort
        oncols_left = onleft
        oncols_right = onright
        right_cols = setdiff(1:length(index(dsr)), oncols_right)
        if !makeunique && !isempty(intersect(_names(dsl), _names(dsr)[right_cols]))
            throw(ArgumentError("duplicate column names, pass `makeunique = true` to make them unique using a suffix automatically." ))
        end
        ranges = Vector{UnitRange{T}}(undef, nrow(dsl))
        if length(oncols_left) == 1 && nrow(dsr)>1
            success, result = _join_left!_dict(dsl, dsr, ranges, oncols_left, oncols_right, right_cols, Val(T); makeunique = makeunique, mapformats = mapformats, check = check, threads = threads, multiple_match = multiple_match, multiple_match_name = multiple_match_name, obs_id = obs_id, obs_id_name = obs_id_name)
            if success
                return result
            end
        end
        idx, uniquemode = _find_permute_and_fill_range_for_join!(ranges, dsr, dsl, oncols_right, oncols_left, stable, alg, mapformats, accelerate, threads = threads)
        for j in 1:length(oncols_left)
            _change_refpool_find_range_for_join!(ranges, dsl, dsr, idx, oncols_left, oncols_right, mapformats[1], mapformats[2], j, threads = threads)
        end
    end
    if !all(x->length(x) <= 1, ranges)
        throw(ArgumentError("`leftjoin!` can only be used when each observation in left data set matches at most one observation from right data set"))
    end

    new_ends = map(x -> max(1, length(x)), ranges)
    cumsum!(new_ends, new_ends)
    total_length = new_ends[end]

    if check
        @assert total_length < 10*nrow(dsl) "the output data set will be very large ($(total_length)×$(ncol(dsl)+length(right_cols))) compared to the left data set size ($(nrow(dsl))×$(ncol(dsl))), make sure that the `on` keyword is selected properly, alternatively, pass `check = false` to ignore this error."
    end

    if multiple_match
        multiple_match_col = _create_multiple_match_col_left(ranges, total_length)
    end

    for j in 1:length(right_cols)
        _res = allocatecol(_columns(dsr)[right_cols[j]], total_length)
        if DataAPI.refpool(_res) !== nothing
            fill_val = DataAPI.invrefpool(_res)[missing]
            _fill_right_cols_table_left!(_res.refs, view(DataAPI.refarray(_columns(dsr)[right_cols[j]]), idx), ranges, new_ends, total_length, fill_val, threads = threads)
        else
            _fill_right_cols_table_left!(_res, view(_columns(dsr)[right_cols[j]], idx), ranges, new_ends, total_length, missing, threads = threads)
        end
        push!(_columns(dsl), _res)
        new_var_name = make_unique([_names(dsl); _names(dsr)[right_cols[j]]], makeunique = makeunique)[end]
        push!(index(dsl), new_var_name)
        setformat!(dsl, index(dsl)[new_var_name], getformat(dsr, _names(dsr)[right_cols[j]]))
    end
    if multiple_match
        insertcols!(dsl, ncol(dsl)+1, multiple_match_name => multiple_match_col, unsupported_copy_cols = false, makeunique = makeunique)
    end
    if obs_id[1]
        obs_id_name1 = Symbol(obs_id_name, "_left")
        obs_id_left = allocatecol(nrow(dsl) < typemax(Int32) ? Int32 : Int64, total_length)
        _fill_oncols_left_table_left!(obs_id_left, 1:nrow(dsl), ranges, new_ends, total_length, missing; threads = threads)
        insertcols!(dsl, ncol(dsl)+1, obs_id_name1 => obs_id_left, unsupported_copy_cols = false, makeunique = makeunique)
    end
    if obs_id[2]
        obs_id_name2 = Symbol(obs_id_name, "_right")
        obs_id_right = allocatecol(T, total_length)
        _fill_right_cols_table_left!(obs_id_right, idx, ranges, new_ends, total_length, missing, threads = threads)
        insertcols!(dsl, ncol(dsl)+1, obs_id_name2 => obs_id_right, unsupported_copy_cols = false, makeunique = makeunique)
    end

    _modified(_attributes(dsl))
    dsl
end

function _join_inner(dsl, dsr::AbstractDataset, ::Val{T}; onleft, onright, onright_range = nothing , makeunique = false, mapformats = [true, true], stable = false, alg = HeapSort, check = true, accelerate = false, droprangecols = true, strict_inequality = [false, false], method = :sort, threads = true, onlyreturnrange = false, multiple_match = false, multiple_match_name = :multiple, obs_id = [false, false], obs_id_name = :obs_id) where T
    (isempty(dsl) || isempty(dsr)) && throw(ArgumentError("in `innerjoin` both left and right tables must be non-empty"))
    oncols_left = onleft
    oncols_right = onright
    type = :both
    right_range_cols = Int[]
    if onright_range !== nothing
        left_range_col = oncols_left[end]

        right_range_cols = index(dsr)[filter!(!isequal(nothing), collect(onright_range))]
        if droprangecols
            right_cols = setdiff(1:length(index(dsr)), [oncols_right; right_range_cols])
        else
            right_cols = setdiff(1:length(index(dsr)), oncols_right)
        end

        oncols_right = [oncols_right; first(right_range_cols)]
        if onright_range[1] !== nothing
            if strict_inequality[1]
                type = :leftstrict
            else
                type = :left
            end
        else
            if strict_inequality[2]
                type = :rightstrict
            else
                type = :right
            end
        end
    else
        right_cols = setdiff(1:length(index(dsr)), oncols_right)
    end
    if !makeunique && !isempty(intersect(_names(dsl), _names(dsr)[right_cols]))
        throw(ArgumentError("duplicate column names, pass `makeunique = true` to make them unique using a suffix automatically." ))
    end

    nsfpaj = true
    # if the columns for inequality like join are PA we cannot use the fast path
    if type != :both
        if any(i-> DataAPI.refpool(_columns(dsr)[i]) !== nothing, right_range_cols)
            nsfpaj = false
        end
    end
    # if (onright_range === nothing || length(onleft) > 1) is false, then we have inequality kind join with no exact match join
    if method == :hash && (onright_range === nothing || length(onleft) > 1)
        if onright_range !== nothing
            ranges, a, idx, minval, reps, sz, right_cols_2 = _find_ranges_for_join_using_hash(dsl, dsr, onleft[1:end-1], oncols_right[1:end-1], mapformats, true, Val(T); threads = threads)
            filter!(!=(0), reps)
            pushfirst!(reps, 1)
            cumsum!(reps, reps)
            pop!(reps)
            grng = GIVENRANGE(idx, reps, Int[], length(reps))
            starts, idx, last_valid_range = _sort_for_join_after_hash(dsr, right_range_cols[1], stable, alg, mapformats, nsfpaj, grng; threads = threads)
            _change_refpool_find_range_for_join!(ranges, dsl, dsr, idx, oncols_left, oncols_right, mapformats[1], mapformats[2], length(oncols_left); type = type, nsfpaj = nsfpaj, threads = threads)
        else
            ranges, a, idx, minval, reps, sz, right_cols = _find_ranges_for_join_using_hash(dsl, dsr, onleft, onright, mapformats, makeunique, Val(T); threads = threads)
        end
    else
        ranges = Vector{UnitRange{T}}(undef, nrow(dsl))
        if length(oncols_left) == 1 && type == :both && nrow(dsr)>1
            success, result =  _join_inner_dict(dsl, dsr, ranges, oncols_left, oncols_right, right_cols, Val(T); makeunique = makeunique, mapformats = mapformats, check = check, threads = threads, multiple_match = multiple_match, multiple_match_name = multiple_match_name, obs_id = obs_id, obs_id_name = obs_id_name)
            if success
                return result
            end
        end
        idx, uniquemode = _find_permute_and_fill_range_for_join!(ranges, dsr, dsl, oncols_right, oncols_left, stable, alg, mapformats, accelerate && (onright_range == nothing || length(oncols_right)>1); nsfpaj = nsfpaj, threads = threads)

        for j in 1:length(oncols_left)-1
            _change_refpool_find_range_for_join!(ranges, dsl, dsr, idx, oncols_left, oncols_right, mapformats[1], mapformats[2], j; nsfpaj = nsfpaj, threads = threads)
        end
        _change_refpool_find_range_for_join!(ranges, dsl, dsr, idx, oncols_left, oncols_right, mapformats[1], mapformats[2], length(oncols_left); type = type, nsfpaj = nsfpaj, threads = threads)
    end

    new_ends = map(length, ranges)
    cumsum!(new_ends, new_ends)
    total_length = new_ends[end]

    inbits = nothing
    revised_ends = nothing
    if length(right_range_cols) == 2
        inbits = zeros(Bool, total_length)
        # TODO any optimisation is needed for pa?
        _fl = identity
        _fr = identity
        if mapformats[1]
            _fl = getformat(dsl, left_range_col)
        end
        if mapformats[2]
            _fr = getformat(dsr, right_range_cols[2])
        end
        revised_ends = _mark_lt_part!(inbits, _columns(dsl)[left_range_col], _columns(dsr)[right_range_cols[2]], _fl, _fr, ranges, idx, new_ends, total_length < typemax(Int32) ? Val(Int32) : Val(Int64); strict = strict_inequality[2], threads = threads)
    end
    if length(right_range_cols) == 2
        total_length = sum(inbits)
    end

    if onlyreturnrange
        return ranges
    end
    if check
        @assert total_length < 10*nrow(dsl) "the output data set will be very large ($(total_length)×$(ncol(dsl)+length(right_cols))) compared to the left data set size ($(nrow(dsl))×$(ncol(dsl))), make sure that the `on` keyword is selected properly, alternatively, pass `check = false` to ignore this error."
    end
    if multiple_match
        multiple_match_col = _create_multiple_match_col_inner(ranges, revised_ends, total_length)
    end

    res = []
    for j in 1:length(index(dsl))
        _res = allocatecol(_columns(dsl)[j], total_length, addmissing = false)
        if DataAPI.refpool(_res) !== nothing
            _fill_oncols_left_table_inner!(_res.refs, DataAPI.refarray(_columns(dsl)[j]), ranges, new_ends, total_length; inbits = inbits, en2 = revised_ends, threads = threads)
        else
            _fill_oncols_left_table_inner!(_res, _columns(dsl)[j], ranges, new_ends, total_length; inbits = inbits, en2 = revised_ends, threads = threads)
        end
        push!(res, _res)
    end
    if dsl isa SubDataset
        newds = Dataset(res, copy(index(dsl)), copycols = false)
    else
        newds = Dataset(res, Index(copy(index(dsl).lookup), copy(index(dsl).names), copy(index(dsl).format)), copycols = false)
    end

    for j in 1:length(right_cols)
        _res = allocatecol(_columns(dsr)[right_cols[j]], total_length, addmissing = false)
        if DataAPI.refpool(_res) !== nothing
            _fill_right_cols_table_inner!(_res.refs, view(DataAPI.refarray(_columns(dsr)[right_cols[j]]), idx), ranges, new_ends, total_length; inbits = inbits, en2 = revised_ends, threads = threads)
        else
            _fill_right_cols_table_inner!(_res, view(_columns(dsr)[right_cols[j]], idx), ranges, new_ends, total_length; inbits = inbits, en2 = revised_ends, threads = threads)
        end
        push!(_columns(newds), _res)

        new_var_name = make_unique([_names(dsl); _names(dsr)[right_cols[j]]], makeunique = makeunique)[end]
        push!(index(newds), new_var_name)
        setformat!(newds, index(newds)[new_var_name], getformat(dsr, _names(dsr)[right_cols[j]]))
    end
    if multiple_match
        insertcols!(newds, ncol(newds)+1, multiple_match_name => multiple_match_col, unsupported_copy_cols = false, makeunique = makeunique)
    end
    if obs_id[1]
        obs_id_name1 = Symbol(obs_id_name, "_left")
        obs_id_left = allocatecol(nrow(dsl) < typemax(Int32) ? Int32 : Int64, total_length)
        _fill_oncols_left_table_inner!(obs_id_left, 1:nrow(dsl), ranges, new_ends, total_length; inbits = inbits, en2 = revised_ends, threads = threads)
        insertcols!(newds, ncol(newds)+1, obs_id_name1 => obs_id_left, unsupported_copy_cols = false, makeunique = makeunique)
    end
    if obs_id[2]
        obs_id_name2 = Symbol(obs_id_name, "_right")
        obs_id_right = allocatecol(T, total_length)
        _fill_right_cols_table_inner!(obs_id_right, idx, ranges, new_ends, total_length; inbits = inbits, en2 = revised_ends, threads = threads)
        insertcols!(newds, ncol(newds)+1, obs_id_name2 => obs_id_right, unsupported_copy_cols = false, makeunique = makeunique)
    end
    newds

end

function _in(dsl::AbstractDataset, dsr::AbstractDataset, ::Val{T}; onleft, onright, mapformats = [true, true], stable = false, alg = HeapSort, accelerate = false, threads = true) where T
    isempty(dsl) && return Bool[]
    oncols_left = onleft
    oncols_right = onright

    # use Set when there is only one column in `on`
    if length(oncols_right) == 1
        if mapformats[1]
            _fl = getformat(dsl, oncols_left[1])
        else
            _fl = identity
        end
        if mapformats[2]
            _fr = getformat(dsr, oncols_right[1])
        else
            _fr = identity
        end
        return _in_use_Set(_columns(dsl)[oncols_left[1]], _columns(dsr)[oncols_right[1]], _fl, _fr, threads = threads)
    end
    ranges = Vector{UnitRange{T}}(undef, nrow(dsl))
    idx, uniquemode = _find_permute_and_fill_range_for_join!(ranges, dsr, dsl, oncols_right, oncols_left, stable, alg, mapformats, accelerate, threads = threads)
    for j in 1:length(oncols_left)
        _change_refpool_find_range_for_join!(ranges, dsl, dsr, idx, oncols_left, oncols_right, mapformats[1], mapformats[2], j, threads = threads)
    end
    map(x -> length(x) == 0 ? false : true, ranges)
end

function _find_right_not_in_left(ranges, n, idx)
    res = trues(n)
    for i in 1:length(ranges)
        view(res, ranges[i]) .= false
    end
    findall(res)
end

function _fill_oncols_left_table_left_outer!(res, x, notinleft, en, total)
    # TODO when x is a pooled array, we cannot use Threads.@threads (maybe we should create pool first and then use Threads)
    for i in en[end]+1:total
        val = x[notinleft[i - en[end]]]
        res[i] = val
    end
end

function _fill_source_for_outer!(res, ranges, notinleft, lval, rval, en, total)
    cnt = 0
    for i in 1:length(ranges)
        if length(ranges[i]) == 0
            cnt += 1
            res[cnt] = lval
        else
            cnt += length(ranges[i])
        end
    end
    for i in en[end]+1:total
        res[i] = rval
    end
end

function _create_source_for_outer(ranges, notinleft, total_length, en)
    res = allowmissing(PooledArray(["left", "right", "both"]))
    resize!(res.refs, total_length)
    fill!(res.refs, get(res.invpool, "both", missing))
    _fill_source_for_outer!(res.refs, ranges, notinleft, get(res.invpool, "left", missing), get(res.invpool, "right", missing), en, total_length)
    res
end

function _create_multiple_match_col_outer(ranges, notinleft, total_length, new_ends)
    res = _missings(Bool, total_length)
    cnt = 0
    for i in 1:length(ranges)
        if length(ranges[i]) == 0
            cnt += 1
            res[cnt] = false
        else
            if length(ranges[i]) == 1
                cnt += 1
                res[cnt] = false
            else
                for j in ranges[i]
                    cnt += 1
                    res[cnt] = true
                end
            end
        end
    end
    res
end

function _join_outer(dsl, dsr::AbstractDataset, ::Val{T}; onleft, onright, makeunique = false, mapformats = [true, true], stable = false, alg = HeapSort, check = true, accelerate = false, method = :sort, threads = true, source::Bool = false, source_col_name = :source, multiple_match = false, multiple_match_name = :multiple, obs_id = [false, false], obs_id_name = :obs_id) where T
    (isempty(dsl) || isempty(dsr)) && throw(ArgumentError("in `outerjoin` both left and right tables must be non-empty"))
    oncols_left = onleft
    oncols_right = onright
    if method == :hash
        ranges, a, idx, minval, reps, sz, right_cols = _find_ranges_for_join_using_hash(dsl, dsr, onleft, onright, mapformats, makeunique, Val(T); threads = threads)
    elseif method == :sort
        right_cols = setdiff(1:length(index(dsr)), oncols_right)
        if !makeunique && !isempty(intersect(_names(dsl), _names(dsr)[right_cols]))
            throw(ArgumentError("duplicate column names, pass `makeunique = true` to make them unique using a suffix automatically." ))
        end
        ranges = Vector{UnitRange{T}}(undef, nrow(dsl))
        if length(oncols_left) == 1 && nrow(dsr)>1
            success, result = _join_outer_dict(dsl, dsr, ranges, oncols_left, oncols_right, oncols_left, oncols_right, right_cols, Val(T); makeunique = makeunique, mapformats = mapformats, check = check, threads = threads, source = source, source_col_name = source_col_name, multiple_match = multiple_match, multiple_match_name = multiple_match_name, obs_id = obs_id, obs_id_name = obs_id_name)
            if success
                return result
            end
        end
        idx, uniquemode = _find_permute_and_fill_range_for_join!(ranges, dsr, dsl, oncols_right, oncols_left, stable, alg, mapformats, accelerate, threads = threads)
        for j in 1:length(oncols_left)
            _change_refpool_find_range_for_join!(ranges, dsl, dsr, idx, oncols_left, oncols_right, mapformats[1], mapformats[2], j, threads = threads)
        end
    end
    new_ends = map(x -> max(1, length(x)), ranges)
    notinleft = _find_right_not_in_left(ranges, nrow(dsr), idx)
    cumsum!(new_ends, new_ends)
    total_length = new_ends[end] + length(notinleft)

    if check
        @assert total_length < 10*nrow(dsl) "the output data set will be very large ($(total_length)×$(ncol(dsl)+length(right_cols))) compared to the left data set size ($(nrow(dsl))×$(ncol(dsl))), make sure that the `on` keyword is selected properly, alternatively, pass `check = false` to ignore this error."
    end

    if source
        source_col = _create_source_for_outer(ranges, notinleft, total_length, new_ends)
    end
    if multiple_match
        multiple_match_col = _create_multiple_match_col_outer(ranges, notinleft, total_length, new_ends)
    end

    res = []
    for j in 1:length(index(dsl))
        _res = allocatecol(_columns(dsl)[j], total_length)
        if DataAPI.refpool(_res) !== nothing
            fill_val = DataAPI.invrefpool(_res)[missing]
            _fill_oncols_left_table_left!(_res.refs, DataAPI.refarray(_columns(dsl)[j]), ranges, new_ends, total_length, fill_val, threads = threads)
        else
            _fill_oncols_left_table_left!(_res, _columns(dsl)[j], ranges, new_ends, total_length, missing, threads = threads)
        end
        push!(res, _res)
    end
    for j in 1:length(oncols_left)
        _fill_oncols_left_table_left_outer!(res[oncols_left[j]], view(_columns(dsr)[oncols_right[j]], idx), notinleft, new_ends, total_length)
    end
    if dsl isa SubDataset
        newds = Dataset(res, copy(index(dsl)), copycols = false)
    else
        newds = Dataset(res, Index(copy(index(dsl).lookup), copy(index(dsl).names), copy(index(dsl).format)), copycols = false)
    end

    for j in 1:length(right_cols)
        _res = allocatecol(_columns(dsr)[right_cols[j]], total_length)
        if DataAPI.refpool(_res) !== nothing
            fill_val = DataAPI.invrefpool(_res)[missing]
            _fill_right_cols_table_left!(_res.refs, view(DataAPI.refarray(_columns(dsr)[right_cols[j]]), idx), ranges, new_ends, total_length, fill_val, threads = threads)
            _fill_oncols_left_table_left_outer!(_res.refs, view(DataAPI.refarray(_columns(dsr)[right_cols[j]]), idx), notinleft, new_ends, total_length)
        else
            _fill_right_cols_table_left!(_res, view(_columns(dsr)[right_cols[j]], idx), ranges, new_ends, total_length, missing, threads = threads)
            _fill_oncols_left_table_left_outer!(_res, view(_columns(dsr)[right_cols[j]], idx), notinleft, new_ends, total_length)
        end
        push!(_columns(newds), _res)
        new_var_name = make_unique([_names(dsl); _names(dsr)[right_cols[j]]], makeunique = makeunique)[end]
        push!(index(newds), new_var_name)
        setformat!(newds, index(newds)[new_var_name], getformat(dsr, _names(dsr)[right_cols[j]]))
    end
    if source
        insertcols!(newds, ncol(newds)+1, source_col_name => source_col, unsupported_copy_cols = false, makeunique = makeunique)
    end
    if multiple_match
        insertcols!(newds, ncol(newds)+1, multiple_match_name => multiple_match_col, unsupported_copy_cols = false, makeunique = makeunique)
    end
    if obs_id[1]
        # Note that the name convention of obs_id_name1 and name2 are used in other places
        obs_id_name1 = Symbol(obs_id_name, "_left")
        obs_id_left = allocatecol(nrow(dsl) < typemax(Int32) ? Int32 : Int64, total_length)
        _fill_oncols_left_table_left!(obs_id_left, 1:nrow(dsl), ranges, new_ends, total_length, missing, threads = threads)
        insertcols!(newds, ncol(newds)+1, obs_id_name1 => obs_id_left, unsupported_copy_cols = false, makeunique = makeunique)
    end
    if obs_id[2]
        obs_id_name2 = Symbol(obs_id_name, "_right")
        obs_id_right = allocatecol(T, total_length)
        _fill_right_cols_table_left!(obs_id_right, idx, ranges, new_ends, total_length, missing, threads = threads)
        _fill_oncols_left_table_left_outer!(obs_id_right, idx, notinleft, new_ends, total_length)
        insertcols!(newds, ncol(newds)+1, obs_id_name2 => obs_id_right, unsupported_copy_cols = false, makeunique = makeunique)
    end
    newds

end


text:
GoTypeEnv = Union{GoEnv,GomokuEnv}

function stone_features(pos)
  features = zeros(UInt8, pos.env.N, pos.env.N, 2pos.env.planes)

  num_deltas_avail = size(pos.board_deltas, 3)
  cumulative_deltas = cumsum(pos.board_deltas, dims=3)
  extra_dim_board = reshape(pos.board, size(pos.board)..., 1)
  last_eight = repeat(extra_dim_board, outer = (1, 1, pos.env.planes))

  # apply deltas to compute previous board states
  last_eight[:, :, 2:num_deltas_avail + 1] .= last_eight[:, :, 2:num_deltas_avail + 1] .- cumulative_deltas
  # if no more deltas are available, just repeat oldest board.
  last_eight[:, :, num_deltas_avail + 2:end] .= last_eight[:, :, num_deltas_avail + 1]

  features[:, :, 1:2:end] .= last_eight .== pos.to_play
  features[:, :, 2:2:end] .= last_eight .== -pos.to_play

  features .+ 0.
end

color_to_play_feature(pos) = pos.to_play * ones(UInt8, pos.env.N, pos.env.N, 1)

get_feats(player::MCTSPlayer) = get_feats(player.root.position)

get_feats(pos) = cat(dims=3, stone_features(pos), color_to_play_feature(pos))

ext: jl


text:
module AMP

using Gumbo
using Cascadia
using Cascadia: matchFirst
using HTTP
using AbstractTrees: PreOrderDFS
using Pkg: project
using FranklinContent: walkfiles, islink, isldjson

const root_dir = Ref("")

function set_root(root=dirname(project().path))
    root_dir[] = root
end

skip_nodes = Set(HTMLElement{sym} for sym in [:iframe, :applet, :audio, :canvas, :embed, :video, :img, :button, :form])

const sel_body = sel"body"
const sel_head = sel"head"
const amp_doc = HTMLDocument("html", HTMLElement(:html))
const amp_head = HTMLElement(:head)
push!(amp_doc.root, amp_head)
const amp_body = HTMLElement(:body)
push!(amp_doc.root, amp_body)

function amp_template()
    html = amp_doc.root
    empty!(html.children)
    empty!(html.attributes)
    empty!(amp_head.children)
    empty!(amp_body.children)

    setattr!(html, "amp", "")
    push!(html, amp_head)
    push!(html, amp_body)
    # amp js
    ampjs = HTMLElement(:script)
    setattr!(ampjs, "async", "")
    setattr!(ampjs, "src", "https://cdn.ampproject.org/v0.js")
    push!(amp_head, ampjs)

    # charset
    charset = HTMLElement(:meta)
    setattr!(charset, "charset", "utf-8")
    push!(amp_head, charset)

    # viewport
    viewport = HTMLElement(:meta)
    setattr!(viewport, "name", "viewport")
    setattr!(viewport, "content", "width=device-width,minimum-scale=1,initial-scale=1")
    push!(amp_head, viewport)

    # amp styles
    # amp-custom goes before boilerplate
    push!(amp_head, style_el)
    style1 = HTMLElement(:style)
    setattr!(style1, "amp-boilerplate", "")
    push!(style1.children, HTMLText("body{-webkit-animation:-amp-start 8s steps(1,end) 0s 1 normal both;-moz-animation:-amp-start 8s steps(1,end) 0s 1 normal both;-ms-animation:-amp-start 8s steps(1,end) 0s 1 normal both;animation:-amp-start 8s steps(1,end) 0s 1 normal both}@-webkit-keyframes -amp-start{from{visibility:hidden}to{visibility:visible}}@-moz-keyframes -amp-start{from{visibility:hidden}to{visibility:visible}}@-ms-keyframes -amp-start{from{visibility:hidden}to{visibility:visible}}@-o-keyframes -amp-start{from{visibility:hidden}to{visibility:visible}}@keyframes -amp-start{from{visibility:hidden}to{visibility:visible}}"))
    push!(amp_head, style1)

    style2 = HTMLElement(:style)
    setattr!(style2, "amp-boilerplate", "")
    push!(style2.children, HTMLText("body{-webkit-animation:none;-moz-animation:none;-ms-animation:none;animation:none}"))
    style2wrapper = HTMLElement(:noscript)
    push!(style2wrapper.children, style2)
    push!(amp_head, style2wrapper)

    (html, amp_head, amp_body)
end

const styles_str = Ref("")
const styles_scripts = Vector{String}()
const style_el = HTMLElement(:style)
setattr!(style_el, "amp-custom", "")
setattr!(style_el, "type", "text/css")

function fetch_style(el, styles)
    style_source = getattr(el, "href")
    if startswith(style_source, "/")
        push!(styles,
              read(joinpath(root_dir[], lstrip(style_source, '/')),
                   String))
    else
        push!(styles,
              String(HTTP.get(style_source).body))
    end
end

function process_head(in_head, out_head, styles)
    canonical_unset = true
    title_unset = true
    for el in in_head.children
        tp = typeof(el)
        if tp ∈ skip_nodes continue end
        if tp === HTMLElement{:link}
            if canonical_unset && islink(el, "canonical")
                push!(out_head, el)
                canonical_unset = false
            elseif islink(el, "icon")
                push!(out_head, el)
            elseif islink(el, "stylesheet")
                fetch_style(el, styles)
            end
        elseif title_unset && tp === HTMLElement{:title}
            push!(out_head, el)
            title_unset = false
        elseif tp === HTMLElement{:script} &&
            isldjson(el)
            push!(out_head, el)
        end
    end
end

function recursive_check(el)
    if typeof(el) <: HTMLElement
        if length(el.children) > 0
            for c in el.children
                recursive_check(c)
            end
        else
            @assert typeof(el) ∉ skip_nodes
        end
    end
end

function process_body(in_el, out_body, out_head, styles, lv=false)
    l = length(in_el.children)
    n = 1
    while n <= l
        el = in_el.children[n]
        tp = typeof(el)
        if tp ∈ skip_nodes
            deleteat!(in_el.children, n)
            l -= 1
            continue
        elseif tp === HTMLElement{:link} &&
            islink(el, "stylesheet")
            fetch_style(el, styles)
            deleteat!(in_el.children, n)
            l -= 1
        elseif tp === HTMLElement{:style}
            # NOTE: this is supposed to be a text element
            push!(styles, el.children[1].text)
            deleteat!(in_el.children, n)
            l -= 1
        elseif tp === HTMLElement{:script}
            # only include ld+json scripts, ignore the rest
            isldjson(el) && push!(out_head, el)
            deleteat!(in_el.children, n)
            l -= 1
        else
            if tp !== HTMLText
                delete!(el.attributes, "onclick")
                length(el.children) !== 0 && process_body(el, out_body, out_head, styles, true)
            end
            # only add top level children
            lv || push!(out_body, el)
            n += 1
        end
    end
end

function amppage(file)
    html = parsehtml(read(file, String))
    in_body, in_head = (matchFirst(sel_body, html.root),
                        matchFirst(sel_head, html.root))
    (out_html, out_head, out_body) = amp_template()
    setattr!(out_html, "amp", "")
    for (a, v) in html.root.attributes
        setattr!(out_html, a, v)
    end

    empty!(styles_scripts)
    styles_str[] = ""

    process_head(in_head, out_head, styles_scripts)

    process_body(in_body, out_body, out_head, styles_scripts)

    # add remaining styles to head
    ss = styles_str[] = join(styles_scripts, "\n") |>
        # NOTE: the replacement should be ordered from most frequent to rarest
        # # remove troublesome animations
        x -> replace(x, r"\s*?@(\-[a-zA-Z]+-)?keyframes\s+?.+?{\s*?.+?({.+?})+?\s*?}"s => "") |>
        # # remove !important hints
        x -> replace(x, r"!important" => "") |>
        # remove charset since not allowed
        x -> replace(x, r"@charset\s+\"utf-8\"\s*;?"i => "")
    # Ensure CSS is less than maximum amp size of 75KB
    # NOTE: this doesn't take into account inline css
    @assert length(ss) < 75000
    empty!(style_el.children)
    push!(style_el, HTMLText(ss))

    string(amp_doc)
end

function ampdir(target; dirs=["posts", "tag", "reads", "_rss"])
    @assert isdir(target) "Path $target is not a valid directory"
    dir = isdirpath(target) ? dirname(target) : target
    rpl = Regex("^$dir/") => ""
    cwd = pwd()
    cd(dir)

    amp_dir = joinpath(dir, "amp")
    if !isdir(amp_dir) mkpath(amp_dir) end
    proc_dirs = Dict{String, Nothing}()

    try
        for file in walkfiles(dir; exts=Set((".html",)),
                              dirs=Set(dirs),
                              ex_dirs=Set(["amp"]),
                              subdir=false)
            html = amppage(file)
            out_file = joinpath(amp_dir, replace(file, rpl))
            out_dir = dirname(out_file)
            if out_dir ∉ keys(proc_dirs) && !isdir(out_dir)
                mkpath(out_dir)
                proc_dirs[out_dir] = nothing
            end
            write(out_file, html)
        end
    finally
        cd(cwd)
    end
end

end


text:
export validate_file

"Note this reads out the streams to their ends"
function validate_equal!(a_stream::IO, b_stream::IO; buffer_len=4*1024)
    a_buf = Vector{UInt8}(buffer_len)
    b_buf = Vector{UInt8}(buffer_len)

    while(true)
        a_len_read = readbytes!(a_stream, a_buf)
        b_len_read = readbytes!(b_stream, b_buf)
        if a_buf[1:a_len_read] != b_buf[1:b_len_read]
            return false
        end
        @assert(a_len_read == b_len_read)
        if a_len_read<length(a_buf)
            return true #we have reached the end
        end
    end

end

"""
Downloads a file file from Swift and checks if it the same as the local file given by `fname`.
Note: this is a comprehense check. Normally you would prefer to just check Hash's.
This method is for testing checking that the hash checking method works.
"""
function validate_file(serv::PyObject, container::AbstractString, objectname::String, fname::String; verbose::Bool=false)
    get_stream(serv::PyObject, container, objectname; verbose=verbose) do rfp
        open(fname, "r") do lfp
            validate_equal!(rfp, lfp)
        end
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule LAMMPS_jll
using Base
using Base: UUID
using LazyArtifacts
include(path) = Base.include(@__MODULE__, path)
include(joinpath("..", ".pkg", "platform_augmentation.jl"))
import JLLWrappers

using Preferences

function known_abis()
    return (:MicrosoftMPI, :MPICH, :MPItrampoline)
end

const abi = @load_preference("abi", Sys.iswindows() ? :MicrosoftMPI : :MPICH)

function set_abi(abi)
    if abi ∉ known_abis()
        error("""
            The MPI ABI $abi is not supported.
            Please set the MPI ABI to one of the following:
            $(known_abis())
        """)
    end
    @set_preferences!("abi" => string(abi))
    @warn "The MPI abi has changed, you will need to restart Julia for the change to take effect" abi
    return abi
end

JLLWrappers.@generate_main_file_header("LAMMPS")
JLLWrappers.@generate_main_file("LAMMPS", UUID("5b3ab26d-9607-527c-88ea-8fe5ba57cafe"))
end  # module LAMMPS_jll

max_stars_repo_name: giordano/JLLWrappers

max_issues_repo_name: giordano/JLLWrappers

max_forks_repo_name: giordano/JLLWrappers


text:
# Tests for MLKernels
function t_mlkernel()

N = 10	# variables
M = 25  # samples

vectordata = rand(M)
matrixdata = rand(N,M)
labels = round.(2*rand(M))
multilabels = round.(2*rand(2,M))

A01 = j4pr.datacell(vectordata)
A02 = j4pr.datacell(vectordata, labels)
A03 = j4pr.datacell(vectordata, multilabels)
A04 = j4pr.datacell(matrixdata)
A05 = j4pr.datacell(matrixdata, labels)
A06 = j4pr.datacell(matrixdata, multilabels)

# Tests
Κ = [MLKernels.PolynomialKernel(),
     MLKernels.LinearKernel()]
for tr in [A01, A02, A03, vectordata]
	for ts in [A01, A02, A03, vectordata]
		for κ in Κ
			Test.@test try
				ts |> (tr |> j4pr.mlkernel(κ));
				true
			catch
				false
			end
		end
	end
end

for tr in [A04, A05, A06, matrixdata]
	for ts in [A04, A05, A06, matrixdata]
		for κ in Κ
			Test.@test try
				ts |> (tr |> j4pr.mlkernel(κ));
				true
			catch
				false
			end
		end
	end
end

# Additional test to cover all methods
κ = MLKernels.PolynomialKernel()
Test.@test try
	A01 |> j4pr.mlkernel(A01,κ)
	true
catch
	false
end

end


text:
function bitCount(n::Int)
	count = 0
	while n>0
		count += n & 1
		n >>= 1
	end
	return count
end

function findFirstBit(n::Int)
	n > 0 || throw(DomainError(" n needs to be larger than 0"))

	breakcondition = 0
	searchposition = 1
	idx = 0
	while breakcondition==0
		breakcondition = n & searchposition
		searchposition <<= 1
		idx += 1
	end

	return idx
end

function getBitPos(n::Int)

	Pos = Vector{Int}([])
	while n > 0

		idx = findFirstBit(n)
		push!(Pos,idx)

		n -= 2^(idx-1)

	end

	return Pos
end

"""
    bitSwap(n, i, j)

swaps bit i with j in n, for e.g.

n = 001001 , i = 1, j = 2 results in bitSwap(n,i,j) = 001010

#Arguments
- n: bit-string represented as an integer
- i: position of the first bit which should be swaped with the second
- j: position of the second bit which should be swaped with the first

"""
function bitSwap(n::Int, i::Int, j::Int)
	i -= 1
	j -= 1

	bit1 = (n >> i) & 1
	bit2 = (n >> j) & 1

	x = bit1 ⊻ bit2

	x = (x << i) | (x << j)

	return (n ⊻ x)
end

function FermionicBitSwap(n::Int, j::Int, k::Int)
	n_new = bitSwap(n,j,k)
	untouchedBit =  n & n_new
	inbetween = j<k ? (j+1:k-1) : (k+1:j-1)
	c = Int(sum(2.0.^[x-1 for x in inbetween]))
	overjumped_fermions = bitCount(untouchedBit & c)
	return n_new, (-1)^overjumped_fermions
end

max_stars_repo_name: Gawatz/exactQmodule

max_issues_repo_name: Gawatz/exactQmodule

max_forks_repo_name: Gawatz/exactQmodule


text:
#=
Created on Tuesday 04 February 2020
Last update: Wedneday 05 February 2020

@author: Michiel Stock
michielfmstock@gmail.com

Dijkstra, A* and Bellman FOrd
=#

module ShortestPath

export reconstruct_path, dijkstra, a_star

using DataStructures, STMO

"""
    reconstruct_path(previous::Dict{T,T}, source::T, sink::T) where {T}

Reconstruct the path from the output of the Dijkstra algorithm.

Inputs:
        - previous : a Dict with the previous node in the path
        - source : the source node
        - sink : the sink node
Ouput:
        - the shortest path from source to sink
"""
function reconstruct_path(previous::Dict{T,T}, source::T, sink::T) where {T}
    path = T[sink]
    current = sink
    while current != source
        current = previous[current]
        push!(path, current)
    end
    return reverse!(path)
end

"""
    dijkstra(graph::AdjList{R,T}, source::T) where {R<:Real,T}

Dijkstra without a specified sink. Give the `graph` and a `source` and this
function uses Dijkstra's algorithm to compute all distances between the nodes
and the source. Returns
    - `distances`: a dictionary with the distances
    - `previous`: a dictionary representing the tree of the shortest paths.
"""
function dijkstra(graph::AdjList{R,T}, source::T) where {R<:Real,T}
    # initialize the tentative distances
    distances = Dict(v => Inf for v in keys(graph))
    distances[source] = 0.0
    previous = Dict{T,T}()
    vertices_to_check = [(0.0, source)]
    while length(vertices_to_check) > 0
        dist, u = heappop!(vertices_to_check)
        for (dist_u_v, v) in graph[u]
            new_dist = dist + dist_u_v
            if distances[v] > new_dist
                distances[v] = new_dist
                previous[v] = u
                heappush!(vertices_to_check, (new_dist, v))
            end
        end
    end
    return distances, previous
end

"""
    dijkstra(graph::AdjList{R,T}, source::T, sink::T) where {R<:Real,T}

Dijkstra's shortest path algorithm.

Inputs:
    - `graph` : adjacency list representing a weighted directed graph
    - `source`
    - `sink`

Outputs:
    - the shortest path
    - the cost of this shortest path
"""
function dijkstra(graph::AdjList{R,T}, source::T, sink::T) where {R<:Real,T}
    # initialize the tentative distances
    distances = Dict(v => Inf for v in keys(graph))
    distances[source] = 0.0
    previous = Dict{T,T}()
    vertices_to_check = [(0.0, source)]
    while length(vertices_to_check) > 0
        dist, u = heappop!(vertices_to_check)
        if u == sink
            return reconstruct_path(previous, source, sink), dist
        end
        for (dist_u_v, v) in graph[u]
            new_dist = dist + dist_u_v
            if distances[v] > new_dist
                distances[v] = new_dist
                previous[v] = u
                heappush!(vertices_to_check, (new_dist, v))
            end
        end
    end
end

"""
    a_star(graph::AdjList{R,T}, source::T, sink::T, heuristic) where {R<:Real,T}

A* shortest path algorithm.

Inputs:
    - `graph` : adjacency list representing a weighted directed graph
    - `source`
    - `sink`
    - `heuristic` : a function that inputs a node and returns an lower bound
            for the distance to the source. Note that a distance can be turned into
            a heuristic using `n -> d(n, sink)`

Outputs:
    - the shortest path
    - the cost of this shortest path
"""
function a_star(graph::AdjList{R,T}, source::T, sink::T, heuristic) where {R<:Real,T}
    # initialize the tentative distances
    distances = Dict(v => Inf for v in keys(graph))
    distances[source] = 0.0
    previous = Dict{T,T}()
    vertices_to_check = [(heuristic(source), source)]
    while length(vertices_to_check) > 0
        dist, u = heappop!(vertices_to_check)
        if u == sink
            return reconstruct_path(previous, source, sink), distances[sink]
        end
        for (dist_u_v, v) in graph[u]
            new_dist = dist + dist_u_v
            if distances[v] > new_dist
                distances[v] = new_dist
                min_dist_to_sink = new_dist + heuristic(v)
                previous[v] = u
                heappush!(vertices_to_check, (min_dist_to_sink, v))
            end
        end
    end
end

end  # module ShortestPaths


text:
using ..CUDASupport
using ..TopOpt: @init_cuda
@init_cuda()

using ..GPUUtils
import ..TopOpt: whichdevice

whichdevice(s::AbstractDisplacementSolver) = whichdevice(s.u)

for T in (:DirectDisplacementSolver, :PCGDisplacementSolver)
    @eval @inline CuArrays.cu(p::$T) = error("$T does not support the GPU.")
end

function update_f!(
    f::CuVector{T},
    values,
    prescribed_dofs,
    applyzero,
    dof_cells,
    cell_dofs,
    black,
    white,
    Kes,
    xmin,
    penalty,
    vars,
    varind,
    M,
) where {T}
    args = (
        f,
        values,
        prescribed_dofs,
        applyzero,
        dof_cells.offsets,
        dof_cells.values,
        cell_dofs,
        black,
        white,
        Kes,
        xmin,
        penalty,
        vars,
        varind,
        M,
    )
    callkernel(dev, bc_kernel, args)
    CUDAdrv.synchronize(ctx)
    return nothing
end

function bc_kernel(
    f::AbstractVector{T},
    values,
    prescribed_dofs,
    applyzero,
    dof_cells_offsets,
    dof_cells_values,
    cell_dofs,
    black,
    white,
    Kes,
    xmin,
    penalty,
    vars,
    varind,
    M,
) where {T}
    ind = @thread_global_index()
    offset = @total_threads()
    while ind <= length(values)
        d = prescribed_dofs[ind]
        v = values[ind]
        m = size(Kes[ind], 1)

        r = dof_cells_offsets[d]:(dof_cells_offsets[d + 1] - 1)
        if !applyzero && v != 0
            for idx in r
                (i, j) = dof_cells_values[idx]
                if PENALTY_BEFORE_INTERPOLATION
                    px = ifelse(
                        black[i],
                        one(T),
                        ifelse(white[i], xmin, density(penalty(vars[varind[i]]), xmin)),
                    )
                else
                    px = ifelse(
                        black[i],
                        one(T),
                        ifelse(white[i], xmin, penalty(density(vars[varind[i]], xmin))),
                    )
                end
                if eltype(Kes) <: Symmetric
                    Ke = Kes[i].data
                else
                    Ke = Kes[i]
                end
                for row in 1:m
                    f[cell_dofs[row, i]] -= px * v * Ke[row, j]
                end
            end
        end
        f[d] = M * v
        ind += offset
    end
    return nothing
end

@define_cu(
    StaticMatrixFreeDisplacementSolver,
    :f,
    :problem,
    :vars,
    :cg_statevars,
    :elementinfo,
    :penalty,
    :prev_penalty,
    :u,
    :fixed_dofs,
    :free_dofs,
    :xes,
    :lhs,
    :rhs
)

whichdevice(m::MatrixFreeOperator) = whichdevice(m.vars)

function mul!(y::TV, A::MatrixFreeOperator, x::TV) where {TV<:CuArrays.CuVector}
    T = eltype(y)
    nels = length(A.elementinfo.Kes)
    ndofs = length(A.elementinfo.fixedload)
    meandiag = A.meandiag

    @unpack Kes, metadata, black, white, varind = A.elementinfo
    @unpack cell_dofs, dof_cells = metadata
    @unpack penalty, xmin, vars, fixed_dofs, free_dofs, xes = A

    args1 = (xes, x, black, white, vars, varind, cell_dofs, Kes, xmin, penalty, nels)
    callkernel(dev, mul_kernel1, args1)
    CUDAdrv.synchronize(ctx)

    args2 = (
        y, x, dof_cells.offsets, dof_cells.values, xes, fixed_dofs, free_dofs, meandiag
    )
    callkernel(dev, mul_kernel2, args2)
    CUDAdrv.synchronize(ctx)

    return y
end

# CUDA kernels
function mul_kernel1(
    xes::AbstractVector{TV},
    x,
    black,
    white,
    vars,
    varind,
    cell_dofs,
    Kes,
    xmin,
    penalty,
    nels,
) where {N,T,TV<:SVector{N,T}}
    i = @thread_global_index()
    offset = @total_threads()
    while i <= nels
        if PENALTY_BEFORE_INTERPOLATION
            px = ifelse(
                black[i],
                one(T),
                ifelse(white[i], xmin, density(penalty(vars[varind[i]]), xmin)),
            )
        else
            px = ifelse(
                black[i],
                one(T),
                ifelse(white[i], xmin, penalty(density(vars[varind[i]], xmin))),
            )
        end
        xe = xes[i]
        for j in 1:N
            xe = @set xe[j] = x[cell_dofs[j, i]]
        end
        if eltype(Kes) <: Symmetric
            xe = SVector{1,T}((px,)) .* (bcmatrix(Kes[i]).data * xe)
        else
            xe = SVector{1,T}((px,)) .* (bcmatrix(Kes[i]) * xe)
        end

        xes[i] = xe
        i += offset
    end

    return nothing
end

function mul_kernel2(
    y, x, dof_cells_offsets, dof_cells_values, xes, fixed_dofs, free_dofs, meandiag
)
    T = eltype(y)
    offset = @total_threads()
    n_fixeddofs = length(fixed_dofs)
    ndofs = length(y)

    i = @thread_global_index()
    while i <= n_fixeddofs
        dof = fixed_dofs[i]
        y[dof] = meandiag * x[dof]
        i += offset
    end

    i = @thread_global_index()
    while n_fixeddofs < i <= ndofs
        dof = free_dofs[i - n_fixeddofs]
        yi = zero(T)
        r = dof_cells_offsets[dof]:(dof_cells_offsets[dof + 1] - 1)
        for ind in r
            k, m = dof_cells_values[ind]
            yi += xes[k][m]
        end
        y[dof] = yi
        i += offset
    end
    return nothing
end

whichdevice(r::LinearElasticityResult) = whichdevice(r.u)

function fill_vars!(vars::GPUArray, topology; round)
    if round
        vars .= round.(typeof(vars)(topology))
    else
        copyto!(vars, topology)
    end
    return vars
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Libksba"
version = v"1.5.1"

# Collection of sources required to build libgcrypt
sources = [
    ArchiveSource("https://gnupg.org/ftp/gcrypt/libksba/libksba-$(version).tar.bz2",
                  "b0f4c65e4e447d9a2349f6b8c0e77a28be9531e4548ba02c545d1f46dc7bf921"),
]

# Bash recipe for building across all platforms

# Tried -no-undefined but still couldn't build for windows
script = raw"""
cd $WORKSPACE/srcdir/libksba-*/
export CPPFLAGS="-I${includedir}"
./configure --prefix=${prefix} --host=${target} --build=${MACHTYPE} --disable-static
make -j${nproc}
make install
if [[ "${target}" == x86_64-*-mingw* ]]; then
    # We have to manually build the shared library for Windows
    cc -shared -fPIC -o "${libdir}/libksba-8.${dlext}" -Wl,$(flagon --whole-archive) "${prefix}/lib/libksba.a" -Wl,$(flagon --no-whole-archive)  -lgpg-error
    rm "${prefix}/lib/libksba.a"
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.  We are manually disabling
# many platforms that do not seem to work.
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libksba", :libksba),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Libgpg_error_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: ctessum/Yggdrasil

max_issues_repo_name: ctessum/Yggdrasil

max_forks_repo_name: ctessum/Yggdrasil


text:
using Random
using LinearAlgebra
using Statistics
using Optim
using DataFrames
using CSV
using HTTP
using GLM
using FreqTables
using ForwardDiff # for bonus at the very end

function wrapall()

url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS3-gev/nlsw88w.csv"
df = CSV.read(HTTP.get(url).body)
X = [df.age df.white df.collgrad]
Z = hcat(df.elnwage1, df.elnwage2, df.elnwage3, df.elnwage4,
         df.elnwage5, df.elnwage6, df.elnwage7, df.elnwage8)
y = df.occupation

function mlogit_with_Z(theta, X, Z, y)

        alpha = theta[1:end-1]
        gamma = theta[end] # Indexing Z variables at the end
        K = size(X,2) # N Variables
        J = length(unique(y))
        N = length(y)
        bigY = zeros(N,J) # Same as PS2.jl ...
        for j=1:J
            bigY[:,j] = y.==j
        end
        bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

        T = promote_type(eltype(X),eltype(theta))
        num   = zeros(T,N,J)
        dem   = zeros(T,N)
        for j=1:J
            num[:,j] = exp.(X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)
            dem .+= num[:,j]
        end

        P = num./repeat(dem,1,J)

        loglike = -sum( bigY.*log.(P) )

        return loglike
    end

startvals = [2*rand(7*size(X,2)).-1; .1] # 7 == J - 1
td = TwiceDifferentiable(theta -> mlogit_with_Z(theta, X, Z, y), startvals; autodiff = :forward)
# run the optimizer
theta_hat_optim_ad = optimize(td, startvals, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
theta_hat_mle_ad = theta_hat_optim_ad.minimizer
# evaluate the Hessian at the estimates
H  = Optim.hessian!(td, theta_hat_mle_ad)
theta_hat_mle_ad_se = sqrt.(diag(inv(H)))
println([theta_hat_mle_ad theta_hat_mle_ad_se]) # these standard errors match Stata

function nested_logit_with_Z(theta, X, Z, y, nesting_structure)

        alpha = theta[1:end-3]
        lambda = theta[end-2:end-1]
        gamma = theta[end]
        K = size(X,2)
        J = length(unique(y))
        N = length(y)
        bigY = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j
        end
        bigAlpha = [repeat(alpha[1:K],1,length(nesting_structure[1])) repeat(alpha[K+1:2K],1,length(nesting_structure[2])) zeros(K)]

        T = promote_type(eltype(X),eltype(theta))
        num   = zeros(T,N,J)
        lidx  = zeros(T,N,J)
        dem   = zeros(T,N)
        for j=1:J
            if j in nesting_structure[1]
                lidx[:,j] = exp.( (X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)./lambda[1] )
            elseif j in nesting_structure[2]
                lidx[:,j] = exp.( (X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)./lambda[2] )
            else
                lidx[:,j] = exp.(zeros(N))
            end
        end
        for j=1:J
            if j in nesting_structure[1]
                num[:,j] = lidx[:,j].*sum(lidx[:,nesting_structure[1][:]];dims=2).^(lambda[1]-1)
            elseif j in nesting_structure[2]
                num[:,j] = lidx[:,j].*sum(lidx[:,nesting_structure[2][:]];dims=2).^(lambda[2]-1)
            else
                num[:,j] = lidx[:,j]
            end
            dem .+= num[:,j]
        end

        P = num./repeat(dem,1,J)

        loglike = -sum( bigY.*log.(P) )

        return loglike
    end
nesting_structure = [[1 2 3], [4 5 6 7]] # This is a list with the mentioned nested structure

startvals = [2*rand(2*size(X,2)).-1;
		1; 1; .1 # Unclear
		]
# Same as before
td = TwiceDifferentiable(theta -> nested_logit_with_Z(theta, X, Z, y, nesting_structure), startvals; autodiff = :forward)
# run the optimizer
nlogit_theta_hat_optim_ad = optimize(td, startvals, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
nlogit_theta_hat_mle_ad = nlogit_theta_hat_optim_ad.minimizer
# evaluate the Hessian at the estimates
H  = Optim.hessian!(td, nlogit_theta_hat_mle_ad)
nlogit_theta_hat_mle_ad_se = sqrt.(diag(inv(H)))
println([nlogit_theta_hat_mle_ad nlogit_theta_hat_mle_ad_se]) # these standard errors match Stata

    return nothing
end
wrapall()


text:
using Distributions
using Random: AbstractRNG, MersenneTwister, default_rng

"""
    f(x::Real)

the function to be integrated.
"""
f(x::Real) = x^10 - 1

"""
    mc_hitormiss(rng::AbstractRNG, N::Int, limits::NTuple{2, Real})

integration function
"""
function mc_hitormiss(f, rng::AbstractRNG, N::Integer, (lb, ub)::NTuple{2, Real})

    hit = 0

    for i = 1:N
        x = rand(rng, Uniform(lb, ub))
        y_guess = rand(rng, Uniform(f(lb), f(ub)))
        y_guess < f(x) && (hit += 1)
    end

    @info "Hit $hit, from $N."

    return (ub-lb) * (f(ub)-f(lb)) * hit / N
end

function mc_hitormiss(f, N::Integer, (lb, ub)::NTuple{2, Real})
    mc_hitormiss(f, default_rng(), N::Integer, (lb, ub)::NTuple{2, Real})
end

function mc_importancesampling(N::Int, limits=[1., 2.])
    """
    integration function
    """
    hit = 0

end

println("Solution of integral: ",
        mc_hitormiss(f, 100000, (1., 2.)))


text:
using DeepcopyModules

using InteractiveUtils
using Test

module Outer end
@testset "deepcopy_module" begin
    m = @eval Outer module M
        const c = 1
        x = 2
    end

    m2 = DeepcopyModules.deepcopy_module(m)

    @test nameof(m2) == nameof(m)
    @test parentmodule(m2) == parentmodule(m)

    @test names(m2, all=true) == names(m, all=true)
end

@testset "mock module" begin
    m1 = @eval Outer module M1
        x = 1
        f1() = x+1
        f2() = 2
        f3() = f1() + f2()
    end

    m2 = DeepcopyModules.deepcopy_module(m1)

    # Override global variable
    @eval m2 x = $m1.x + 1
    @test m2.x != m1.x
    @test m2.f1() == m1.f1() + 1

    # Override function
    m2.f2() = m1.f1() + 1
    @test m2.f2() == m1.f1() + 1

    # Test overridden references
    @test m2.f3() == m1.f3() + 2  # (+1 for each change)
end

@testset "deepcopy_function" begin
    m1 = @eval module M1 end
    @eval m1 begin
        f() = 2
        foo() = f() + f(); m = methods(foo, ()).ms[1];
    end

    m2 = @eval module M2 end
    DeepcopyModules.deepcopy_function(M2, M1.foo)

    @test M2.foo != M1.foo
    # Successfully pointing at a _new version of f_!
    @test_throws UndefVarError(:f) M2.foo()

    # So copy over f as well
    DeepcopyModules.deepcopy_function(M2, M1.f)
    # And now it works
    @test M2.foo() == 4

    # Redefine f
    M2.f() = 3
    @test M2.foo() == 6

    # And the original is unchanged!
    @test M1.foo() == 4
end

@testset "functions with typevar" begin
    @testset "simple UnionAll" begin
        m1 = @eval module M1
            foo(x::T) where {T} = T
        end
        @assert m1.foo(5) == Int

        m2 = @eval module M2 end
        DeepcopyModules.deepcopy_function(M2, M1.foo)

        @test M2.foo != M1.foo
        @test M2.foo(5) == Int
    end

    @testset "nested UnionAll" begin
        m1 = @eval module M1
            _dims(a::Array{T, N}) where T where N = N
            _eltype(a::Array{T, N}) where T where N = T
        end
        @assert m1._dims([1,2,3]) == 1
        @assert m1._eltype([1,2,3]) == Int

        m2 = @eval module M2 end
        DeepcopyModules.deepcopy_function(m2, m1._dims)
        DeepcopyModules.deepcopy_function(m2, m1._eltype)

        @test m2._dims != m1._dims
        @test m2._eltype != m1._eltype
        @test m2._dims([1,2,3]) == 1
        @test m2._eltype([1,2,3]) == Int
    end
end

@testset "Modules with type definitions" begin
    # TODO: Currently struct definitions are not deep-copied, so this
    # test fails.
    m1 = @eval module M1
        struct S1 x end
    end
    m2 = deepcopy_module(M1)
    @test_broken @eval m2 S1() = S1(1)

    # Make sure that the new method in m2 isn't part of m1
    @test_throws MethodError m1.S1()
end

# ----- Implementation functions: --------

@testset "replace modules" begin
    m2 = @eval module M2 end
    m1 = @eval module M1
        f() = 3
        foo() = f() + f()
    end
    ci = @code_lowered(m1.foo())
    ci2 = DeepcopyModules.code_replace_module(ci, m1=>m2)
    @test ci2 != ci
    @test ci2.code[1].args[1].mod == m2
end

@testset "setparent!" begin
    m = @eval module M end
    @assert parentmodule(m) == Main

    DeepcopyModules.setparent!(m, m)
    @test parentmodule(m) == m

    DeepcopyModules.setparent!(m, Base)
    @test parentmodule(m) == Base

    @testset "create false submodule" begin
        m2 = @eval module M2 end
        DeepcopyModules.setparent!(m2, m)

        @test parentmodule(m2) == m
    end
end


text:
export getInterpolation

function getInterpolation(AT::SparseMatrixCSC, S::SparseMatrixCSC, coloring::Array, n::Int64)
	S.nzval .= 1.0
    S .= AT .* S
	## Pp - P.colptr
	## P.j - P.rowval
	## Here we build PT.
	Pp = getInterpolation1(AT, S, coloring, n) ## Essentially counts non-zeros per row.
	Pp .= Pp .+ 1
	#(Pp, Pj, Px) = getInterpolation2(AT, S, Pp, coloring, n)
	(Pp, Pj, Px) = getDirectInterpolation2(AT, S, Pp, coloring, n)
	Pj .= Pj .+ 1
	R = SparseMatrixCSC(maximum(Pj), n, Pp, Pj, Px)
	P = sparse(R')
	return P, R
end

function getInterpolation1(AT::SparseMatrixCSC, S::SparseMatrixCSC, coloring::Array, n::Int64)
	Pp = zeros(Int, size(AT.colptr))
	nz = 0
	for i=1:n
		if (coloring[i] == 1)
			nz += 1
		else
			for j in nzrange(S, i)
				if S.rowval[j] != i && coloring[S.rowval[j]] == 1
					nz += 1
				end
			end
		end
		Pp[i+1] = nz
	end
	return Pp
end

# Direct interpolation.
# A little faster
# Discussed here : https://computation.llnl.gov/projects/hypre-scalable-linear-solvers-multigrid-methods/yang1.pdf
# and on "Multigrid" book.
# First one I implemented for testing. Implementation from PyAMG (originally in C++).

function getDirectInterpolation2(AT::SparseMatrixCSC, S::SparseMatrixCSC, Pp::Array, coloring::Array, n::Int64)
	Px = zeros(Float64, Pp[end] - 1)
	Pj = zeros(Int, Pp[end] - 1)

	for i = 1:n
		# Make coarse node the i row of identity
		if coloring[i] == 1
			Pj[Pp[i]] = i
			Px[Pp[i]] = 1
			continue
		end

		# get positive and negative sums of nodes that are strongly connected to i
		(sum_strong_pos, sum_strong_neg) = getStrongSum(S, coloring, i)
		# get sum a_ij for j in 1:n , split to positive and negative
		(sum_all_pos, sum_all_neg, a_ii) = getAllSum(AT, coloring, i)

		i_alpha = sum_all_neg / sum_strong_neg
		i_beta = sum_all_pos / sum_strong_pos

		if sum_strong_pos == 0
			a_ii += sum_all_pos
			i_beta = 0.0
		end

		neg = -1 * i_alpha / a_ii
		pos = -1 * i_beta / a_ii

		nz = Pp[i]

		for j in nzrange(S, i)
			if coloring[S.rowval[j]] == 1 && S.rowval[j] != i
				Pj[nz] = S.rowval[j]
				if S.nzval[j] > 0
					Px[nz] = pos * S.nzval[j]
				else
					Px[nz] = neg * S.nzval[j]
				end
				nz += 1
			end
		end
	end

	# for node i in F : Pi = sum(w_ij*e_j). Here we have the transpose
	sum_map = zeros(Int, n)
	p_sum = 0
	for i = 1:n
		sum_map[i] = p_sum
		p_sum += coloring[i]
	end
	Pj .= sum_map[Pj]

   return Pp, Pj, Px
end

# get sum(a_ii + sum(a_in)) for n in weakly connected set
function getDenominator(AT::SparseMatrixCSC, S::SparseMatrixCSC, i::Int64)
	denominator = 0
	for j in nzrange(AT, i)
		denominator += AT.nzval[j]
	end

	for j in nzrange(S,i)
		if S.rowval[j] != i
			denominator -= S.nzval[j]
		end
	end
	return denominator
end

# Sum for strong coarse neighbours, returns tow sums : one for positive elements, and the other for negative elements
function getStrongSum(S::SparseMatrixCSC, coloring::Array, i::Int64)
	sum_positive = 0.0
	sum_negative = 0.0
	for j in nzrange(S, i)
		if coloring[S.rowval[j]] == 1 && S.rowval[j] != i
			if S.nzval[j] > 0
				sum_positive += S.nzval[j]
			else
				sum_negative += S.nzval[j]
			end
		end
	end
	sum_positive, sum_negative
end

function getAllSum(AT::SparseMatrixCSC, coloring::Array, i::Int64)
	sum_all_pos = 0.0
	sum_all_neg = 0.0
	diag = 0.0
	for j in nzrange(AT, i)
		if AT.rowval[j] == i
			diag += AT.nzval[j]
		else
			if AT.nzval[j] < 0
				sum_all_neg += AT.nzval[j]
			else
				sum_all_pos += AT.nzval[j]
			end
		end
	end
	sum_all_pos, sum_all_neg, diag
end

function getInnerDenominator(S::SparseMatrixCSC, coloring::Array, i::Int64, m::Int64)
	inner_denominator = 0
	# inner_denominator is sum of a_mk such that:
	# k is in C_i (coarse strong connection) and m is in F_i (fine strong connection)
	for k in nzrange(S,S.rowval[m])
		for kk in nzrange(S,i)
			if S.rowval[k] == S.rowval[kk] && coloring[S.rowval[k]] == 1
				inner_denominator+=S.nzval[k]
			end
		end
	end
	return inner_denominator
end

## This is the textbook version, discussed in "A Multigrid Tutorial".
# We calculate the i component of the operator :
# For every coarse node, the row is unity (e_i)
# For every free node, we calculate the weight and put sum(w_ij * e_j)

function getInterpolation2(AT::SparseMatrixCSC, S::SparseMatrixCSC, Pp::Array, coloring::Array, n::Int64)

	Px = zeros(Float64, Pp[end] - 1)
	Pj = zeros(Int, Pp[end] - 1)

	for i=1:n
		# Make coarse node the i row of identity
		## coloring = 1: C node, 0: F node
		if coloring[i] == 1
			Pj[Pp[i]] = i
			Px[Pp[i]] = 1
			continue
		end
		nz = Pp[i]
		# denominator holds all a_ij for weakly connected i-j.
		# with sparse arrays it's easier to take all, and subtract strong connections.
		i_denominator = getDenominator(AT, S, i)

		assert(i_denominator != 0)
		#j_numerator : AT[j][i] + ( AT[m][i] * AT[j][m] /  ( AT[k][m] for k in nz(S,i) if coloring(k) == 1) for m in nz(S,i) if coloring (k == 0))

		for j in nzrange(S,i)
			if S.rowval[j] == i || coloring[S.rowval[j]] == 0
				continue
			end

			Pj[nz] = S.rowval[j]

			j_numerator = S.nzval[j]
			for m in nzrange(S,i)
				# Ignore strong coarse connections
				if S.rowval[m] == i || coloring[S.rowval[m]] == 1
					continue
				end
				# check that AT[m][i], AT[j][m] != 0, otherwise there's no contribution to the numerator.
				inner_denominator = getInnerDenominator(S, coloring, i, m)
				assert(inner_denominator != 0)
				for l in nzrange(AT, S.rowval[m])
					if AT.rowval[l] == S.rowval[j]
						j_numerator += (S.nzval[m] * AT.nzval[l]) / inner_denominator
					end
				end
			end
			Px[nz] = - (j_numerator / i_denominator)
			nz += 1
		end
	end

	# for node i in F : Pi = sum(w_ij*e_j). Here we have the transpose.
	sum_map = zeros(Int, n)
	p_sum = 0
	for i = 1:n
		sum_map[i] = p_sum
		p_sum += coloring[i]
	end
	Pj .= sum_map[Pj]
	return Pp, Pj, Px
end


text:
immutable Gamma <: ContinuousUnivariateDistribution
    shape::Float64
    scale::Float64
    function Gamma(sh::Real, sc::Real)
        sh > zero(sh) && sc > zero(sc) ||
            error("Both shape and scale must be positive")
        new(float64(sh), float64(sc))
    end
end

Gamma(sh::Real) = Gamma(sh, 1.0)
Gamma() = Gamma(1.0, 1.0) # Standard exponential distribution

scale(d::Gamma) = d.scale
rate(d::Gamma) = 1.0 / d.scale

@_jl_dist_2p Gamma gamma

function entropy(d::Gamma)
    x = (1.0 - d.shape) * digamma(d.shape)
    x + lgamma(d.shape) + log(d.scale) + d.shape
end

insupport(::Gamma, x::Real) = zero(x) <= x < Inf
insupport(::Type{Gamma}, x::Real) = zero(x) <= x < Inf

kurtosis(d::Gamma) = 6.0 / d.shape

mean(d::Gamma) = d.shape * d.scale

median(d::Gamma) = quantile(d, 0.5)

mgf(d::Gamma, t::Real) = (1.0 - t * d.scale)^(-d.shape)

cf(d::Gamma, t::Real) = (1.0 - im * t * d.scale)^(-d.shape)

function mode(d::Gamma)
    d.shape >= 1.0 ? d.scale * (d.shape - 1.0) : error("Gamma has no mode when shape < 1.0")
end

modes(d::Gamma) = [mode(d)]

# rand()
#
#  A simple method for generating gamma variables - Marsaglia and Tsang (2000)
#  http://www.cparity.com/projects/AcmClassification/samples/358414.pdf
#  Page 369
#  basic simulation loop for pre-computed d and c
function randg2(d::Float64, c::Float64)
    while true
        x = v = 0.0
        while v <= 0.0
            x = randn()
            v = 1.0 + c * x
        end
        v = v^3
        U = rand()
        x2 = x^2
        if U < 1.0 - 0.331 * x2^2 ||
           log(U) < 0.5 * x2 + d * (1.0 - v + log(v))
            return d * v
        end
    end
end

# sampling from Gamma(α, 1)
function randg(α::Float64)
    dpar = (α <= 1.0 ? α + 1.0 : α) - 1.0 / 3.0
    cpar = 1.0 / sqrt(9.0 * dpar)
    randg2(dpar, cpar) * (α > 1.0 ? 1.0 : rand()^(1.0 / α))
end

rand(d::Gamma) = d.scale * randg(d.shape)

function rand!(d::Gamma, A::Array{Float64})
    α = d.shape
    dpar = (α <= 1.0 ? α + 1.0 : α) - 1.0 / 3.0
    cpar = 1.0 / sqrt(9.0 * dpar)
    n = length(A)
    for i in 1:n
        A[i] = randg2(dpar, cpar)
    end
    if α <= 1.0
        ainv = 1.0 / α
        for i in 1:n
            A[i] *= rand()^ainv
        end
    end
    multiply!(A, d.scale)
end

skewness(d::Gamma) = 2.0 / sqrt(d.shape)

var(d::Gamma) = d.shape * d.scale * d.scale

## Fit model

immutable GammaStats
    sx::Float64      # (weighted) sum of x
    slogx::Float64   # (weighted) sum of log(x)
    tw::Float64      # total sample weight

    GammaStats(sx::Real, slogx::Real, tw::Real) = new(float64(sx), float64(slogx), float64(tw))
end

function suffstats(::Type{Gamma}, x::Array)
    sx = 0.
    slogx = 0.
    for xi = x
        sx += xi
        slogx += log(xi)
    end
    GammaStats(sx, slogx, length(x))
end

function suffstats(::Type{Gamma}, x::Array, w::Array{Float64})
    n = length(x)
    if length(w) != n
        throw(ArgumentError("Inconsistent argument dimensions."))
    end

    sx = 0.
    slogx = 0.
    tw = 0.
    for i = 1:n
        xi = x[i]
        wi = w[i]
        sx += wi * xi
        slogx += wi * log(xi)
        tw += wi
    end
    GammaStats(sx, slogx, tw)
end

function gamma_mle_update(logmx::Float64, mlogx::Float64, a::Float64)
    ia = 1.0 / a
    z = ia + (mlogx - logmx + log(a) - digamma(a)) / (abs2(a) * (ia - trigamma(a)))
    1.0 / z
end

function fit_mle(::Type{Gamma}, ss::GammaStats;
    alpha0::Float64=NaN, maxiter::Int=1000, tol::Float64=1.0e-16)

    mx = ss.sx / ss.tw
    logmx = log(mx)
    mlogx = ss.slogx / ss.tw

    a::Float64 = isnan(alpha0) ? 0.5 / (logmx - mlogx) : alpha0
    converged = false

    t = 0
    while !converged && t < maxiter
        t += 1
        a_old = a
        a = gamma_mle_update(logmx, mlogx, a)
        converged = abs(a - a_old) <= tol
    end

    Gamma(a, mx / a)
end


text:
@testset "$TEST $G" begin

g3 = PathGraph(5, G)
@test adjacency_matrix(g3)[3,2] == 1
@test adjacency_matrix(g3)[2,4] == 0
@test laplacian_matrix(g3)[3,2] == -1
@test laplacian_matrix(g3)[1,3] == 0

# check adjacency matrices with self loops
g = copy(g3)
add_edge!(g,1,1)
@test adjacency_matrix(g)[1,1] == 2

g10 = CompleteGraph(10, G)
B, em = nonbacktracking_matrix(g10)
@test length(em) == 2*ne(g10)
@test size(B) == (2*ne(g10),2*ne(g10))
for i=1:10
    @test sum(B[:,i]) == 8
    @test sum(B[i,:]) == 8
end
@test !issymmetric(B)

@test adjacency_matrix(g3) ==
    adjacency_matrix(g3, :out) ==
    adjacency_matrix(g3, :in) ==
    adjacency_matrix(g3, :all)

@test_throws ErrorException adjacency_matrix(g3, :purple)

g5 = DG(4)
add_edge!(g5,1,2); add_edge!(g5,2,3); add_edge!(g5,1,3); add_edge!(g5,3,4)
#that call signature works
inmat   = adjacency_matrix(g5, :in, Int)
outmat  = adjacency_matrix(g5, :out, Int)
bothmat = adjacency_matrix(g5, :all, Int)

#relations that should be true
@test inmat' == outmat
@test all((bothmat - outmat) .>= 0)
@test all((bothmat - inmat)  .>= 0)

#check properties of the undirected laplacian carry over.
for dir in [:in, :out, :all]
    amat = adjacency_matrix(g5, dir, Float64)
    lmat = laplacian_matrix(g5, dir, Float64)
    @test isa(amat, SparseMatrixCSC{Float64, Int64})
    @test isa(lmat, SparseMatrixCSC{Float64, Int64})
    evals = eigvals(Matrix(lmat))
    @test all(evals .>= -1e-15) # positive semidefinite
    @test isapprox(minimum(evals),0, atol=1e-13)
end

g4 = PathDiGraph(5, DG)
# testing incidence_matrix, first directed graph
@test size(incidence_matrix(g4)) == (5,4)
@test incidence_matrix(g4)[1,1] == -1
@test incidence_matrix(g4)[2,1] == 1
@test incidence_matrix(g4)[3,1] == 0

g3 = PathGraph(5, G)
# now undirected graph
@test size(incidence_matrix(g3)) == (5,4)
@test incidence_matrix(g3)[1,1] == 1
@test incidence_matrix(g3)[2,1] == 1
@test incidence_matrix(g3)[3,1] == 0

i3o = incidence_matrix(g3; oriented=true)
@test i3o == incidence_matrix(g4)
@test laplacian_matrix(g3) == i3o * i3o'

n = 10; k = 5
pg = CompleteGraph(n, G)
# ϕ1 = nonbacktrack_embedding(pg, k)'
B, emap = nonbacktracking_matrix(pg)
## TODO add test

# spectral distance checks
for n=3:10
  polygon = random_regular_graph(n, 2)
  @test isapprox(spectral_distance(polygon, polygon), 0, atol=1e-8)
  @test isapprox(spectral_distance(polygon, polygon, 1), 0, atol=1e-8)
end

end # testset


text:
# K3.5 = D K3 D⁻¹ ,
#      [ I   0   0       0  ]
# D² = [ 0   I   0       0  ]
#      [ 0   0  S_l⁻¹    0  ]
#      [ 0   0   0     S_u⁻¹]
#
# [-Q - ρI    Aᵀ   √S_l  -√S_u ][ Δx  ]     [    -rc        ]
# [   A       δI    0      0   ][ Δy  ]     [    -rb        ]
# [  √S_l     0    X-L     0   ][Δs_l2] = D [σμe - (X-L)S_le]
# [ -√S_u     0     0     U-X  ][Δs_u2]     [σμe + (U-X)S_ue]
export K3_5StructuredParams

"""
Type to use the K3.5 formulation with a Krylov method, using the package
[`Krylov.jl`](https://github.com/JuliaSmoothOptimizers/Krylov.jl).
The outer constructor

    K3_5StructuredParams(; uplo = :U,   kmethod::Symbol = :trimr,
                         atol0 = 1.0e-4, rtol0 = 1.0e-4,
                         atol_min = 1.0e-10, rtol_min = 1.0e-10,
                         ρ_min = 1e4 * sqrt(eps()), δ_min = 1e4 * sqrt(eps()))

creates a [`RipQP.SolverParams`](@ref) that should be used to create a [`RipQP.InputConfig`](@ref).
The available methods are:
- `:tricg`
- `:trimr`

"""
struct K3_5StructuredParams <: SolverParams
  uplo::Symbol
  kmethod::Symbol
  atol0::Float64
  rtol0::Float64
  atol_min::Float64
  rtol_min::Float64
  ρ_min::Float64
  δ_min::Float64
end

function K3_5StructuredParams(;
  uplo::Symbol = :U,
  kmethod::Symbol = :trimr,
  atol0::T = 1.0e-4,
  rtol0::T = 1.0e-4,
  atol_min::T = 1.0e-10,
  rtol_min::T = 1.0e-10,
  ρ_min::T = 1e4 * sqrt(eps()),
  δ_min::T = 1e4 * sqrt(eps()),
) where {T <: Real}
  return K3_5StructuredParams(uplo, kmethod, atol0, rtol0, atol_min, rtol_min, ρ_min, δ_min)
end

mutable struct PreallocatedDataK3_5Structured{
  T <: Real,
  S,
  L1 <: LinearOperator,
  L2 <: LinearOperator,
  L3 <: LinearOperator,
  Ksol <: KrylovSolver,
} <: PreallocatedDataNewtonStructured{T, S}
  rhs1::S
  rhs2::S
  regu::Regularization{T}
  δv::Vector{T}
  As::L1
  Qreg::AbstractMatrix{T} # regularized Q
  QregF::LDLFactorizations.LDLFactorization{T, Int, Int, Int}
  Qregop::L2 # factorized matrix Qreg
  opBR::L3
  KS::Ksol
  atol::T
  rtol::T
  atol_min::T
  rtol_min::T
end

function opAsprod!(
  res::AbstractVector{T},
  nvar::Int,
  ncon::Int,
  ilow::Vector{Int},
  iupp::Vector{Int},
  nlow::Int,
  s_l::AbstractVector{T},
  s_u::AbstractVector{T},
  A::AbstractMatrix{T},
  v::AbstractVector{T},
  α::T,
  β::T,
  uplo::Symbol,
) where {T}
  if β == 0
    res[(ncon + 1):(ncon + nlow)] .= @views α .* sqrt.(s_l) .* v[ilow]
    res[(ncon + nlow + 1):end] .= @views (-α) .* sqrt.(s_u) .* v[iupp]
  else
    res[(ncon + 1):(ncon + nlow)] .=
      @views α .* sqrt.(s_l) .* v[ilow] .+ β .* res[(ncon + 1):(ncon + nlow)]
    res[(ncon + nlow + 1):end] .=
      @views (-α) .* sqrt.(s_u) .* v[iupp] .+ β .* res[(ncon + nlow + 1):end]
  end
  if uplo == :U
    @views mul!(res[1:ncon], A', v, α, β)
  else
    @views mul!(res[1:ncon], A, v, α, β)
  end
end

function opAstprod!(
  res::AbstractVector{T},
  nvar::Int,
  ncon::Int,
  ilow::Vector{Int},
  iupp::Vector{Int},
  nlow::Int,
  s_l::AbstractVector{T},
  s_u::AbstractVector{T},
  A::AbstractMatrix{T},
  v::AbstractVector{T},
  α::T,
  β::T,
  uplo::Symbol,
) where {T}
  if uplo == :U
    @views mul!(res, A, v[1:ncon], α, β)
  else
    @views mul!(res, A', v[1:ncon], α, β)
  end
  res[ilow] .+= @views α .* sqrt.(s_l) .* v[(ncon + 1):(ncon + nlow)]
  res[iupp] .-= @views α .* sqrt.(s_u) .* v[(ncon + nlow + 1):end]
end

function opBRprod!(
  res::AbstractVector{T},
  ncon::Int,
  nlow::Int,
  x_m_lvar::AbstractVector{T},
  uvar_m_x::AbstractVector{T},
  δv::AbstractVector{T},
  v::AbstractVector{T},
  α::T,
  β::T,
) where {T <: Real}
  if β == zero(T)
    res[1:ncon] .= @views (α / δv[1]) .* v[1:ncon]
    res[(ncon + 1):(ncon + nlow)] .= @views α ./ x_m_lvar .* v[(ncon + 1):(ncon + nlow)]
    res[(ncon + nlow + 1):end] .= @views α ./ uvar_m_x .* v[(ncon + nlow + 1):end]
  else
    res[1:ncon] .= @views (α / δv[1]) .* v[1:ncon] .+ β .* res[1:ncon]
    res[(ncon + 1):(ncon + nlow)] .=
      @views α ./ x_m_lvar .* v[(ncon + 1):(ncon + nlow)] .+ β .* res[(ncon + 1):(ncon + nlow)]
    res[(ncon + nlow + 1):end] .=
      @views α ./ uvar_m_x .* v[(ncon + nlow + 1):end] .+ β .* res[(ncon + nlow + 1):end]
  end
end

function update_kresiduals_history!(
  res::AbstractResiduals{T},
  Qreg::AbstractMatrix{T},
  As::AbstractLinearOperator{T},
  δ::T,
  solx::AbstractVector{T},
  soly::AbstractVector{T},
  rhs1::AbstractVector{T},
  rhs2::AbstractVector{T},
  s_l::AbstractVector{T},
  s_u::AbstractVector{T},
  x_m_lvar::AbstractVector{T},
  uvar_m_x::AbstractVector{T},
  nvar::Int,
  ncon::Int,
  nlow::Int,
  ilow::Vector{Int},
  iupp::Vector{Int},
) where {T <: Real}
  if typeof(res) <: ResidualsHistory
    @views mul!(res.Kres[1:nvar], Symmetric(Qreg, :U), solx)
    @views mul!(res.Kres[1:nvar], As', soly, one(T), one(T))
    @views mul!(res.Kres[(nvar + 1):end], As, solx)
    res.Kres[(nvar + 1):(nvar + ncon)] .+= δ .* soly[1:ncon]
    res.Kres[(nvar + ncon + 1):(nvar + ncon + nlow)] .+=
      @views s_l .* solx[ilow] .+ x_m_lvar .* soly[(ncon + 1):(ncon + nlow)]
    res.Kres[(nvar + ncon + nlow + 1):end] .+=
      @views .-s_u .* solx[iupp] .+ uvar_m_x .* soly[(ncon + nlow + 1):end]
    res.Kres[1:nvar] .-= rhs1
    res.Kres[(nvar + 1):end] .-= rhs2
  end
end

function PreallocatedData(
  sp::K3_5StructuredParams,
  fd::QM_FloatData{T},
  id::QM_IntData,
  itd::IterData{T},
  pt::Point{T},
  iconf::InputConfig{Tconf},
) where {T <: Real, Tconf <: Real}

  # init Regularization values
  if iconf.mode == :mono
    regu =
      Regularization(T(sqrt(eps()) * 1e3), T(sqrt(eps()) * 1e4), T(sp.ρ_min), T(sp.δ_min), :classic)
  else
    regu = Regularization(
      T(sqrt(eps()) * 1e5),
      T(sqrt(eps()) * 1e5),
      T(sqrt(eps(T)) * 1e0),
      T(sqrt(eps(T)) * 1e0),
      :classic,
    )
  end
  δv = [regu.δ] # put it in a Vector so that we can modify it without modifying opK2prod!
  As = LinearOperator(
    T,
    id.ncon + id.nlow + id.nupp,
    id.nvar,
    false,
    false,
    (res, v, α, β) -> opAsprod!(
      res,
      id.nvar,
      id.ncon,
      id.ilow,
      id.iupp,
      id.nlow,
      pt.s_l,
      pt.s_u,
      fd.A,
      v,
      α,
      β,
      fd.uplo,
    ),
    (res, v, α, β) -> opAstprod!(
      res,
      id.nvar,
      id.ncon,
      id.ilow,
      id.iupp,
      id.nlow,
      pt.s_l,
      pt.s_u,
      fd.A,
      v,
      α,
      β,
      fd.uplo,
    ),
  )

  Qreg = fd.Q + regu.ρ_min * I
  QregF = ldl(Symmetric(Qreg, :U))
  Qregop = LinearOperator(T, id.nvar, id.nvar, true, true, (res, v) -> ldiv!(res, QregF, v))
  Rd = δv[1] * I(id.ncon)
  opBR = LinearOperator(
    T,
    id.ncon + id.nlow + id.nupp,
    id.ncon + id.nlow + id.nupp,
    true,
    true,
    (res, v, α, β) -> opBRprod!(res, id.ncon, id.nlow, itd.x_m_lvar, itd.uvar_m_x, δv, v, α, β),
  )
  rhs1 = similar(fd.c, id.nvar)
  rhs2 = similar(fd.c, id.ncon + id.nlow + id.nupp)
  kstring = string(sp.kmethod)
  KS = eval(KSolver(sp.kmethod))(As', rhs1)

  return PreallocatedDataK3_5Structured(
    rhs1,
    rhs2,
    regu,
    δv,
    As,
    Qreg,
    QregF,
    Qregop,
    opBR,
    KS,
    sp.atol0,
    sp.rtol0,
    sp.atol_min,
    sp.rtol_min,
  )
end

function solver!(
  dd::AbstractVector{T},
  pad::PreallocatedDataK3_5Structured{T},
  dda::DescentDirectionAllocs{T},
  pt::Point{T},
  itd::IterData{T},
  fd::Abstract_QM_FloatData{T},
  id::QM_IntData,
  res::AbstractResiduals{T},
  cnts::Counters,
  T0::DataType,
  step::Symbol,
) where {T <: Real}
  Δs_l = itd.Δs_l
  Δs_u = itd.Δs_u
  pad.rhs1 .= @views step == :init ? fd.c : dd[1:(id.nvar)]
  if step == :init && all(pad.rhs1 .== zero(T))
    pad.rhs1 .= one(T)
  end
  pad.rhs2[1:(id.ncon)] .=
    @views (step == :init && all(dd[(id.nvar + 1):end] .== zero(T))) ? one(T) :
           dd[(id.nvar + 1):end]
  pad.rhs2[(id.ncon + 1):(id.ncon + id.nlow)] .= (step == :init) ? one(T) : Δs_l ./ sqrt.(pt.s_l)
  pad.rhs2[(id.ncon + id.nlow + 1):end] .= (step == :init) ? one(T) : Δs_u ./ sqrt.(pt.s_u)
  ksolve!(
    pad.KS,
    pad.As',
    pad.rhs1,
    pad.rhs2,
    pad.Qregop,
    pad.opBR,
    verbose = 0,
    atol = pad.atol,
    rtol = pad.rtol,
  )
  update_kresiduals_history!(
    res,
    pad.Qreg,
    pad.As,
    pad.regu.δ,
    pad.KS.x,
    pad.KS.y,
    pad.rhs1,
    pad.rhs2,
    pt.s_l,
    pt.s_u,
    itd.x_m_lvar,
    itd.uvar_m_x,
    id.nvar,
    id.ncon,
    id.nlow,
    id.ilow,
    id.iupp,
  )

  dd[1:(id.nvar)] .= @views pad.KS.x
  dd[(id.nvar + 1):end] .= @views pad.KS.y[1:(id.ncon)]
  Δs_l .= @views pad.KS.y[(id.ncon + 1):(id.ncon + id.nlow)] .* sqrt.(pt.s_l)
  Δs_u .= @views pad.KS.y[(id.ncon + id.nlow + 1):end] .* sqrt.(pt.s_u)

  return 0
end

# function update_upper_Qreg!(Qreg, Q, ρ)
#   n = size(Qreg, 1)
#   nnzQ = nnz(Q)
#   if nnzQ > 0
#     for i=1:n
#       diag_idx = Qreg.colptr[i+1] - 1
#       ptrQ = Q.colptr[i+1] - 1
#       Qreg.nzval[diag_idx] = (ptrQ ≤ nnzQ && Q.rowval[ptrQ] == i) ? Q.nzval[ptrQ] + ρ : ρ
#     end
#   end
# end

function update_pad!(
  pad::PreallocatedDataK3_5Structured{T},
  dda::DescentDirectionAllocs{T},
  pt::Point{T},
  itd::IterData{T},
  fd::Abstract_QM_FloatData{T},
  id::QM_IntData,
  res::AbstractResiduals{T},
  cnts::Counters,
  T0::DataType,
) where {T <: Real}
  if cnts.k != 0
    update_regu!(pad.regu)
  end
  # if cnts.k == 4
  # update_upper_Qreg!(pad.Qreg, fd.Q, pad.regu.ρ)
  # pad.Qreg[diagind(pad.Qreg)] .= fd.Q[diagind(fd.Q)] .+ pad.regu.ρ
  # ldl_factorize!(Symmetric(pad.Qreg, :U), pad.QregF)
  # end

  if pad.atol > pad.atol_min
    pad.atol /= 10
  end
  if pad.rtol > pad.rtol_min
    pad.rtol /= 10
  end

  pad.δv[1] = pad.regu.δ

  return 0
end


text:
# estimate observed time by bayes theorem
# multiple observation at 11 o'clock

module BayesEstAt11
    using DataFrames, CSV, FreqTables, Plots

    function bayes_estimation(sensor_value, cond_z_t, current_estimation)
        new_estimation = zeros(24)
        for i in Array(1:24)
            new_estimation[i] = cond_z_t[Name(sensor_value), Name(i - 1)] * current_estimation[i]
        end
        return new_estimation / sum(new_estimation)
    end

    function main(is_test=false)
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_600.txt")
        df_org = CSV.read(data_path, DataFrame,
                      header=["date", "time", "ir", "lidar"],
                      delim=' ')

        hours = [Int64(floor(e/10000)) for e in df_org.time]

        # add hour array to df
        df_h = DataFrame(hour=hours)
        df_new = hcat(df_org, df_h)

        # P(t)
        freqs = freqtable(df_new, :lidar, :hour)
        probs = freqs / length(df_org.lidar)
        keys = names(probs, 2)
        p_t = [sum(probs[begin:end, Name(key)]) for key in keys]
        cond_z_t = probs / p_t[findfirst(==(0), keys)]

        # estimate
        observations_11 = [617, 624, 619]
        estimation = p_t
        for observ in observations_11
            estimation = bayes_estimation(observ, cond_z_t, estimation)
        end
        plot(estimation, label="P(t=11|z=617,624,619)")

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/complex_dist/bayes_est_at_11/bayes_est_at_11.png")
            savefig(save_path)
        end
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy