text:
const INTEGERS = Union{Int8, Int16, Int32, Int64, UInt8, UInt16, UInt32, UInt64, Bool}
const FLOATS = Union{Float16, Float32, Float64}

# work around slow allocation of type union in julia
function _our_vect_alloc(T, len)
	if len > 0
		res = DataAPI.defaultarray(T, 1)(undef, 1)
		resize!(res, len)
	else
		DataAPI.defaultarray(T, 1)(undef, len)
	end
end

_missings(::Type{T}, len) where {T} = fill!(_our_vect_alloc(Union{T, Missing}, len), missing)

#macro for using Threaded for if needed
macro _threadsfor(threads, exp)
	esc(:(
		if $threads
			Threads.@threads $exp
		else
			$exp
		end
	))
end

# modified return_type to suit for our purpose
function return_type(f::Function, x)
	eltype(x) == Missing && return Missing
    CT = nonmissingtype(eltype(x))
    if CT <: AbstractVector
        return return_type_tuple(f, x)
    end
    T = Core.Compiler.return_type(f, (Vector{CT}, ))
    # workaround for SubArray type
    if T <: SubArray
        return Core.Compiler.return_type(f, (typeof(x), ))
    elseif T <: AbstractVector
        T = AbstractVector{Union{Missing, eltype(T)}}
    elseif T <: Tuple
        T = Union{Missing, Core.Compiler.return_type(f, (Vector{eltype(x)}, ))}
    else
        T = Union{Missing, T}
    end
    T
end

function return_type_tuple(f::Function, x)
    CT = ntuple(i -> nonmissingtype(eltype(x[i])), length(x))
    T = Core.Compiler.return_type(f, ntuple(i->Vector{CT[i]}, length(x)))
    # workaround for SubArray type
    if T <: SubArray
        return Core.Compiler.return_type(f, typeof.(x))
    elseif T <: AbstractVector
        T = AbstractVector{Union{Missing, eltype(T)}}
    elseif T <: Tuple
        T = Union{Missing, Core.Compiler.return_type(f, ntuple(i->Vector{eltype(x[i])}, length(x)))}
    else
        T = Union{Missing, T}
    end
    T
end

Missings.allowmissing(a::PooledArray) = convert(PooledArray{Union{Missing, eltype(a)}}, a)

function allocatecol(x::AbstractVector, len; addmissing = true)
    if addmissing
        @assert len > 0 "cannot allocate a column with length zero"
    end
    if DataAPI.refpool(x) !== nothing
        if x isa PooledArray
            _res = PooledArray(PooledArrays.RefArray(x.refs[1:1]), DataAPI.invrefpool(x), DataAPI.refpool(x), PooledArrays.refcount(x))
        else
            # TODO not optimised for Categorical arrays
            _res = copy(x)
        end
        resize!(_res, len)
        if addmissing
            _res = allowmissing(_res)
            _res[1] = missing
        end
    else
        _res = _our_vect_alloc(Union{Missing, eltype(x)}, len)
    end
    return _res
end

function allocatecol(T, len)
    _our_vect_alloc(Union{Missing, T}, len)
end

function _generate_inverted_dict_pool(x)
    invp = DataAPI.invrefpool(x)
    if invp isa Dict
        return Dict{valtype(invp), Union{Missing, keytype(invp)}}(values(invp) .=> keys(invp))
    elseif invp.invpool isa Dict
        a = Dict{valtype(invp.invpool), Union{Missing, keytype(invp.invpool)}}(values(invp.invpool) .=> keys(invp.invpool))
        push!(a, 0 => missing)
    else
        throw(ArgumentError("$(typeof(x)) is not supported, used PooledArray or Categorical Array"))
    end
end

function _hp_map_a_function!(fv, f, v)
    Threads.@threads for i in 1:length(v)
        fv[i] = f(v[i])
    end
end

function _hp_map!_a_function!(x, f)
    Threads.@threads for i in 1:length(x)
        x[i] = f(x[i])
    end
end

function _first_nonmiss(x)
    for i in 1:length(x)
        res = x[i]
        !ismissing(res) && return res
    end
    res
end

# define a structure for gathered data
mutable struct START_END
	start::Bool
	sz::Int
	where
end

function _f_barrier_give_end!(y, sz)
	for i in 1:length(y)-1
		y[i] = y[i+1] - 1
	end
	y[end] = sz
end
function _f_barrier_give_start!(y)
	for i in length(y):-1:2
		y[i] = y[i-1] + 1
	end
	y[1] = 1
end

function Base.reverse!(x::START_END)
	if x.start
		_f_barrier_give_end!(x.where, x.sz)
		x.start = false
		return x
	else
		_f_barrier_give_start!(x.where)
		x.start = true
		return x
	end
end

struct GIVENRANGE
    idx
    starts
	starts_loc
	lastvalid
end

function _sortitout!(res, starts, x)
    fill!(starts, 0)
    starts[1] = 1
    for i in 1:length(x)
        starts[x[i] + 1] += 1
    end
	starts_normalised = map(>(0), starts)
    cumsum!(starts, starts)
    for i in 1:length(x)
        label = x[i]
        res[starts[label]] = i
        starts[label] += 1
    end
	starts .-= 1
	reverse!(START_END(false, length(x), starts))
	return starts_normalised[2:end]
end

function _divide_for_fast_join_barrier!(res, starts, x, f, ::Val{T}; threads = true) where T
    nc = length(starts) - 1
    _hashed = Vector{T}(undef, length(x))
    @_threadsfor threads for i in 1:length(x)
        _hashed[i] = hash(f(x[i])) % nc + 1
    end
    starts_normalised = _sortitout!(res, starts, _hashed)
    return starts_normalised
end

function _remove_unwantedstarts!(starts, sz)
    curloc=2
    i = 1
    while true
        if starts[curloc] == starts[i]
            curloc += 1
        else
            i += 1
            starts[i] = starts[curloc]
        end
        starts[i]>sz && break
    end
    return resize!(starts, i-1)
end

function _divide_for_fast_join(x, f, chunk; threads = true) # chunk = 2^10 then data are divided to 1024 pieces
    T = length(x) < typemax(Int32) ? Int32 : Int64
    res = Vector{T}(undef, length(x))
    starts = Vector{T}(undef, chunk + 1)
    starts_loc = _divide_for_fast_join_barrier!(res, starts, x, f, chunk < typemax(UInt8) ? Val(UInt8) : chunk < typemax(UInt16) ? Val(UInt16) : Val(UInt32); threads = threads)
	starts = _remove_unwantedstarts!(starts, length(x))
	GIVENRANGE(res, starts, starts_loc, length(starts))
end

function _calculate_ends(groups, ngroups, ::Val{T}) where T
    where = zeros(T, ngroups)
    @inbounds for i = 1:length(groups)
        where[groups[i]] += 1
    end
    START_END(false, length(groups), cumsum!(where, where))
end

# From DataFrames.jl

function do_call(f::Base.Callable, incols::NTuple{2, AbstractVector}, r)
    return f(view(incols[1], r), view(incols[2], r))
end

function do_call(f::Base.Callable, incols::NTuple{3, AbstractVector}, r)
    return f(view(incols[1], r), view(incols[2], r),  view(incols[3], r))
end

function do_call(f::Base.Callable, incols::NTuple{4, AbstractVector}, r)
    return f(view(incols[1], r), view(incols[2], r),  view(incols[3], r), view(incols[4], r))
end

function do_call(f::Base.Callable, incols::Tuple, r)
    return f(map(c -> view(c, r), incols)...)
end

# Date & Time should be treated as integer
_date_value(::Missing) = missing
_date_value(x::TimeType) = Dates.value(x)::Int
_date_value(x::Period) = Dates.value(x)::Int
_date_value(x) = x

function _create_dictionary_unstable!(prev_groups, groups, gslots, rhashes, f, v, prev_max_group, ::Val{T}; threads = true) where T
	_which_to_process = _find_groups_with_more_than_one_observation(prev_groups, prev_max_group)[1]
	offsets = findall(_which_to_process)
	_sum_of_which_to_process = sum(_which_to_process)
	@_threadsfor threads for i in 1:length(v)
		@inbounds if _which_to_process[prev_groups[i]]
			rhashes[i] = hash(f(v[i]), hash(prev_groups[i]))
		end
	end

    n = length(v)
    # sz = 2 ^ ceil(Int, log2(n)+1)
    sz = length(gslots)
    # fill!(gslots, 0)
    @_threadsfor threads for i in 1:sz
        @inbounds gslots[i] = 0
    end
    szm1 = sz - 1
    ngroups = prev_max_group - _sum_of_which_to_process
    flag = true
	@inbounds for i in eachindex(rhashes)
		if _which_to_process[prev_groups[i]]
        	slotix = rhashes[i] & szm1 + 1
        	gix = -1
        	probe = 0
        	while true
            	g_row = gslots[slotix]
            	if g_row == 0
                	gslots[slotix] = i
                	gix = ngroups += 1
                	break
            #check hash collision
            	elseif rhashes[i] == rhashes[g_row]
                	if isequal(prev_groups[i],prev_groups[g_row]) && isequal(f(v[i]), f(v[g_row]))
                    	gix = groups[g_row]
                    	break
                	end
            	end
            	slotix = slotix & szm1 + 1
            	probe += 1
            	@assert probe < sz
        	end
        	groups[i] = gix
    	end
	end

	@_threadsfor threads for i in 1:length(rhashes)
		@inbounds if !_which_to_process[prev_groups[i]]
			pos = searchsortedlast(offsets, prev_groups[i])
			prev_groups[i] -= pos
		else
			@inbounds _which_to_process[prev_groups[i]] ? prev_groups[i] = groups[i] : nothing
		end
    end
	if threads
		ngroups = hp_maximum(prev_groups)
	else
		ngroups = stat_maximum(prev_groups)
	end

    if ngroups == n
        flag = false
        return flag, ngroups
    end

    # copy!(prev_groups, rhashes)
    return flag, ngroups
end

function _create_dictionary!(prev_groups, groups, gslots, rhashes, f, v, prev_max_group, stable, ::Val{T}; threads = true) where T
	if !stable
		return _create_dictionary_unstable!(prev_groups, groups, gslots, rhashes, f, v, prev_max_group, Val(T); threads = threads)
	end
    @_threadsfor threads for i in 1:length(v)
        @inbounds rhashes[i] = hash(f(v[i]), hash(prev_groups[i])) #hash(prev_groups[i]) is used to prevent reduce probe, the question is: is it working?
    end
    n = length(v)
    # sz = 2 ^ ceil(Int, log2(n)+1)
    sz = length(gslots)
    # fill!(gslots, 0)
    @_threadsfor threads for i in 1:sz
        @inbounds gslots[i] = 0
    end
    szm1 = sz - 1
    ngroups = 0
    flag = true
    @inbounds for i in eachindex(rhashes)
        slotix = rhashes[i] & szm1 + 1
        gix = -1
        probe = 0
        while true
            g_row = gslots[slotix]
            if g_row == 0
                gslots[slotix] = i
                gix = ngroups += 1
                break
            #check hash collision
            elseif rhashes[i] == rhashes[g_row]
                if isequal(prev_groups[i],prev_groups[g_row]) && isequal(f(v[i]), f(v[g_row]))
                    gix = groups[g_row]
                    break
                end
            end
            slotix = slotix & szm1 + 1
            probe += 1
            @assert probe < sz
        end
        groups[i] = gix
    end
    @_threadsfor threads for i in 1:length(rhashes)
        @inbounds prev_groups[i] = groups[i]
    end
    if ngroups == n
        flag = false
        return flag, ngroups
    end

    # copy!(prev_groups, rhashes)
    return flag, ngroups
end

function _create_dictionary_int_fast!(prev_groups, f, v, prev_max_group, minval, rangelen, ::Val{T}) where T
    offset = 1 - minval
    n = length(v)
    ngroups = 0
    flag = true

    remap = zeros(T, prev_max_group, rangelen + 1)

    @inbounds for i in 1:length(v)
        slotix = f(v[i]) + offset
        if ismissing(slotix)
            slotix = rangelen + 1
        end
		prv_grp = prev_groups[i]
        if remap[prv_grp, slotix] == 0
            ngroups += 1
            remap[prv_grp, slotix] = ngroups
            prev_groups[i] = remap[prv_grp, slotix]
        else
            prev_groups[i] = remap[prv_grp, slotix]
        end
    end
    if ngroups == n
        flag = false
    end

    return flag, ngroups
end

function _gather_groups(ds, cols, ::Val{T}; mapformats = false, stable = true, threads = true) where T
    colidx = index(ds)[cols]
    _max_level = nrow(ds)
    prev_max_group = UInt(1)
    prev_groups = ones(T, nrow(ds))
    groups = T[]
    # rhashes = Vector{UInt}(undef, nrow(ds))
    rhashes = UInt[]
    seen_nonint = false
    sz = max(1 + ((5 * _max_level) >> 2), 16)
    sz = 1 << (8 * sizeof(sz) - leading_zeros(sz - 1))
    @assert 4 * sz >= 5 * _max_level
    gslots = T[]

    for j in 1:length(colidx)
        _f = _date_value
        if mapformats
            _f = _date_value∘getformat(ds, colidx[j])
        end

        if DataAPI.refpool(_columns(ds)[colidx[j]]) !== nothing
            if _f == _date_value∘identity || !mapformats
                v = DataAPI.refarray(_columns(ds)[colidx[j]])
            else
                v = DataAPI.refarray(map(_f, _columns(ds)[colidx[j]]))
            end
            _f = identity
        else
            v = _columns(ds)[colidx[j]]
        end
        if nonmissingtype(Core.Compiler.return_type(_f, (nonmissingtype(eltype(v)),))) <: Union{Missing, INTEGERS}
			if threads
				_minval = hp_minimum(_f, v)
			else
				_minval = stat_minimum(_f, v)
			end
            if ismissing(_minval)
                continue
            else
                minval::Integer = _minval
            end
			if threads
				maxval::Integer = hp_maximum(_f, v)
			else
				maxval = stat_maximum(_f, v)
			end
			rnglen = BigInt(maxval) - BigInt(minval) + 1
			o1 = false
			if rnglen < typemax(Int)
				o1 = true
				rangelen = Int(rnglen)
			end
			o2 = false
			if o1 && BigInt(prev_max_group)*rangelen < 2*length(v)
				o2 = true
			end
            if o1 && o2 && maxval < typemax(Int)
                flag, prev_max_group = _create_dictionary_int_fast!(prev_groups, _f, v, prev_max_group, minval, rangelen, Val(T))
            else
                if !seen_nonint
                    seen_nonint = true
                    resize!(rhashes, nrow(ds))
                    resize!(gslots, sz)
                    resize!(groups, nrow(ds))
                end
				prev_max_group > nrow(ds)/100 ? stablegather = stable : stablegather = false
                flag, prev_max_group = _create_dictionary!(prev_groups, groups, gslots, rhashes, _f, v, prev_max_group, stablegather, Val(T); threads = threads)

            end
        else
            if !seen_nonint
                seen_nonint = true
                resize!(rhashes, nrow(ds))
                resize!(gslots, sz)
                resize!(groups, nrow(ds))
            end
			prev_max_group > nrow(ds)/100 ? stablegather = stable : stablegather = false
            flag, prev_max_group = _create_dictionary!(prev_groups, groups, gslots, rhashes, _f, v, prev_max_group, stablegather, Val(T); threads = threads)
        end
        !flag && break
    end
    return prev_groups, gslots, prev_max_group
end

function _find_groups_with_more_than_one_observation(groups, ngroups)
    res = trues(length(groups))
    seen_groups = falses(ngroups)

    _nonunique_barrier!(res, groups, seen_groups)

	fill!(seen_groups, false)

    _find_groups_with_more_than_one_observation_barrier!(res, groups, seen_groups)
	seen_groups, res

end

function _find_groups_with_more_than_one_observation_barrier!(res, groups, seen_groups)
    @inbounds for i in 1:length(res)
        res[i] && !seen_groups[groups[i]] ? seen_groups[groups[i]] = true : nothing
    end
    nothing
end

function _gather_groups_old_version(ds, cols, ::Val{T}; mapformats = false) where T
    colidx = index(ds)[cols]
    _max_level = nrow(ds)
    prev_max_group = UInt(1)
    prev_groups = ones(T, nrow(ds))
    groups = Vector{T}(undef, nrow(ds))
    rhashes = Vector{UInt}(undef, nrow(ds))
    sz = max(1 + ((5 * _max_level) >> 2), 16)
    sz = 1 << (8 * sizeof(sz) - leading_zeros(sz - 1))
    @assert 4 * sz >= 5 * _max_level
    gslots = Vector{T}(undef, sz)

    for j in 1:length(colidx)
        _f = identity
        if mapformats
            _f = getformat(ds, colidx[j])
        end

        if DataAPI.refpool(_columns(ds)[colidx[j]]) !== nothing
            if _f == identity
                v = DataAPI.refarray(_columns(ds)[colidx[j]])
            else
                v = DataAPI.refarray(map(_f, _columns(ds)[colidx[j]]))
            end
            _f = identity
        else
            v = _columns(ds)[colidx[j]]
        end
        flag, prev_max_group = InMemoryDatasets._create_dictionary!(prev_groups, groups, gslots, rhashes, _f, v, prev_max_group)

        !flag && break
    end
    return groups, gslots, prev_max_group
end

# ds assumes is grouped based on cols and groups are gathered togther
function _find_starts_of_groups(ds, cols::MultiColumnIndex, ::Val{T}; mapformats = true, threads = true) where T
    colsidx = index(ds)[cols]
    ranges = Vector{T}(undef, nrow(ds))
    inbits = zeros(Bool, nrow(ds))
    inbits[1] = true
    last_valid_index = 1
    for j in 1:length(colsidx)
        if mapformats
            _f = getformat(ds, colsidx[j])
        else
            _f = identity
        end
        _find_starts_of_groups!(_columns(ds)[colsidx[j]], _get_perms(ds; threads = threads), _f , inbits; threads = threads)
	all(inbits) && break
    end
    @inbounds for i in 1:length(inbits)
        if inbits[i]
            ranges[last_valid_index] = i
            last_valid_index += 1
        end
    end
    return collect(colsidx), ranges, (last_valid_index - 1)
end

_find_starts_of_groups(ds, col::ColumnIndex, ::Val{T}; mapformats = true, threads = true) where T = _find_starts_of_groups(ds, [col], Val(T), mapformats = mapformats, threads = threads)

function _find_starts_of_groups!(x, perm, f, inbits; threads = true)
    @_threadsfor threads for i in 2:length(inbits)
        @inbounds if !inbits[i]
		inbits[i] = !isequal(f(x[perm[i]]), f(x[perm[i-1]]))
	end
    end
end
# function _find_starts_of_groups!(x, perm, f, inbits, starts, ngroups)
# 	Threads.@threads for j in 1:ngroups
# 		i = starts[j]
# 		@inbounds inbits[i] = inbits[i]==1 ? 1 : !isequal(f(x[perm[i]]), f(x[perm[i-1]]))
# 	end
# end

function make_unique!(names::Vector{Symbol}, src::AbstractVector{Symbol};
                      makeunique::Bool=false)
    if length(names) != length(src)
        throw(DimensionMismatch("Length of src doesn't match length of names."))
    end
    seen = Set{Symbol}()
    dups = Int[]
	dups_dict = Dict{Symbol, Int}()
    for i in 1:length(names)
        name = src[i]
        if in(name, seen)
            push!(dups, i)
			if ismissing(get(dups_dict, src[i], missing))
				dups_dict[src[i]] = 1
			end
        else
            names[i] = src[i]
            push!(seen, name)
        end
    end

    if length(dups) > 0
        if !makeunique
            dupstr = join(string.(':', unique(src[dups])), ", ", " and ")
            msg = "Duplicate variable names: $dupstr. Pass makeunique=true " *
                  "to make them unique using a suffix automatically."
            throw(ArgumentError(msg))
        end
    end
    for i in dups
        nm = src[i]
		dup_info = get(dups_dict, src[i], missing)
		if ismissing(dup_info)
			throw(ErrorException("Something is wrong"))
		else
			k = dup_info
			cnt = 1
			while true
				newnm = Symbol("$(nm)_$(k)")
				if !in(newnm, seen)
	                names[i] = newnm
	                push!(seen, newnm)
					break
	            end
				k += 1
				cnt += 1
			end
			dups_dict[src[i]] += cnt
		end
    end

    return names
end

function make_unique(names::AbstractVector{Symbol}; makeunique::Bool=false)
    make_unique!(similar(names), names, makeunique=makeunique)
end

"""
    gennames(n::Integer)

Generate standardized names for columns of a DataFrame.
The first name will be `:x1`, the second `:x2`, etc.
"""
function gennames(n::Integer)
    res = Vector{Symbol}(undef, n)
    for i in 1:n
        res[i] = Symbol(@sprintf "x%d" i)
    end
    return res
end

function funname(f)
    local n
    try
        n = nameof(f)
    catch
        return :function
    end
    if String(n)[1] == '#'
        :function
    elseif String(n) == "Fix2"
        nameof(f.f)
    elseif String(n) == "Fix1"
        nameof(f.f)
    else
        n
    end
end

if isdefined(Base, :ComposedFunction) # Julia >= 1.6.0-DEV.85
    using Base: ComposedFunction
else
    using Compat: ComposedFunction
end

funname(c::ComposedFunction) = Symbol(funname(c.outer), :_, funname(c.inner))

_findall(B) = findall(B)

function _findall(B::AbstractVector{Bool})
    @assert firstindex(B) == 1
    nnzB = count(B)

    # fast path returning range
    nnzB == 0 && return 1:0
    len = length(B)
    nnzB == len && return 1:len
    start::Int = findfirst(B)
    nnzB == 1 && return start:start
    start + nnzB - 1 == len && return start:len
    stop::Int = findnext(!, B, start + 1) - 1
    start + nnzB == stop + 1 && return start:stop

    # slow path returning Vector{Int}
    I = Vector{Int}(undef, nnzB)
    @inbounds for i in 1:stop - start + 1
        I[i] = start + i - 1
    end
    cnt = stop - start + 2
    @inbounds for i in stop+1:len
        if B[i]
            I[cnt] = i
            cnt += 1
        end
    end
    @assert cnt == nnzB + 1
    return I
end

@inline _blsr(x) = x & (x-1)

# findall returning a range when possible (all true indices are contiguous), and optimized for B::BitVector
# the main idea is taken from Base.findall(B::BitArray)
function _findall(B::BitVector)::Union{UnitRange{Int}, Vector{Int}}
    nnzB = count(B)
    nnzB == 0 && return 1:0
    nnzB == length(B) && return 1:length(B)
    local I
    Bc = B.chunks
    Bi = 1 # block index
    i1 = 1 # index of current block beginng in B
    i = 1  # index of the _next_ one in I
    c = Bc[1] # current block

    start = -1 # the begining of ones block
    stop = -1  # the end of ones block

    @inbounds while true # I not materialized
        if i > nnzB # all ones in B found
            Ir = start:start + i - 2
            @assert length(Ir) == nnzB
            return Ir
        end

        if c == 0
            if start != -1 && stop == -1
                stop = i1 - 1
            end
            while c == 0 # no need to return here as we returned above
                i1 += 64
                Bi += 1
                c = Bc[Bi]
            end
        end
        if c == ~UInt64(0)
            if stop != -1
                I = Vector{Int}(undef, nnzB)
                for j in 1:i-1
                    I[j] = start + j - 1
                end
                break
            end
            if start == -1
                start = i1
            end
            while c == ~UInt64(0)
                if Bi == length(Bc)
                    Ir = start:length(B)
                    @assert length(Ir) == nnzB
                    return Ir
                end

                i += 64
                i1 += 64
                Bi += 1
                c = Bc[Bi]
            end
        end
        if c != 0 # mixed ones and zeros in block
            tz = trailing_zeros(c)
            lz = leading_zeros(c)
            co = c >> tz == (one(UInt64) << (64 - lz - tz)) - 1 # block of countinous ones in c
            if stop != -1  # already found block of ones and zeros, just not materialized
                I = Vector{Int}(undef, nnzB)
                for j in 1:i-1
                    I[j] = start + j - 1
                end
                break
            elseif !co # not countinous ones
                I = Vector{Int}(undef, nnzB)
                if start != -1
                    for j in 1:i-1
                        I[j] = start + j - 1
                    end
                end
                break
            else # countinous block of ones
                if start != -1
                    if tz > 0 # like __1111__ or 111111__
                        I = Vector{Int}(undef, nnzB)
                        for j in 1:i-1
                            I[j] = start + j - 1
                        end
                        break
                    else # lz > 0, like __111111
                        stop = i1 + (64 - lz) - 1
                        i += 64 - lz

                        # return if last block
                        if Bi == length(Bc)
                            Ir = start:stop
                            @assert length(Ir) == nnzB
                            return Ir
                        end

                        i1 += 64
                        Bi += 1
                        c = Bc[Bi]
                    end
                else # start == -1
                    start = i1 + tz

                    if lz > 0 # like __111111 or like __1111__
                        stop = i1 + (64 - lz) - 1
                        i += stop - start + 1
                    else # like 111111__
                        i += 64 - tz
                    end

                    # return if last block
                    if Bi == length(Bc)
                        Ir = start:start + i - 2
                        @assert length(Ir) == nnzB
                        return Ir
                    end

                    i1 += 64
                    Bi += 1
                    c = Bc[Bi]
                end
            end
        end
    end
    @inbounds while true # I materialized, process like in Base.findall
        if i > nnzB # all ones in B found
            @assert nnzB == i - 1
            return I
        end

        while c == 0 # no need to return here as we returned above
            i1 += 64
            Bi += 1
            c = Bc[Bi]
        end

        while c == ~UInt64(0)
            for j in 0:64-1
                I[i + j] = i1 + j
            end
            i += 64
            if Bi == length(Bc)
                @assert nnzB == i - 1
                return I
            end
            i1 += 64
            Bi += 1
            c = Bc[Bi]
        end

        while c != 0
            tz = trailing_zeros(c)
            c = _blsr(c) # zeros last nonzero bit
            I[i] = i1 + tz
            i += 1
        end
    end
    @assert false "should not be reached"
end

# this function is needed as == does not allow for comparison between tuples and vectors
function _equal_names(r1, r2)
    n1 = _getnames(r1)
    n2 = _getnames(r2)
    length(n1) == length(n2) || return false
    for (a, b) in zip(n1, n2)
        a == b || return false
    end
    return true
end

# a structure for vcat to vector without allocation
struct Cat2Vec{F1, F2, CT, T, S, A, B} <: AbstractVector{Union{T, S}}
	vec1::A
	vec2::B
	f1::F1
	f2::F2
	len1::Int
	len2::Int
	function Cat2Vec(x , y, f1::F1, f2::F2) where {F1, F2}
		if length(x) > length(y)
			if DataAPI.invrefpool(x) !== nothing
				if f1 != identity
					v = map(f1, x)
				else
					v = x
				end
				dict = DataAPI.invrefpool(v)
				# It is workaround for Categorical data
				if hasproperty(dict, :invpool)
					vtype = valtype(dict.invpool)
				else
					vtype = valtype(dict)
				end
				res = Vector{Union{Missing, vtype}}(undef, length(y))
				_rev_map_invrefpool!(res, dict, y, f2)
				new{typeof(identity), typeof(identity), Union{Missing, vtype}, Union{vtype, Missing}, Union{vtype, Missing}, typeof(res), typeof(res)}(DataAPI.refarray(v), res, identity, identity, length(x), length(y))
			elseif DataAPI.invrefpool(y) !== nothing
				if f2 != identity
					v = map(f2, y)
				else
					v = y
				end
				dict = DataAPI.invrefpool(v)
				if hasproperty(dict, :invpool)
					vtype = valtype(dict.invpool)
				else
					vtype = valtype(dict)
				end
				res =  Vector{Union{Missing, vtype}}(undef, length(x))
				_rev_map_invrefpool!(res, dict, x, f1)
				new{typeof(identity), typeof(identity), Union{Missing, vtype}, Union{vtype, Missing}, Union{vtype, Missing}, typeof(res), typeof(res)}(res, DataAPI.refarray(v), identity, identity, length(x), length(y))
			else
				new{F1, F2, promote_type(Core.Compiler.return_type(f1, (eltype(x), )), Core.Compiler.return_type(f2, (eltype(y), ))), eltype(x), eltype(y), typeof(x), typeof(y)}(x, y, f1, f2, length(x), length(y))
			end
		else
			if DataAPI.invrefpool(y) !== nothing
				if f2 != identity
					v = map(f2, y)
				else
					v = y
				end
				dict = DataAPI.invrefpool(v)
				if hasproperty(dict, :invpool)
					vtype = valtype(dict.invpool)
				else
					vtype = valtype(dict)
				end
				res =  Vector{Union{Missing, vtype}}(undef, length(x))
				_rev_map_invrefpool!(res, dict, x, f1)
				new{typeof(identity), typeof(identity), Union{Missing, vtype}, Union{vtype, Missing}, Union{vtype, Missing}, typeof(res), typeof(res)}(res, DataAPI.refarray(v), identity, identity, length(x), length(y))
			elseif DataAPI.invrefpool(x) !== nothing
				if f1 != identity
					v = map(f1, x)
				else
					v = x
				end
				dict = DataAPI.invrefpool(v)
				if hasproperty(dict, :invpool)
					vtype = valtype(dict.invpool)
				else
					vtype = valtype(dict)
				end
				res = Vector{Union{Missing, vtype}}(undef, length(y))
				_rev_map_invrefpool!(res, dict, y, f2)
				new{typeof(identity), typeof(identity), Union{Missing ,vtype}, Union{vtype, Missing}, Union{vtype, Missing}, typeof(res), typeof(res)}(DataAPI.refarray(v), res, identity, identity, length(x), length(y))
			else
				new{F1, F2, promote_type(Core.Compiler.return_type(f1, (eltype(x), )), Core.Compiler.return_type(f2, (eltype(y), ))), eltype(x), eltype(y), typeof(x), typeof(y)}(x, y, f1, f2, length(x), length(y))
			end
		end
	end
	#
	# Cat2Vec(x,y,f1::F1,f2::F2) where {F1, F2}= new{F1, F2, promote_type(Core.Compiler.return_type(f1, (eltype(x), )), Core.Compiler.return_type(f2, (eltype(y), ))), eltype(x), eltype(y), typeof(x), typeof(y)}(x, y, f1, f2, length(x), length(y))
end

function _rev_map_invrefpool!(res, dict, y, f)
	Threads.@threads for i in 1:length(res)
		res[i] = get(dict, DataAPI.unwrap(f(y[i])), missing)

	end
end

function __getindex(v::Cat2Vec{F1, F2, CT, T,S, A, B}, i::Int, f1, f2 )::CT where {F1, F2, CT, T, S, A, B}
    if i <= v.len1
		f1(v.vec1[i]::T)
	else
		f2(v.vec2[i-v.len1]::S)
	end
end

function Base.getindex(v::Cat2Vec{F1, F2, CT, T, S, A, B}, i::Int)::CT where {F1, F2, CT, T, S, A, B}
	__getindex(v, i, v.f1, v.f2)
end

Base.IndexStyle(::Type{<:Cat2Vec}) = Base.IndexLinear()
Base.size(v::Cat2Vec) = (length(v),)
Base.length(v::Cat2Vec) = v.len1 + v.len2
Base.eltype(v::Cat2Vec{F1, F2, CT, T, S, A, B}) where {F1, F2, CT, T, S, A, B} = CT


text:
module HijriConverter
# Main module of the hijri-converter package.
include("ummalqura.jl")
include("helpers.jl")
include("locales.jl")
using .ummalqura
using .helpers
using .locales
using Parameters
using Dates
import Base: isless, isgreater, isequal, !=, <=, >=
import Dates: year, month, day

export Hijri, hijri_today, to_hijri, to_gregorian
export isless, isgreater, iseequal, !=, <=, >=

#-----------------PART-I--------------------
# ----------------HIJRI-------------------

@doc  """
      A Hijri object represents a date (year, month and day) in Hijri
      calendar.
      """
@with_kw struct Hijri
	year::Int64
	month::Int64
	day::Int64
	validate::Bool
end

Hijri(year, month, day) = Hijri(year, month, day, true)

@doc     """Return month name.

        :param language: Language tag for localized month name. Full locale name can
            be used, e.g. ``en-US`` or ``en_US.UTF-8``. Supported languages are
            ``en``, ``ar`` and ``bn``.
        :type language: str
        """
function month_name(date::Hijri, language::String = "en")
        return locales.month_name(locales.get_locale(language), month(date))
end

@doc    """Return month’s index in ummalqura month starts"""
function _month_index(date::Hijri)
        prior_months = ((year(date) - 1) * 12) + month(date) - 1
        index = prior_months - ummalqura.HIJRI_OFFSET
        return index
end

@doc    """Check date values if within valid range."""
function _check_date_(date::Hijri)
        # check year
        min_year, max_year = [d[1] for d in ummalqura.HIJRI_RANGE]
        if min_year > year(date) ||  year(date) > max_year
            throw(DomainError("date out of range"))
        end
        # check month
        if 1 > month(date) || month(date) > 12
            throw(DomainError("month must be in 1..12"))
        end
        # check day
        month_len = month_length(date)
        if 1 > day(date) || day(date) > month_len
            throw(DomainError(string("day must be in 1..", month_len, " for month")))
        end
end

@doc    """Return calendar era notation.

        :param language: Language tag for localized notation. Full locale name can
            be used, e.g. ``en-US`` or ``en_US.UTF-8``. Supported languages are
            ``en``, ``ar`` and ``bn``.
        :type language: str
        """
function notation(date::Hijri, language::String = "en")

        return locales.get_locale(language).notation
end

@doc    """Return Gregorian object for the corresponding Hijri date.

        :rtype: Gregorian
        """
function to_gregorian(date::Hijri)
	_check_date_(date)
        jdn = to_julian(date)
        n = helpers.jdn_to_ordinal(jdn)
        return fromordinal(n)
end

@doc    """Return corresponding Julian day number (JDN)."""
function to_julian(date::Hijri)
        month_starts = ummalqura.MONTH_STARTS
        index = _month_index(date)
	rjd = month_starts[index + 1] + day(date) - 1
        jdn = helpers.rjd_to_jdn(rjd)
        return jdn
end

@doc    """Return day of week, where Monday is 0 and Sunday is 6."""
function weekday(date::Hijri)
        jdn = to_julian(date)
        return Int(jdn % 7)
end

@doc    """Return day of week, where Monday is 1 and Sunday is 7."""
function isoweekday(date::Hijri)
        jdn = to_julian(date)
        return Int(jdn % 7) + 1
end

@doc     """Return number of days in month."""
function month_length(date::Hijri)
        month_starts = ummalqura.MONTH_STARTS
        index = _month_index(date)
        len = [month_starts[index] - month_starts[mod1(index - 1, end)] for index in 1:length(month_starts)][1]
        return len
end

#-----------------PART-II--------------------
# ----------------GREGORIAN-------------------

@doc    """Check date values if within valid range."""
function _check_date_(date::Union{Date, DateTime})
        # check year
        min_year, max_year = [d[1] for d in ummalqura.GREGORIAN_RANGE]
        if min_year > year(date) ||  year(date) > max_year
            throw(error("date out of range")) # replace with overflow error
        end
        # check month
        if 1 > month(date) || month(date) > 12
            throw(ValueError("month must be in 1..12"))
        end
        # check day
        month_len = month_length(date) #TODO: fix this
        if 1 > day(date) || day(date) > month_len
            throw(ValueError(string("day must be in 1..", month_len, " for month")))
        end
end

@doc    """Return Hijri object for the corresponding Gregorian date.

        :rtype: Hijri
        """
function to_hijri(date::Union{Date, DateTime})
        _check_range(date)
        jdn = Int(round(datetime2julian(date)))
        rjd = helpers.jdn_to_rjd(jdn)
        month_starts = ummalqura.MONTH_STARTS
        index = helpers.bisect_left(month_starts, rjd) - 1
        months = index + ummalqura.HIJRI_OFFSET - 1

        years = floor(months / 12)
        year = years + 1
        month = months - (years * 12) + 1
        day = rjd - month_starts[index] + 2 #TODO: I don't why this is "+2" rather than "+1"

        @assert day > 0 && day < 32
        @assert month > 0 && month < 13
        return Hijri(year, month, day, false)
end

@doc    """Return proleptic Gregorian ordinal for the year, month and day.
        January 1 of year 1 is day 1.  Only the year, month and day values
        contribute to the result.
        """
function toordinal(date::Union{Date, DateTime})
        return helpers._ymd2ord(year(date), month(date), day(date))
end

@doc    """Construct a date from a proleptic Gregorian ordinal.
        January 1 of year 1 is day 1.  Only the year, month and day are
        non-zero in the result.
        """
function fromordinal(n)
        y, m, d = helpers._ord2ymd(n)
        return DateTime(y, m, d)
end

@doc    """Check if Gregorian date is within valid range."""
function _check_range(date::Union{Date, DateTime})
        min_date, max_date = ummalqura.GREGORIAN_RANGE
        if min_date > (year(date), month(date), day(date)) || (year(date), month(date), day(date)) > max_date
            throw(error("date out of range")) # replace with OverFlowError
        end
end

@doc    """Return calendar era notation.

        :param language: Language tag for localized notation. Full locale name can
            be used, e.g. ``en-US`` or ``en_US.UTF-8``. Supported languages are
            ``en``, ``ar`` and ``bn``.
        :type language: str
        """
function notation(date::Union{Date,DateTime}, language::String = "en")

        return locales.get_locale(language).gregorian_notation
end

@doc    """Return month name.

        :param language: Language tag for localized month name. Full locale name can
            be used, e.g. ``en-US`` or ``en_US.UTF-8``. Supported languages are
            ``en``, ``ar`` and ``bn``.
        :type language: str
        """
function month_name(date::Union{Date, DateTime}, language::String = "en")

        return gregorian_month_name(locales.get_locale(language), month(date))
end

#-----------------PART-III--------------------
# ----------------COMMON METHODS-------------------

@doc    """Construct Hijri object from an ISO formatted Hijri date
        ``YYYY-MM-DD``.

        :param date_string: Hijri date in ISO format ``YYYY-MM-DD``.
        :type date_string: str
        :rtype: Hijri
        """
function fromisoformat(date_string::String)
        year = parse(Int, date_string[1:4])
        month = parse(Int, date_string[6:7])
        day = parse(Int, date_string[9:10])
        return Hijri(year, month, day)
end

@doc    """Construct Hijri object from today's date.

        :rtype: Hijri
        """
function hijri_today() # Gregorian struct is incomplete
        return to_hijri(DateTime(today()))
end

@doc     """Return year as an integer."""
function year(date::Hijri)
        return date.year
end

@doc     """Return month as an integer."""
function month(date::Hijri)
        return date.month
end

@doc     """Return day as an integer."""
function day(date::Hijri)
        return date.day
end

@doc     """Return date in ISO format i.e. ``YYYY-MM-DD``."""
function isoformat(date::T)  where {T<:Union{Hijri, Date, DateTime}}
	year = lpad(year(date), 4, "0")
	month = lpad(month(date), 2, "0")
	day = lpad(day(date), 2, "0")

	return string(year, "-", month, "-", day)
end

@doc    """Return date in day-month-year format (``DD/MM/YYYY`` by default).

        :param separator: String that separates the day, month, and year values.
        :type separator: str
        :param padding: Whether to add a leading zero as a padding character to fill
            day and month values when less than 10.
        :type padding: bool
        """
function dmyformat(date::T; separator::String = "/", padding::Bool = true) where {T <: Union{Hijri, Date, DateTime}}

        day = padding ? lpad(day(date), 2, "0") : day(date)
        month = padding ? lpad(month(date), 2, "0") : month(date)
        return string(day, separator, month, separator, year(date))
end

@doc    """Return day name.

        :param language: Language tag for localized day name. Full locale name can
            be used, e.g. ``en-US`` or ``en_US.UTF-8``. Supported languages are
            ``en``, ``ar`` and ``bn``.
        :type language: str
        """
function day_name(date::T, language::String = "en") where {T <: Union{Hijri, Date, DateTime}} #TODO: implement isoweekday for gregorian

        return locales.day_name(locales.get_locale(language), isoweekday(date))
end

function datetuple(date::T) where {T <: Union{Hijri, Date, DateTime}}
	return (year(date), month(date), day(date))
end

# NUMERIC COMPARISONS
function isless(a::Hijri, b::Hijri)
	tuple_a = datetuple(a)
	tuple_b = datetuple(b)

	return tuple_a < tuple_b
end

function isgreater(a::Hijri, b::Hijri)
	tuple_a = datetuple(a)
	tuple_b = datetuple(b)

	return tuple_a > tuple_b
end

function isequal(a::Hijri, b::Hijri)
	tuple_a = datetuple(a)
	tuple_b = datetuple(b)

	return tuple_a == tuple_b
end

function !=(a::Hijri, b::Hijri)
	tuple_a = datetuple(a)
	tuple_b = datetuple(b)

	return tuple_a != tuple_b
end

function <=(a::Hijri, b::Hijri)
	tuple_a = datetuple(a)
	tuple_b = datetuple(b)

	return tuple_a <= tuple_b
end

function >=(a::Hijri, b::Hijri)
	tuple_a = datetuple(a)
	tuple_b = datetuple(b)

	return tuple_a >= tuple_b
end
end

ext: jl

max_stars_repo_name: ceferisbarov/hijri-converter

max_issues_repo_name: ceferisbarov/hijri-converter


text:
module ILDG_format
    using CLIME_jll
    using EzXML

    import ..IOmodule:IOFormat
    import ..Gaugefields:GaugeFields,SU3GaugeFields,SU2GaugeFields,set_wing!

    struct LIME_header
        doc::EzXML.Document
        function LIME_header(L,field,version,precision)
            doc = XMLDocument()
            elm = ElementNode("ildgFormat")
            addelement!(elm, "version", "$version")
            addelement!(elm, "field", "$field")
            addelement!(elm, "precision", "$precision")
            addelement!(elm, "lx", "$(L[1])")
            addelement!(elm, "ly", "$(L[2])")
            addelement!(elm, "lz", "$(L[3])")
            addelement!(elm, "lt", "$(L[4])")
            setroot!(doc, elm)
            return new(doc)
        end
    end

    struct ILDG <: IOFormat
        header::Array{Dict,1}
        filename::String
        ILDG(filename) = new(read_header(filename),filename)
    end

    function Base.length(ildg::ILDG)
        return length(ildg.header)
    end

    function Base.getindex(ildg::ILDG,i)
        return ildg.header[i]
    end

    function save_binarydata(U,filename)

        NX = U[1].NX
        NY = U[1].NY
        NZ = U[1].NZ
        NT = U[1].NT
        NC = U[1].NC

        #li = LIME_header((NX,NY,NZ,NT),"su3gauge",1,64)
        #print(li.doc)
        #write("test.xml", li.doc)

        fp = open("testbin.dat","w")
        i = 0
        i = 0
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    for ix=1:NX
                        for μ=1:4
                            for ic2 = 1:NC
                                for ic1 = 1:NC
                                    i+= 1
                                    #rvalue = read(fp, floattype)
                                    rvalue = real(U[μ][ic2,ic1,ix,iy,iz,it])
                                    ivalue = imag(U[μ][ic2,ic1,ix,iy,iz,it])
                                    write(fp,hton(rvalue))
                                    write(fp,hton(ivalue))
                                end
                            end
                        end
                    end
                end
            end
        end
        close(fp)

        fp = open("filelist.dat","w")
        #println(fp,"test.xml ","ildg-format")
        println(fp,"testbin.dat ","ildg-binary-data")
        close(fp)

        lime_pack() do exe
            run(`$exe filelist.dat $filename`)
        end

        return

    end

    function load_binarydata!(U,NX,NY,NZ,NT,NC,filename,precision)
        if precision == 32
            floattype = Float32
        else
            floattype = Float64
        end
        fp = open(filename,"r")
        totalnum = NX*NY*NZ*NT*NC*NC*2*4
        #println(countlines(filename))

        i = 0
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    for ix=1:NX
                        for μ=1:4
                            for ic2 = 1:NC
                                for ic1 = 1:NC
                                    i+= 1
                                    #rvalue = read(fp, floattype)
                                    rvalue = ntoh(read(fp, floattype))
                                    #rvalue = ntoh(read(fp, ComplexF64))
                                    i+= 1
                                    ivalue = ntoh(read(fp, floattype))
                                    #ivalue = read(fp, floattype)
                                    U[μ][ic2,ic1,ix,iy,iz,it] = rvalue + im*ivalue
                                    #U[μ][ic1,ic2,ix,iy,iz,it] = rvalue# + im*ivalue
                                    #println(i,"/",totalnum,"= ",i/totalnum)
                                    #println(rvalue,"\t",ivalue)
                                end
                            end
                        end
                    end
                end
            end
        end

        set_wing!(U)

        close(fp)
    end

    function load_gaugefield!(U,i,ildg::ILDG,L,NC;NDW = 1)
        NX = L[1]
        NY = L[2]
        NZ = L[3]
        NT = L[4]
        data = ildg[i]
        filename = ildg.filename
        @assert U[1].NX == NX "NX mismatch"
        @assert U[1].NY == NY "NY mismatch"
        @assert U[1].NZ == NZ "NZ mismatch"
        @assert U[1].NT == NT "NT mismatch"
        @assert U[1].NC == NC "NC mismatch"

        message_no = data["message_no"]
        reccord_no = data["reccord_no"]
        if haskey(data,"precision")
            precision = data["precision"]
        else
            precision = 64
        end

        lime_extract_record() do exe
            run(`$exe $filename $message_no $reccord_no tempconf.dat`)
        end

        load_binarydata!(U,NX,NY,NZ,NT,NC,"tempconf.dat",precision)

        return
    end

    function load_gaugefield(i,ildg::ILDG,L,NC;NDW = 1)
        NX = L[1]
        NY = L[2]
        NZ = L[3]
        NT = L[4]
        data = ildg[i]
        filename = ildg.filename

        if NC == 3
            U = Array{SU3GaugeFields,1}(undef,4)
        elseif NC == 2
            U = Array{SU2GaugeFields,1}(undef,4)
        end

        for μ=1:4
            U[μ] = GaugeFields(NC,NDW,NX,NY,NZ,NT)
        end

        load_gaugefield!(U,i,ildg::ILDG,L,NC;NDW = 1)
        return U

    end

    function load_gaugefield(i,ildg::ILDG;NDW = 1)
        #@assert length(ildg) != 0 "the header file is not found"
        data = ildg[i]
        filename = ildg.filename
        if haskey(data,"L")
            L = data["L"]
            NX = L[1]
            NY = L[2]
            NZ = L[3]
            NT = L[4]
        else
            error("header file is not found. Please put lattice size")
        end
        if haskey(data,"L")
            error("header file is not found. Please put NC")
            NC = data["NC"]
        end
        load_gaugefield(i,ildg::ILDG,L,NC;NDW = NDW)

    end

    function read_header(filename)
        contents_data = ""
        header = nothing
        lime_contents() do exe
            contents_data = read(`$exe $filename`,String)
            #println(contents_data)
            contents_data = split(string(contents_data),"\n")
            #println(contents_data)
            header = extract_info(contents_data)

        end
        #println(header)
        return header
    end

    function extract_info(contents_data)
        #println(typeof(contents_data))
        i = 0
        message_no = 0
        reccord_no = 0
        datatype = ""
        header = Dict[]
        NX = 0
        NY = 0
        NZ = 0
        NT = 0
        NC = 3
        precision = 32
        headerfound = false
        headerdic = Dict()
        for data in contents_data
            u = split(data)
            #println(u)
            if length(u) ≥ 2
                if u[1] == "Message:"
                    message_no = parse(Int64,u[2])
                    #println(message_no)
                end
                if u[1] == "Record:"
                    reccord_no = parse(Int64,u[2])
                end
                if u[1] == "Type:"
                    datatype = u[2]

                end
                if u[1] == "Data:"
                    #println("message_no = $(message_no)")
                    #println("reccord_no = $reccord_no")
                    #println("datatype  = $datatype ")
                    if datatype == "ildg-format"
                        headerdic = Dict()
                        #println(data[2:end])
                        ist =  findfirst('\"',data)
                        ien =  findlast('\"',data)

                        #ien =  findlast(""\",data)
                        doc = parsexml(data[ist+1:ien-1])
                        #elm_lx = ElementNode("lx")
                        #println("lx = ",elm_lx.content)
                        ildgFormat  = root(doc)
                        #systemdata = elements(ildgFormat)
                        #println(systemdata["version"])

                        for d in eachelement(ildgFormat)
                            if d.name == "lx"
                                NX = parse(Int64,d.content)
                            elseif d.name == "ly"
                                NY = parse(Int64,d.content)
                            elseif d.name == "lz"
                                NZ = parse(Int64,d.content)
                            elseif    d.name == "lt"
                                NT = parse(Int64,d.content)
                            elseif d.name == "field"
                                gauge = d.content
                                if findfirst("su3",gauge) != nothing
                                    NC = 3
                                elseif findfirst("su2",gauge) != nothing
                                    NC = 2
                                else
                                    error("not supported. gauge is ",gauge)
                                end
                            elseif d.name == "precision"
                                precision = parse(Int64,d.content)
                            end

                        end
                        headerdic["L"] = (NX,NY,NZ,NT)
                        headerdic["NC"] = NC
                        headerdic["precision"] = precision
                        headerdic["headertype"] = "ildg-format"

                        headerfound = true

                    elseif datatype == "scidac-private-file-xml"
                        headerdic = Dict()
                        #println(data[2:end])

                        ist =  findfirst('\"',data)
                        ien =  findlast('\"',data)

                        doc = parsexml(data[ist+1:ien-1])

                        scidacFile  = root(doc)
                        #systemdata = elements(ildgFormat)
                        #println(systemdata["version"])

                        for d in eachelement(scidacFile)
                            #println(d)
                            if d.name == "dims"
                                L = parse.(Int64,split(d.content))
                                NX = L[1]
                                NY = L[2]
                                NZ = L[3]
                                NT = L[4]
                                #=
                            elseif d.name == "colors"
                                NC = parse(Int64,d.content)
                                println(NC)
                            elseif d.name == "precision"
                                if d.content == "F"
                                    precision = 32
                                elseif d.content == "D"
                                    precision = 64
                                end
                                =#
                            end

                        end
                        headerdic["L"] = (NX,NY,NZ,NT)
                        headerdic["NC"] = NC
                        headerdic["headertype"] = "scidac-private-file-xml"
                        headerdic["precision"] = precision

                        headerfound = true

                    end

                    if datatype == "ildg-binary-data" #&& headerfound
                        #println("message_no = $(message_no)")
                        #println("reccord_no = $reccord_no")
                        headerdic["message_no"] = message_no
                        headerdic["reccord_no"] = reccord_no
                        push!(header,headerdic)
                        headerfound = false

                    end

                end

            end

        end

        return header
    end

    function test()
        contents_data = ""
        header = nothing
        lime_contents() do exe
            contents_data = read(`$exe $(ARGS[1])`,String)
            #println(contents_data)
            contents_data = split(string(contents_data),"\n")
            #println(contents_data)
            header = extract_info(contents_data)

        end
        println(header)

        #println(contents_data)

        #println(text)
        #=
        lime_extract_record() do exe
            run(`$exe $(ARGS[1]) 1 1 out11.dat`)
            run(`$exe $(ARGS[1]) 1 2 out12.dat`)
            run(`$exe $(ARGS[1]) 1 3 out13.dat`)
            run(`$exe $(ARGS[1]) 2 1 out21.dat`)
            run(`$exe $(ARGS[1]) 2 2 out22.dat`)
            run(`$exe $(ARGS[1]) 3 3 out23.dat`)
        end
        =#
    end
end


text:
"""
    coil_combination_method(x)
    coil_combination_method!(x, val)

The coil combination method used for image acquisiton. Almost all fMRI studies
using phased-array coils use root-sum-of-squares (rSOS) combination, but other
methods exist. The image reconstruction is changed by the coil combination method
(as for the matrix coil mode above), so anything non-standard should be reported.
"""
@defprop CoilCombinationMethod{:coil_combination_method}::String

"""
    gradient_set_type(x)
    gradient_set_type!(x, val)

The gradient set type. It should be possible to infer the gradient coil from the
scanner model. If not, e.g. because of a custom upgrade or use of a gradient insert
set, then the specifications of the actual gradient coil should be reported independently.
"""
@defprop GradientSetType{:gradient_set_type}::String

"""
    magnetic_field_strength(x)
    magnetic_field_strength!(x, val)

The nominal field strength of MR magnet in Tesla.
"""
@defprop MagneticFieldStrength{:magnetic_field_strength}::F64Tesla

"""
    matrix_coil_mode(x)
    matrix_coil_mode!(x, val)

Returns the matrix coil mode. A method for reducing the number of independent
channels by combining in analog the signals from multiple coil elements. There
are typically different default modes when using un-accelerated or accelerated
(e.g. GRAPPA, SENSE) imaging.
"""
@defprop MatrixCoilMode{:matrix_coil_mode}::String

"""
    mr_transmit_coil_sequence(x)
    mr_transmit_coil_sequence!(x, val)

A sequence that provides information about the transmit coil used. This is a relevant
field if a non-standard transmit coil is used.
"""
@defprop MRTransmitCoilSequence{:mr_transmit_coil_sequence}::String

"""
    receive_coil_name(x)
    receive_coil_name!(x, val)

The information describing the receiver coil. Corresponds to DICOM Tag
Receive Coil Name, although not all vendors populate that DICOM Tag, in which
case this field can be derived from an appropriate private DICOM field.
"""
@defprop ReceiveCoilName{:receive_coil_name}::String

"""
    receive_coil_active_elements(x)
    receive_coil_active_elements!(x, val)

Information describing the active/selected elements of the receiver coil. This
doesn’t correspond to a tag in the DICOM ontology. The vendor-defined terminology
for active coil elements can go in this field. As an example, for
Siemens, coil channels are typically not activated/selected individually, but
rather in pre-defined selectable "groups" of individual channels, and the list
of the "groups" of elements that are active/selected in any given scan
populates the Coil String entry in Siemens’ private DICOM fields (e.g., HEA;HEP
for the Siemens standard 32 ch coil when both the anterior and posterior groups
are activated). This is a flexible field that can be used as most appropriate
for a given vendor and coil to define the "active" coil elements. Since
individual scans can sometimes not have the intended coil elements selected, it
is preferable for this field to be populated directly from the DICOM for each
individual scan, so that it can be used as a mechanism for checking that a
given scan was collected with the intended coil elements selected.
"""
@defprop ReceiveCoilActiveElements{:receive_coil_active_elements}::String

"""
    flip_angle(x)
    flip_angle!(x, val)

Returns the flip angle for the acquisition in degrees.
"""
@defprop FlipAngle{:flip_angle}::IntDegree

"Returns the multiband factor, for multiband acquisitions."
@defprop MultibandAccelerationFactor{:multiband_acceleration_factor}::String

"""
    negative_contrast(x)
    negative_contrast!(x, val)

Specifies whether increasing voxel intensity (within sample voxels)
denotes a decreased value with respect to the contrast suffix. This is commonly
the case when Cerebral Blood Volume is estimated via usage of a contrast agent
in conjunction with a T2* weighted acquisition protocol.
"""
@defprop NegativeContrast{:negative_contrast}::Bool

"""
    contrast_bolus_ingredient(x)
    contrast_bolus_ingredient!(x, val)

Return active ingredient of constrast agent. See [`ContrastIngredient`](@ref) for more details.
"""
@defprop ContrastBolusIngredient{:contrast_bolus_ingredient}::ContrastIngredient

"""
    nvol_discarded_by_scanner(x)
    nvol_discarded_by_scanner!(x, val)

The number of volumes ("dummy scans") discarded by the scanner (as opposed
to those discarded by the user post hoc) before saving the imaging file. For example,
a sequence that automatically discards the first 4 volumes before saving would
have this field as 4. A sequence that doesn't discard dummy scans would have
this set to 0. Please note that the onsets recorded in the _event.tsv file
should always refer to the beginning of the acquisition of the first volume in
the corresponding imaging file - independent of the value of
number_of_volumes_discarded_by_scanner field.
"""
@defprop NumberOfVolumesDiscardedByScanner{:nvol_discarded_by_scanner}::Int

"""
    nvol_discarded_by_user(x)
    nvol_discarded_by_user!(x, val)

Returns the  number of volumes ("dummy scans") discarded by the user before
including the file in the dataset. If possible, including all of the volumes is
strongly recommended. Please note that the onsets recorded in the _event.tsv file
should always refer to the beginning of the acquisition of the first volume in
the corresponding imaging file - independent of the value of
number_of_volumes_discarded_by_user field.
"""
@defprop NumberOfVolumesDiscardedByUser{:nvol_discarded_by_user}::Int


text:
@symbol_func function K_KF(cur_reactor::AbstractReactor)
  cur_K = cur_reactor.max_q_95

  cur_K /= 5

  cur_K /= cur_reactor.epsilon ^ 2

  cur_K /= _f_q(cur_reactor)

  cur_K
end

@symbol_func function _f_q(cur_reactor::AbstractReactor)
  cur_epsilon = cur_reactor.epsilon

  cur_kappa = cur_reactor.kappa_95

  cur_delta = cur_reactor.delta_95

  cur_left_term = 1.17 - 0.65 * cur_epsilon

  cur_left_term /= 2 * ( 1 - cur_epsilon ^ 2 ) ^ 2

  cur_right_term = 1 + 2 * cur_delta ^ 2 - 1.2 * cur_delta ^ 3

  cur_right_term = 1 + cur_kappa ^ 2 * cur_right_term

  cur_f_q = cur_left_term

  cur_f_q *= cur_right_term

  cur_f_q
end


text:
int_rules_7_1_6 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 7.1.6*Miscellaneous*inverse*hyperbolic*sine =#
    @apply_utils Antiderivative((~(a') + ~(b') * asinh(~c + ~(d') * ~x)) ^ ~(n'), ~x) => (1 / ~d) * Subst(Antiderivative((~a + ~b * asinh(~x)) ^ ~n, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(a') + ~(b') * asinh(~c + ~(d') * ~x)) ^ ~(n'), ~x) => (1 / ~d) * Subst(Antiderivative(((~d * ~e - ~c * ~f) / ~d + (~f * ~x) / ~d) ^ ~m * (~a + ~b * asinh(~x)) ^ ~n, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x)
    @apply_utils Antiderivative((~(A') + ~(B') * ~x + ~(C') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asinh(~c + ~(d') * ~x)) ^ ~(n'), ~x) => (1 / ~d) * Subst(Antiderivative((~C / (~d) ^ 2 + (~C / (~d) ^ 2) * (~x) ^ 2) ^ ~p * (~a + ~b * asinh(~x)) ^ ~n, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~A, ~B, ~C, ~n, ~p], ~x) && (EqQ(~B * (1 + (~c) ^ 2) - 2 * ~A * ~c * ~d, 0) && EqQ(2 * ~c * ~C - ~B * ~d, 0))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(A') + ~(B') * ~x + ~(C') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asinh(~c + ~(d') * ~x)) ^ ~(n'), ~x) => (1 / ~d) * Subst(Antiderivative(((~d * ~e - ~c * ~f) / ~d + (~f * ~x) / ~d) ^ ~m * (~C / (~d) ^ 2 + (~C / (~d) ^ 2) * (~x) ^ 2) ^ ~p * (~a + ~b * asinh(~x)) ^ ~n, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~C, ~m, ~n, ~p], ~x) && (EqQ(~B * (1 + (~c) ^ 2) - 2 * ~A * ~c * ~d, 0) && EqQ(2 * ~c * ~C - ~B * ~d, 0))
    @apply_utils Antiderivative(sqrt(~(a') + ~(b') * asinh(~c + ~(d') * (~x) ^ 2)), ~x) => (~x * sqrt(~a + ~b * asinh(~c + ~d * (~x) ^ 2)) - (sqrt(pi) * ~x * (cosh(~a / (2 * ~b)) - ~c * sinh(~a / (2 * ~b))) * FresnelC(sqrt(-(~c) / (pi * ~b)) * sqrt(~a + ~b * asinh(~c + ~d * (~x) ^ 2)))) / (sqrt(-(~c / ~b)) * (cosh(asinh(~c + ~d * (~x) ^ 2) / 2) + ~c * sinh(asinh(~c + ~d * (~x) ^ 2) / 2)))) + (sqrt(pi) * ~x * (cosh(~a / (2 * ~b)) + ~c * sinh(~a / (2 * ~b))) * FresnelS(sqrt(-(~c) / (pi * ~b)) * sqrt(~a + ~b * asinh(~c + ~d * (~x) ^ 2)))) / (sqrt(-(~c / ~b)) * (cosh(asinh(~c + ~d * (~x) ^ 2) / 2) + ~c * sinh(asinh(~c + ~d * (~x) ^ 2) / 2))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~c) ^ 2, -1)
    @apply_utils Antiderivative((~(a') + ~(b') * asinh(~c + ~(d') * (~x) ^ 2)) ^ ~n, ~x) => (~x * (~a + ~b * asinh(~c + ~d * (~x) ^ 2)) ^ ~n - (2 * ~b * ~n * sqrt(2 * ~c * ~d * (~x) ^ 2 + (~d) ^ 2 * (~x) ^ 4) * (~a + ~b * asinh(~c + ~d * (~x) ^ 2)) ^ (~n - 1)) / (~d * ~x)) + 4 * (~b) ^ 2 * ~n * (~n - 1) * Antiderivative((~a + ~b * asinh(~c + ~d * (~x) ^ 2)) ^ (~n - 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ((~c) ^ 2, -1) && GtQ(~n, 1))
    @apply_utils Antiderivative(1 / (~(a') + ~(b') * asinh(~c + ~(d') * (~x) ^ 2)), ~x) => (~x * (~c * cosh(~a / (2 * ~b)) - sinh(~a / (2 * ~b))) * CoshIntegral((~a + ~b * asinh(~c + ~d * (~x) ^ 2)) / (2 * ~b))) / (2 * ~b * (cosh(asinh(~c + ~d * (~x) ^ 2) / 2) + ~c * sinh((1 / 2) * asinh(~c + ~d * (~x) ^ 2)))) + (~x * (cosh(~a / (2 * ~b)) - ~c * sinh(~a / (2 * ~b))) * SinhIntegral((~a + ~b * asinh(~c + ~d * (~x) ^ 2)) / (2 * ~b))) / (2 * ~b * (cosh(asinh(~c + ~d * (~x) ^ 2) / 2) + ~c * sinh((1 / 2) * asinh(~c + ~d * (~x) ^ 2)))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~c) ^ 2, -1)
    @apply_utils Antiderivative(1 / sqrt(~(a') + ~(b') * asinh(~c + ~(d') * (~x) ^ 2)), ~x) => ((~c + 1) * sqrt(pi / 2) * ~x * (cosh(~a / (2 * ~b)) - sinh(~a / (2 * ~b))) * erfi(sqrt(~a + ~b * asinh(~c + ~d * (~x) ^ 2)) / sqrt(2 * ~b))) / (2 * sqrt(~b) * (cosh(asinh(~c + ~d * (~x) ^ 2) / 2) + ~c * sinh(asinh(~c + ~d * (~x) ^ 2) / 2))) + ((~c - 1) * sqrt(pi / 2) * ~x * (cosh(~a / (2 * ~b)) + sinh(~a / (2 * ~b))) * erf(sqrt(~a + ~b * asinh(~c + ~d * (~x) ^ 2)) / sqrt(2 * ~b))) / (2 * sqrt(~b) * (cosh(asinh(~c + ~d * (~x) ^ 2) / 2) + ~c * sinh(asinh(~c + ~d * (~x) ^ 2) / 2))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~c) ^ 2, -1)
    @apply_utils Antiderivative(1 / (~(a') + ~(b') * asinh(~c + ~(d') * (~x) ^ 2)) ^ (3 / 2), ~x) => (-(sqrt(2 * ~c * ~d * (~x) ^ 2 + (~d) ^ 2 * (~x) ^ 4)) / (~b * ~d * ~x * sqrt(~a + ~b * asinh(~c + ~d * (~x) ^ 2))) - ((-(~c) / ~b) ^ (3 / 2) * sqrt(pi) * ~x * (cosh(~a / (2 * ~b)) - ~c * sinh(~a / (2 * ~b))) * FresnelC(sqrt(-(~c) / (pi * ~b)) * sqrt(~a + ~b * asinh(~c + ~d * (~x) ^ 2)))) / (cosh(asinh(~c + ~d * (~x) ^ 2) / 2) + ~c * sinh(asinh(~c + ~d * (~x) ^ 2) / 2))) + ((-(~c) / ~b) ^ (3 / 2) * sqrt(pi) * ~x * (cosh(~a / (2 * ~b)) + ~c * sinh(~a / (2 * ~b))) * FresnelS(sqrt(-(~c) / (pi * ~b)) * sqrt(~a + ~b * asinh(~c + ~d * (~x) ^ 2)))) / (cosh(asinh(~c + ~d * (~x) ^ 2) / 2) + ~c * sinh(asinh(~c + ~d * (~x) ^ 2) / 2)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~c) ^ 2, -1)
    @apply_utils Antiderivative(1 / (~(a') + ~(b') * asinh(~c + ~(d') * (~x) ^ 2)) ^ 2, ~x) => -(sqrt(2 * ~c * ~d * (~x) ^ 2 + (~d) ^ 2 * (~x) ^ 4)) / (2 * ~b * ~d * ~x * (~a + ~b * asinh(~c + ~d * (~x) ^ 2))) + (~x * (cosh(~a / (2 * ~b)) - ~c * sinh(~a / (2 * ~b))) * CoshIntegral((~a + ~b * asinh(~c + ~d * (~x) ^ 2)) / (2 * ~b))) / (4 * (~b) ^ 2 * (cosh(asinh(~c + ~d * (~x) ^ 2) / 2) + ~c * sinh(asinh(~c + ~d * (~x) ^ 2) / 2))) + (~x * (~c * cosh(~a / (2 * ~b)) - sinh(~a / (2 * ~b))) * SinhIntegral((~a + ~b * asinh(~c + ~d * (~x) ^ 2)) / (2 * ~b))) / (4 * (~b) ^ 2 * (cosh(asinh(~c + ~d * (~x) ^ 2) / 2) + ~c * sinh(asinh(~c + ~d * (~x) ^ 2) / 2))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~c) ^ 2, -1)
    @apply_utils Antiderivative((~(a') + ~(b') * asinh(~c + ~(d') * (~x) ^ 2)) ^ ~n, ~x) => (-(~x) * (~a + ~b * asinh(~c + ~d * (~x) ^ 2)) ^ (~n + 2)) / (4 * (~b) ^ 2 * (~n + 1) * (~n + 2)) + (sqrt(2 * ~c * ~d * (~x) ^ 2 + (~d) ^ 2 * (~x) ^ 4) * (~a + ~b * asinh(~c + ~d * (~x) ^ 2)) ^ (~n + 1)) / (2 * ~b * ~d * (~n + 1) * ~x) + (1 / (4 * (~b) ^ 2 * (~n + 1) * (~n + 2))) * Antiderivative((~a + ~b * asinh(~c + ~d * (~x) ^ 2)) ^ (~n + 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ((~c) ^ 2, -1) && (LtQ(~n, -1) && NeQ(~n, -2)))
    @apply_utils Antiderivative(asinh(~(a') * (~x) ^ ~p) ^ ~(n') / ~x, ~x) => (1 / ~p) * Subst(Antiderivative((~x) ^ ~n * coth(~x), ~x), ~x, asinh(~a * (~x) ^ ~p)) <-- FreeQ([~a, ~p], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative(~(u') * asinh(~(c') / (~(a') + ~(b') * (~x) ^ ~(n'))) ^ ~(m'), ~x) => Antiderivative(~u * acsch(~a / ~c + (~b * (~x) ^ ~n) / ~c) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~m], ~x)
    @apply_utils Antiderivative(asinh(sqrt(-1 + ~(b') * (~x) ^ 2)) ^ ~(n') / sqrt(-1 + ~(b') * (~x) ^ 2), ~x) => (sqrt(~b * (~x) ^ 2) / (~b * ~x)) * Subst(Antiderivative(asinh(~x) ^ ~n / sqrt(1 + (~x) ^ 2), ~x), ~x, sqrt(-1 + ~b * (~x) ^ 2)) <-- FreeQ([~b, ~n], ~x)
    @apply_utils Antiderivative((~f) ^ (~(c') * asinh(~(a') + ~(b') * ~x) ^ ~(n')), ~x) => (1 / ~b) * Subst(Antiderivative((~f) ^ (~c * (~x) ^ ~n) * cosh(~x), ~x), ~x, asinh(~a + ~b * ~x)) <-- FreeQ([~a, ~b, ~c, ~f], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative((~x) ^ ~(m') * (~f) ^ (~(c') * asinh(~(a') + ~(b') * ~x) ^ ~(n')), ~x) => (1 / ~b) * Subst(Antiderivative((-(~a) / ~b + sinh(~x) / ~b) ^ ~m * (~f) ^ (~c * (~x) ^ ~n) * cosh(~x), ~x), ~x, asinh(~a + ~b * ~x)) <-- FreeQ([~a, ~b, ~c, ~f], ~x) && (IGtQ(~m, 0) && IGtQ(~n, 0))
    @apply_utils Antiderivative(asinh(~u), ~x) => ~x * asinh(~u) - Antiderivative(SimplifyIntegrand((~x * D(~u, ~x)) / sqrt(1 + (~u) ^ 2), ~x), ~x) <-- InverseFunctionFreeQ(~u, ~x) && Not(FunctionOfExponentialQ(~u, ~x))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~(a') + ~(b') * asinh(~u)), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * (~a + ~b * asinh(~u))) / (~d * (~m + 1)) - (~b / (~d * (~m + 1))) * Antiderivative(SimplifyIntegrand(((~c + ~d * ~x) ^ (~m + 1) * D(~u, ~x)) / sqrt(1 + (~u) ^ 2), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (NeQ(~m, -1) && (InverseFunctionFreeQ(~u, ~x) && (Not(FunctionOfQ((~c + ~d * ~x) ^ (~m + 1), ~u, ~x)) && Not(FunctionOfExponentialQ(~u, ~x)))))
    @apply_utils Antiderivative(~v * (~(a') + ~(b') * asinh(~u)), ~x) => With([w = IntHide(~v, ~x)], Dist(~a + ~b * asinh(~u), w, ~x) - ~b * Antiderivative(SimplifyIntegrand((w * D(~u, ~x)) / sqrt(1 + (~u) ^ 2), ~x), ~x) <-- InverseFunctionFreeQ(w, ~x)) <-- FreeQ([~a, ~b], ~x) && (InverseFunctionFreeQ(~u, ~x) && Not(MatchQ(~v, (~(c') + ~(d') * ~x) ^ ~(m') <-- FreeQ([c, d, m], ~x))))
    @apply_utils Antiderivative(ℯ ^ (~(n') * asinh(~u)), ~x) => Antiderivative((~u + sqrt(1 + (~u) ^ 2)) ^ ~n, ~x) <-- IntegerQ(~n) && PolyQ(~u, ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * ℯ ^ (~(n') * asinh(~u)), ~x) => Antiderivative((~x) ^ ~m * (~u + sqrt(1 + (~u) ^ 2)) ^ ~n, ~x) <-- RationalQ(~m) && (IntegerQ(~n) && PolyQ(~u, ~x))
end

max_stars_repo_path: src/rules/7 Inverse hyperbolic functions/7.1 Inverse hyperbolic sine/7.1.6 Miscellaneous inverse hyperbolic sine.jl


text:
include("tsne_imports.jl")

get_output_path(path_to_output) = let
  joinpath(splitdir(path_to_output)[1], "output.bson")
end

results_dir = "$(ENV["SCRATCH_DIR"])/cytof/results/repfam/test-sim-6-8-7"
csvdir = "viz/csv"
mkpath(csvdir)

function get_best_lam(samples)
  lams = [s[:theta__lam] for s in samples]
  Ws = [s[:theta__W] for s in samples]
  Zs = [s[:theta__Z] for s in samples]
  I = length(lams[1])

  best_idx = [PlotUtils.estimate_ZWi_index(Zs, Ws, i) for i in 1:I]
  best_lam = [lams[best_idx[i]][i] for i in 1:I]

  return best_lam
end

for zind in 1:3
  for pmiss in (0.0, 0.2)
    for phi in (0, 1, 10, 25, 100)
      # Get simulation name
      simname = "pmiss$(pmiss)-phi$(phi)-zind$(zind)"
      println("Processing: ", simname)

      # Get path to simuldate data.
      path_to_output = joinpath(results_dir, simname, "output.bson")

      # Get output
      output = BSON.load(path_to_output)

      # Get relevant parameters
      samples =  [s[1] for s in output[:samples][1]]

      # Get best lam (for each sample)
      best_lam = vcat(get_best_lam(samples)...)

      # Write output to file.
      open("$(csvdir)/rfam-$(simname).csv", "w") do io
        writedlm(io, best_lam)
      end
    end
  end
end


text:
using BinaryBuilder

# To distinguish between upstream changes and changes to the JLL, we use:
#     version = 100 * upstream version + offset
# See C/Coin-OR/coin-or-common.jl for more details.
version = v"200.0.800"

# Collection of sources required to build ECOSBuilder
sources = [
    GitSource("https://github.com/embotech/ecos.git", "3b98fe0376ceeeb8310a06694b0a84ac59920f3f")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/ecos*
make shared

mkdir -p ${libdir}
cp libecos.${dlext} ${libdir}
cp -r include ${prefix}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(;experimental=true)

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

# Dependencies that must be installed before this package can be built
dependencies = Dependency[]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using CollisionDetection
using StaticArrays
using JLD
using Base.Test

fn = normpath(joinpath(dirname(@__FILE__),"center_sizes.jld"))
d = load(fn)

tmp = d["ctrs"]
ctrs = [SVector(q...) for q in tmp]
rads = d["rads"]

tree = Octree(ctrs, rads)

# extract all the triangles that (potentially) intersect octant (+,+,+)
pred(i) = all(ctrs[i].+rads[i] .> 0)
bb = SVector(0.5, 0.5, 0.5), 0.5
ids = collect(searchtree(pred, tree, bb))
@test length(ids) == 154

# using MATLAB
# A = [c[i] for i in 1:3, c in ctrs].'
# @mput A rads
# @matlab begin
#     plot3(A[:,1],A[:,2],A[:,3],".")
# end
#
# B = [ctrs[j][i] for i in 1:3, j in ids].'
# @mput B
# @matlab begin
#     hold("on")
#     plot3(B[:,1],B[:,2],B[:,3],"or")
# end

max_stars_repo_name: shakisparki/CollisionDetection

max_issues_repo_name: shakisparki/CollisionDetection

max_forks_repo_name: shakisparki/CollisionDetection


text:
pywavelet= false
try
    global pywt = pyimport("pywt")
    global swt = pywt.swt
    global pywavelet = true
catch
    @warn "Skipping joSWTtests - PyCall clouldn't import PyWavelets"
end

families = ("haar", "db", "sym", "coif")
pywavelet ? wavelets = vcat([pywt.wavelist(name) for name in families]...) : wavelets = []

tsname="joSWT"
@testset "$tsname" begin
for w in wavelets
    for m in [124, 125]
        Swt = joSWT(m, w)
        u = rand(m)
        v = Swt*u
        verbose && println("$tsname, $w, $m")
        @testset "$w, $m" begin
            @test isadjoint(Swt)[1]
            @test islinear(Swt)[1]
            # Test orthogonality
            @test norm((adjoint(Swt)*Swt)* u - u) < joTol
            @test norm((Swt*adjoint(Swt))* v - v) < joTol
        end
    end
end # end test loop
end


text:
function test_Group_interface(G::Group)
    @testset "Group interface" begin
        @testset "Iteration protocol" begin
            IS = Base.IteratorSize(typeof(G))
            if IS isa Base.HasLength || IS isa Base.HasShape
                @test isfinite(G) == true
                @test length(G) isa Integer
                @test length(G) > 0

                @test eltype(G) <: GroupElement
                @test GroupsCore.elem_type(typeof(G)) == eltype(G)
                @test one(G) isa eltype(G)

                @test first(iterate(G)) isa eltype(G)
                _, s = iterate(G)
                @test first(iterate(G, s)) isa eltype(G)
                @test isone(first(G))
            else
                @test isfinite(G) == false
            end
        end

        @testset "Group generators" begin
            @test GroupsCore.hasgens(G) isa Bool

            if GroupsCore.hasgens(G)
                @test ngens(G) isa Int
                @test gens(G) isa AbstractVector{eltype(G)}
                @test length(gens(G)) == ngens(G)
                @test first(gens(G)) == gens(G, 1)
                @test last(gens(G)) == gens(G, ngens(G))
            else
                # TODO: throw something more specific
                @test_throws ErrorException gens(G)
                @test_throws ErrorException ngens(G)
            end
        end

        @testset "order, rand" begin
            if isfinite(G)
                @test order(Int16, G) isa Int16
                @test order(G) isa Integer
                @test order(G) >= 1
            end
            @test rand(G) isa GroupElement
            @test rand(G, 2) isa AbstractVector{eltype(G)}
            g, h = rand(G, 2)
            @test parent(g) === parent(h) === G

            @test GroupsCore.pseudo_rand(G) isa eltype(G)
            @test GroupsCore.pseudo_rand(G, 2, 2) isa AbstractMatrix{eltype(G)}

            g, h = GroupsCore.pseudo_rand(G, 2)
            @test parent(g) === parent(h) === G
        end
    end
end

function test_GroupElement_interface(g::GEl, h::GEl) where {GEl<:GroupElement}

    @assert parent(g) === parent(h)

    @testset "GroupElement interface" begin

        @testset "Parent methods" begin
            @test parent(g) isa Group
            G = parent(g)

            @test eltype(G) == typeof(g)
            @test GroupsCore.parent_type(typeof(g)) == typeof(G)
            @test GroupsCore.elem_type(GroupsCore.parent_type(g)) == typeof(g)

            @test one(g) isa eltype(G)

            @test one(G) == one(g) == one(h)

            if !isbits(g)
                @test one(G) !== one(g)
            end

            @test isone(one(G))
        end

        @testset "Equality, deepcopy && hash" begin
            @test (g == h) isa Bool

            @test istrulyequal(g, h) isa Bool
            @test g == g
            @test istrulyequal(h, h)

            if g != h
                @test !istrulyequal(g, h)
            end

            @test deepcopy(g) isa typeof(g)
            @test deepcopy(g) == g
            if !isbits(g)
                @test deepcopy(g) !== g
            end
            k = deepcopy(g)
            @test parent(k) === parent(g)
            @test hash(g) isa UInt
            @test hash(g) == hash(k)
        end

        @testset "Group operations" begin
            old_g, old_h = deepcopy(g), deepcopy(h)

            # check that the default operations don't mutate their arguments
            @test inv(g) isa typeof(g)
            @test (g, h) == (old_g, old_h)

            @test g * h isa typeof(g)
            @test (g, h) == (old_g, old_h)

            @test g^2 == g * g
            @test (g, h) == (old_g, old_h)

            @test g^-3 == inv(g) * inv(g) * inv(g)
            @test (g, h) == (old_g, old_h)

            @test (g * h)^-1 == inv(h) * inv(g)
            @test (g, h) == (old_g, old_h)

            @test conj(g, h) == inv(h) * g * h
            @test (g, h) == (old_g, old_h)

            @test ^(g, h) == inv(h) * g * h
            @test (g, h) == (old_g, old_h)

            @test comm(g, h) == g^-1 * h^-1 * g * h
            @test (g, h) == (old_g, old_h)

            @test comm(g, h, g) == conj(inv(g), h) * conj(conj(g, h), g)
            @test (g, h) == (old_g, old_h)

            @test isone(g * inv(g)) && isone(inv(g) * g)
            @test (g, h) == (old_g, old_h)

            @test g / h == g * inv(h)
            @test (g, h) == (old_g, old_h)
        end

        @testset "Misc GroupElement methods" begin
            @test one(g) isa typeof(g)
            @test isone(g) isa Bool
            @test isone(one(g))

            @test hasorder(g) isa Bool

            if hasorder(g)
                @test order(g) isa Integer
                @test order(Int16, g) isa Int16
                @test order(g) >= 1
                @test iszero(rem(order(parent(g)), order(g)))

                if g^2 != one(g)
                    @test order(g) > 2
                end
                @test order(inv(g)) == order(g)
                @test order(one(g)) == 1
            end

            @test similar(g) isa typeof(g)
        end

        one!, inv!, mul!, conj!, comm!, div_left!, div_right! = (
            GroupsCore.one!,
            GroupsCore.inv!,
            GroupsCore.mul!,
            GroupsCore.conj!,
            GroupsCore.comm!,
            GroupsCore.div_left!,
            GroupsCore.div_right!,
        )

        @testset "In-place operations" begin
            old_g, old_h = deepcopy(g), deepcopy(h)
            out = similar(g)

            @test isone(one!(g))
            g = deepcopy(old_g)

            @test inv!(out, g) == inv(old_g)
            @test g == old_g
            @test inv!(out, g) == inv(old_g)
            g = deepcopy(old_g)

            @testset "mul!" begin
                @test mul!(out, g, h) == old_g * old_h
                @test (g, h) == (old_g, old_h)

                @test mul!(out, g, h) == old_g * old_h
                @test (g, h) == (old_g, old_h)

                @test mul!(g, g, h) == old_g * old_h
                @test h == old_h
                g = deepcopy(old_g)

                @test mul!(h, g, h) == old_g * old_h
                @test g == old_g
                h = deepcopy(old_h)

                @test mul!(g, g, g) == old_g * old_g
                g = deepcopy(old_g)
            end

            @testset "conj!" begin
                res = old_h^-1 * old_g * old_h
                @test conj!(out, g, h) == res
                @test (g, h) == (old_g, old_h)

                @test conj!(g, g, h) == res
                @test h == old_h
                g = deepcopy(old_g)

                @test conj!(h, g, h) == res
                @test g == old_g
                h = deepcopy(old_h)

                @test conj!(g, g, g) == old_g
                g = deepcopy(old_g)
            end

            @testset "comm!" begin
                res = old_g^-1 * old_h^-1 * old_g * old_h

                @test comm!(out, g, h) == res
                @test (g, h) == (old_g, old_h)

                @test comm!(out, g, h) == res
                @test (g, h) == (old_g, old_h)

                @test comm!(g, g, h) == res
                @test h == old_h
                g = deepcopy(old_g)

                @test comm!(h, g, h) == res
                @test g == old_g
                h = deepcopy(old_h)
            end

            @testset "div_[left|right]!" begin
                res = g * h^-1
                @test div_right!(out, g, h) == res
                @test (g, h) == (old_g, old_h)

                @test div_right!(g, g, h) == res
                @test h == old_h
                g = deepcopy(old_g)

                @test div_right!(h, g, h) == res
                @test g == old_g
                h = deepcopy(old_h)

                @test div_right!(g, g, g) == one(g)
                g = deepcopy(old_g)

                res = h^-1 * g
                @test div_left!(out, g, h) == res
                @test (g, h) == (old_g, old_h)

                @test div_left!(g, g, h) == res
                @test h == old_h
                g = deepcopy(old_g)

                @test div_left!(h, g, h) == res
                @test g == old_g
                h = deepcopy(old_h)

                @test div_left!(g, g, g) == one(g)
                g = deepcopy(old_g)
            end
        end
    end
end