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, Tuple{Vector{CT}})
    # workaround for SubArray type
    if T <: SubArray
        return Core.Compiler.return_type(f, Tuple{typeof(x)})
    elseif T <: AbstractVector
        T = AbstractVector{Union{Missing, eltype(T)}}
    elseif T <: Tuple
        T = Union{Missing, Core.Compiler.return_type(f, Tuple{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, Tuple{ntuple(i->Vector{CT[i]}, length(x))...})
    # workaround for SubArray type
    if T <: SubArray
        return Core.Compiler.return_type(f, Tuple{(typeof.(x))...})
    elseif T <: AbstractVector
        T = AbstractVector{Union{Missing, eltype(T)}}
    elseif T <: Tuple
        T = Union{Missing, Core.Compiler.return_type(f, Tuple{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)

	if nrow(ds) > 2^23 && !stable && 5<length(colidx)<16 # the result is stable anyway
		if !mapformats || all(==(identity), getformat.(Ref(ds), colidx))
			return _gather_groups_hugeds_multicols(ds, cols, Val(T); threads = threads)
		end
	end

    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, Tuple{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

### Special path for huge ds and multiple cols - trade off between compilation and performance
# table columns are passed as a tuple of vectors to ensure type specialization - From DataFrames.jl
isequal_row(cols::Tuple{AbstractVector}, r1::Int, r2::Int) =
    isequal(cols[1][r1], cols[1][r2])
isequal_row(cols::Tuple{Vararg{AbstractVector}}, r1::Int, r2::Int) =
    isequal(cols[1][r1], cols[1][r2]) && isequal_row(Base.tail(cols), r1, r2)

isequal_row(cols1::Tuple{AbstractVector}, r1::Int, cols2::Tuple{AbstractVector}, r2::Int) =
    isequal(cols1[1][r1], cols2[1][r2])
isequal_row(cols1::Tuple{Vararg{AbstractVector}}, r1::Int,
            cols2::Tuple{Vararg{AbstractVector}}, r2::Int) =
    isequal(cols1[1][r1], cols2[1][r2]) &&
        isequal_row(Base.tail(cols1), r1, Base.tail(cols2), r2)

_grabrefs(x) = DataAPI.refpool(x) == nothing ? x : DataAPI.refarray(x)
function _gather_groups_hugeds_multicols(ds, cols, ::Val{T}; threads = true) where T
	colidx = index(ds)[cols]
	rhashes = byrow(ds, hash, cols, threads = threads)
	colsvals = ntuple(i->_grabrefs(_columns(ds)[colidx[i]]), length(colidx))
	create_dict_hugeds_multicols(colsvals, rhashes, Val(T))
end

function create_dict_hugeds_multicols(colvals, rhashes, ::Val{T}) where T
	sz = max(1 + ((5 * length(rhashes)) >> 2), 16)
    sz = 1 << (8 * sizeof(sz) - leading_zeros(sz - 1))
    @assert 4 * sz >= 5 * length(rhashes)
    szm1 = sz-1
    gslots = zeros(T, sz)
	groups = Vector{T}(undef, length(rhashes))
    ngroups = 0
    @inbounds for i in eachindex(rhashes)
        # find the slot and group index for a row
        slotix = rhashes[i] & szm1 + 1
        gix = -1
        probe = 0
        while true
            g_row = gslots[slotix]
            if g_row == 0 # unoccupied slot, current row starts a new group
                gslots[slotix] = i
                gix = ngroups += 1
                break
            elseif rhashes[i] == rhashes[g_row] # occupied slot, check if miss or hit
                if isequal_row(colvals, i, Int(g_row)) # hit
                    gix = groups[g_row]
                    break
                end
            end
            slotix = slotix & szm1 + 1 # check the next slot
            probe += 1
            @assert probe < sz
        end
        groups[i] = gix
    end
    return groups, gslots, ngroups
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)

_findall(B::AbstractVector{Bool}) = findall(convert(BitVector, B))

# 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 two vectors 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, Tuple{eltype(x)}), Core.Compiler.return_type(f2, Tuple{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, Tuple{eltype(x)}), Core.Compiler.return_type(f2, Tuple{eltype(y)})), eltype(x), eltype(y), typeof(x), typeof(y)}(x, y, f1, f2, length(x), length(y))
			end
		end
	end

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:
function layersGroundTruth1D(gt::DataGroundTruth)
	PL = []
	push!(PL, gt.plt_lyr_cluster[1])
	push!(PL, gt.plt_lyr_cluster[2])
	PL
end

function layersGroundTruth2D(id::Int, gt::DataGroundTruth)
	PL = Gadfly.Context[]
	PL = union(PL, gt.plt_lyr_cluster[id])
	PL
end
function layersGroundTruth2D(ids::Vector{Int}, gt::DataGroundTruth)
	PL = Gadfly.Context[]
	for id in ids
		PL = union(PL, gt.plt_lyr_cluster[id])
	end
	PL
end

function layersExpertBelief1D(cs::ClassificationSystem;
															margdim=1,
															colors::NothingUnion{Dict{Int, String}}=nothing )
	# use defaults
	colors = colors == nothing ? cs.colors : colors

  # prep data
	pr, cl = BallTreeDensity[], String[]
	for p in cs.expertBelief
		push!(pr, marginal(p[2],[margdim]))
		push!(cl, colors[p[1]]) # get accompanying color
	end

	# plot and return layers
	pl_init = plotKDE(pr,c=cl);
	pl_init.layers
end

function layersExpertBelief2D(id::Int, cs::ClassificationSystem;
															margdim=1,
															colors::NothingUnion{Dict{Int, String}}=nothing )
	# use defaults
	colors = colors == nothing ? cs.colors : colors

  # prep data
	pr, cl = BallTreeDensity[], String[]
	for p in cs.expertBelief
		push!(pr, marginal(p[2],[margdim]))
		push!(cl, colors[p[1]]) # get accompanying color
	end

	# plot and return layers
	pl_init = plotKDE(pr,c=cl);
	pl_init.layers
end

function layersCurrentBelief1D(cs::ClassificationSystem;
															margdim=1,
															colors::NothingUnion{Dict{Int, String}}=nothing )
	# use defaults
	colors = colors == nothing ? cs.colors : colors

  # prep data
	pr, cl = BallTreeDensity[], String[]
	for p in cs.currentBelief
		push!(pr, marginal(p[2],[margdim]))
		push!(cl, colors[p[1]]) # get accompanying color
	end

	# plot and return layers
	pl_init = plotKDE(pr,c=cl);
	pl_init.layers
end

function plotUtil1D(;
		sampledata::NothingUnion{SampleData}=nothing,
		groundtruth::NothingUnion{DataGroundTruth}=nothing,
		cs::NothingUnion{ClassificationSystem}=nothing,
		expertcolor::NothingUnion{Dict{Int,String}}=nothing,
		drawcurrent::Bool=false,
		currentcolor::NothingUnion{Dict{Int,String}}=nothing
		)

		PL = []
		PL = groundtruth == nothing ? PL : union(PL, layersGroundTruth1D(groundtruth))
		# push!(PL, Guide.title("True classification for intersecting clusters in 1 dimension"))

		PL = cs == nothing ? PL : union(PL, layersExpertBelief1D(cs, colors=expertcolor))

		PL = (cs == nothing || !drawcurrent) ? PL : union(PL, layersCurrentBelief1D(cs, colors=currentcolor))

		# layer(x=data.samples,y=zeros(size(data.samples)),Geom.point,Theme(default_color=colorant"gray")),

		plot(PL...)
end

function plotUtil2D(id::Union{Int, Vector{Int}};
		sampledata::NothingUnion{SampleData}=nothing,
		groundtruth::NothingUnion{DataGroundTruth}=nothing,
		cs::NothingUnion{ClassificationSystem}=nothing,
		expertcolor::NothingUnion{Dict{Int,String}}=nothing,
		drawcurrent::Bool=false,
		currentcolor::NothingUnion{Dict{Int,String}}=nothing
		)

		PL = []
		PL = groundtruth == nothing ? PL : union(PL, layersGroundTruth2D(id, groundtruth))
		# push!(PL, Guide.title("True classification for intersecting clusters in 1 dimension"))

		PL = cs == nothing ? PL : union(PL, layersExpertBelief2D(id, cs, colors=expertcolor))
		#
		# PL = (cs == nothing || !drawcurrent) ? PL : union(PL, layersCurrentBelief1D(cs, colors=currentcolor))

		plot(PL...)
end

function plotEMStatus(params, stats)
  plfrac = plot(
  layer(x=1:params.EMiters, y=stats.POPFRAC[1,:], Geom.line, Theme(default_color=colorant"blue" )),
  layer(x=1:params.EMiters, y=stats.POPFRAC[2,:], Geom.line, Theme(default_color=colorant"red"  )),
  # layer(x=1:params.EMiters, y=abs(dbg.INDV_MISASSIGN_C[1] - dbg.INDV_MISASSIGN_C[2]), Geom.line, Theme(default_color=colorant"magenta"  )),
  Guide.title("Population fraction estimates"),
  Guide.ylabel("%")
  )

	plkl = plot(
	layer(x=1:params.EMiters, y=abs.(stats.SEQKLDIVERG[1][:]), Geom.line, Theme(default_color=colorant"blue" )),
	layer(x=1:params.EMiters, y=abs.(stats.SEQKLDIVERG[2][:]), Geom.line, Theme(default_color=colorant"red"  )),
	# layer(x=1:params.EMiters, y=abs(dbg.INDV_MISASSIGN_C[1] - dbg.INDV_MISASSIGN_C[2]), Geom.line, Theme(default_color=colorant"magenta"  )),
	Guide.title("Abs sequential KL divergence, | D(p_{k+1} || p_{k}) |"),
	)

	vstack(plfrac, plkl)
end

function plotClassificationStats(params, dbg)
  pl_accur = plot(
  layer(x=1:params.EMiters, y=dbg.ACCUR_C[1], Geom.line, Theme(default_color=colorant"blue" )),
  layer(x=1:params.EMiters, y=dbg.ACCUR_C[2], Geom.line, Theme(default_color=colorant"red"  )),
  Guide.title("Absolute percentage error in sample count among classifications"),
  Guide.ylabel("%")
  );

  pl_rel_accur = plot(
  layer(x=1:params.EMiters, y=dbg.REL_ACCUR_C[1], Geom.line, Theme(default_color=colorant"blue" )),
  layer(x=1:params.EMiters, y=dbg.REL_ACCUR_C[2], Geom.line, Theme(default_color=colorant"red"  )),
  Guide.title("Relative percentage error in sample count among classifications"),
  Guide.ylabel("%")
  );

  # pl_indvmissassign = plot(
  # layer(x=1:params.EMiters, y=dbg.INDV_MISASSIGN_C[1], Geom.line, Theme(default_color=colorant"blue" )),
  # layer(x=1:params.EMiters, y=dbg.INDV_MISASSIGN_C[2], Geom.line, Theme(default_color=colorant"red"  )),
  # # layer(x=1:params.EMiters, y=abs(dbg.INDV_MISASSIGN_C[1] - dbg.INDV_MISASSIGN_C[2]), Geom.line, Theme(default_color=colorant"magenta"  )),
  # Guide.title("Common percentage error count among classifications (unweighted)"),
  # Guide.ylabel("%")
  # );

  vstack(pl_accur, pl_rel_accur)
end

ext: jl


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# ╔═╡ e660c4b8-6999-11eb-3f49-f171c07a7062
begin
	using StatsPlots
	using Random
	using TimeSeries
	using Dates
	using Statistics
	using DataFrames
end

# ╔═╡ ba9e0f02-6999-11eb-003c-ff442833cf0b
md"""
# Introduction to Time Series III

> Author: Edwin Bedolla
>
> Data: Original, May 6th, 2020. This version, 7th February 2021.

This document closes the introduction to Time Series by exploring the last concepts which
are the *estimation* of the **autocovariance** and **autocorrelation** functions, with a
simple example.

As always, we import the needed packages.

"""

# ╔═╡ f6a2939c-6999-11eb-0bd5-855cff3d8171
gr();

# ╔═╡ f92da192-6999-11eb-386d-b185e1e75364
# Ensure reproducibility of the results
rng = MersenneTwister(5489);

# ╔═╡ fc8ba370-6999-11eb-3438-1fe2a079401e
md"""
## Estimation of correlation

In the last part we studied the **theoretical** definitions of the autocorrelation (ACF)
and autocovariance functions (AVF), where we did some extensive examples on how to compute the
theoretical values.

But in reality we only have *sampled data* and we usually just have a subset or
*sample* of the original *population.* When that is the case, we can only estimate
the respective statistics by using **estimated** methods. In this document we will
explore the definitions and algorithms to estimate the autocorrelation and
autocovariance functions, we close the document exploring the **lag operator** and the
importance of it.

We start with a simple example.

"""

# ╔═╡ 0720775c-699a-11eb-3fdf-aba950119a6f
# Create a range of time for a year, spaced evenly every 1 minute
dates = DateTime(2018, 1, 1, 1):Dates.Minute(1):DateTime(2018, 12, 31, 24);

# ╔═╡ 0fc5a500-699a-11eb-3703-5d2262a3ea77
# Create a white noise array
xt = randn(rng, length(dates));

# ╔═╡ 14348eb0-699a-11eb-2427-b31244499c48
md"""
### Autoregressive time series

In this section we will build an **autoregressive time series**; this is a big topic within
Time Series analysis so we will have a special chapter for that, but here we want
to introduce the concept.

An autoregressive time series takes information from past intervals to build
future information. We will build an autoregressive time series now and plot it.
"""

# ╔═╡ 1f184cea-699a-11eb-180f-93c070a28df9
# This is the autoregressive time series
ar = @. 5.0 + xt[1:(end-1)] - 0.7 * xt[2:end];

# ╔═╡ 2d913bb0-699a-11eb-268b-4f142ec0c262
# Note that we need to remove the last element from the dates
df_ar = TimeArray(dates[1:(end-1)], ar) |> DataFrame;

# ╔═╡ 3fd5eda2-699a-11eb-0dc3-1dd47b2abd60
# Plotting indices
idxs = 1:200;

# ╔═╡ 327b97a6-699a-11eb-2a3b-bd798f261a50
# We can now plot the time series, only 100 elements though
@df df_ar plot(:timestamp[idxs], :A[idxs])

# ╔═╡ 4cb0bb06-699a-11eb-21e3-f9dc016d4ef3
md"""
As we can see, this time series is not really different from the ones that we
have studied so far, but one of the fundamental properties of the autoregressive
model or time series is that is shows **periodicity**; we will explore all of this
later.

"""

# ╔═╡ 5a5dfe3a-699a-11eb-079e-33d9f0391bc1
md"""
### Correlation of an autoregressive time series

As we said before, we are concerned with estimating the ACF and AVF of a time
series when we only have sampled data at hand. In the case of the autoregressive
time series we are studying, I have obtained the true value of the AVF,
``\gamma`` which is the following

"""

# ╔═╡ 6e9d7434-699a-11eb-3686-d95a3b2eb686
true_γ = 1.49;

# ╔═╡ 781bce66-699a-11eb-3298-69d8c28bd143
md"""
We compare it to `Julia`'s built-in function to compute the covariance
value of the time series.
"""

# ╔═╡ 80747da6-699a-11eb-05b7-2d59cfa8f535
julia_γ = cov(df_ar[:, :A], df_ar[:, :A])

# ╔═╡ 8aa25e42-699a-11eb-300a-eda4c5ae4ef9
md"""

#### Autocovariance function

We now turn our attention to defining the estimated, *unbiased* AVF, which
is defined as follows

```math
\hat{\gamma}(h)=\frac{1}{n} \sum_{t=1}^{n-h} (x_{t+h}-\bar{x}) (x_{t}-\bar{x})
```

where the **sample mean**, ``\bar{x}``, is estimated as

```math
\bar{x}=\frac{1}{N} \sum_{i}^N x_i
```

We will now implement the estimated AVF in `Julia`.
"""

# ╔═╡ a7178d88-699a-11eb-35f6-a921abd1b7dd
function avf(x, h::Integer)
    # Get the total number of elements in the time series
    n = length(x)
    # Compute sample mean
    μ = sum(x) / n
    # Initialize the sample autocovariance
    γ = 0.0

    # Get the final index according to the lag
    endidx = Int(n - h)
    # Use the definition for the sample autocovariance
    @inbounds for i in 1:endidx
        γ += (x[i + h] - μ) * (x[i] - μ)
    end

    # And normalize
    return γ / n
end;

# ╔═╡ b452ff9e-699a-11eb-1c30-e58c82796b47
md"""
Let us now use this new function and compute the estimated AVF.
"""

# ╔═╡ bb0c3172-699a-11eb-1956-fdba69b12e4f
γ_by_hand = avf(df_ar[:, :A], 0)

# ╔═╡ c21299a8-699a-11eb-1cc6-fbefa6e663b1
md"""
We can see that the results are very similar to those obtained with
the built-in `cov`, and very close to the true theoretical value. We can obtain
the true theoretical value with increasing samples from the time series.
"""

# ╔═╡ c9c05052-699a-11eb-1cf9-13bbde5a16c4
md"""
#### Autocorrelation function

A similar estimator for the ACF can be used, and in particular we have

```math
\hat{\rho}(h)=\frac{\hat{\gamma}(h)}{\hat{\gamma}(0)}=\frac{\sum_{t=1}^{N-h} (x_{t+h}-\bar{x}) (x_{t}-\bar{x})}{\sum_{t=1}^{N} (x_{t}-\bar{x})^2}
```

so we can use this to our advantage and employ the AVF to compute the ACF.

In our implementation we exploit this property bellow.
"""

# ╔═╡ e2d0731a-699a-11eb-0b57-5d8d2f14aa05
function acf(x, h::Integer)
    # Get the total number of elements in the time series
    n = length(x)
    # Obtain the sample mean
    μ = sum(x) / n

    # Compute the autocovariance for a time series
    # without lag
    autocov = sum((x[i] - μ)^2 for i = eachindex(x))
    # Normalize
    autocov /= n

    # Get the autocovariance with lag and divide it
    # by the one without lag
    return avf(x, h) / autocov
end;

# ╔═╡ eb28f34a-699a-11eb-2a8e-05a9d4c85a39
md"""
We now compute the estimated ACF. Notice that we need to specify the value of
``h=1`` which adds a *lag* to the time series.
"""

# ╔═╡ f3f50a02-699a-11eb-2260-5389e475279c
ρ_by_hand = acf(df_ar[:, :A], 1)

# ╔═╡ fb2c920e-699a-11eb-1981-d3610f2bec4e
md"""
We now compare it to the one obtained by the Pearson correlation function with
the built-in `cor` function in `Julia`.
"""

# ╔═╡ 045e3a6e-699b-11eb-175e-3db7a6a91aa0
ρ_julia = cor(df_ar[1:(end - 1), :A], df_ar[2:end, :A])

# ╔═╡ 0a0c7a5a-699b-11eb-3407-1b91a137ddeb
md"""
The results are quite similar, just as we expected.
"""

# ╔═╡ 1286d234-699b-11eb-02f9-87c97528f072
md"""
## Lag in Time Series

Up until now we have mentioned several times the word **lag** in different scenarios.
It is now time to address the meaning of this word in the context of Time Series.
Time Series are a process in different time intervals, and we are sometimes interested
in what happens between two time intervals ``s`` and ``t``
within the same Time Series.

A special type of operation, called the **lag operator** produces the value of a previous
time interval with respect to the current time interval.
The lag operator is denoted by ``L`` and when applied to a Time Series
``X=\{x_1,x_2,x_3,\dots\}`` we get the following

```math
Lx_k = x_{k-1}
```

When doing something like this `df_ar[1:(end-1), :A]` we are essentially applying the lag
operator by hand, were the lag is one in this case. The `TimeSeries` package has support for this. Let us do a simple example computing the ACF like before.

"""

# ╔═╡ 65f6025a-699b-11eb-3ade-79ceaaab9d13
# We need to use a TimeArray in order to apply the lag operator
ts = TimeArray(dates[1:(end-1)], ar);

# ╔═╡ 6e93168c-699b-11eb-019b-e14265dbc1dd
# We use lag to put the past value in the current value of the time series
ts_lag = lag(ts) |> DataFrame;

# ╔═╡ 72f6ffea-699b-11eb-1a40-0fcd4586faf2
# Convert to DataFrame for easier handling
df_ts = DataFrame(ts);

# ╔═╡ 7a25936c-699b-11eb-1814-5348ae852ed3
# Compute the ACF, note that instead of using indexing we just use the
# lagged time series first
lag_value = cor(ts_lag[:, :A], df_ts[2:end, :A]);

# ╔═╡ 7e4ab922-699b-11eb-0940-31f4a6eef3b5
# And it should be the same as before
lag_value == ρ_julia

# ╔═╡ 91586a8c-699b-11eb-33e9-c1a053cb23c2
md"""
We can go back to an arbitrary time interval if we *raise* the lag operator to a given
*power*, like so

```math
L^n x_k = x_{k-n}
```

for example ``L^2 x_k = x_{x-2}`` so we actually got the
previous-to-last value from the time series.

Furthermore, because the lag operator can be raised to arbitrary powers we can build
**polynomials** out of the lag operator, like so

```math
a(L) = a_0 + a_1 L + a_2 L^2 + \cdots \\
a(L)x_t = a_0 x_t + a_1 x_{t-1} + a_2 x_{t-2} + a_3 x_{t-3} \cdots
```

which resembles an autoregressive time series. In fact it is the use in autoregressive
models and moving averages where the lag operator is of the utmost importance.

For example, in this document we used the autoregressive model

```math
y_t = 5+x_t-0.7x_{t-1}
```

which could be written in terms of the lag operator as

```math
y_t = 5+x_t-0.7Lx_t = 5+x_t(1-0.7L) .
```

We will come back to this when we study the
[ARIMA processes](https://en.wikipedia.org/wiki/Autoregressive_integrated_moving_average) later on.

"""

# ╔═╡ Cell order:
# ╟─ba9e0f02-6999-11eb-003c-ff442833cf0b
# ╠═e660c4b8-6999-11eb-3f49-f171c07a7062
# ╠═f6a2939c-6999-11eb-0bd5-855cff3d8171
# ╠═f92da192-6999-11eb-386d-b185e1e75364
# ╟─fc8ba370-6999-11eb-3438-1fe2a079401e
# ╠═0720775c-699a-11eb-3fdf-aba950119a6f
# ╠═0fc5a500-699a-11eb-3703-5d2262a3ea77
# ╟─14348eb0-699a-11eb-2427-b31244499c48
# ╠═1f184cea-699a-11eb-180f-93c070a28df9
# ╠═2d913bb0-699a-11eb-268b-4f142ec0c262
# ╠═3fd5eda2-699a-11eb-0dc3-1dd47b2abd60
# ╠═327b97a6-699a-11eb-2a3b-bd798f261a50
# ╟─4cb0bb06-699a-11eb-21e3-f9dc016d4ef3
# ╟─5a5dfe3a-699a-11eb-079e-33d9f0391bc1
# ╠═6e9d7434-699a-11eb-3686-d95a3b2eb686
# ╟─781bce66-699a-11eb-3298-69d8c28bd143
# ╠═80747da6-699a-11eb-05b7-2d59cfa8f535
# ╟─8aa25e42-699a-11eb-300a-eda4c5ae4ef9
# ╠═a7178d88-699a-11eb-35f6-a921abd1b7dd
# ╟─b452ff9e-699a-11eb-1c30-e58c82796b47
# ╠═bb0c3172-699a-11eb-1956-fdba69b12e4f
# ╟─c21299a8-699a-11eb-1cc6-fbefa6e663b1
# ╟─c9c05052-699a-11eb-1cf9-13bbde5a16c4
# ╠═e2d0731a-699a-11eb-0b57-5d8d2f14aa05
# ╟─eb28f34a-699a-11eb-2a8e-05a9d4c85a39
# ╠═f3f50a02-699a-11eb-2260-5389e475279c
# ╟─fb2c920e-699a-11eb-1981-d3610f2bec4e
# ╠═045e3a6e-699b-11eb-175e-3db7a6a91aa0
# ╟─0a0c7a5a-699b-11eb-3407-1b91a137ddeb
# ╟─1286d234-699b-11eb-02f9-87c97528f072
# ╠═65f6025a-699b-11eb-3ade-79ceaaab9d13
# ╠═6e93168c-699b-11eb-019b-e14265dbc1dd
# ╠═72f6ffea-699b-11eb-1a40-0fcd4586faf2
# ╠═7a25936c-699b-11eb-1814-5348ae852ed3
# ╠═7e4ab922-699b-11eb-0940-31f4a6eef3b5
# ╟─91586a8c-699b-11eb-33e9-c1a053cb23c2

max_issues_repo_name: BioLogs/TimeSeriesConcepts


text:
"""
Hypercube iterates over a hypercube of radius γ around xᶜ, bounded inside the space {1..sz}ᴺ.

# Examples
```jldoctest
julia> hc= Hypercube((3,3),1,(10,10));

julia> foreach(println, hc)
(2, 2)
(3, 2)
(4, 2)
(2, 3)
(3, 3)
(4, 3)
(2, 4)
(3, 4)
(4, 4)
```
"""
struct Hypercube{N}
  "xᶜ are the center coordinates"
  xᶜ::NTuple{N,Int}
  "γ is the hypercube radius. Points < γ away from xᶜ belong to the hypercube"
  γ::Int
  "sz is the upper bound of the enclosing space"
  sz::NTuple{N,Int}
  indices::CartesianIndices{N}
end
Hypercube(xᶜ,γ,sz)= Hypercube(xᶜ,γ,sz, start(xᶜ,γ,sz))

# Hypercube implements Iterator interface
start(xᶜ,γ,sz)= CartesianIndices(map( (a,b)-> a:b,
                    max.(xᶜ .- γ, 1),
                    min.(xᶜ .+ γ, sz) ))
Base.iterate(hc::Hypercube)= begin
  i= iterate(hc.indices)
  !isnothing(i) ? (i[1].I,i[2]) : nothing
end
Base.iterate(hc::Hypercube{N}, state) where N= begin
  i= Base.iterate(hc.indices, state)::Union{Nothing, Tuple{CartesianIndex{N}, CartesianIndex{N}}}
  !isnothing(i) ? (i[1].I,i[2]) : nothing
end
#Base.eltype(hc::Hypercube)= eltype(hc.indices)
Base.length(hc::Hypercube)= length(hc.indices)
Base.size(hc::Hypercube)= size(hc.indices)

# Good first approximation!
"""
Hypersphere is approximated with a Hypercube for simplicity

See also: [`Hypercube`](@ref)
"""
const Hypersphere{N}= Hypercube{N}

# Create a hypercube that wraps around the edges of the enclosing space
function wrapping_hypercube(xᶜ::NTuple{N,T},radius::Int,dims::NTuple{N,T}) where {N,T}
  xᶜ= convert(NTuple{N,UIntSP},xᶜ)
  dims= convert(NTuple{N,UIntSP},dims)
  wrapping_hypercube(xᶜ,radius,dims, startWrapping(xᶜ,radius,dims))
end

# BUG Wrapping can't work with a single CartesianIndices! It needs 2/4 of them!
#   from UP -> DIM, 0 -> DOWN !!!
startWrapping(xᶜ::NTuple{N,T},radius::T,dims::NTuple{N,T}) where {N,T}= begin
    CartesianIndices( map((a,b)-> a:b,
      wrapUnder.(xᶜ .- Int(radius), dims),   # BUG wrapUnder does smt funky, should return 10
      wrapOver.(xᶜ .+ Int(radius), dims)))
    end

wrapOver(a::T, lim::T) where T= a>lim ? a-lim : a
wrapUnder(a::T, lim::T) where T= a>lim ? lim-(typemax(T)-a)-T(1) : a


text:
using SparseArrays

abstract type AbstractSolution  end
abstract type AbstractLinearOperator end

############################
###### Solution Types ######
############################

mutable struct FluidSolution <: AbstractSolution
  u::Vector{Float64}
  v::Vector{Float64}
  p::Vector{Float64}
end

mutable struct ScalarSolution <: AbstractSolution
  u::Vector{Float64}
end

#####################################
###### LinearOperator Type ##########
#####################################

mutable struct LinearOperator{T1,T2} <: AbstractLinearOperator
  Op::SparseMatrixCSC{T1,T2}
  rhs::Vector{T1}
end

##### Constructor for LinearOperator ######
function LinearOperator(mesh::AbstractMesh,prob::AbstractProblem)
  Op  = spzeros(length(mesh.xy),length(mesh.xy))
  rhs = zeros(length(mesh.xy))

  return LinearOperator(Op,rhs)
end

function LinearOperator(mesh::FluidMesh,prob::AbstractProblem)
  nu = length(mesh.xy)
  np = length(mesh.xyp)
  nTotal = 2*nu+np
  Op  = spzeros(nTotal,nTotal)
  rhs = zeros(nTotal)

  return LinearOperator(Op,rhs)
end


text:
include("phylo.jl")
include("external.jl")
include("phylo_extra.jl")

max_stars_repo_name: RuneBlaze/clade-phylo

max_issues_repo_name: RuneBlaze/clade-phylo

max_forks_repo_name: RuneBlaze/clade-phylo


text:
module Transit

import JSON
import Base.getindex
import Base

using DataStructures
using URIParser
using Decimals
using Dates

export parse, write,  # intended public use
       Decoder, Encoder, decode, encode  # used by tests

include("tagged_value.jl")
include("tsymbol.jl")
include("tset.jl")
include("turi.jl")
include("link.jl")
include("constants.jl")
include("cache.jl")
include("decoder.jl")
include("utilities.jl")
include("emitter.jl")
include("encoder.jl")
include("writer.jl")
include("reader.jl")

end


text:
struct Sinus{T} <: Distances.SemiMetric
# struct Sinus{T} <: Distances.UnionSemiMetric
    r::Vector{T}
end

Distances.parameters(d::Sinus) = d.r
@inline Distances.eval_op(::Sinus, a::Real, b::Real, p::Real) = abs2(sinpi(a - b) / p)
@inline (dist::Sinus)(a::AbstractArray, b::AbstractArray) = Distances._evaluate(dist, a, b)
@inline (dist::Sinus)(a::Number, b::Number) = abs2(sinpi(a - b) / first(dist.r))

Distances.result_type(::Sinus{T}, Ta::Type, Tb::Type) where {T} = promote_type(T, Ta, Tb)

@inline function Distances._evaluate(d::Sinus, a::AbstractVector, b::AbstractVector) where {T}
    @boundscheck if (length(a) != length(b)) || length(a) != length(d.r)
        throw(DimensionMismatch("Dimensions of the inputs are not matching : a = $(length(a)), b = $(length(b)), r = $(length(d.r))"))
    end
    return sum(abs2, sinpi.(a - b) ./ d.r)
end


text:
# test that PPE values are update with solve, see issue #548

using Test
using IncrementalInference
using DistributedFactorGraphs

@testset "test PPE update during solve" begin

fg = generateCanonicalFG_Kaess(graphinit=true)
ensureAllInitialized!(fg)

# list of variables to check
vars = listVariables(fg)

# fetch values before solve
before = Dict()
for vs in vars
  before[vs] = getVariablePPE(getVariable(fg, vs)) |> getPPESuggested
end

# do the solve
# getSolverParams(fg).dbg = true

# tree = buildTreeReset!(fg)
# drawTree(tree, show=true)

# solveCliqUp!(fg, tree, :l2)
# solveCliqUp!(fg, tree, :x3)
# solveCliqUp!(fg, tree, :x2)

solveTree!(fg)

after = Dict()
for vs in vars
  after[vs] = getVariablePPE(getVariable(fg, vs)) |> getPPESuggested
end

# before and after should be noticably different, because first inferred values have been found
for vs in vars
  errd = norm(before[vs] - after[vs])
  # @show vs, errd
  @test 1e-5 < errd
end

# force recalc and update each PPE
force = Dict()
for vs in vars
  setVariablePosteriorEstimates!(fg, vs)
  force[vs] = getVariablePPE(getVariable(fg, vs)) |> getPPESuggested
  # these need to be close to the same as after
  errd = norm(force[vs] - after[vs])
  # @show vs, errd
  @test errd < 0.1
end

## suspect cliqSubFg updated, but not back to main dfg object... test via load graph

end


text:
data = Data[Data(X), Data(y)]
elements = ModelSanitizer._elements(data)

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

m = LinearModel{Float64}()
fit!(m, X[training_rows, :], y[training_rows])

Test.@test predict(m) isa AbstractVector
Test.@test predict(m, X[training_rows, :]) isa AbstractVector
Test.@test predict(m, X[testing_rows, :]) isa AbstractVector

Test.@test predict(m) == predict(m, X[training_rows, :])
Test.@test mse(m) == mse(m, X[training_rows, :], y[training_rows])
Test.@test rmse(m) == rmse(m, X[training_rows, :], y[training_rows])
Test.@test r2(m) == r2(m, X[training_rows, :], y[training_rows])

training_r2_before_sanitization = r2(m, X[training_rows, :], y[training_rows])
testing_r2_before_sanitization = r2(m, X[testing_rows, :], y[testing_rows])

Test.@test m.X == X[training_rows, :]
Test.@test m.y == y[training_rows]
Test.@test all(m.X .== X[training_rows, :])
Test.@test all(m.y .== y[training_rows])
Test.@test !all(m.X .== 0)
Test.@test !all(m.y .== 0)

sanitize!(Model(m), Data(X), Data(y))

Test.@test predict(m, X[training_rows, :]) isa AbstractVector
Test.@test predict(m, X[testing_rows, :]) isa AbstractVector

training_r2_after_sanitization = r2(m, X[training_rows, :], y[training_rows])
testing_r2_after_sanitization = r2(m, X[testing_rows, :], y[testing_rows])
Test.@test training_r2_before_sanitization == training_r2_after_sanitization
Test.@test testing_r2_before_sanitization == testing_r2_after_sanitization

Test.@test m.X != X[training_rows, :]
Test.@test m.y != y[training_rows]
Test.@test !all(m.X .== X[training_rows, :])
Test.@test !all(m.y .== y[training_rows])
Test.@test all(m.X .== 0)
Test.@test all(m.y .== 0)

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

m = LinearModel{Float64}()
fit!(m, X[training_rows, :], y[training_rows])

Test.@test m.X == X[training_rows, :]
Test.@test m.y == y[training_rows]
Test.@test all(m.X .== X[training_rows, :])
Test.@test all(m.y .== y[training_rows])
Test.@test !all(m.X .== 0)
Test.@test !all(m.y .== 0)

sanitize!(ForceSanitize(m.X))

Test.@test m.y == y[training_rows]
Test.@test all(m.y .== y[training_rows])
Test.@test !all(m.y .== 0)

Test.@test m.X != X[training_rows, :]
Test.@test !all(m.X .== X[training_rows, :])
Test.@test all(m.X .== 0)

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

m = LinearModel{Float64}()
fit!(m, X[training_rows, :], y[training_rows])

Test.@test m.X == X[training_rows, :]
Test.@test m.y == y[training_rows]
Test.@test all(m.X .== X[training_rows, :])
Test.@test all(m.y .== y[training_rows])
Test.@test !all(m.X .== 0)
Test.@test !all(m.y .== 0)

sanitize!(ForceSanitize(m.X), ForceSanitize(m.y))

Test.@test m.X != X[training_rows, :]
Test.@test m.y != y[training_rows]
Test.@test !all(m.X .== X[training_rows, :])
Test.@test !all(m.y .== y[training_rows])
Test.@test all(m.X .== 0)
Test.@test all(m.y .== 0)

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

m = LinearModel{Float64}()
fit!(m, X[training_rows, :], y[training_rows])

Test.@test m.X == X[training_rows, :]
Test.@test m.y == y[training_rows]
Test.@test all(m.X .== X[training_rows, :])
Test.@test all(m.y .== y[training_rows])
Test.@test !all(m.X .== 0)
Test.@test !all(m.y .== 0)

sanitize!(Model(m))

Test.@test m.X == X[training_rows, :]
Test.@test m.y == y[training_rows]
Test.@test all(m.X .== X[training_rows, :])
Test.@test all(m.y .== y[training_rows])
Test.@test !all(m.X .== 0)
Test.@test !all(m.y .== 0)

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

m = LinearModel{Float64}()
fit!(m, X[training_rows, :], y[training_rows])

Test.@test m.X == X[training_rows, :]
Test.@test m.y == y[training_rows]
Test.@test all(m.X .== X[training_rows, :])
Test.@test all(m.y .== y[training_rows])
Test.@test !all(m.X .== 0)
Test.@test !all(m.y .== 0)

sanitize!(Model(m), Data(X))

Test.@test m.y == y[training_rows]
Test.@test all(m.y .== y[training_rows])
Test.@test !all(m.y .== 0)

Test.@test m.X != X[training_rows, :]
Test.@test !all(m.X .== X[training_rows, :])
Test.@test all(m.X .== 0)

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

m = LinearModel{Float64}()
fit!(m, X[training_rows, :], y[training_rows])

Test.@test m.X == X[training_rows, :]
Test.@test m.y == y[training_rows]
Test.@test all(m.X .== X[training_rows, :])
Test.@test all(m.y .== y[training_rows])
Test.@test !all(m.X .== 0)
Test.@test !all(m.y .== 0)

sanitize!(Model(m), Data(X), Data(y))

Test.@test m.X != X[training_rows, :]
Test.@test m.y != y[training_rows]
Test.@test !all(m.X .== X[training_rows, :])
Test.@test !all(m.y .== y[training_rows])
Test.@test all(m.X .== 0)
Test.@test all(m.y .== 0)

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

m = LinearModel{Float64}()
fit!(m, X[training_rows, :], y[training_rows])

Test.@test m.X == X[training_rows, :]
Test.@test m.y == y[training_rows]
Test.@test all(m.X .== X[training_rows, :])
Test.@test all(m.y .== y[training_rows])
Test.@test !all(m.X .== 0)
Test.@test !all(m.y .== 0)

sanitize!(Model(m), Data[Data(X), Data(y)])

Test.@test m.X != X[training_rows, :]
Test.@test m.y != y[training_rows]
Test.@test !all(m.X .== X[training_rows, :])
Test.@test !all(m.y .== y[training_rows])
Test.@test all(m.X .== 0)
Test.@test all(m.y .== 0)

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

m = LinearModel{Float64}()
fit!(m, X[training_rows, :], y[training_rows])

Test.@test m.X == X[training_rows, :]
Test.@test m.y == y[training_rows]
Test.@test all(m.X .== X[training_rows, :])
Test.@test all(m.y .== y[training_rows])
Test.@test !all(m.X .== 0)
Test.@test !all(m.y .== 0)

sanitize!(Model(m), Data([[[[[[], [[[[X], [y]]]], []]]]]]))

Test.@test m.X != X[training_rows, :]
Test.@test m.y != y[training_rows]
Test.@test !all(m.X .== X[training_rows, :])
Test.@test !all(m.y .== y[training_rows])
Test.@test all(m.X .== 0)
Test.@test all(m.y .== 0)

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

m = LinearModel{Float64}()
fit!(m, X[training_rows, :], y[training_rows])

Test.@test m.X == X[training_rows, :]
Test.@test m.y == y[training_rows]
Test.@test all(m.X .== X[training_rows, :])
Test.@test all(m.y .== y[training_rows])
Test.@test !all(m.X .== 0)
Test.@test !all(m.y .== 0)

sanitize!(Model(m), Data([[[[[[], [[[[[X, y]]]]], []]]]]]))

Test.@test m.X != X[training_rows, :]
Test.@test m.y != y[training_rows]
Test.@test !all(m.X .== X[training_rows, :])
Test.@test !all(m.y .== y[training_rows])
Test.@test all(m.X .== 0)
Test.@test all(m.y .== 0)

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

m = LinearModel{Float64}()
fit!(m, X[training_rows, :], y[training_rows])

Test.@test m.X == X[training_rows, :]
Test.@test m.y == y[training_rows]
Test.@test all(m.X .== X[training_rows, :])
Test.@test all(m.y .== y[training_rows])
Test.@test !all(m.X .== 0)
Test.@test !all(m.y .== 0)

sanitize!(Model(m), Data[Data(X[1:20]), Data(y[1:20])])

Test.@test m.X != X[training_rows, :]
Test.@test m.y != y[training_rows]
Test.@test !all(m.X .== X[training_rows, :])
Test.@test !all(m.y .== y[training_rows])
Test.@test all(m.X .== 0)
Test.@test all(m.y .== 0)


text:
module TestDPO
using Test
using Catlab.Graphs
using Catlab.Present
using Catlab.WiringDiagrams
using Catlab.CategoricalAlgebra.FinSets
using Catlab, Catlab.Theories, Catlab.CategoricalAlgebra

# Wiring diagrams
#################

A, B, C, D = [:A], [:B], [:C], [:D];
f,g,h,i = Box(:f, A, B), Box(:g, B, C),Box(:h, C, D), Box(:i, [:D,:B], [:D,:B]);

GWD = WiringDiagram([:A,:A],[:D,:B])
fv, gv, hv, fv2 = add_box!(GWD, f), add_box!(GWD, g), add_box!(GWD, h), add_box!(GWD, f);
add_wires!(GWD, Pair[
  (input_id(GWD),1) => (fv,1),
  (fv,1) => (gv,1),
  (gv,1) => (hv,1),
  (hv,1) => (output_id(GWD),1),
  (input_id(GWD),2) => (fv2,1),
  (fv2,1) => (output_id(GWD),2)]);

RWD = WiringDiagram([], []);
fv, iv, gv = add_box!(RWD, Box(:f,[],B)), add_box!(RWD, i), add_box!(RWD, Box(:g,B,[]));
add_wires!(RWD, Pair[(iv,1) => (iv, 1),
                     (fv, 1) => (iv, 2),
                     (iv, 2)=>(gv,1)]);

LWD = WiringDiagram([], []);
fv, gv = add_box!(LWD, Box(:f,[],B)), add_box!(LWD, Box(:g,B,[]));
add_wires!(LWD, Pair[
  (fv,1) => (gv,1),
]);

IWD = WiringDiagram([], []);
fv, gv = add_box!(IWD, Box(:f,[],[])), add_box!(IWD, Box(:g,[],[]));

XWD = WiringDiagram([:A,:A],[:D,:B])
fv, gv, hv = add_box!(XWD, f), add_box!(XWD, g), add_box!(XWD, h)
fv2, iv = add_box!(XWD, f), add_box!(XWD, i);
add_wires!(XWD, Pair[
  (input_id(XWD),1) => (fv,1),
  (iv,1) => (iv, 1), (fv, 1) => (iv, 2), (iv, 2)=>(gv,1),
  (gv,1) => (hv,1),
  (hv,1) => (output_id(XWD),1),
  (input_id(XWD),2) => (fv2,1),
  (fv2,1) => (output_id(XWD),2)]);

L=ACSetTransformation(IWD.diagram, LWD.diagram, Box=[1,2]);
R=ACSetTransformation(IWD.diagram, RWD.diagram, Box=[1,3]);
m=ACSetTransformation(LWD.diagram, GWD.diagram, Box=[1,2],InPort=[2],OutPort=[1],Wire=[1]);
@test is_isomorphic(XWD.diagram, rewrite_match(L,R,m))

# Graphs with attributes
########################

@present TheoryDecGraph(FreeSchema) begin
  E::Ob
  V::Ob
  src::Hom(E,V)
  tgt::Hom(E,V)

  X::AttrType
  dec::Attr(E,X)
end

@present TheoryLabelledDecGraph <: TheoryDecGraph begin
  label::Attr(V,X)
end

@acset_type LabelledDecGraph(TheoryLabelledDecGraph, index=[:src,:tgt])

aI2 = @acset LabelledDecGraph{String} begin
  V = 2
  E = 0
  label = ["a","b"]
end

aarr = @acset LabelledDecGraph{String} begin
  V = 2
  E = 1
  src = [1]
  tgt = [2]
  dec = ["e1"]
  label = ["a","b"]
end

abiarr = @acset LabelledDecGraph{String} begin
  V = 2
  E = 2
  src = [1,2]
  tgt = [2,1]
  dec = ["e1","rev_e1"]
  label = ["a","b"]
end

aspan = @acset LabelledDecGraph{String} begin
  V = 3
  E = 2
  src = [1,1]
  tgt = [2,3]
  dec = ["e1","e2"]
  label = ["a","b","c"]
end

expected = @acset LabelledDecGraph{String} begin
  V = 3
  E = 3
  src = [1,1,2]
  tgt = [2,3,1]
  dec = ["e1","e2","rev_e1"]
  label = ["a","b","c"]
end

L = ACSetTransformation(aI2, aarr, V=[1,2]);
R = ACSetTransformation(aI2, abiarr, V=[1,2]);
m = ACSetTransformation(aarr, aspan, V=[2,1], E=[1]);  # sends 'a'->'b' and 'b'->'a'

@test_throws ErrorException("ACSet colimit does not exist: label attributes a != b") rewrite_match(L,R,m)

m = ACSetTransformation(aarr, aspan, V=[1,2], E=[1]);

@test is_isomorphic(expected, rewrite_match(L, R, m))

# Graphs
########

# Example graphs
I2 = Graph(2)
I3 = Graph(3)
#   e1   e2
# 1 <- 2 -> 3
span = Graph(3)
add_edges!(span, [2,2],[1,3])
# 1 -> 2
arr = path_graph(Graph, 2)
# 1 <-> 2
biarr = Graph(2)
add_edges!(biarr, [1,2],[2,1])
# 1 -> 2 -> 3 -> 1
tri = cycle_graph(Graph, 3)
# 4 <- 1 -> 2 and 2 <- 3 -> 4
dispan = Graph(4)
add_edges!(dispan, [1,1,3,3],[2,4,2,4])

#      e1
#    1 -> 2
# e2 |  / ^
#    vv   | e4
#    4 -> 3
#      e5
squarediag = Graph(4)
add_edges!(squarediag, [1,1,2,3,4],[2,4,4,2,3])

# Add a reverse arrow to span
span_w_arrow = Graph(3)
add_edges!(span_w_arrow,[1,1,2],[2,3,1])

L = CSetTransformation(I2, arr, V=[1,2])
R = CSetTransformation(I2, biarr, V=[1,2])
m = CSetTransformation(arr, span, V=[2,1], E=[1])
@test is_isomorphic(span_w_arrow, rewrite_match(L, R, m))

# Remove apex of a subspan (top left corner of squarediag, leaves the triangle behind)
L = CSetTransformation(I2, span, V=[1,3])
m = CSetTransformation(span, squarediag, V=[2,1,4], E=[1,2])
@test is_isomorphic(tri, rewrite_match(L,id(I2),m))

# Remove self-edge using a *non-monic* match morphism
two_loops = Graph(2)
add_edges!(two_loops,[1,2],[1,2]) # ↻1   2↺
one_loop = Graph(2)
add_edges!(one_loop,[2],[2]) # 1   2↺

L = CSetTransformation(I2, arr, V=[1,2])
m = CSetTransformation(arr, two_loops, V=[1, 1], E=[1])
@test is_isomorphic(one_loop, rewrite_match(L,id(I2),m))

# Simplest non-trivial, non-monic exmaple
@present TheoryFinSet(FreeSchema) begin
  X::Ob
end
@acset_type FinSetType(TheoryFinSet)

I, L, G = [FinSetType() for _ in 1:3]
add_parts!(I,:X,2)
add_parts!(L,:X,1)
add_parts!(G,:X,1)
l = CSetTransformation(I,L,X=[1,1])
m = CSetTransformation(L,G,X=[1])
@test can_pushout_complement(l,m)
ik, kg = pushout_complement(l,m)
# There are 3 functions `ik` that make this a valid P.C.
# codom=1 with [1,1], codom=2 with [1,2] or [2,1]
K = codom(ik)
@test nparts(K, :X) == 1 # algorithm currently picks the first option

# Non-discrete interface graph. Non-monic matching
arr_loop= Graph(2)
add_edges!(arr_loop,[1,2],[2,2])  # 1->2↺
arrarr = Graph(2)
add_edges!(arrarr, [1,1],[2,2]) #  1⇉2
arrarr_loop = Graph(2)
add_edges!(arrarr_loop,[1,1,2],[2,2,2]) # 1⇉2↺
arr_looploop = Graph(2)
add_edges!(arr_looploop, [1,2,2],[2,2,2]) # 1-> ↻2↺

L = CSetTransformation(arr, arr, V=[1,2],E=[1]) # identity
R = CSetTransformation(arr, arrarr, V=[1,2], E=[1])
m = CSetTransformation(arr, arr_loop, V=[2,2], E=[2]) # NOT MONIC
@test is_isomorphic(arr_looploop, rewrite_match(L,R,m))

# only one monic match
@test is_isomorphic(arrarr_loop, rewrite(L, R, arr_loop, monic=true))

# two possible morphisms L -> squarediag, but both violate dangling condition
L = CSetTransformation(arr, span, V=[1,2], E=[1]);
m = CSetTransformation(span, squarediag, V=[2,1,4], E=[1,2]);
@test (:src, 5, 4) in dangling_condition(ComposablePair(L,m))

# violate id condition because two orphans map to same point
L = CSetTransformation(I2, biarr, V=[1,2]); # delete both arrows
m = CSetTransformation(biarr, arr_loop, V=[2,2], E=[2,2]);
@test (1, 2) in id_condition(ComposablePair(L[:E],m[:E]))[2]
L = CSetTransformation(arr, biarr, V=[1,2], E=[1]); # delete one arrow
@test 1 in id_condition(ComposablePair(L[:E],m[:E]))[1]

span_triangle = Graph(3); # 2 <- 1 -> 3 (with edge 2->3)
add_edges!(span_triangle,[1,1,2], [2,3,3]);

L = CSetTransformation(arr, tri, V=[1,2], E=[1]);
m = CSetTransformation(tri, squarediag, V=[2,4,3], E=[3,5,4]);
@test is_isomorphic(span_triangle, rewrite_match(L,id(arr),m))

k, g = pushout_complement(L, m); # get PO complement to do further tests

# the graph interface is equal to the final graph b/c we only delete things
@test is_isomorphic(span_triangle, codom(k))

# Check pushout properties 1: apex is the original graph
@test is_isomorphic(squarediag, ob(pushout(L, k))) # recover original graph

# Check pushout properties 2: the diagram commutes
Lm = compose(L,m);
kg = compose(k,g);
for I_node in 1:2
  @test Lm[:V](I_node) == kg[:V](I_node)
end
@test Lm[:E](1) == kg[:E](1)

# Check pushout properties 3: for every pair of unmatched things between K and L, they are NOT equal
for K_node in 1:3
  @test m[:V](3) != g[:V](K_node)
end

for K_edge in 2:3
  @test m[:E](3) != g[:E](K_edge)
end

# Undirected bipartite graphs
#############################

# 1 --- 1
#    /
# 2 --- 2

z_ = UndirectedBipartiteGraph()
add_vertices₁!(z_, 2)
add_vertices₂!(z_, 2)
add_edges!(z_, [1,2,2], [1,1,2])

line = UndirectedBipartiteGraph()
add_vertices₁!(line, 1)
add_vertices₂!(line, 2)
add_edges!(line, [1], [1])

parallel = UndirectedBipartiteGraph()
add_vertices₁!(parallel, 2)
add_vertices₂!(parallel, 2)
add_edges!(parallel, [1,2], [1,2])

merge = UndirectedBipartiteGraph()
add_vertices₁!(merge, 2)
add_vertices₂!(merge, 2)
add_edges!(merge, [1,2], [1,1])

Lspan = UndirectedBipartiteGraph()
add_vertices₁!(Lspan, 1)
add_vertices₂!(Lspan, 2)
add_edges!(Lspan, [1,1],[1,2])

I = UndirectedBipartiteGraph()
add_vertices₁!(I, 1)
add_vertices₂!(I, 2)

L = CSetTransformation(I, Lspan, V₁=[1], V₂=[1,2])
R = CSetTransformation(I, line, V₁=[1], V₂=[1,2])
m1 = CSetTransformation(Lspan, z_, V₁=[1], V₂=[1,2], E=[1, 2])
m2 = CSetTransformation(Lspan, z_, V₁=[1], V₂=[2,1], E=[2, 1])

@test is_isomorphic(parallel, rewrite_match(L, R, m1))
@test is_isomorphic(merge, rewrite_match(L, R, m2))

end


text:
# Code written entirely by Eric Alcaide.

# An implementation of Uninformed search methods algorithms
# in julia for a future package: SearchMethods.jl

""" Breadth First Search method algorithm. Expands always the shallowest node.
    Inputs:
        * root: the initial node. Any type.
        * tree: the search space. Any type.
        * goal_check: checks if the goal is met. Function.
        * expand: returns an array of new nodes. Function.
    Returns:
        * the goal if it's found in tree.
        * false (boolean value) otherwise.
"""
function bfs(root, tree, goal_check, expand)
    explored = Set([])
    frontier = [root]
    # While frontier is active, pop from left (FIFO)
    while !isempty(frontier)
    	node = shift!(frontier)
    	push!(explored, node)
    	# Check if goal found. Else => continue.
    	if goal_check(node) return node end
    	# Explore node
    	for neighbor in expand(node, tree)
    		if !in(neighbor, union(explored, frontier))
    			push!(frontier, neighbor)
    		end
    	end
    end
    # Return false if goal not found.
    return false
end

""" Depth First Search method algorithm. Expands always the deepest node.
    Inputs:
        * root: the initial node. Any type.
        * tree: the search space. Any type.
        * goal_check: checks if the goal is met. Function.
        * expand: returns an array of new nodes. Function.
    Returns:
        * the goal if it's found in tree.
        * false (boolean value) otherwise.
"""
function dfs(root, tree, goal_check, expand)
    explored = Set([])
    frontier = [root]
    # While frontier is active, pop from right (LIFO)
    while !isempty(frontier)
    	node = pop!(frontier)
    	push!(explored, node)
    	# Check if goal found. Else => continue.
    	if goal_check(node) return node end
    	# Explore node
    	for neighbor in expand(node, tree)
    		if !in(neighbor, union(explored, frontier))
    			push!(frontier, neighbor)
    		end
    	end
    end
    # Return false if goal not found.
    return false
end

""" Depth Limited Search method algorithm. Expands always the deepest node
    until maximum depth is reached.
    Inputs:
        * root: the initial node. Any type.
        * tree: the search space. Any type.
        * goal_check: checks if the goal is met. Function.
        * expand: returns an array of new nodes. Function.
        * depth: maximum depth. Integer.
    Returns:
        * the goal if it's found in tree.
        * false (boolean value) otherwise.
"""
function dls(root, tree, goal_check, expand, depth)
	explored = Set([])
	# Keep track of node, depth while node in frontier
    frontier = [[root, 0]]
    # While frontier is active, pop from right (LIFO)
    while !isempty(frontier)
    	node = pop!(frontier)
    	push!(explored, node[1])
    	# Check if goal found. Else => continue.
    	if goal_check(node[1]) return node[1] end
    	# Explore node only if depth-admissible
    	if node[end] < depth
        	for neighbor in expand(node[1], tree)
        		if !in(neighbor, union(explored, [f[1] for f in frontier]))
        			push!(frontier, [neighbor, node[end]+1])
        		end
        	end
        end
    end
    # Return false if goal not found.
    return false
end

""" Iterative Deepening Search method algorithm. Expands always the deepest node
    until maximum depth is reached. Repeat the process iteratively increasing
    maximum depth.
    Inputs:
        * root: the initial node. Any type.
        * tree: the search space. Any type.
        * goal_check: checks if the goal is met. Function.
        * expand: returns an array of new nodes. Function.
        * depth: maximum depth. Integer.
        * it: iteratively increase depth by it. Integer.
    Returns:
        * the goal if it's found in tree.
        * false (boolean value) otherwise.
"""
function ids(root, tree, goal_check, expand, depth; it = 1)
	final = false
	# Normal dfs but check if bottom of tree. If not, increase depth by it.
	while final == false
		final = true
		explored = Set([])
    	# Keep track of node, depth while node in frontier
        frontier = [[root, 0]]
        # While frontier is active, pop from right (LIFO)
        while !isempty(frontier)
        	node = pop!(frontier)
        	push!(explored, node[1])
        	# Check if goal found. Else => continue.
        	if goal_check(node[1]) return node[1] end
        	# Explore node only if depth-admissible
        	if node[end] < depth
	        	for neighbor in expand(node[1], tree)
	        		if !in(neighbor, union(explored, [f[1] for f in frontier]))
	        			push!(frontier, [neighbor, node[end]+1])
	        		end
	        	end
	        # Check for tree end. Continue if not end yet.
	        elseif node[end] == depth && final == true
	        	if !isempty(expand(node[1], tree))
	        		final = false
	        	end
	        end
        end
        # Update depth parameter
        depth += it
	end
	# if not goal found
	return false
end

""" Uniform Cost Search method algorithm. Expands always the least cost node.
    Inputs:
        * root: the initial node. Any type.
        * tree: the search space. Any type.
        * goal_check: checks if the goal is met. Function.
        * expand: returns an array of new nodes. Function.
        * cost: returns the cost of expanding a node. Function.
    Returns:
        * the goal if it's found in tree.
        * false (boolean value) otherwise.
"""
function ucs(root, tree, goal_check, expand, cost)
	explored = Set([])
	# Keep track of node, depth while node in frontier
    frontier = PriorityQueue()
    enqueue!(frontier, root, 0)
    # While frontier is active, pop least cost node.
    while !isempty(frontier)
    	# Peek least cost value and cost. Eliminate from frontier.
    	node, prev_cost = map(Any, peek(frontier))
    	dequeue!(frontier)
    	# Check if goal found. Else => continue.
    	if goal_check(node) return node end
    	# Explore node
    	for neighbor in expand(node, tree)
    		if !in(neighbor, union(explored, keys(frontier)))
    			enqueue!(frontier, neighbor, cost(node,neighbor) + prev_cost)
    		end
    	end
    end
    # Return false if goal not found.
    return false
end


text:
export SharedProbCircuit, SharedProbLeafNode, SharedProbInnerNode, SharedProbLiteralNode,
SharedMulNode, SharedSumNode, num_components

#####################
# Probabilistic circuits which share the same structure
#####################

"""
Root of the shared probabilistic circuit node hierarchy
"""
abstract type SharedProbCircuit <: ProbCircuit end

"""
A shared probabilistic leaf node
"""
abstract type SharedProbLeafNode <: SharedProbCircuit end

"""
A shared probabilistic inner node
"""
abstract type SharedProbInnerNode <: SharedProbCircuit end

"""
A shared probabilistic literal node
"""
mutable struct SharedProbLiteralNode <: SharedProbLeafNode
    literal::Lit
    data
    counter::UInt32
    SharedProbLiteralNode(l) = new(l, nothing, 0)
end

"""
A shared probabilistic multiplcation node
"""
mutable struct SharedMulNode <: SharedProbInnerNode
    children::Vector{<:SharedProbCircuit}
    data
    counter::UInt32
    SharedMulNode(children) = new(children, nothing, 0)
end

"""
A shared probabilistic summation node
"""
mutable struct SharedSumNode <: SharedProbInnerNode
    children::Vector{<:SharedProbCircuit}
    log_probs::Matrix{Float64}
    data
    counter::UInt32
    SharedSumNode(children, n_mixture) = begin
        new(children, log.(ones(Float64, length(children), n_mixture) / length(children)), nothing, 0)
    end
end

#####################
# traits
#####################

import LogicCircuits.GateType # make available for extension
@inline GateType(::Type{<:SharedProbLiteralNode}) = LiteralGate()
@inline GateType(::Type{<:SharedMulNode}) = ⋀Gate()
@inline GateType(::Type{<:SharedSumNode}) = ⋁Gate()

#####################
# methods
#####################

import LogicCircuits: children # make available for extension
@inline children(n::SharedProbInnerNode) = n.children

@inline num_parameters_node(n::SharedSumNode) = length(n.log_probs)

"How many components are mixed together in this shared circuit?"
@inline num_components(n::SharedSumNode) = size(n.log_probs,2)
@inline num_components(n::SharedMulNode) = begin
    @assert length(n.children) > 0 "Invalid SharedMulNode (it has no child node)."
    num_components(n.children[1])
end

#####################
# constructors and conversions
#####################

function multiply(arguments::Vector{<:SharedProbCircuit};
    reuse=nothing)
    @assert length(arguments) > 0
    reuse isa SharedMulNode && children(reuse) == arguments && return reuse
    return SharedMulNode(arguments)
end

function summate(arguments::Vector{<:SharedProbCircuit}, num_components=0;
       reuse=nothing)
    @assert length(arguments) > 0
    reuse isa SharedSumNode && children(reuse) == arguments && return reuse
    return SharedSumNode(arguments, num_components) # unknwown number of components; resize later
end

compile(::Type{<:SharedProbCircuit}, l::Lit) =
    SharedProbLiteralNode(l)

function compile(::Type{<:SharedProbCircuit}, circuit::LogicCircuit, num_components::Int)
    f_con(n) = error("Cannot construct a probabilistic circuit from constant leafs: first smooth and remove unsatisfiable branches.")
    f_lit(n) = compile(SharedProbCircuit, literal(n))
    f_a(_, cns) = multiply(cns)
    f_o(n, cns) = begin
        prod = summate(cns, num_components)
        if n isa ProbCircuit
            prod.log_probs .= n.log_probs
        end
        prod
    end
    foldup_aggregate(circuit, f_con, f_lit, f_a, f_o, SharedProbCircuit)
end

import LogicCircuits: fully_factorized_circuit #extend

function fully_factorized_circuit(::Type{<:SharedProbCircuit}, n::Int, num_mix::Int)
    ff_logic_circuit = fully_factorized_circuit(PlainLogicCircuit, n)
    compile(SharedProbCircuit, ff_logic_circuit, num_mix)
end