While working on a major rewrite of the Datetime Julia library for Dates and DateTimes, increased performance was a key goal. And though the original implementation was already fast, Julia developers tend to be greedy, so naturally I wanted more.

One of the known bottlenecks at the time, was how leap seconds were being dealt with. The DateTime type accounted for leap seconds like any other second, and naturally this complicated the internals a bit. See, the internal type is just stored as a machine instant, a monotonically increasing number from zero.

julia> dt = datetime(2014,6,19)  
2014-06-19T00:00:00 UTC

julia> int(dt)  
63538819225000  

So the tricky part, comes in the algorithm to calculate this machine instant given the parts of a date (year, month, day, etc.). In a leap-seconds-don’t-exist kind of world, this is trivial:

totalms = ms + 1000*(s + 60mi + 3600h + 86400*totaldays(y,m,d))  

The hard part, is when trying to convert the above into leap-seconds-do-exist world.

Without going into too much detail, the solution basically involves using two cached arrays, one in the leap second timeline specifying the instant a leap second occurred; while the other is on a non-leap second timeline specifying the instant right before a leap second would occur (yes, it totally makes me think of alternate timelines).

The bottleneck problem boils down to the fact that every time you create or show a DateTime, you have to do a lookup into these leap second arrays to figure out how many leap seconds need to be corrected in either setting the machine instant, or displaying the right date parts.

The natural thing here is to employ binary search, which Base Julia provides an excellent implmentation of in the searchsorted family of functions. Basically, we calculate our totalms from above, then do binary search to figure out when the previous leap second instant occured and how many leap seconds it represents and need to be corrected for our calculated instant (in this case, we could use the searchsortedlast function to get the last leap second instant before our calculated instant). Simple enough, right? End of story? Move on?

But let’s step back a minute and think about the nature of the problem. A major insight is the fact that our leap second arrays are static–nothing is being added or removed apart from library releases and when the IERS decides things are getting too out of whack with earth’s rotation. So how could we use this knowledge to improve performance? Well, if we take our searchsortedlast binary search to the extreme, it would ideally look like:

const LEAPSECONDS = [...]  
function setleaps(totalms)  
    if totalms < 62624707200000  # our middle value of LEAPSECONDS
        if totalms < 62388144000000 # cut 1st half in half again
            if totalms < 62230377600000 # last cut needed
                return 1  # only one leap second has occured
            else
                return 2  # 2 leap seconds
            end
        else
        ...
    else
    ...
end  

Basically, a hand-written searchsortedlast binary search tree using our static array values that elimates all recursive/function call overhead and is close to the metal. But who wants to write all that out by hand? And when our leap second arrays do get updated, we’ll have to manually rebalance our tree, and hand-write it all over again…..yuck. But wait, I seem to recall a famous PR by Steven Johnson where he used a macro for hand-unrolling array coefficients that other libraries wouldn’t dream of trying to maintain. Perhaps we too can leverage Julia’s meta-programming capabilities to do our dirty work.

# Recursively build binary search tree w/ known lookup values
# A is sorted array of lookup values
# R is return values for each index of A
# i.e. R[1] is returned for values < A[1], R[2] for < A[2], etc.
function searchsortedlasttree(A,R,var_name)  
    l = length(A)
    mid = iseven(l) ? l>>>1 : (l>>>1)+1
    # Handle base case
    if mid == 1
        if l == 1
            return :($var_name < $(A[1]) ? $(R[1]) : $(R[2]))
        else # l == 2
            return :($var_name < $(A[1]) ? $(R[1]) : 
                     $var_name < $(A[2]) ? $(R[2]) : $(R[3]))
        end
    end
    iftree = Expr(:if)
    iftree.args = Array(Any,3)
    iftree.args[1] = :($var_name < $(A[mid])) # condition
    iftree.args[2] = searchsortedlasttree(A[1:mid-1],R[1:mid],var_name)
    iftree.args[3] = searchsortedlasttree(A[mid+1:end],R[mid+1:end],var_name)
    return iftree
end

macro leapstree(a,r,var_name)  
    # var_name is the variable name that will be passed in thru the parent function
    A = eval(a)  # safe to eval here because 'a' is known at compile time
    R = eval(r)  # same here
    ret = Expr(:block)  # we'll store everything in a big block expression
    # First we manually handle before and after the ends of our leaps array
    push!(ret.args,:($var_name < $(A[1]) && return 0))
    push!(ret.args,:($var_name >= $(A[end]) && return $(endof(A)*1000)))
    # Now we call the recursive searchsortedlasttree to get the rest
    push!(ret.args,searchsortedlasttree(A[2:(endof(A)-1)],R,var_name))
    return ret
end

const SETLEAPS = [62214480000000,62230377600000,62261913600000, 62293449600000,62324985600000,62356608000000,62388144000000, 62419680000000,62451216000000,62498476800000,62530012800000, 62561548800000,62624707200000,62703676800000,62766835200000, 62798371200000,62845632000000,62877168000000,62908704000000, 62956137600000,63003398400000,63050832000000,63271756800000, 63366451200000,63476784000000]

function setleaps(ms)  
    # @leapstree as a macro ensures everything gets expanded
    # at compile time
    @leapstree(SETLEAPS,[1000:1000:((endof(SETLEAPS)-1)*1000)],ms)
end  

I tried to add useful comments so you can follow what’s going on above, but basically we have our setleaps function that will determine how many leap seconds need to be corrected given a totalms value. setleaps is built at compile time by expanding the macro leapstree which manually builds a binary search tree given the SETLEAPS array of leap second instants and the return values for each slot. The binary search tree is built primarily by searchsortedlasttree, which is similar to the searchsortedlast implementation, except instead of returning values, it returns expressions. It is also built recursively to handle any depth of tree.

Hmmm….seems like a bit of hassle, does it really improve things much? Well, for one thing, the above code is robust, I’ve written it once and don’t have to do anything more except tack on the next leap second to SETLEAPS the next time one is announced. Our manually built binary search tree will be rebalanced automatically. What about performance?

julia> @time [searchsortedlast(SETLEAPS, 63366451200000) for i = 1:10000000];  
elapsed time: 0.50351136 seconds (80000048 bytes allocated)

julia> @time [setleaps(63366451200000) for i = 1:10000000];  
elapsed time: 0.067922885 seconds (80000048 bytes allocated)  

Boom! Almost 10x!

Anyway, I thought this was an interesting application of meta-programming in Julia where we’re basically taking static data and generating an operation on a data structure in the code itself. A recent thread in the Julia-Users forum asked about the long-term hopes of Julia achieving true C/Fortran performance (currently usually 1.2-2x). I think another strong point to the argument of Julia vs. C/Fortran is what I’ve shown above. Powerful language design choices like strong meta-programming functionality allows you to do things in Julia that would be almost impossible in many other languages, but that also allow for unique performance gains.

Stay tuned for more adventures in Julia land and maybe for my next post, I’ll talk a little more about Julia’s code introspection tools and how they can be so darn handy while developing. As a teaser, here’s the generated code for our setleaps method above:

In  [62]: @code_lowered setleaps(62561548800000)

Out [62]: 1-element Array{Any,1}:  
 :($(Expr(:lambda, {:ms}, {{},{{:ms,:Any,0}},{}}, :(begin  # In[30], line 43:
        unless ms < 62214480000000 goto 0
        return 0
        0: 
        unless ms >= 63476784000000 goto 1
        return 25000
        1: 
        unless ms < 62624707200000 goto 15
        unless ms < 62388144000000 goto 8
        unless ms < 62293449600000 goto 4
        unless ms < 62230377600000 goto 2
        return 1000
        2: 
        unless ms < 62261913600000 goto 3
        return 2000
        3: 
        return 3000
        goto 7
        4: 
        unless ms < 62324985600000 goto 5
        return 4000
        5: 
        unless ms < 62356608000000 goto 6
        return 5000
        6: 
        return 6000
        7: 
        goto 14
        8: 
        unless ms < 62498476800000 goto 11
        unless ms < 62419680000000 goto 9
        return 7000
        9: 
        unless ms < 62451216000000 goto 10
        return 8000
        10: 
        return 9000
        goto 14
        11: 
        unless ms < 62530012800000 goto 12
        return 10000
        12: 
        unless ms < 62561548800000 goto 13
        return 11000
        13: 
        return 12000
        14: 
        goto 28
        15: 
        unless ms < 62908704000000 goto 22
        unless ms < 62798371200000 goto 18
        unless ms < 62703676800000 goto 16
        return 13000
        16: 
        unless ms < 62766835200000 goto 17
        return 14000
        17: 
        return 15000
        goto 21
        18: 
        unless ms < 62845632000000 goto 19
        return 16000
        19: 
        unless ms < 62877168000000 goto 20
        return 17000
        20: 
        return 18000
        21: 
        goto 28
        22: 
        unless ms < 63050832000000 goto 25
        unless ms < 62956137600000 goto 23
        return 19000
        23: 
        unless ms < 63003398400000 goto 24
        return 20000
        24: 
        return 21000
        goto 28
        25: 
        unless ms < 63271756800000 goto 26
        return 22000
        26: 
        unless ms < 63366451200000 goto 27
        return 23000
        27: 
        return 24000
        28: 
    end))))