genomic-intersection.jl

#!/usr/bin/env julia

using IntervalTrees
using Bio.Intervals: Interval

function random_interval(seqname, minlen, maxlen, minstart, maxstart)
    first = rand(minstart:maxstart)
    last = first + rand(minlen:maxlen) - 1
    return Interval(seqname[1:end], first, last)
end


function run()
    srand(1234)
    seqname = "1"
    minstart = 1
    maxstart = 1000000
    minlen = 1
    maxlen = 10000
    N = 100000
    TreeType = IntervalTree{Int64,Interval{Void}}
    xs = [random_interval(seqname, minlen, maxlen, minstart, maxstart) for _ in 1:N]
    ys = [random_interval(seqname, minlen, maxlen, minstart, maxstart) for _ in 1:N]
    sort!(xs)
    sort!(ys)

    us = TreeType(xs)
    vs = TreeType(ys)

    #method = :iterative
    method = :successive

    for (u, v) in intersect(us, vs, method=method)
    end

    #count = 0
    #@profile for (u, v) in intersect(us, vs, method=method)
        #count += 1
    #end
    #@show count
    #Profile.print()

    count = 0
    @time for (u, v) in intersect(us, vs, method=method)
        count += 1
    end
    @show count
end

run()