text:
struct QueryException <: Exception
	msg::String
	context::Any
	QueryException(msg::AbstractString, context=nothing) = new(msg, context)
end

function Base.showerror(io::IO, ex::QueryException)
	print(io, "QueryException: $(ex.msg)")
	if ex.context !== nothing
		print(io, " at $(ex.context)")
	end
end

function helper_namedtuples_replacement(ex)
	return postwalk(ex) do x
		if x isa Expr && x.head==:braces
			new_ex = Expr(:tuple, x.args...)

			for (j,field_in_NT) in enumerate(new_ex.args)
				if isa(field_in_NT, Expr) && field_in_NT.head==:.
					name_to_use = field_in_NT.args[2].value
					new_ex.args[j] = Expr(:(=), name_to_use, field_in_NT)
				elseif isa(field_in_NT, Symbol)
					new_ex.args[j] = Expr(:(=), field_in_NT, field_in_NT)
				end
			end

			return new_ex
		else
			return x
		end
	end
end

function helper_replace_anon_func_syntax(ex)
	if !(isa(ex, Expr) && ex.head==:->)
		new_symb = gensym()
		new_symb2 = gensym()
		two_args = false
		new_ex = postwalk(ex) do x
			if isa(x, Symbol)
				if x==:_
					return new_symb
				elseif x==:__
					two_args = true
					return new_symb2
				else
					return x
				end
			else
				return x
			end
		end

		if two_args
			return :(($new_symb, $new_symb2) -> $(new_ex) )
		else
			return :($new_symb -> $(new_ex) )
		end
	else
		return ex
	end
end

function query_expression_translation_phase_A(qe)
	i = 1
	while i<=length(qe)
		clause = qe[i]
		if ismacro(clause, "@left_outer_join")
			clause.args[1] = Symbol("@join")
			temp_name = gensym()
			x1 = clause.args[3].args[2]
			push!(clause.args, :into)
			push!(clause.args, temp_name)
			nested_from = :(@from $x1 in $(QueryOperators.default_if_empty)($temp_name))
			insert!(qe,i+1,nested_from)
		end
		i+=1
	end
end

function query_expression_translation_phase_B(qe)
	i = 1
	while i<=length(qe)
		qe[i] = helper_namedtuples_replacement(qe[i])
		clause = qe[i]

		# for l=length(clause.args):-1:1
		# 	if clause.args[l] isa LineNumberNode
		# 		deleteat!(clause.args,l)
		# 	end
		# end

		if i==1 && ismacro(clause, "@from")
			# Handle the case of a nested query. We are essentially detecting
			# here that the subquery starts with convert2nullable
			# and then we don't escape things.
			subq = clause.args[3].args[3]
			if isa(subq, Expr) && subq.head==:call && isa(subq.args[1],Expr) && subq.args[1].head==:. && subq.args[1].args[1]==:QueryOperators
				clause.args[3].args[3] = :($(QueryOperators.query)($(subq)))
			elseif !(isa(subq, Expr) && subq.head==:macrocall && isa(subq.args[1],Expr) && subq.args[1].head==:. && subq.args[1].args[1]==:QueryOperators)
				clause.args[3].args[3] = :($(QueryOperators.query)($(esc(subq))))
			end
		elseif ismacro(clause, "@from")
			clause.args[3].args[3] = :($(QueryOperators.query)($(clause.args[3].args[3])))
		elseif ismacro(clause, "@join")
			clause.args[3].args[3] = :($(QueryOperators.query)($(esc(clause.args[3].args[3]))))
		end
		i+=1
	end
end

function query_expression_translation_phase_1(qe)
	done = false
	while !done
		group_into_index = findfirst(i->ismacro(i, "@group", 6) && i.args[6]==:into, qe)
		select_into_index = findfirst(i->ismacro(i, "@select", 4) && i.args[4]==:into, qe)
		if length(qe)>=2 && ismacro(qe[1], "@from") && group_into_index!==nothing
			x = qe[group_into_index].args[7]

			sub_query = Expr(:block, qe[1:group_into_index]...)

			deleteat!(sub_query.args[end].args,7)
			deleteat!(sub_query.args[end].args,6)

			translate_query(sub_query)
			length(sub_query.args)==1 || throw(QueryException("@group ... into subquery too long", sub_query))

			qe[1] = :( @from $x in $(sub_query.args[1]) )
			deleteat!(qe, 2:group_into_index)
		elseif length(qe)>=2 && ismacro(qe[1], "@from") && select_into_index!==nothing
			x = qe[select_into_index].args[5]

			sub_query = Expr(:block, qe[1:select_into_index]...)
			deleteat!(sub_query.args[end].args,5)
			deleteat!(sub_query.args[end].args,4)

			translate_query(sub_query)

			length(sub_query.args)==1 || throw(QueryException("@select ... into subquery too long", sub_query))

			qe[1] = :( @from $x in $(sub_query.args[1]) )
			deleteat!(qe, 2:select_into_index)
		else
			done = true
		end
	end
end

function query_expression_translation_phase_3(qe)
	done = false
	while !done
		if length(qe)>=2 && ismacro(qe[1], "@from") && ismacro(qe[2], "@select") && qe[1].args[3].args[2]==qe[2].args[3]
			x = qe[1].args[3].args[2]
			e = qe[1].args[3].args[3]

			qe[1] = :( QueryOperators.@map($e,x->x) )
			deleteat!(qe,2)
		else
			done = true
		end
	end
end

function query_expression_translation_phase_4(qe)
	done = false
	while !done
		if length(qe)>=3 && ismacro(qe[1], "@from") && ismacro(qe[2], "@from") && ismacro(qe[3], "@select")
			x1 = qe[1].args[3].args[2]
			x2 = qe[2].args[3].args[2]
			e1 = qe[1].args[3].args[3]
			e2 = qe[2].args[3].args[3]
			v = qe[3].args[3]

			f_collection_selector = Expr(:->, x1, e2)
			f_result_selector = Expr(:->, Expr(:tuple,x1,x2), v)

			qe[1] = :( QueryOperators.@mapmany($e1, $(esc(f_collection_selector)), $(esc(f_result_selector))) )
			deleteat!(qe,3)
			deleteat!(qe,2)
		elseif length(qe)>=3 && ismacro(qe[1], "@from") && ismacro(qe[2], "@from")
			x1 = qe[1].args[3].args[2]
			x2 = qe[2].args[3].args[2]
			e1 = qe[1].args[3].args[3]
			e2 = qe[2].args[3].args[3]

			f_collection_selector = Expr(:->, x1, e2)
			f_result_selector = Expr(:->, Expr(:tuple,x1,x2), :(($x1=$x1,$x2=$x2)))

			qe[1].args[3].args[2] = Expr(:transparentidentifier, gensym(:t), x1, x2)
			qe[1].args[3].args[3] = :( QueryOperators.@mapmany($e1, $(esc(f_collection_selector)), $(esc(f_result_selector))) )
			deleteat!(qe,2)
		elseif length(qe)>=3 && ismacro(qe[1], "@from") && ismacro(qe[2], "@let")
			x = qe[1].args[3].args[2]
			e = qe[1].args[3].args[3]
			y = qe[2].args[3].args[1]
			f = qe[2].args[3].args[2]

			f_selector = Expr(:->, x, :(($x=$x,$y=$f)))

			qe[1].args[3].args[2] = Expr(:transparentidentifier, gensym(:t), x, y)
			qe[1].args[3].args[3] = :( QueryOperators.@map($e,$(esc(f_selector))) )
			deleteat!(qe,2)
		elseif length(qe)>=3 && ismacro(qe[1], "@from") && ismacro(qe[2], "@where")
			x = qe[1].args[3].args[2]
			e = qe[1].args[3].args[3]
			f = qe[2].args[3]

			f_condition = Expr(:->, x, f)

			qe[1].args[3].args[3] = :( QueryOperators.@filter($e,$(esc(f_condition))) )
			deleteat!(qe,2)
		elseif length(qe)>=3 && ismacro(qe[1], "@from") && ismacro(qe[2], "@join", 6) && ismacro(qe[3], "@select")
			outer = qe[1].args[3].args[3]
			inner = qe[2].args[3].args[3]
			outer_range_var = qe[1].args[3].args[2]
			inner_range_var = qe[2].args[3].args[2]
			f_outer_key = Expr(:->, outer_range_var, qe[2].args[5])
			f_inner_key = Expr(:->, inner_range_var, qe[2].args[7])
			f_result = Expr(:->, Expr(:tuple,outer_range_var,inner_range_var), qe[3].args[3])
			qe[1] = :(
				QueryOperators.@join($outer, $inner, $(esc(f_outer_key)), $(esc(f_inner_key)), $(esc(f_result)))
				)

			deleteat!(qe,3)
			deleteat!(qe,2)
		elseif length(qe)>=3 && ismacro(qe[1], "@from") && ismacro(qe[2], "@join", 6)
			e1 = qe[1].args[3].args[3]
			e2 = qe[2].args[3].args[3]
			x1 = qe[1].args[3].args[2]
			x2 = qe[2].args[3].args[2]
			k1 = qe[2].args[5]
			k2 = qe[2].args[7]
			f_outer_key = Expr(:->, x1, k1)
			f_inner_key = Expr(:->, x2, k2)
			f_result = Expr(:->, Expr(:tuple,x1,x2), :(($x1=$x1,$x2=$x2)) )

			qe[1].args[3].args[2] = Expr(:transparentidentifier, gensym(:t), x1, x2)
			qe[1].args[3].args[3] = :( QueryOperators.@join($e1,$e2,$(esc(f_outer_key)), $(esc(f_inner_key)), $(esc(f_result))) )
			deleteat!(qe,2)
		elseif length(qe)>=3 && ismacro(qe[1], "@from") && ismacro(qe[2], "@join", 8) && ismacro(qe[3], "@select")
			e1 = qe[1].args[3].args[3]
			e2 = qe[2].args[3].args[3]
			x1 = qe[1].args[3].args[2]
			x2 = qe[2].args[3].args[2]
			k1 = qe[2].args[5]
			k2 = qe[2].args[7]
			g = qe[2].args[9]
			v = qe[3].args[3]
			f_outer_key = Expr(:->, x1, k1)
			f_inner_key = Expr(:->, x2, k2)
			f_result = Expr(:->, Expr(:tuple,x1,g), v)
			qe[1] = :( QueryOperators.@groupjoin($e1, $e2, $(esc(f_outer_key)), $(esc(f_inner_key)), $(esc(f_result))) )

			deleteat!(qe,3)
			deleteat!(qe,2)
		elseif length(qe)>=3 && ismacro(qe[1], "@from") && ismacro(qe[2], "@join", 8)
			e1 = qe[1].args[3].args[3]
			e2 = qe[2].args[3].args[3]
			x1 = qe[1].args[3].args[2]
			x2 = qe[2].args[3].args[2]
			k1 = qe[2].args[5]
			k2 = qe[2].args[7]
			g = qe[2].args[9]
			f_outer_key = Expr(:->, x1, k1)
			f_inner_key = Expr(:->, x2, k2)
			f_result = Expr(:->, Expr(:tuple,x1,g), :(($x1=$x1,$g=$g)) )

			qe[1].args[3].args[2] = Expr(:transparentidentifier, gensym(:t), x1, g)
			qe[1].args[3].args[3] = :( QueryOperators.@groupjoin($e1,$e2,$(esc(f_outer_key)), $(esc(f_inner_key)), $(esc(f_result))) )
			deleteat!(qe,2)
		elseif length(qe)>=3 && ismacro(qe[1], "@from") && ismacro(qe[2], "@orderby")
			e = qe[1].args[3].args[3]
			x = qe[1].args[3].args[2]
			ks = []
			if isa(qe[2].args[3], Expr) && qe[2].args[3].head==:tuple
				for sort_clause in qe[2].args[3].args
					if isa(sort_clause, Expr) && iscall(sort_clause, :descending)
						k = sort_clause.args[2]
						direction = :descending
					elseif isa(sort_clause, Expr) && iscall(sort_clause, :ascending)
						k = sort_clause.args[2]
						direction = :ascending
					else
						k = sort_clause
						direction = :ascending
					end
					push!(ks, (k, direction))
				end
			else
				if isa(qe[2].args[3], Expr) && iscall(qe[2].args[3], :descending)
					k = qe[2].args[3].args[2]
					direction = :descending
				elseif isa(qe[2].args[3], Expr) && iscall(qe[2].args[3], :ascending)
					k = qe[2].args[3].args[2]
					direction = :ascending
				else
					k = qe[2].args[3]
					direction = :ascending
				end
				push!(ks, (k, direction))
			end

			for (i,sort_clause) in enumerate(ks)
				f_condition = Expr(:->, x, sort_clause[1])

				if sort_clause[2]==:ascending
					if i==1
						qe[1].args[3].args[3] = :( QueryOperators.@orderby($e,$(esc(f_condition))) )
					else
						qe[1].args[3].args[3] = :( QueryOperators.@thenby($(qe[1].args[3].args[3]),$(esc(f_condition))) )
					end
				elseif sort_clause[2]==:descending
					if i==1
						qe[1].args[3].args[3] = :( QueryOperators.@orderby_descending($e,$(esc(f_condition))) )
					else
						qe[1].args[3].args[3] = :( QueryOperators.@thenby_descending($(qe[1].args[3].args[3]),$(esc(f_condition))) )
					end
				end
			end
			deleteat!(qe,2)
		else
			done = true
		end
	end
end

function query_expression_translation_phase_5(qe)
	i = 1
	while i<=length(qe)
		clause = qe[i]
		if ismacro(clause, "@select")
			from_clause = qe[i-1]
			ismacro(from_clause, "@from") || throw(QueryException("Phase 5: expected @from before @select", from_clause))
			range_var = from_clause.args[3].args[2]
			source = from_clause.args[3].args[3]
			if clause.args[3]==range_var
				qe[i-1] = source
			else
				func_call = Expr(:->, range_var, clause.args[3])
				qe[i-1] = :( QueryOperators.@map($source, $(esc(func_call))) )
			end
			deleteat!(qe,i)
		else
			i+=1
		end
	end
end

function query_expression_translation_phase_6(qe)
	done = false
	while !done
		if length(qe)>=2 && ismacro(qe[1], "@from") && ismacro(qe[2], "@group")
			e = qe[1].args[3].args[3]
			x = qe[1].args[3].args[2]
			v = qe[2].args[3]
			k = qe[2].args[5]

			f_elementSelector = Expr(:->, x, k)
			f_resultSelector = Expr(:->, x, v)

			if v==x
				qe[1] = :( QueryOperators.@groupby_simple($e, $(esc(f_elementSelector))) )
			else
				qe[1] = :( QueryOperators.@groupby($e, $(esc(f_elementSelector)), $(esc(f_resultSelector))) )
			end
			deleteat!(qe,2)
		else
			done = true
		end
	end
end

# Phase 7

function replace_transparent_identifier_in_anonym_func(ex::Expr, names_to_put_in_scope)
	for (i,child_ex) in enumerate(ex.args)
		if isa(child_ex, Expr)
			replace_transparent_identifier_in_anonym_func(child_ex, names_to_put_in_scope)
		elseif isa(child_ex, Symbol) && !(ex.head == :(=) && i==1) # Don't replace names if they appear to the left of an assignment
			index_of_name = findfirst(j->child_ex==j[2], names_to_put_in_scope)
			if index_of_name!==nothing && !(ex.head==Symbol("=>") && i==1)
				ex.args[i] = Expr(:., names_to_put_in_scope[index_of_name][1], QuoteNode(child_ex))
			end
		end
	end
end

# translates (:a, :((b.c).d)) to :(((a.b).c).d)
function shift_access!(sym::Symbol, ex::Expr)
	while ex.args[1] isa Expr
		ex = ex.args[1]
	end
	ex.args[1] = Expr(:., sym, QuoteNode(ex.args[1]))
	return
end

function find_names_to_put_in_scope(ex::Expr)
	names = []
	for child_ex in ex.args[2:end]
		if isa(child_ex,Expr) && child_ex.head==:transparentidentifier
			child_names = find_names_to_put_in_scope(child_ex)
			for child_name in child_names
				c = child_name[1]
				if c isa Symbol
					xx = (Expr(:., ex.args[1], QuoteNode(c)), child_name[2])
				elseif c.args[1] isa Symbol
					xx = (Expr(:., Expr(:., ex.args[1], QuoteNode(c.args[1])), c.args[2]), child_name[2])
				else
					shift_access!(ex.args[1], c)
					xx = (c, child_name[2])
				end
				push!(names, xx)
			end
		elseif isa(child_ex, Symbol)
			push!(names,(ex.args[1],child_ex))
		elseif isa(child_ex, Expr) && child_ex.head==:.
			push!(names,(ex.args[1],child_ex.args[2].value))
		else
			throw(QueryException("identifier expected", child_ex))
		end
	end
	return names
end

function find_and_translate_transparent_identifier(ex::Expr)
	# First expand any transparent identifiers in lambdas
	if ex.head==:-> && isa(ex.args[1], Expr) && ex.args[1].head==:transparentidentifier
		names_to_put_in_scope = find_names_to_put_in_scope(ex.args[1])
		ex.args[1] = ex.args[1].args[1]
		replace_transparent_identifier_in_anonym_func(ex, names_to_put_in_scope)
	elseif ex.head==:-> && isa(ex.args[1], Expr) && ex.args[1].head==:tuple
		names_to_put_in_scope = []
		for (i, child_ex) in enumerate(ex.args[1].args)
			if isa(child_ex, Expr) && child_ex.head==:transparentidentifier
				append!(names_to_put_in_scope, find_names_to_put_in_scope(child_ex))
				ex.args[1].args[i] = child_ex.args[1]
			end
		end
		replace_transparent_identifier_in_anonym_func(ex, names_to_put_in_scope)
	end

	for (i,child_ex) in enumerate(ex.args)
		if isa(child_ex, Expr) && child_ex.head==:transparentidentifier
			ex.args[i] = child_ex.args[1]
		elseif isa(child_ex, Expr)
			find_and_translate_transparent_identifier(child_ex)
		end
	end
end

function query_expression_translation_phase_7(qe)
	for clause in qe
		isa(clause, Expr) && find_and_translate_transparent_identifier(clause)
	end
end

function query_expression_translation_phase_D(qe)
	i = 1
	while i<=length(qe)
		clause = qe[i]
		if ismacro(clause, "@collect")
			previous_clause = qe[i-1]
			if length(clause.args)==2
				qe[i-1] = :( collect($previous_clause) )
			else
				qe[i-1] = :( $(esc(clause.args[3]))($previous_clause) )
			end
			deleteat!(qe,i)
		else
			i+=1
		end
	end
end

function translate_query(body)
	debug_output = false

	debug_output && println("AT START")
	debug_output && println(body)

	query_expression_translation_phase_1(body.args)
	debug_output && println("AFTER 1")
	debug_output && println(body)

	query_expression_translation_phase_A(body.args)
	debug_output && println("AFTER A")
	debug_output && println(body)

	query_expression_translation_phase_B(body.args)
	debug_output && println("AFTER B")
	debug_output && println(body)

	query_expression_translation_phase_3(body.args)
	debug_output && println("AFTER 3")
	debug_output && println(body)

	query_expression_translation_phase_4(body.args)
	debug_output && println("AFTER 4")
	debug_output && println(body)

	query_expression_translation_phase_5(body.args)
	debug_output && println("AFTER 5")
	debug_output && println(body)

	query_expression_translation_phase_6(body.args)
	debug_output && println("AFTER 6")
	debug_output && println(body)

	query_expression_translation_phase_7(body.args)
	debug_output && println("AFTER 7")
	debug_output && println(body)

	query_expression_translation_phase_D(body.args)
	debug_output && println("AFTER D")
	debug_output && println(body)
end


text:
using BandedMatrices
using LinearAlgebra
using Test

include("test_banded.jl")
include("test_subarray.jl")
include("test_gbmm.jl")
include("test_broadcasting.jl")
include("test_indexing.jl")
include("test_bandedlu.jl")
include("test_bandedqr.jl")
include("test_symbanded.jl")
include("test_tribanded.jl")
include("test_interface.jl")
include("test_miscs.jl")

ext: jl


text:
using Gen
import MacroTools

normalize(ex) = MacroTools.prewalk(MacroTools.rmlines, ex)

# dynamic
@testset "tilde syntax smoke test (dynamic)" begin
    @gen function foo(a, b)
        p ~ beta(a, b)
        return ({:coin => 1} ~ bernoulli(p)) + ({:coin => 2} ~ bernoulli(p))
    end

    trace = simulate(foo, (1, 1))
    @test (trace[:coin => 1] + trace[:coin => 2]) == get_retval(trace)
    @test get_score(trace) == sum([(trace[:coin => i] ? log(trace[:p]) : log(1-trace[:p])) for i=1:2])
end

@testset "tilde syntax desugars as expected (dynamic)" begin
    expected = normalize(:(
    @gen function foo()
        x = @trace(normal(0, 1), :x)
        y = @trace(normal(0, 1), :y => :z)
        z = @trace(bar())
    end))

    actual = normalize(Gen.desugar_tildes(:(
    @gen function foo()
        x ~ normal(0, 1)
        y = ({:y => :z} ~ normal(0, 1))
        z = ({*} ~ bar())
    end)))

    @test actual == expected
end

# static

@testset "tilde syntax smoke test (static)" begin
    @gen (static) function bar(r)
        a ~ normal(0, 1)
        return r
    end

    @gen (static) function foo()
        x ~ bar(1)
        yz = {:y => :z} ~ bar(2)
        u ~ normal(0, 1)
        {:v => :w} ~ normal(0, 1)
        retrate = 7
        return {:ret} ~ poisson(retrate)
    end

    Gen.load_generated_functions()

    trace = simulate(foo, ())

    # random choices
    @test trace[:u] isa Float64
    @test trace[:v => :w] isa Float64
    @test trace[:x => :a] isa Float64
    @test trace[:y => :z => :a] isa Float64

    # auxiliary state
    @test trace[:x] == 1
    @test trace[:y => :z] == 2

    # return value
    @test trace[] == trace[:ret]
end

@testset "tilde syntax desugars as expected (static)" begin
expected = normalize(:(
@gen (static) function foo()
    x = @trace(normal(0, 1), :x)
    y = @trace(normal(0, 1), :y)
end))

actual = normalize(Gen.desugar_tildes(:(
@gen (static) function foo()
    x ~ normal(0, 1)
    y = ({:y} ~ normal(0, 1))
end)))

@test actual == expected
end

max_stars_repo_name: ali-ramadhan/Gen

max_issues_repo_name: ali-ramadhan/Gen

max_forks_repo_name: ali-ramadhan/Gen


text:
using Test
using BED

using Distributions
using Documenter
using FormatSpecimens
using GenomicFeatures

import Random
import ColorTypes: RGB
import FixedPointNumbers: N0f8
import BGZFStreams

# Generate an array of n random Interval{Int} object. With sequence names
# samples from seqnames, and intervals drawn to lie in [1, maxpos].
function random_intervals(seqnames, maxpos::Int, n::Int)
    seq_dist = Categorical(length(seqnames))
    strand_dist = Categorical(2)
    length_dist = Normal(1000, 1000)
    intervals = Vector{Interval{Int}}(undef, n)
    for i in 1:n
        intlen = maxpos
        while intlen >= maxpos || intlen <= 0
            intlen = ceil(Int, rand(length_dist))
        end
        first = rand(1:maxpos-intlen)
        last = first + intlen - 1
        strand = rand(strand_dist) == 1 ? STRAND_POS : STRAND_NEG
        intervals[i] = Interval{Int}(seqnames[rand(seq_dist)],
                                     first, last, strand, i)
    end
    return intervals
end

@testset "BED" begin

    @testset "Record" begin
        record = BED.Record("chr1\t17368\t17436")
        @test BED.chrom(record) == "chr1"
        @test BED.chromstart(record) === 17369
        @test BED.chromend(record) === 17436
        @test !BED.hasname(record)
        @test !BED.hasscore(record)
        @test !BED.hasstrand(record)
        @test !BED.hasthickstart(record)
        @test !BED.hasthickend(record)
        @test !BED.hasitemrgb(record)
        @test !BED.hasblockcount(record)
        @test !BED.hasblocksizes(record)
        @test !BED.hasblockstarts(record)
        @test startswith(repr(record), "BED.Record:\n")
        @test string(record) == "chr1\t17368\t17436"

        record = BED.Record("chrXIII\t854794\t855293\tYMR292W\t0\t+\t854794\t855293\t0\t2\t22,395,\t0,104,")
        @test BED.chrom(record) == "chrXIII"
        @test BED.chromstart(record) === 854795
        @test BED.chromend(record) === 855293
        @test BED.name(record) == "YMR292W"
        @test BED.score(record) === 0
        @test BED.strand(record) === STRAND_POS
        @test BED.thickstart(record) === 854795
        @test BED.thickend(record) === 855293
        @test BED.itemrgb(record) === RGB(0, 0, 0)
        @test BED.blockcount(record) === 2
        @test BED.blocksizes(record) == [22, 395]
        @test BED.blockstarts(record) == [1, 105]
        @test startswith(repr(record), "BED.Record:\n")
        @test string(record) == "chrXIII\t854794\t855293\tYMR292W\t0\t+\t854794\t855293\t0\t2\t22,395,\t0,104,"

        record = BED.Record("chrX\t151080532\t151081699\tCHOCOLATE1\t0\t-\t151080532\t151081699\t255,127,36")
        @test BED.chrom(record) == "chrX"
        @test BED.chromstart(record) === 151080533
        @test BED.chromend(record) === 151081699
        @test BED.name(record) == "CHOCOLATE1"
        @test BED.score(record) === 0
        @test BED.strand(record) === STRAND_NEG
        @test BED.thickstart(record) === 151080533
        @test BED.thickend(record) === 151081699
        @test BED.itemrgb(record) === RGB(map(x -> reinterpret(N0f8, UInt8(x)), (255, 127, 36))...)
        @test !BED.hasblockcount(record)
        @test !BED.hasblocksizes(record)
        @test !BED.hasblockstarts(record)
        @test startswith(repr(record), "BED.Record:\n")
        @test string(record) == "chrX\t151080532\t151081699\tCHOCOLATE1\t0\t-\t151080532\t151081699\t255,127,36"
    end

    function check_bed_parse(filename)
        # Reading from a stream
        for interval in BED.Reader(open(filename))
        end

        # Reading from a regular file
        for interval in open(BED.Reader, filename)
        end

        # in-place parsing
        stream = open(BED.Reader, filename)
        entry = eltype(stream)()
        while !eof(stream)
            read!(stream, entry)
        end
        close(stream)

        # Check round trip
        output = IOBuffer()
        writer = BED.Writer(output)
        expected_entries = BED.Record[]
        for interval in open(BED.Reader, filename)
            write(writer, interval)
            push!(expected_entries, interval)
        end
        flush(writer)

        seekstart(output)
        read_entries = BED.Record[]
        for interval in BED.Reader(output)
            push!(read_entries, interval)
        end

        return expected_entries == read_entries
    end

    dir_bed = path_of_format("BED")

    for specimen in list_valid_specimens("BED")

        if hastag(specimen, "gzip")
            # skip compressed files
            continue
        end

        filepath = joinpath(dir_bed, filename(specimen))

        @test check_bed_parse(filepath)

    end

    for specimen in list_invalid_specimens("BED")

        if hastag(specimen, "gzip")
            # skip compressed files
            continue
        end

        filepath = joinpath(dir_bed, filename(specimen))

        @test_throws Exception check_bed_parse(filepath)
    end

    #=
    Testing strategy: there are two entirely separate intersection algorithms for IntervalCollection and IntervalStream.
    Here we test them both by checking that they agree by generating and intersecting random BED files.
    =#

    function check_intersection(filename_a, filename_b)
        ic_a = IntervalCollection{BED.Record}()
        open(BED.Reader, filename_a) do reader
            for record in reader
                push!(ic_a, Interval(record))
            end
        end

        ic_b = IntervalCollection{BED.Record}()
        open(BED.Reader, filename_b) do reader
            for record in reader
                push!(ic_b, Interval(record))
            end
        end

        # This is refactored out to close streams
        fa = open(BED.Reader, filename_a)
        fb = open(BED.Reader, filename_b)
        xs = sort(collect(eachoverlap(fa, fb)))
        close(fa)
        close(fb)

        ys = sort(collect(eachoverlap(ic_a, ic_b)))

        return xs == ys
    end

    function write_intervals(filename, intervals)
        open(filename, "w") do out
            for interval in sort(intervals)
                println(out, interval.seqname, "\t", interval.first - 1, "\t", interval.last, "\t", interval.metadata, "\t", 1000, "\t", interval.strand)
            end
        end

    end

    n = 10000
    Random.seed!(1234)
    intervals_a = random_intervals(["one", "two", "three", "four", "five"], 1000000, n)
    intervals_b = random_intervals(["one", "two", "three", "four", "five"], 1000000, n)

    mktempdir() do dir
        filename_a = joinpath(dir, "test_a.bed")
        filename_b = joinpath(dir, "test_b.bed")
        write_intervals(filename_a, intervals_a)
        write_intervals(filename_b, intervals_b)
        @test check_intersection(filename_a, filename_b)
    end

    @testset "eachoverlap" begin
        path = joinpath(path_of_format("BED"), "ws245Genes.WBGene.bed.bgz")
        stream = BGZFStreams.BGZFStream(path)
        reader = BED.Reader(stream, index=string(path, ".tbi"))
        for (interval, n_records) in [
                (Interval("chrII", 500_000:10_000_000), 38),
                (Interval("chrII", 1:70_000), 0),
                (Interval("chrIV", 500_000:10_000_000), 43),
                (Interval("chrI",  1_000_000:12_000_000), 27),]
            n = 0
            for record in eachoverlap(reader, interval)
                n += 1
            end
            @test n == n_records
        end
        @test isa(BED.Reader(path), BED.Reader)
    end

    # Include doctests.
    DocMeta.setdocmeta!(BED, :DocTestSetup, :(using BED); recursive=true)
    doctest(BED; manual = false)

end


text:
__precompile__()

module RadioTools

#Constants
include("constants.jl")
#Convertion
include("convert.jl")
#Antenna
include("antenna.jl")
#Physics
include("physics.jl")
#Flux
include("flux.jl")

end


text:
using XCB
using Base.Test

scr = Ref{Cint}()
# connect to xserver
connection = XCB.xcb_connect(C_NULL, scr)
println(connection)
#
setup = XCB.xcb_get_setup(connection)
println(setup)
iter = XCB.xcb_setup_roots_iterator(setup)
screen = unsafe_load(iter.data, 1)

window = XCB.xcb_generate_id(connection)
mask = XCB.XCB_CW_BACK_PIXEL | XCB.XCB_CW_EVENT_MASK
value_list = zeros(UInt32, 32)
value_list[1] = screen.white_pixel
value_list[2] = XCB.XCB_EVENT_MASK_EXPOSURE | XCB.XCB_EVENT_MASK_KEY_PRESS

XCB.xcb_create_window(
	connection, UInt8(XCB.XCB_COPY_FROM_PARENT), window,
	screen.root, Int16(0), Int16(0), UInt16(512), UInt16(512), UInt16(1),
	UInt16(XCB.XCB_WINDOW_CLASS_INPUT_OUTPUT), screen.root_visual,
	mask, value_list
)
mask = XCB.XCB_GC_FOREGROUND | XCB.XCB_GC_GRAPHICS_EXPOSURES
value_list[1] = screen.black_pixel
value_list[2] = 0
g = XCB.xcb_generate_id(connection)
XCB.xcb_create_gc(connection, g, window, mask, value_list)

XCB.xcb_map_window(connection, window)
XCB.xcb_flush(connection)
r = Ref(XCB.xcb_rectangle_t(20,20,60,60))
# event loop
done = false
while !done
    e = unsafe_load(XCB.xcb_wait_for_event(connection))
    if e.response_type == XCB.XCB_EXPOSE
        println("EXPOSE")
        XCB.xcb_poly_fill_rectangle(connection, window, g, UInt32(1), r)
        XCB.xcb_flush(connection)
    elseif e.response_type == XCB.XCB_KEY_PRESS
        # exit on keypress */
        done = true
    end
end
XCB.xcb_disconnect(connection)


text:
module OutputSizes

using Flux
using LinearAlgebra

"""
    Nil <: Number

Nil is a singleton type with a single instance `nil`. Unlike
`Nothing` and `Missing` it subtypes `Number`.

"""
struct Nil <: Number end

const nil = Nil()

Nil(::T) where T<:Number = nil

Base.copy(::Nil) = nil

Base.:+(::Nil) = nil
Base.:-(::Nil) = nil

Base.:+(::Nil, ::Nil) = nil
Base.:+(::Nil, ::Number) = nil
Base.:+(::Number, ::Nil) = nil

Base.:-(::Nil, ::Nil) = nil
Base.:-(::Nil, ::Number) = nil
Base.:-(::Number, ::Nil) = nil

Base.:*(::Nil, ::Nil) = nil
Base.:*(::Nil, ::Number) = nil
Base.:*(::Number, ::Nil) = nil

Base.:/(::Nil, ::Nil) = nil
Base.:/(::Nil, ::Number) = nil
Base.:/(::Number, ::Nil) = nil

Base.inv(::Nil) = nil

Base.isless(::Nil, ::Nil) = true
Base.isless(::Nil, ::Number) = true
Base.isless(::Number, ::Nil) = true

Base.abs(::Nil) = nil
Base.exp(::Nil) = nil

Base.typemin(::Type{Nil}) = nil
Base.typemax(::Type{Nil}) = nil
Base.:^(::Nil, ::Nil) = nil

Base.promote(x::Nil, y::Nil) = (nil, nil)
Base.promote(x::Nil, y) = (nil, nil)
Base.promote(x, y::Nil) = (nil, nil)
Base.promote(x::Nil, y, z) = (nil, nil, nil)
Base.promote(x, y::Nil, z) = (nil, nil, nil)
Base.promote(x, y, z::Nil) = (nil, nil, nil)
Base.promote(x::Nil, y, z::Nil) = (nil, nil, nil)
Base.promote(x::Nil, y::Nil, z::Nil) = (nil, nil, nil)
Base.promote(x::Nil, y::Nil, z) = (nil, nil, nil)

LinearAlgebra.adjoint(::Nil) = nil
LinearAlgebra.transpose(::Nil) = nil

# Since some Flux layers don't work out of the box,
# we're "cheating" here by implementing the layers
# for `Nil` with the correct output size.

function (c::Flux.Conv)(x::AbstractArray{Nil})
    return fill(nil, Flux.outdims(c, size(x))..., size(c.weight, 4), size(x, 4))
end

function (c::Flux.DepthwiseConv)(x::AbstractArray{Nil})
    return fill(nil, Flux.outdims(c, size(x))..., size(c.weight, 4), size(x, 4))
end

function (c::Flux.ConvTranspose)(x::AbstractArray{Nil})
    return fill(nil, Flux.outdims(c, size(x))..., size(c.weight, 4), size(x, 4))
end

# Interface

function outputsize(model, sz; addbatch = true)
    return obssize(model(fill(nil, sz..., 1)))
end

# Utils

obssize(a::AbstractArray) = size(a)[1:end-1]
obssize(t::Union{Tuple, NamedTuple}) = map(obssize, t)

export nil, Nil, outputsize

end  # module


text:
using MKtest
using Documenter

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

makedocs(;
    modules=[MKtest],
    authors="Jesús Murga-Moreno, Lawrence Uricchio, David Enard",
    repo="https://github.com/jmurga/MKtest.jl/blob/{commit}{path}#{line}",
    sitename="MKtest.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jmurga.github.io/MKtest.jl",
        assets=String[],
    ),
        pages     = [
                "Home" => "index.md",
                "Package Overview" => [
                        "Analytical estimations" => "analytical.md",
                        "Processing data" => "data.md",
                        "MK approaches" => "mk.md"
                ],
                "Infering the rate and strength of adaptation" =>[
                        "Empirical estimation" => "empirical.md",
                        "Rates" => "rates.md",
                        "Input data" => "input.md",
                        "Summary statistics" => "summstat.md",
                        "ABC inference" => "abc.md",
                        "Multiple dataset" => "multiple.md",
                ],
                "Command-Line interface" => "cli.md",
                "Reference" => "reference.md"
        ]

)

deploydocs(;
    repo="github.com/jmurga/MKtest.jl",
    devbranch="main",
)


text:
#CMP FILE
#Every parameters and fonctions about analysed components databases
#---------------------------------------------------------------------
# CMP Structure : Need 3 files x_DB.txt , xxxx_PATTERN.bmp & xxxx_LB.txt
#	----DB---- From x_DB.txt
#	ToAnalyse? ( -(or other) No, + Yes)
# 	Name
#	heightAnalysisZone (xPx)
#	widthAnalysisZone (yPx)
#	AnalysisFile (.txt)
#	Others Parameters -> Height|Temperature|etc..
#  	----ANALYSIS---- From xxxx_PATTERN.bmp & xxxx_LB.txt
#	IdealParam
#	FactorWeight
#	CategorizingLimit
#	----COORD----
#	VFC (array of type FootprintCoord (X,Y,statParam,difference)) (see PROCESSING.jl)
#		(!) the VFC array is first filled with FootprintCoord that had most lowest sifference score (btw other footprints)
#		then the array is filtered again with gaussian analysis so it only remains the "ready to extract" coordinates.

#		Only the DB section is processed at the start to avoid waiting time.

#	x_DB.txt
#	Name1
#	xPx|yPx
#	AnalysisFile1 (.txt)
#	Others Parameters1 -> Height|Temperature|etc..
#	Name2
#	xPx|yPx
#	Etc...

#	xxxx_LB.txt
# 	IdealParam
#	FactorWeight
#	CategorizingLimit
#	statParam (with spaces as separator)


text:
function psd_projection!(v::Vector{Float64}, a::AuxiliaryData, cones::ConicSets, opt::Options, p::Params, arc_list, iter::Int64)::Nothing

    p.min_eig, current_rank, sqrt_2 = zeros(length(cones.sdpcone)), 0, sqrt(2.)

    # Build symmetric matrix(es) X
    @timeit "reshape1" begin
        cont = 1
        @inbounds for (idx, sdp) in enumerate(cones.sdpcone), j in 1:sdp.sq_side, i in 1:j#j:sdp.sq_side
            if i != j
                a.m[idx].data[i,j] = v[cont] / sqrt_2
            else
                a.m[idx].data[i,j] = v[cont]
            end
            cont += 1
        end
    end

    # Project onto the p.s.d. cone
    for (idx, sdp) in enumerate(cones.sdpcone)
        p.current_rank[idx] = 0

        if sdp.sq_side == 1

            a.m[idx][1] = max(0., a.m[idx][1])
            p.min_eig[idx] = a.m[idx][1]

        elseif !opt.full_eig_decomp && p.target_rank[idx] <= opt.max_target_rank_krylov_eigs && sdp.sq_side > opt.min_size_krylov_eigs
            @timeit "eigs" begin
                eig!(arc_list[idx], a.m[idx], p.target_rank[idx], iter, opt.warm_start_eig)
                if hasconverged(arc_list[idx])
                    fill!(a.m[idx].data, 0.)
                    for i in 1:p.target_rank[idx]
                        if unsafe_getvalues(arc_list[idx])[i] > 0.
                            p.current_rank[idx] += 1
                            vec = unsafe_getvectors(arc_list[idx])[:, i]
                            LinearAlgebra.BLAS.gemm!('N', 'T', unsafe_getvalues(arc_list[idx])[i], vec, vec, 1., a.m[idx].data)
                        end
                    end
                end
            end
            if hasconverged(arc_list[idx])
                @timeit "get min eig" p.min_eig[idx] = minimum(unsafe_getvalues(arc_list[idx]))
            else
                @timeit "eigfact" full_eig!(a, idx, opt, p)
            end

        else
            p.min_eig[idx] = 0.
            @timeit "eigfact" full_eig!(a, idx, opt, p)
        end
    end

    @timeit "reshape2" begin
        cont = 1
        @inbounds for (idx, sdp) in enumerate(cones.sdpcone), j in 1:sdp.sq_side, i in 1:j#j:sdp.sq_side
            if i != j
                v[cont] = a.m[idx].data[i, j] * sqrt_2
            else
                v[cont] = a.m[idx].data[i, j]
            end
            cont += 1
        end
    end

    return nothing
end

function full_eig!(a::AuxiliaryData, idx::Int, opt::Options, p::Params)::Nothing
    p.current_rank[idx] = 0
    fact = eigen!(a.m[idx])
    fill!(a.m[idx].data, 0.)
    for i in 1:length(fact.values)
        if fact.values[i] > 0.
            LinearAlgebra.BLAS.gemm!('N', 'T', fact.values[i], fact.vectors[:, i], fact.vectors[:, i], 1., a.m[idx].data)
            if fact.values[i] > opt.tol_psd
                p.current_rank[idx] += 1
            end
        end
    end

    return nothing
end

function soc_projection!(v::Vector{Float64}, a::AuxiliaryData, cones::ConicSets, opt::Options, p::Params)::Nothing
    for (idx, soc) in enumerate(cones.socone)
        soc_projection!(a.soc_v[idx], a.soc_s[idx])
    end

    return nothing
end

function soc_projection!(v::ViewVector, s::ViewScalar)::Nothing
    nv = norm(v, 2)
    if nv <= -s[]
        s[] = 0.
        v .= 0.
    elseif nv <= s[]
        #do nothing
    else
        val = .5 * (1. + s[] / nv)
        v .*= val
        s[] = val * nv
    end

    return nothing
end

function box_projection!(v::Array{Float64,1}, aff::AffineSets, step::Float64)::Nothing
    # Projection onto = b
    @inbounds @simd for i in 1:length(aff.b)
        v[i] = aff.b[i]
    end
    # Projection onto <= h
    @inbounds @simd for i in 1:length(aff.h)
        v[aff.p+i] = min(v[aff.p+i] / step, aff.h[i])
    end

    return nothing
end


text:
struct Nilpotent{Ti<:Integer}
    size::Ti
    degree::Ti
    offset::Ti
    nilpVec::AbstractVector
end

function factorial(num::Tv, depth::Ti) where {Tv <: Real, Ti <: Integer}
    fact = num
    for i = 1:depth-1
        fact *= (num-i)
    end

    return fact
end

"""
    Nilp(nbOne::Ti, size::Ti) where {Ti<:Integer}

Create a nilpotent matrix object with parameters `nbOne` and `size`. `size` refers to the size of
nilpotent matrix to be generated, and `nbOne` refers to the number of continuous `1` on the non-zero
diagonal of generated nilptent matrix.

This is the default type of nilptent matrix used by SMG2S, in which the non-zero diagonal is selected
as the one of offset `1`. This diagonal starts with `nbOne` number of `1`, then a `0`, then `nbOne` number
of `1`, ... until to the end.

## Examples
```jldoctest; setup = :(using SMG2S)
julia> SMG2S.Nilp(2,8)
Nilpotent{Int64}(2, 8, 3,
  ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ )
```
"""
function Nilp(nbOne::Ti, size::Ti) where {Ti<:Integer}

    degree::Ti = nbOne + 1

    if nbOne > size - 1
        error("for constructing niloptent matrix, nbOne should ≦ size - 1 ")
    end

    nilpVec = zeros(size-1)

    for i = 1:size-1
        if (i % (nbOne + 1)) != 0
            nilpVec[i] = 1.0
        end
    end

    offset = 1
    return Nilpotent(size, degree, offset, nilpVec)
end

"""
    Nilp(vector::AbstractVector, size::Ti) where {Ti<:Integer}

Create a nilpotent matrix whose dimension is `size` from user-provided vector. The
non-zero diagonal of nilpotent matrix is fixed as the one with offset `1`. Therefore, the size
of given `vector` should at least be `size-1`.

## Examples
```jldoctest; setup = :(using SMG2S)
julia> vec=[1; 1; 0; 1; 1; 1; 0]
7-element Vector{Int64}:
 1
 1
 0
 1
 1
 1
 0

julia> SMG2S.Nilp(vec, 8)
Nilpotent{Int64}(3, 8, 4,
  ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ )
```
"""
function Nilp(vector::AbstractVector, size::Ti) where {Ti<:Integer}

    #check
    for i = 1:size-1
        if vector[i] != 1.0 && vector[i] != 0
            vector[i] /= vector[i]
        end
    end

    nbOne::Ti = 0
    cnt::Ti = 0

    for i = 1:size-1
        if vector[i] == 1.0
            cnt += 1
        elseif vector[i] == 0.0
            if cnt > nbOne
                nbOne = cnt
            end
            cnt = 0
        end
    end

    degree::Ti = nbOne + 1

    offset = 1

    return Nilpotent(size, degree, offset, vector[1:size-1])

end

"""
    Nilp(vec::AbstractVector, diag::Ti, size::Ti; maxdegree::Ti=80) where{Ti <: Integer}

Create a nilpotent matrix whose dimension is `size` from user-provided vector. The
non-zero diagonal of nilpotent matrix is set as the one with offset `diag`. Therefore, the size
of given `vector` should at least be `size-diag`. This matrix should be nilpotent, and its nilpotency
should not be larger than `80`. Before generation, SMG2S.jl will check if the user-provided matrix
is nilpotent matrix with nilpotency no larger than `80`. If it doesn't satisfy any of the two,
an error message will appear.

## Examples
```jldoctest; setup = :(using SMG2S; using SparseArrays)
julia> vec=[1; 1; 0; 1; 1; 1; 0]
7-element Vector{Int64}:
 1
 1
 0
 1
 1
 1
 0

julia> SMG2S.Nilp(vec, 3, 8)
┌ Info: the degree of given nilpotent matrix is:
└   degree = 3
Nilpotent{Int64}(1, 8, 3,
 ⋅  ⋅  ⋅  1  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  1  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅)

```
"""
function Nilp(vec::AbstractVector, diag::Ti, size::Ti; maxdegree::Ti=80) where{Ti <: Integer}
    length = size-diag

    nilpvec = vec[1:length]

    vectmps = nilpvec

    diag_offset = diag
    degree = 1
    while degree <= maxdegree
        for i in 1:(size-diag - diag_offset)
            vectmps[i] = nilpvec[i] * vectmps[i+diag]
        end
        vectmps = vectmps[1:(size-diag - diag_offset)]
        degree = degree + 1
        if all(vectmps .== 0)
            break
        end
        diag_offset = diag_offset + diag
    end

    if degree > maxdegree
        error("the given nilpotent matrix is invalid")
    end

    @info "the degree of given nilpotent matrix is: " degree

    return Nilpotent(size, degree, diag, nilpvec)
end

"""
    Nilp(nbOne::Ti, diag::Ti, size::Ti) where{Ti <: Integer}

Create a nilpotent matrix object with parameters `nbOne`, `diag` and `size`. `size` refers to the size of
nilpotent matrix to be generated, and `nbOne` refers to the maximum number of continuous `1` on the non-zero
diagonal of generated nilptent matrix. The offset of non-zero diagonal is determined by
the parameter `diag`.

The non-zero diagonal is generated as:
1. randomly generating an interger number in the interval `[1, nbOne]`

2. this number determines the number of continous `1` for this time of sampling

3. then generating another within the same interval, which determines the number of continuous `0`

4. then re-samling for a new sequence of continuous `1`

5. ...

## Examples
```julia
julia> SMG2S.Nilp(4, 3, 15)
┌ Info: the degree of given nilpotent matrix is:
└   degree = 3
Nilpotent{Int64}(1, 15, 3,
 ⋅  ⋅  ⋅  1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅)
```
"""
function Nilp(nbOne::Ti, diag::Ti, size::Ti) where{Ti <: Integer}
    length = size-diag

    nilpvec = zeros(Ti, length)

    cnt = 0

    while cnt < length
        l₁ = rand(1:nbOne)
        l₀ = rand(1:nbOne)

        for i = cnt+1:l₁+cnt
            if i <= length
                nilpvec[i] = 1
            end
        end

        cnt += l₀ + l₁

    end

    return Nilp(nilpvec, diag, size)
end

"""
    NilpMat(nilp::Nilpotent)
    Explicitly construct the nilpotent matrix
"""
function NilpMat(nilp::Nilpotent)
    return spdiagm(nilp.offset => sparse(nilp.nilpVec))
end

"""
    NilpxM(M::AbstractSparseMatrix, nilp::Nilpotent)
    Nilpotent x M
"""
function NilpxM(M::AbstractSparseMatrix, nilp::Nilpotent)
    rowind, colind, val = findnz(M)
    offset = nilp.offset
    size = nilp.size
    zerosMask = findall(x -> x == 0 , nilp.nilpVec)
    mask =  findall(x -> (x > offset) && !(x-offset in zerosMask), rowind)

    prod = sparse(rowind[mask] .- offset, colind[mask], val[mask], size, size)
    return prod
end

"""
    NilpxM(M::AbstractSparseMatrix, nilp::Nilpotent)
    M x Nilpotent
"""
function MxNilp(M::AbstractSparseMatrix, nilp::Nilpotent)
    rowind, colind, val = findnz(M)
    offset = nilp.offset
    size = nilp.size

    zerosMask = findall(x -> x == 0 , nilp.nilpVec)
    mask =  findall(x -> (x <= size - offset) && !(x in zerosMask) , colind)

    prod = sparse(rowind[mask], colind[mask].+ offset, val[mask], size, size)
    return prod
end

function MNilpM(M::AbstractSparseMatrix, nilp::Nilpotent)
    rowind, colind, val = findnz(M)
    offset = nilp.offset
    size = nilp.size
    zerosMask = findall(x -> x == 0 , nilp.nilpVec)
    rmask =  findall(x -> (x > offset) && !(x-offset in zerosMask), rowind)
    cmask =  findall(x -> (x <= size - offset) && !(x in zerosMask) , colind)

    p1 = sparse(rowind[cmask], colind[cmask].+ offset, val[cmask], size, size)
    p2 = sparse(rowind[rmask] .- offset, colind[rmask], val[rmask], size, size)

    return p1-p2
end


text:
using MIToS.MSA
using MIToS.Information

const msa_long = read("../data/PF00089_aligned.fasta", FASTA)
const msa_wide = read("../data/PF16957_aligned.fasta", FASTA)

percentidentity(msa_long)
percentidentity(msa_wide)

println("[BENCH] Percent Identity PF00089: ", @elapsed percentidentity(msa_long))
println("[BENCH] Percent Identity PF16957: ", @elapsed percentidentity(msa_wide))

max_stars_repo_name: diegozea/mitos-benchmarks

max_issues_repo_name: diegozea/mitos-benchmarks

max_forks_repo_name: diegozea/mitos-benchmarks


text:
mutable struct Rule{W<:AbstractWord}
    lhs::W
    rhs::W
    id::UInt
    active::Bool
end

deactivate!(r::Rule) = r.active = false
isactive(r::Rule) = r.active

function update_rhs!(r::Rule, new_rhs)
    store!(r.rhs, new_rhs)
    r.id = hash(r.lhs, hash(r.rhs))
    return r
end

function Rule{W}(l::AbstractWord, r::AbstractWord, o::Ordering) where W
    lhs, rhs = lt(o, l, r) ? (r, l) : (l, r)
    @assert !lt(o, lhs, rhs) "$lhs should be larger than $rhs"
    return Rule{W}(lhs, rhs, hash(lhs, hash(rhs)), true)
end
Rule(l::W, r::W, o::Ordering) where W = Rule{W}(l, r, o)

function Rule{W}(p::Pair) where W
    lhs, rhs = p
    return Rule{W}(lhs, rhs, hash(lhs, hash(rhs)), true)
end

Rule(p::Pair{W,W}) where W = Rule{W}(p)

function Base.:(==)(rule1::Rule{W}, rule2::Rule{W}) where W
    rule1.id == rule2.id || return false
    res = (rule1.lhs == rule2.lhs) && (rule1.rhs == rule2.rhs)
    res || @warn "hash collision between $rule1 and $rule2"
    return res
end
Base.hash(r::Rule, h::UInt) = hash(r.id, h)

Base.iterate(r::Rule) = r.lhs, 1
Base.iterate(r::Rule, ::Any) = r.rhs, nothing
Base.iterate(r::Rule, ::Nothing) = nothing
Base.length(r::Rule) = 2
Base.last(r::Rule) = first(iterate(r, 1))
Base.eltype(r::Rule{W}) where W = W

Base.show(io::IO, r::Rule) = ((a,b) = r; print(io, a, " ⇒ ", b))

function rules(::Type{W}, o::WordOrdering) where {W<:AbstractWord}
    A = alphabet(o)
    res = Rule{W}[]

    for l in letters(A)
        hasinverse(l, A) || continue
        L = inv(A, l);
        x = W([A[l], A[L]])
        push!(res, Rule(x, one(x), o))
    end
    res
end

mutable struct RulesIter{R, S}
    rules::R
    inner_state::S
end

Base.eltype(::Type{RulesIter{R}}) where R = eltype(R)
IteratorEltype(::Type{RulesIter{R}}) where {R} = IteratorEltype(I)
Base.IteratorSize(::Type{<:RulesIter}) = Base.SizeUnknown()

function Base.iterate(itr::RulesIter, state...)
    y = Base.iterate(itr.rules, itr.inner_state)
    while y !== nothing
        itr.inner_state = last(y)
        isactive(first(y)) && return y
        y = Base.iterate(itr.rules, itr.inner_state)
    end
    return nothing
end

function remove_inactive!(rws, RI, RJ)
    i = 0
    j = 0
    for (idx, r) in enumerate(rws.rwrules)
        if !isactive(r)
            if idx ≤ RI.inner_state
                i += 1
            end
            if idx ≤ RI.inner_state
                j += 1
            end
        end
        idx ≥ max(RI.inner_state, RJ.inner_state) && break
    end
    RI.inner_state -= i
    RJ.inner_state -= j
    remove_inactive!(rws)
    return rws
end

function remove_inactive!(rws, RI)
    i = 0
    for (idx, r) in enumerate(rws.rwrules)
        if !isactive(r)
            if idx ≤ RI.inner_state
                i += 1
            end
        end
        idx ≥ RI.inner_state && break
    end
    RI.inner_state -= i
    remove_inactive!(rws)
    return rws
end

remove_inactive!(rws) = (filter!(isactive, rws.rwrules); rws)


text:
@testset "ReservoirArrayBuffer" begin
    b = ReservoirArrayBuffer{Int}(3, 2)
    @assert size(b) == (3, 0)

    push!(b, [1, 1, 1])
    @assert size(b) == (3, 1)
    @test all(b .== [1; 1; 1])

    push!(b, [2, 2, 2])
    @assert size(b) == (3, 2)
    @test all(b .== [1 2; 1 2; 1 2])

    push!(b, [0, 0, 0])

    @test size(b) == (3, 2)
end


text:
push!(LOAD_PATH,"../../src/")

#=
    bornbem.jl

    Comparison of local and nonlocal exterior potentials of Born ions.
=#

using NESSie
using NESSie.BEM
using NESSie.TestModel
using NESSie.Format
using PyPlot: figure, plot, show, legend, xlabel, ylabel, title, axvline

function plotexterior(name::String, maxdist::Float64=10., resolution::Int=100)
    name          = lowercase(name)
    born          = bornion(name)
    model         = readoff("../../data/born/$name.off")
    model.charges = readpqr("../../data/born/$name.pqr")
    model.params  = Option(1., 78., 1., 23.)
    Ξ             = collect(obspoints_line(
                        [0., 0., born.radius],
                        [0., 0., maxdist],
                        resolution)
                    )
    figure()
    title("$name ion (radius=$(born.radius) Å)")
    axvline(born.radius, color=(.5, .5, .5), linestyle="--", linewidth=3)

    xvals = [ξ[3] for ξ in Ξ]
    print("Local Born... ")
    plot(
        xvals,
        [TestModel.φΣ(LocalES, ξ, born, model.params) for ξ in Ξ],
        label="Born (local)"
    )

    print("done!\nNonlocal Born... ")
    plot(
        xvals,
        [TestModel.φΣ(NonlocalES, ξ, born, model.params) for ξ in Ξ],
        label="Born (nonlocal)"
    )

    print("done!\nLocal BEM... ")
    plot(
        xvals,
        BEM.φΣ(Ξ, solve(LocalES, model)),
        label="BEM (local)"
    )

    print("done!\nNonlocal BEM... ")
    plot(
        xvals,
        BEM.φΣ(Ξ, solve(NonlocalES, model)),
        label="BEM (nonlocal)"
    )
    println("done!")

    legend()
    xlabel("Distance from point charge [Å]")
    ylabel("Exterior potential [V]")
    show()
end

if length(ARGS) < 1
    println("\n\e[1mUsage\e[0m: julia bornbem.jl ION_NAME")
    println("\nwhere ION_NAME is a valid Born ion name (see documentation).")
    exit(1)
end

plotexterior(ARGS[1])


text:
# module Bukdu

immutable Resource{AC<:ApplicationController}
    path::String
    param::String
    controller::Type{AC}
    actions::Vector
    route::Dict
    member::Dict
    collection::Dict
    singleton::Bool
end

module RouterResource

import ..Bukdu: ApplicationController
import ..Bukdu: Keyword, Naming
import ..Bukdu: RouterScope
import ..Bukdu: Resource
import ..Bukdu: index, edit, new, show, create, update, delete

const default_param_key = "id"
const controller_actions = [index, edit, new, show, create, update, delete]

function build{AC<:ApplicationController}(path::String, ::Type{AC}, options::Dict)
    path = RouterScope.validate_path(path)
    alias = Keyword.get(options, :alias, "")
    param = Keyword.get(options, :param, default_param_key)
    name = Keyword.get(options, :name, Naming.resource_name(AC, "Controller"))
    as      = Keyword.get(options, :as, name)
    private = Keyword.get(options, :private, Dict{Symbol,Any}())
    assigns = Keyword.get(options, :assigns, Dict{Symbol,Any}())

    singleton = Keyword.get(options, :singleton, false)
    actions   = extract_actions(options, singleton)

    route       = Dict(:as => as, :private => private, :assigns => assigns)
    collection  = Dict(:path => path, :as => as, :private => private, :assigns => assigns)
    member_path = singleton ? path : string(path, "/:", name, "_", param)
    member      = Dict(:path => member_path, :as => as, :alias => alias, :private => private, :assigns => assigns)

    Resource(path, param, AC, actions, route, member, collection, singleton)
end

function extract_actions(options::Dict, singleton::Bool)::Vector{Function}
    if haskey(options, :only)
        setdiff(controller_actions, setdiff(controller_actions, options[:only]))
    else
        except = Keyword.get(options, :except, [])
        setdiff(default_actions(singleton), except)
    end
end

function default_actions(singleton::Bool)::Vector{Function}
    singleton ? setdiff(controller_actions, [index]) : controller_actions
end

end # module Bukdu.RouterResource


text:
export DSOSPoly, SDSOSPoly, SOSPoly

function PolyJuMP.bridges(::Type{<:PositiveSemidefinite2x2ConeTriangle})
    return [Bridges.Variable.PositiveSemidefinite2x2Bridge]
end
function PolyJuMP.bridges(::Type{<:ScaledDiagonallyDominantConeTriangle})
    return [Bridges.Variable.ScaledDiagonallyDominantBridge]
end
function PolyJuMP.bridges(::Type{<:CopositiveInnerCone})
    return [Bridges.Variable.CopositiveInnerBridge]
end

function JuMP.value(p::GramMatrix{<:JuMP.AbstractJuMPScalar})
    GramMatrix(map(JuMP.value, p.Q), p.x)
end

for poly in (:DSOSPoly, :SDSOSPoly, :SOSPoly)
    @eval begin
        struct $poly{PB<:PolyJuMP.AbstractPolynomialBasis} <: PolyJuMP.AbstractPoly
            polynomial_basis::PB
        end
        $poly(x::AbstractVector{<:MP.APL}) = $poly(MonomialBasis(x))
    end
end

matrix_cone_type(::SOSPoly) = MOI.PositiveSemidefiniteConeTriangle
matrix_cone_type(::DSOSPoly) = DiagonallyDominantConeTriangle
matrix_cone_type(::SDSOSPoly) = ScaledDiagonallyDominantConeTriangle

const PosPoly{PB} = Union{DSOSPoly{PB}, SDSOSPoly{PB}, SOSPoly{PB}}

JuMP.variable_type(m::JuMP.Model, p::PosPoly) = PolyJuMP.polytype(m, p, p.polynomial_basis)
function PolyJuMP.polytype(m::JuMP.Model, cone::PosPoly,
                           basis::PolyJuMP.MonomialBasis{MT, MV}) where {MT<:MP.AbstractMonomial, MV<:AbstractVector{MT}}
    return GramMatrix{JuMP.VariableRef, MT, MV}
end

# Sum-of-Squares polynomial

_polytype(m::JuMP.Model, ::PosPoly, x::MVT) where {MT<:MP.AbstractMonomial, MVT<:AbstractVector{MT}} = GramMatrix{JuMP.VariableRef, MT, MVT}

function moi_add_variable(model::MOI.ModelLike, set::MOI.AbstractVectorSet,
                          binary::Bool, integer::Bool)
    Q, con_Q = MOI.add_constrained_variables(model, set)
    if binary
        for q in Q
            MOI.add_constraint(model, q, MOI.ZeroOne())
        end
    end
    if integer
        for q in Q
            MOI.add_constraint(model, q, MOI.Integer())
        end
    end
    return Q
end

function JuMP.add_variable(model::JuMP.AbstractModel,
                           v::PolyJuMP.Variable{<:PosPoly{<:PolyJuMP.MonomialBasis}},
                           name::String="")
    monos = v.p.polynomial_basis.monomials
    set = matrix_cone(matrix_cone_type(v.p), length(monos))
    # FIXME There is no variable bridge mechanism yet:
    #       https://github.com/JuliaOpt/MathOptInterface.jl/issues/710
    #       so there is not equivalent to `BridgeableConstraint`.
    #       Yet, we may need constraint bridges here if it goes through
    #       the `generic_variable_bridge`.
    #       We don't need the `ScaledDiagonallyDominantBridge` since it does
    #       not use the `generic_variable_bridge`.
    JuMP.add_bridge(model, Bridges.Variable.PositiveSemidefinite2x2Bridge)
    JuMP.add_bridge(model, Bridges.Variable.ScaledDiagonallyDominantBridge)
    JuMP.add_bridge(model, Bridges.Variable.CopositiveInnerBridge)
    JuMP.add_bridge(model, Bridges.Constraint.EmptyBridge)
    JuMP.add_bridge(model, Bridges.Constraint.PositiveSemidefinite2x2Bridge)
    JuMP.add_bridge(model, Bridges.Constraint.DiagonallyDominantBridge)
    Q = moi_add_variable(backend(model), set, v.binary, v.integer)
    return build_gram_matrix(
        JuMP.VariableRef[JuMP.VariableRef(model, vi) for vi in Q], monos)
end


text:
const Detection2D = Tuple{Float64, Float64, Int64}

"""
    get_line_segment(camera_params::Camera_Params, params::Video_Params, detection::Detection2D)

Returns a line segment of the possible 3D positions of the object that could have caused a 2D detection
"""
function get_line_segment(camera_params::Camera_Params, params::Video_Params, detection::Detection2D)

    x = detection[1]
    y = detection[2]

    #now (0,0) is the center of the image
    x = x - params.image_dim_x/2
    y = y - params.image_dim_y/2

    angle_from_vertical = x / (params.image_dim_x/2) * deg2rad(params.horizontal_FoV/2)
    angle_from_horizontal = -1 * y / (params.image_dim_y/2) * deg2rad(params.vertical_FoV/2)

    #println("angle_from_vertical ", rad2deg(angle_from_vertical))
    #println("angle_from_horizontal ", rad2deg(angle_from_horizontal))

    #something wrong

    (a_vertical, b_vertical, c_vertical) = get_vertical_plane(camera_params)
    # println((a_vertical, b_vertical, c_vertical))
    (a_horizontal, b_horizontal, c_horizontal) = get_horizontal_plane(camera_params, a_vertical, b_vertical, c_vertical)
    # println((a_horizontal, b_horizontal, c_horizontal))
    # normal_to_image = (camera_params.camera_focus.x - camera_params.camera_location.x,
    # camera_params.camera_focus.y - camera_params.camera_location.y,
    # camera_params.camera_focus.z - camera_params.camera_location.z)
    # #get one more point on the ray. don't care which point, but could be on same plane as camera's focus
    # l = get_distance(camera_params.camera_location, camera_params.camera_focus)
    # println("l ", l)
    # w = 2*l*tan(deg2rad(params.horizontal_FoV/2)) #total image width in meters (non-pixel space)
    # #x = angle_from_vertical / (deg2rad(params.horizontal_FoV/2)) * w / 2
    # x = x * (w / 2) / (params.image_dim_x/2)
    # h = 2*l*tan(deg2rad(params.vertical_FoV/2)) #total image heigth in meters (non-pixel space)
    # #y = angle_from_horizontal / (deg2rad(params.vertical_FoV/2)) * h / 2
    # y = -1 * y * (h / 2) / (params.image_dim_y/2)

    l = get_distance(camera_params.camera_location, camera_params.camera_focus)
    x = l*tan(angle_from_vertical)
    y = l*tan(angle_from_horizontal)

    #println("x ", x)
    #println("y ", y)

    normalized_v = [a_vertical, b_vertical, c_vertical]./norm([a_vertical, b_vertical, c_vertical])
    normalized_h = [a_horizontal, b_horizontal, c_horizontal]./norm([a_horizontal, b_horizontal, c_horizontal])
    #println("normalized_v ", normalized_v)
    #println("normalized_h ", normalized_h)
    focus = [camera_params.camera_focus.x, camera_params.camera_focus.y, camera_params.camera_focus.z]
    #println("focus ", focus)
    #println("dot that should be 0 ", dot(x .* normalized_v + y .* normalized_h, normal_to_image))
    point = focus + x .* normalized_v + y .* normalized_h
    point = Coordinate(point[1], point[2], point[3])

    #println("point ", point)

    #check each wall for an intersection within limits
    (endpoint_x, endpoint_y, endpoint_z, a, b, c) = check_walls(point, camera_params.camera_location, params)

    #return Line_Segment(camera_params.camera_location, Coordinate(endpoint_x, endpoint_y, endpoint_z), a, b, c)
    #Flip line segment to LHR coordinate system
    l = Line_Segment(Coordinate(camera_params.camera_location.x, camera_params.camera_location.y, -camera_params.camera_location.z),
        Coordinate(endpoint_x, endpoint_y, -endpoint_z), a, b, -c)
    return l
end

function get_distance(p1::Coordinate, p2::Coordinate)
    return sqrt((p1.x-p2.x)^2 + (p1.y-p2.y)^2 + (p1.z-p2.z)^2)
end

"""
    check_walls(point::Coordinate, camera::Coordinate, params::Video_Params)

Calculates parameters for the line along which the object could appear, alongside
the coordinates corresponding to the farthest the object could be from the
camera while still being within the boundary of the scene.

Arguments:
- point: possible 3d location of the object
- camera: location of the camera
- params:: Video_Params

Returns:
- (x, y, z): farthest the location from the camera the object could be while
still being within the boundary of the scene.
- (a, b, c): parametrization of the line along which the object could appear.
when t = 1, it gives the farthest point (x, y, z)
"""
function check_walls(point::Coordinate, camera::Coordinate, params::Video_Params)
    a = point.x - camera.x
    b = point.y - camera.y
    c = point.z - camera.z

    #get intersection between line and plane

    #check x = x_min
    x = params.x_min
    t = (x - camera.x)/a
    y = camera.y + b*t
    z = camera.z + c*t
    if y < params.y_max && y > params.y_min && z < params.z_max && z > params.z_min && t > 0 #t>0 checks right side
        #redoing a,b,c so when t=1, at this fartest point
        a = x - camera.x
        b = y - camera.y
        c = z - camera.z
        return (x, y, z, a, b, c)
    end

    #check x = x_max
    x = params.x_max
    t = (x - camera.x)/a
    y = camera.y + b*t
    z = camera.z + c*t
    if y < params.y_max && y > params.y_min && z < params.z_max && z > params.z_min && t > 0
        a = x - camera.x
        b = y - camera.y
        c = z - camera.z
        return (x, y, z, a, b, c)
    end

    #check y = y_min
    y = params.y_min
    t = (y - camera.y)/b
    x = camera.x + a*t
    z = camera.z + c*t
    if x < params.x_max && x > params.x_min && z < params.z_max && z > params.z_min && t > 0
        a = x - camera.x
        b = y - camera.y
        c = z - camera.z
        return (x, y, z, a, b, c)
    end

    #check y = y_max
    y = params.y_max
    t = (y - camera.y)/b
    x = camera.x + a*t
    z = camera.z + c*t
    if x < params.x_max && x > params.x_min && z < params.z_max && z > params.z_min && t > 0
        a = x - camera.x
        b = y - camera.y
        c = z - camera.z
        return (x, y, z, a, b, c)
    end

    #check z = z_min
    z = params.z_min
    t = (z - camera.z)/c
    x = camera.x + a*t
    y = camera.y + b*t
    if x < params.x_max && x > params.x_min && y < params.y_max && y > params.y_min && t > 0
        a = x - camera.x
        b = y - camera.y
        c = z - camera.z
        return (x, y, z, a, b, c)
    end

    #check z = z_max
    z = params.z_max
    t = (z - camera.z)/c
    x = camera.x + a*t
    y = camera.y + b*t
    if x < params.x_max && x > params.x_min && y < params.y_max && y > params.y_min && t > 0
        a = x - camera.x
        b = y - camera.y
        c = z - camera.z
        return (x, y, z, a, b, c)
    end
end

export get_line_segment

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
struct PolicyTestResult
    test::PolicyTest
    prediction::Vector{Float32}
end

struct PolicyTestResultSVG
    result::PolicyTestResult
end

svg(x::PolicyTestResult) = PolicyTestResultSVG(x)
policypredsvg(state::AbstractState, ::PolicyTestResult) = @assert false OVERRIDE_MSG

function Base.show(io::IO, m::MIME"image/svg+xml", x::PolicyTestResultSVG)
    test = x.result.test
    pred = x.result.prediction
    s = test.state
    player_name = s.to_play.name
    actions = test.actions
    state_grid = getgrid(s; header="to play: $player_name", margin_left=0)
    next_states = Rect[]
    from = 1
    len = length(actions)
    to = min(len, 3)
    for (idx,a) in enumerate(actions[from:to])
        push!(
            next_states,
            getgrid(s, a, 0f0, getnextstate(s, a); header="action: $a", margin_left=0)
        )
    end
    if length(actions) > 3
        push!(next_states, Grid(1,1, (r,c)->"..."; line_thickness=0))
    end
    input_grid = Grid(
        [state_grid],
        :,
        1,
        header="INPUT",
        margin_left=0,
        line_thickness=0
    )
    policy_grid_header = "ACTUAL π*"
    policy_grid = Grid(
        next_states,
        :,
        1,
        header=(text=policy_grid_header, height=40),
        margin_left=0,
        line_thickness=0,
    )
    pred_states = Rect[]
    limit = 3
    from = 1
    to = min(limit, length(pred))
    pred_actions = reverse(sortperm(pred))[from:to]
    prs = pred[pred_actions]
    for (a,pr) in zip(pred_actions, prs)
        pr_str = @sprintf "%.2f" pr
        push!(
            pred_states,
            getgrid(
                s,
                a,
                0f0,
                getnextstate(s, a);
                margin_left=0,
                header="action: $a",
                subheader=(text=pr_str, color="#666"),
            )
        )
    end
    pred_grid = Grid(
        pred_states,
        :,
        1;
        header="PRED π*",
        margin_left=0,
        line_thickness=0,
    )
    g = Grid([input_grid,pred_grid,policy_grid]; cell_valign=:top, header="POLICY TEST RESULT #$(test.idx)", line_thickness=0)
    show(io, m, g)
end


text:
BOOLEAN_EQUALS(s::State) = if length(s.boolean) >= 2
  push!(s.boolean, pop!(s.boolean) == pop!(s.boolean))
end

# Don't bother actually popping and pushing. Edit the top of the stack.
BOOLEAN_NOT(s::State) = if !isempty(s.boolean)
  s.boolean[1] = !s.boolean[1]
end

BOOLEAN_AND(s::State) = if length(s.boolean) >= 2
  push!(s.boolean, pop!(s.boolean) & pop!(s.boolean))
end

BOOLEAN_OR(s::State) = if length(s.boolean) >= 2
  push!(s.boolean, pop!(s.boolean) | pop!(s.boolean))
end

BOOLEAN_DEFINE(s::State) = if !isempty(s.boolean) && !isempty(s.name)
  s.definitions[pop!(s.name)] = pop!(s.boolean)
end

BOOLEAN_DUP(s::State) = if !isempty(s.boolean)
  push!(s.boolean, peek(s.boolean))
end

BOOLEAN_FLUSH(s::State) = if !isempty(s.boolean)
  empty!(s.boolean)
end

BOOLEAN_FROM_FLOAT(s::State) = if !isempty(s.float)
  push!(s.boolean, pop!(s.float) != zero(Float32))
end

BOOLEAN_FROM_INT(s::State) = if !isempty(s.integer)
  push!(s.boolean, pop!(s.integer) != zero(Int32))
end

BOOLEAN_POP(s::State) = pop!(s.boolean)

BOOLEAN_ROT(s::State) = if length(s.boolean) > 2
  s.boolean[end], s.boolean[end-1], s.boolean[end-2] =
    s.boolean[end-1], s.boolean[end-2], s.boolean[end]
end

BOOLEAN_SHOVE(s::State) = if !isempty(s.integer) && !isempty(s.boolean)
  shove!(s.boolean, pop!(s.integer))
end

BOOLEAN_RAND(s::State) = push!(s.boolean, random_bool(s))

BOOLEAN_STACK_DEPTH(s::State) = push!(s.integer, length(s.boolean))

BOOLEAN_SWAP(s::State) = if length(s.boolean) >= 2
  s.boolean[end], s.boolean[end-1] = s.boolean[end-1], s.boolean[end]
end

BOOLEAN_YANK(s::State) = if !isempty(s.integer) && !isempty(s.boolean)
  yank!(s.boolean, pop!(s.integer))
end

BOOLEAN_YANKDUP(s::State) = if !isempty(s.integer) && !isempty(s.boolean)
  yankdup!(s.boolean, pop!(s.integer))
end

BOOLEAN_ERC(s::State) = return

# could tell it which stacks we use?
Push.register("BOOLEAN.=",           BOOLEAN_EQUALS)
Push.register("BOOLEAN.NOT",         BOOLEAN_NOT)
Push.register("BOOLEAN.AND",         BOOLEAN_AND)
Push.register("BOOLEAN.OR",          BOOLEAN_OR)
Push.register("BOOLEAN.DEFINE",      BOOLEAN_DEFINE)
Push.register("BOOLEAN.DUP",         BOOLEAN_DUP)
Push.register("BOOLEAN.FLUSH",       BOOLEAN_FLUSH)
Push.register("BOOLEAN.FROMFLOAT",   BOOLEAN_FROM_FLOAT)
Push.register("BOOLEAN.FROMINT",     BOOLEAN_FROM_INT)
Push.register("BOOLEAN.POP",         BOOLEAN_POP)
Push.register("BOOLEAN.ROT",         BOOLEAN_ROT)
Push.register("BOOLEAN.SHOVE",       BOOLEAN_SHOVE)
Push.register("BOOLEAN.RAND",        BOOLEAN_RAND)
Push.register("BOOLEAN.STACKDEPTH",  BOOLEAN_STACK_DEPTH)
Push.register("BOOLEAN.SWAP",        BOOLEAN_SWAP)
Push.register("BOOLEAN.YANK",        BOOLEAN_YANK)
Push.register("BOOLEAN.YANKDUP",     BOOLEAN_YANKDUP)
Push.register("BOOLEAN.ERC",         BOOLEAN_ERC)


text:
module GridViewer

include("main.jl")

export D3Grid, run, update
export AbstractGrid

end # module


text:
module TestingBinaryHeaps

using Test

using Base: Ordering, ForwardOrdering, ReverseOrdering, Forward, Reverse, lt

using QuickHeaps
using QuickHeaps:
    AbstractBinaryHeap,
    FastForwardOrdering, FastForward, FastReverse,
    FastMin, FastMax, SafeMin, SafeMax,
    isheap, heapify, heapify!,
    ordering

is_max_ordering(x) = !is_min_ordering(x)
is_min_ordering(o::Ordering) = is_min_ordering(typeof(o))
is_min_ordering(::AbstractBinaryHeap{<:Any,O}) where {O} = is_min_ordering(O)
is_min_ordering(::Type{<:FastForwardOrdering}) = true
is_min_ordering(::Type{<:ForwardOrdering}) = true
is_min_ordering(::Type{<:ReverseOrdering{O}}) where {O} = is_max_ordering(O)

function is_sorted(o::Base.Ordering, x::AbstractVector)
    flag = false
    for i in 2:length(x)
        flag |= lt(o, x[i], x[i-1])
    end
    return !flag
end

function heap_test(T::Type, h::AbstractBinaryHeap, n::Integer = 7)
    # Build array of values.  For n ≥ 3, it should not be a heap in any order.
    n ≥ 3 || error("number of elements must be at least 3")
    A = T[(2:2:n)..., (1:2:n)...]
    @test !isheap(A, Forward)
    @test !isheap(A, Reverse)

    @test eltype(h) === T
    @test size(h) == (length(h),)
    @test IndexStyle(h) == IndexLinear()
    empty!(h)
    @test isempty(h)
    @test isheap(h; check=true)
    @test_throws ArgumentError peek(h)
    @test_throws ArgumentError pop!(h)
    for i in 1:n
        push!(h, rand(T))
        @test isheap(h; check=true)
        @test !isempty(h)
        @test length(h) == i
        v1 = (is_min_ordering(h) ? minimum(h) : maximum(h))
        @test peek(h) == v1
        @test first(h) == v1
    end
    for i in 1:n
        h[i] = rand(T)
        @test isheap(h; check=true)
        @test length(h) == n
        v1 = (is_min_ordering(h) ? minimum(h) : maximum(h))
        @test peek(h) == v1
        @test first(h) == v1
    end
    x = T[]
    while !isempty(h)
        push!(x, pop!(h))
    end
    @test length(x) == n
    @test is_sorted(ordering(h), x)
end

@testset "Arrays as binary heaps" begin
    # Min-heap.
    let A = [1, 4, 2, 6, 5, 3, 8, 7, 11, 13, 10, 14, 9, 12]
        n = length(A)
        @test isheap(A)
        @test isheap(A, n)
        @test isheap(A, n - 1) # last element can be discarded
        for o in (Forward, FastForward)
            @test isheap(A, o)
            @test isheap(o, A)
            @test isheap(A, n, o)
            @test isheap(A, o, Int16(n))
            @test isheap(o, A, n)
            @test isheap(o, A, n - 1) # last element can be discarded
            @test !isheap(A, ReverseOrdering(o))
        end
    end
    # Max-heap.
    let A = [14, 11, 13, 9, 10, 12, 7, 8, 4, 2, 5, 1, 6, 3]
        n = length(A)
        @test !isheap(A)
        for o in (Reverse, FastReverse)
            @test isheap(A, o)
            @test isheap(o, A)
            @test isheap(A, n, o)
            @test isheap(A, o, n)
            @test isheap(o, A, n)
            @test isheap(o, A, n - 1) # last element can be discarded
            @test !isheap(A, ReverseOrdering(o))
        end
    end
    # Starting with an array which is not a heap in any order.
    let A = [13, 1, 14, 8, 10, 9, 7, 2, 5, 4, 12, 3, 6, 11]
        n = length(A)
        @test !isheap(A)
        @test !isheap(A, n)
        @test !isheap(A, Reverse)
        @test !isheap(A, FastForward)
        @test !isheap(A, FastReverse)
        B = heapify(A)
        @test isheap(B)
        @test !isheap(B, Reverse)
        B = heapify(A, Reverse, Int16(n))
        @test !isheap(B)
        @test isheap(B, Reverse)
        C = copy(A)
        B = heapify!(C)
        @test B === C
        @test isheap(B)
        B = heapify!(C, Reverse, Int16(n))
        @test B === C
        @test isheap(B, Reverse)
    end
end

@testset "Binary heaps          " begin
    n = 8
    for T in (Float64, Int)
        heap_test(T, BinaryHeap{T}(), n)
        heap_test(T, BinaryHeap{T}(FastMin), n)
        heap_test(T, BinaryHeap{T}(FastMax), n)
        heap_test(T, BinaryHeap{T}(SafeMin), n)
        heap_test(T, BinaryHeap{T}(SafeMax), n)
    end
end

@testset "Fast binary heaps     " begin
    n = 7
    for T in (Float32, Int16)
        heap_test(T, FastBinaryHeap{T}(), n)
        heap_test(T, FastBinaryHeap{T}(FastMin), n)
        heap_test(T, FastBinaryHeap{T}(FastMax), n)
        heap_test(T, FastBinaryHeap{T}(SafeMin), n)
        heap_test(T, FastBinaryHeap{T}(SafeMax), n)
    end
end

end # module


text:
using StructuresKit

shape_name = "W14X90"
#define discretization along member
zcoords = 0:1:20
n=(4, 4, 4, 4, 4)

shape_info = StructuresKit.CrossSection.AISC(shape_name)

#A Ix Iy J Cw xc yc xs ys
section_properties = [(shape_info.A, shape_info.Ix, shape_info.Iy, shape_info.J, shape_info.Cw, shape_info.d/2, 0.0, shape_info.d/2, 0.0)]

#E  ν
material_properties = [(29000,0.30)]

#kx ky kϕ hx hy
springs = [0.0, 0.0, 0.0, 0.0, 0.0]

#member information
#L dL SectionProperties MaterialProperties Springs
member_definitions = [(13.0*12, (13*12)/12,1,1,1)]

num_nodes = Int(member_definitions[1][1]/member_definitions[1][2]) + 1

#P
loads = 4000.0 * ones(num_nodes)

#end boundary conditions
#type=1 u''=v''=ϕ''=0 (simply supported), type=2 u'=v'=ϕ'=0  (fixed), type=3 u''=v''=ϕ''=u'''=v'''=ϕ'''=0 (free end, e.g., a cantilever)
end_boundary_conditions = [1 1]

#supports
#location where u=v=ϕ=0
supports = [0.0 13.0*12]

#imperfection
L = member_definitions[1][1]
dL = member_definitions[1][2]
z = 0: dL :L
Δo = L / 1000 * sin.(π*z/L)

uo = Δo  #in the weak axis
vo = zeros(Float64, num_nodes)
ϕo = zeros(Float64, num_nodes)

imperfections = [uo, vo, ϕo]

P = 0:100:4000

u_total_predicted_max = zeros(Float64, length(P))
u_total_theor_max = zeros(Float64, length(P))

#Calculate weak axis flexural buckling load
E = material_properties[1][1]
Iy = shape_info.Iy
Pcry = π^2*E*Iy/L^2   #weak axis buckling loads

for i=1:length(P)

    loads = P[i] * ones(num_nodes)
    u, v, ϕ, properties = Column.solve(member_definitions, section_properties, material_properties, loads, springs, end_boundary_conditions, supports, imperfections)

    u_total_predicted_max[i] = maximum(u) + L/1000

    #displacement at midheight
    #Chajes Eq. 1.70
    u_total_theor_max[i] = L/1000*(1/(1-P[i]/Pcry))

end

# using Plots
# plot(u_total_predicted_max, P, markershape = :o)
# plot!(u_total_theor_max, P, markershape = :hexagon)

#calculate percent error at every load compared to Mathematica solution
u_difference= abs.((u_total_predicted_max .- u_total_theor_max))./ u_total_theor_max

MaxError=maximum(u_difference[2:end])
#accept 7% error from numerical solution
#solution becomes less consistent with Chajes as P approaches Pcry
#this could be because the Chajes solution is an approximation, hard to tell if computed solution or Chajes is more correct

@test MaxError <= 0.07


text:
function initQuadTree(img, var, branch)
  qTree = QuadTree(Rectangle(Int(size(img)[1]/2),
    Int(size(img)[1]/2),
    Int(size(img)[1]/2),
    Int(size(img)[2]/2)), var, branch)
  return qTree
end

max_stars_repo_name: neuro-myoung/CellSegmenter

max_issues_repo_name: neuro-myoung/CellSegmenter

max_forks_repo_name: neuro-myoung/CellSegmenter


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1ServicePort <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    nodePort::Any # spec type: Union{ Nothing, Int32 } # spec name: nodePort
    port::Any # spec type: Union{ Nothing, Int32 } # spec name: port
    protocol::Any # spec type: Union{ Nothing, String } # spec name: protocol
    targetPort::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgUtilIntstrIntOrString } # spec name: targetPort

    function IoK8sApiCoreV1ServicePort(;name=nothing, nodePort=nothing, port=nothing, protocol=nothing, targetPort=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ServicePort, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiCoreV1ServicePort, Symbol("nodePort"), nodePort)
        setfield!(o, Symbol("nodePort"), nodePort)
        validate_property(IoK8sApiCoreV1ServicePort, Symbol("port"), port)
        setfield!(o, Symbol("port"), port)
        validate_property(IoK8sApiCoreV1ServicePort, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        validate_property(IoK8sApiCoreV1ServicePort, Symbol("targetPort"), targetPort)
        setfield!(o, Symbol("targetPort"), targetPort)
        o
    end
end # type IoK8sApiCoreV1ServicePort

const _property_map_IoK8sApiCoreV1ServicePort = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("nodePort")=>Symbol("nodePort"), Symbol("port")=>Symbol("port"), Symbol("protocol")=>Symbol("protocol"), Symbol("targetPort")=>Symbol("targetPort"))
const _property_types_IoK8sApiCoreV1ServicePort = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("nodePort")=>"Int32", Symbol("port")=>"Int32", Symbol("protocol")=>"String", Symbol("targetPort")=>"IoK8sApimachineryPkgUtilIntstrIntOrString")
Base.propertynames(::Type{ IoK8sApiCoreV1ServicePort }) = collect(keys(_property_map_IoK8sApiCoreV1ServicePort))
Swagger.property_type(::Type{ IoK8sApiCoreV1ServicePort }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1ServicePort[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ServicePort }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ServicePort[property_name]

function check_required(o::IoK8sApiCoreV1ServicePort)
    (getproperty(o, Symbol("port")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ServicePort }, name::Symbol, val)
end


text:
@testset "types.jl: Confinement" begin
    @test isconfined(BulkConfinement()) == false
    @test isconfined(BoxConfinement()) == true
    @test isconfined(SlabConfinement()) == true
    @test isconfined(DiskConfinement()) == true
    @test isconfined(SphereConfinement()) == true
    @test isconfined(CylinderConfinement()) == true
end

@testset "types.jl: Charge" begin
    @test ischarged(Neutral()) == false
    @test ischarged(SmearedCharge()) == true
    @test ischarged(DiscreteCharge()) == true
end

@testset "types.jl: Specie" begin
    segment = KuhnSegment(:A)
    @test segment.b == 1.0
    @test segment.M == 1.0
end

@testset "types.jl: Block" begin
    sA = KuhnSegment(:A)
    e1A = FreeEnd(:A1)
    e2A = FreeEnd(:A2)
    A = PolymerBlock(:A, sA, 1.0, e1A, e2A)
    @test A.f == 1.0
end

@testset "types.jl: Component" begin
    sA = KuhnSegment(:A)
    sB = KuhnSegment(:B)
    e1A = FreeEnd(:A1)
    e2A = BranchPoint(:A2)
    e1B = FreeEnd(:B1)
    e2B = BranchPoint(:B2)
    A = PolymerBlock(:A, sA, 0.4, e1A, e2A)
    B = PolymerBlock(:B, sB, 0.6, e1B, e2B)
    c = BlockCopolymer(:AB, [A,B])
    pc = Component(c)
    @test pc.α == 1.0
    @test pc.ϕ == 1.0

    smol = SmallMolecule(:S)
    @test smol.b == 1.0
    @test smol.M == 1.0
    sc = Component(smol)
    @test sc.α == 1.0
    @test sc.ϕ == 1.0
end

@testset "types.jl: Polymer System AB+S" begin
    sA = KuhnSegment(:A)
    sB = KuhnSegment(:B)
    e1A = FreeEnd(:A1)
    e2A = BranchPoint(:A2)
    e1B = FreeEnd(:B1)
    e2B = BranchPoint(:B2)
    A = PolymerBlock(:A, sA, 0.4, e1A, e2A)
    B = PolymerBlock(:B, sB, 0.6, e1B, e2B)
    chain = BlockCopolymer(:AB, [A,B])
    smol = SmallMolecule(:S)
    polymer = Component(chain, 1.0, 0.1)
    solvent = Component(smol, 0.01, 0.9)
    χN_map = Dict(
        Set([:A,:B])=>20.0,
        Set([:A,:S])=>80.0,
        Set([:B,:S])=>120.0)
    m = χNMatrix(χN_map)
    ABS = PolymerSystem([polymer, solvent], m)
    @test isconfined(ABS.confinement) == false
    @test ABS.C == 1.0

    @test multicomponent(ABS) == true
    @test ncomponents(ABS) == 2
    @test species(chain) == [:A, :B]
    @test nspecies(chain) == 2
    @test species(smol) == [:S]
    @test nspecies(smol) == 1
    @test species(polymer) == species(chain)
    @test nspecies(polymer) == nspecies(chain)
    @test species(solvent) == species(smol)
    @test nspecies(solvent) == nspecies(smol)
    @test species(ABS) == [:A, :B, :S]
    @test nspecies(ABS) == 3

    @test systemtype(ABS) == PolymerSolution()
end

@testset "types.jl: Polymer System AB+A" begin
    sA = KuhnSegment(:A)
    sB = KuhnSegment(:B)
    eA = FreeEnd(:A1)
    eAB = BranchPoint(:AB)
    eB = FreeEnd(:B1)
    A = PolymerBlock(:A, sA, 0.4, eA, eAB)
    B = PolymerBlock(:B, sB, 0.6, eB, eAB)
    chainAB = BlockCopolymer(:AB, [A,B])
    hA = PolymerBlock(:hA, sA, 1.0, FreeEnd(:hA1), FreeEnd(:hA2))
    chainA = BlockCopolymer(:hA, [hA])
    polymerAB = Component(chainAB, 1.0, 0.5)
    polymerA = Component(chainA, 0.5, 0.5)
    χN_map = Dict(Set([:A,:B])=>20.0)
    AB_A = PolymerSystem([polymerAB, polymerA], χN_map)

    @test multicomponent(AB_A) == true
    @test ncomponents(AB_A) == 2
    @test species(chainAB) == [:A, :B]
    @test nspecies(chainAB) == 2
    @test species(chainA) == [:A]
    @test nspecies(chainA) == 1
    @test species(polymerAB) == species(chainAB)
    @test nspecies(polymerAB) == nspecies(chainAB)
    @test species(polymerA) == species(chainA)
    @test nspecies(polymerA) == nspecies(chainA)
    @test species(AB_A) == [:A, :B]
    @test nspecies(AB_A) == 2

    @test systemtype(AB_A) == PolymerBlend()
end


text:
type Timedata{T} <: AbstractTimedata
    vals::DataFrame
    idx::Array{T, 1}

    function Timedata(vals::DataFrame, idx::Array{T, 1})
        chkIdx(idx)
        if(size(vals, 1) != length(idx))
            if (length(idx) == 0) | (size(vals, 1) == 0)
                df = convert(DataFrame, rand(2, 2))
                return new(df[1:2, [false, false]], Array{T, 1}[])
            end
            error(length(idx), " idx entries, but ", size(vals, 1), " rows of data")
        end
        return new(vals, idx)
    end
end

function Timedata{T}(vals::DataFrame, idx::Array{T, 1})
    return Timedata{T}(vals, idx)
end

#######################
## find boolean true ##
#######################

import Base.find
function find(td::TimeData.Timedata)
    ## apply find to column DataArrays
    (allRowInds, allColInds) = find2sub(td)
    return sub2ind(size(td), allRowInds, allColInds)
end

function find2sub(td::TimeData.Timedata)
    ## apply find to column DataArrays and return tuple of subscripts
    allRowInds = Int[]
    allColInds = Int[]
    for (name, col) in eachcol(td.vals)
        rowInds = find(col)
        colInd = td.vals.colindex[name]
        colInds = repmat([colInd], length(rowInds))
        append!(allRowInds, [rowInds])
        append!(allColInds, [colInds])
    end
    return (allRowInds, allColInds)
end

function find(f::Function, td::TimeData.AbstractTimedata)
    ## apply function elementwise to td
    ## apply find to column DataArrays

    (allRowInds, allColInds) = find2sub(f, td)
    return sub2ind(size(td), allRowInds, allColInds)
end

function find2sub(f::Function, td::TimeData.AbstractTimedata)
    ## apply function elementwise to td
    ## apply find to column DataArrays

    (nObs, nVars) = size(td)

    allRowInds = Int[]
    allColInds = Int[]

    for ii=1:nObs
        for jj=1:nVars
            fval = f(get(td, ii, jj))
            if !isa(fval, Bool) & !isna(fval)
                error("function must return boolean value or NA")
            end
            if isequal(fval, true)
                push!(allRowInds, ii)
                push!(allColInds, jj)
            end
        end
    end
    return (allRowInds, allColInds)
end


text:
@testset "Block" begin
    @testset "Parse" begin
        block_raw = hex2bytes("020000208ec39428b17323fa0ddec8e887b4a7c53b8c0a0a220cfd0000000000000000005b0750fce0a889502d40508d39576821155e9c9e3f5c3157f961db38fd8b25be1e77a759e93c0118a4ffd71d")
        stream = IOBuffer(block_raw)
        block = Bitcoin.io2blockheader(stream)
        @test block.version == 0x20000002
        want = hex2bytes("8ec39428b17323fa0ddec8e887b4a7c53b8c0a0a220cfd000000000000000000")
        @test block.prev_block == want
        want = hex2bytes("5b0750fce0a889502d40508d39576821155e9c9e3f5c3157f961db38fd8b25be")
        @test block.merkle_root == want
        @test block.timestamp == 0x59a7771e
        @test block.bits == 0x18013ce9
        @test block.nonce == 0x1dd7ffa4
    end
    @testset "Serialize" begin
        block_raw = hex2bytes("020000208ec39428b17323fa0ddec8e887b4a7c53b8c0a0a220cfd0000000000000000005b0750fce0a889502d40508d39576821155e9c9e3f5c3157f961db38fd8b25be1e77a759e93c0118a4ffd71d")
        stream = IOBuffer(block_raw)
        block = Bitcoin.io2blockheader(stream)
        @test Bitcoin.serialize(block) == block_raw
    end
    @testset "Hash" begin
        block_raw = hex2bytes("020000208ec39428b17323fa0ddec8e887b4a7c53b8c0a0a220cfd0000000000000000005b0750fce0a889502d40508d39576821155e9c9e3f5c3157f961db38fd8b25be1e77a759e93c0118a4ffd71d")
        stream = IOBuffer(block_raw)
        block = Bitcoin.io2blockheader(stream)
        Bitcoin.hash(block) == hex2bytes("0000000000000000007e9e4c586439b0cdbe13b1370bdd9435d76a644d047523")
    end
    @testset "BIP9" begin
        block_raw = hex2bytes("020000208ec39428b17323fa0ddec8e887b4a7c53b8c0a0a220cfd0000000000000000005b0750fce0a889502d40508d39576821155e9c9e3f5c3157f961db38fd8b25be1e77a759e93c0118a4ffd71d")
        stream = IOBuffer(block_raw)
        block = Bitcoin.io2blockheader(stream)
        @test Bitcoin.bip9(block)
        block_raw = hex2bytes("0400000039fa821848781f027a2e6dfabbf6bda920d9ae61b63400030000000000000000ecae536a304042e3154be0e3e9a8220e5568c3433a9ab49ac4cbb74f8df8e8b0cc2acf569fb9061806652c27")
        stream = IOBuffer(block_raw)
        block = Bitcoin.io2blockheader(stream)
        @test !Bitcoin.bip9(block)
    end
    @testset "BIP91" begin
        block_raw = hex2bytes("1200002028856ec5bca29cf76980d368b0a163a0bb81fc192951270100000000000000003288f32a2831833c31a25401c52093eb545d28157e200a64b21b3ae8f21c507401877b5935470118144dbfd1")
        stream = IOBuffer(block_raw)
        block = Bitcoin.io2blockheader(stream)
        @test Bitcoin.bip91(block)
        block_raw = hex2bytes("020000208ec39428b17323fa0ddec8e887b4a7c53b8c0a0a220cfd0000000000000000005b0750fce0a889502d40508d39576821155e9c9e3f5c3157f961db38fd8b25be1e77a759e93c0118a4ffd71d")
        stream = IOBuffer(block_raw)
        block = Bitcoin.io2blockheader(stream)
        @test !Bitcoin.bip91(block)
    end
    @testset "BIP141" begin
        block_raw = hex2bytes("020000208ec39428b17323fa0ddec8e887b4a7c53b8c0a0a220cfd0000000000000000005b0750fce0a889502d40508d39576821155e9c9e3f5c3157f961db38fd8b25be1e77a759e93c0118a4ffd71d")
        stream = IOBuffer(block_raw)
        block = Bitcoin.io2blockheader(stream)
        @test Bitcoin.bip141(block)
        block_raw = hex2bytes("0000002066f09203c1cf5ef1531f24ed21b1915ae9abeb691f0d2e0100000000000000003de0976428ce56125351bae62c5b8b8c79d8297c702ea05d60feabb4ed188b59c36fa759e93c0118b74b2618")
        stream = IOBuffer(block_raw)
        block = Bitcoin.io2blockheader(stream)
        @test !Bitcoin.bip141(block)
    end
    @testset "Target" begin
        block_raw = hex2bytes("020000208ec39428b17323fa0ddec8e887b4a7c53b8c0a0a220cfd0000000000000000005b0750fce0a889502d40508d39576821155e9c9e3f5c3157f961db38fd8b25be1e77a759e93c0118a4ffd71d")
        stream = IOBuffer(block_raw)
        block = Bitcoin.io2blockheader(stream)
        @test target(block) == parse(BigInt, "13ce9000000000000000000000000000000000000000000", base=16)
        @test difficulty(block) == 888171856257
    end
    @testset "Check POW" begin
        block_raw = hex2bytes("04000000fbedbbf0cfdaf278c094f187f2eb987c86a199da22bbb20400000000000000007b7697b29129648fa08b4bcd13c9d5e60abb973a1efac9c8d573c71c807c56c3d6213557faa80518c3737ec1")
        stream = IOBuffer(block_raw)
        block = Bitcoin.io2blockheader(stream)
        @test check_pow(block) == true
        block_raw = hex2bytes("04000000fbedbbf0cfdaf278c094f187f2eb987c86a199da22bbb20400000000000000007b7697b29129648fa08b4bcd13c9d5e60abb973a1efac9c8d573c71c807c56c3d6213557faa80518c3737ec0")
        stream = IOBuffer(block_raw)
        block = Bitcoin.io2blockheader(stream)
        @test !check_pow(block)
    end
    @testset "Validate merkle root" begin
        hashes_hex = [
           "f54cb69e5dc1bd38ee6901e4ec2007a5030e14bdd60afb4d2f3428c88eea17c1",
           "c57c2d678da0a7ee8cfa058f1cf49bfcb00ae21eda966640e312b464414731c1",
           "b027077c94668a84a5d0e72ac0020bae3838cb7f9ee3fa4e81d1eecf6eda91f3",
           "8131a1b8ec3a815b4800b43dff6c6963c75193c4190ec946b93245a9928a233d",
           "ae7d63ffcb3ae2bc0681eca0df10dda3ca36dedb9dbf49e33c5fbe33262f0910",
           "61a14b1bbdcdda8a22e61036839e8b110913832efd4b086948a6a64fd5b3377d",
           "fc7051c8b536ac87344c5497595d5d2ffdaba471c73fae15fe9228547ea71881",
           "77386a46e26f69b3cd435aa4faac932027f58d0b7252e62fb6c9c2489887f6df",
           "59cbc055ccd26a2c4c4df2770382c7fea135c56d9e75d3f758ac465f74c025b8",
           "7c2bf5687f19785a61be9f46e031ba041c7f93e2b7e9212799d84ba052395195",
           "08598eebd94c18b0d59ac921e9ba99e2b8ab7d9fccde7d44f2bd4d5e2e726d2e",
           "f0bb99ef46b029dd6f714e4b12a7d796258c48fee57324ebdc0bbc4700753ab1"]
       hashes = [hex2bytes(x) for x in hashes_hex]
       stream = IOBuffer(hex2bytes("00000020fcb19f7895db08cadc9573e7915e3919fb76d59868a51d995201000000000000acbcab8bcc1af95d8d563b77d24c3d19b18f1486383d75a5085c4e86c86beed691cfa85916ca061a00000000"))
       block = Bitcoin.io2blockheader(stream)
       @test Bitcoin.validate_merkle_root(block, hashes)
    end
end