Idempotent + Moving Window is simply a reduction

I want to calculate the moving window max or min in a data parallel way.
Lets start with the two argument function max

max(x,y) returns the greater of x and y
max(x,x) returns x (idempotent)

max(5,3) return 5

max(max(5,3),max(5,3)) reduces to max(5,5) which just returns 5

If we want max of an entire list we can simply think of it as a reduction in the lisp/APL sense

max(head (x), (max (head (tail (x), max( head(tail(tail(x))), …. )

or in a more readable way, replace max with |, and insert it between every element in the list

x[0] | x[1] | x[2] | x[3] this is a standard reduction/over.

Here is a concrete example:

x:5 4 3 2 7 2 9 1

Max over(x) –> 9

|/x –> 9

We can look at the intermediate results

Max scan(x) –> 5 5 5 5 7 7 9 9

|\x –> 5 5 5 5 7 7 9 9

Now let’s say that we want to look at the 3 slice moving window.

Let’s take advantage of the fact that max(x,x) yields x, idempotent

We can calculate the max between each pair in our list. Read max each-prior (‘:)

(|’:)x –> 5 5 4 3 7 7 9 9

applying this function n-1 times gives us a moving max

so (|’:) (|’:) x is 3 slice moving window, which we can rewrite as
2 |’:/ x

5 5 5 4 7 7 9 9

We can see that |\x is equivalent to (count[x]-1) |’:/ x which is data parallel by construction. In other words, we are doing an adjacent transform.

To make this a bit clearer, we can show intermediate results by using a scan(\) instead of an over(/)

(count[x]-1)|’:\x

5 4 3 2 7 2 9 1
5 5 4 3 7 7 9 9
5 5 5 4 7 7 9 9
5 5 5 5 7 7 9 9
5 5 5 5 7 7 9 9
5 5 5 5 7 7 9 9
5 5 5 5 7 7 9 9
5 5 5 5 7 7 9 9

One thing we can notice is that we could terminate early, since we know that if an adjacent element did not change there is no way for the max value to propagate further.

Which means we can rewrite (count[x]-1)|’:\x to simply (|’:\)x

and indeed this gives us:

(|’:\)x

5 4 3 2 7 2 9 1
5 5 4 3 7 7 9 9
5 5 5 4 7 7 9 9
5 5 5 5 7 7 9 9

Technically this is still slower than maxs(x) but if we had gpu support for each-prior(‘:) we could calculate maxs in n calls to max prior which has a parallelization factor of n/2. Depending on the range of x, the n calls might be bound significantly such that for small ranges, for example you are dealing with a short int (8 bytes), you can effectively compute the maxs in O(c*n/2) time with parallelism. Where c is a function of the effective range of your inputs — max[x]-min[x].

Now let’s get back to the problem of max in a sliding window, like this leetcode problem. This problem is classified as hard and the reason is that leetcode doesn’t want you to use parallelism to solve sliding window, it wants you to do it in O(n) time directly.

We know we can solve the problem in O(n) time when the window is equal to 1,2, and n.
identity({x}), max each prior ({|’:x}), and maxs ({|\x}) respectively.

The key to solving it for all window sizes is to recognize that each entry only depends on itself and its neighbors and that the result is the same if you duplicate neighbors – again idempotentcy comes to the rescue. max(1 2 3) is equal to max(1 1 2 2 3 3). Another good explanation can be found here. This is technically a dynamic programming approach with two pointers, but I have never seen it explained in what I think is the most straightforward way.
Let’s borrow the example from leetcode and work through what should happen.

Suppose list of numbers n=1 4 3 0 5 2 6 7 and window size k=3. Then reusing the prior code to show intermediate steps, we get:

1 4 3 0 5 2 6 7
1 4 4 3 5 5 6 7
1 4 4 4 5 5 6 7

If we look at a particular slot, 0 for example, it gets overtaken by 4, but the 5 slot only needs to compete with 3. In other words, each element is looking at at most n-1 elements to the right and n-1 elements to the left. So we can rearrange n to do this in a natural way. Let’s reshape n, so that there are k columns. q supports ragged arrays, so our last row is not the same length which is nice.

q)(0N;k)#n
1 4 3
0 5 2
6 7

We can then compute the maxs of each row

q)maxs each (0N;k)#n
1 4 4
0 5 5
6 7

And flatten it back:

q)l:raze maxs each (0N;k)#n
q)l
1 4 4 0 5 5 6 7

This gives us our left window. Now we will repeat the steps but instead of going from the left we will go from the right.
The simplest way go from the right is to reverse the list apply the function as normal from the left and reverse the list again. In J this operation is called under, as in going under anesthesia, perform the operation and then wake you from the anesthesia. This is equivalent to looking at the cummulative max from the right or walking backward through the array.

q)under:{[f;g](f g f ::)}
q)(reverse maxs reverse)~under[reverse;maxs] /this is the same
1b 

q)(reverse maxs reverse ::) each (0N;k)#n
4 4 3
5 5 2
7 7

Now we just need to flatten the list and call it r.

q)r:raze (reverse maxs reverse ::) each (0N;k)#n

Because r is going from the right, we need to shift it by k-1 elements to the right so that we are not using data from the future when looking at what the current max should be. Printing l on top of r, we get:

q)(l;(k-1) xprev r)
1 4 4 0 5 5 6 7
0N 0N 4 4 3 5 5 2

As we can see, a simple max along the columns will give us the correct answer. (Note: the 0N in the beginning correspond to nulls)

q)max (l;(k-1) xprev r)
1 4 4 4 5 5 6 7

This idea can be generalized, so let’s write a generic function that can take an idempotent operation and create a moving window version.

fmx:{[f;g;h;m;x]l:raze (f')w:(0N;m)#x;r:raze (g')w;h[l;(m-1) xprev r]}

fmx takes a function f that will be applied to generate the left window, a function g that will be used to generate the right window and h a function that will combine left and right windows shifting r based on the window.

We can now generate (mmax, mmin, mfill):
fmmax:fmx[maxs;(reverse maxs reverse ::);|]
fmmin:fmx[mins;(reverse mins reverse ::);{&[x;x^y]}] /fill the nulls with x
fmfill:fmx[fills;(reverse ({y^x}) reverse ::);{y^x}]

We can also generate a function that will allow us to inspect what elements we have access to in the right and left windows so that we can debug/make new functions, with some small modifications.

inspect:fmx[,\;(reverse (,\) reverse ::);{([]l:x;r:y)}]

Here we define f to concatenate the elements seen so far, g is the reverse concatenate, and h is return a table of the l and r.

When we run inspect on our original n, we see that every row has the information from the appropriate sliding window in it, though sometimes more than once.

q)inspect[k;n]
l      r       
--------------
,1     `long$()
1 4    `long$()
1 4 3  1 4 3   
,0     4 3     
0 5    3       
0 5 2  0 5 2   
,6     5 2     
6 7    2       

We can see that the combination of left and right windows will always have at least the three elements required.

To summarize, we saw that idempotent functions, like max and min, allow for parallelization and allow us to use the dynamic programming two pointer technique to solve sliding window calculation in O(n).



Below is the python code for sliding maximum window, written in numpy/APL style, I’m not sure it makes the concept clearer, but more people read python than q. numpy doesn’t like ragged arrays, so there is a bit of extra code to handle the cases where count of n is not evenly divisible by k.

def maxSlidingWindow(k,nums)
  cmax=np.maximum.accumulate
  n=len(nums)
  z=np.zeros(k*np.ceil(n/k).astype(int))
  z[:n]=nums
  z=z.reshape(-1,k)
  l=np.resize(cmax(z,1).reshape(z.size),n)
  r=np.resize(np.flip(cmax(np.flip(z,1),1),1).reshape(z.size),n)
  return list(np.max(np.stack([l[k-1:],r[:r.size-(k-1)]]),0).astype(int))

all code can be found here:
https://github.com/pindash/q_misc/blob/master/sliding_window_max.q

Advertisement

Terse Code

According to this nature article additive solutions are preferred to subtractive solutions. This heuristic may encode a form of Chesterton’s fence, but it blinds people to finding solutions that are better by taking for granted what is. The article says people don’t dismiss the subtractive solutions, rather they never consider them in the first place. The idea behind terse code, may simply reverse this heuristic thus opening the mind to novel connections that would be hidden by a strictly additive approach. This might explain why constraints are so often a boon to creativity. They take the place of the mind’s general purpose heuristics in narrowing the search space, but leave open areas pertinent to the problem. So for example, trying to fit a piece of code into n chars, forces the user to rethink complicated methodologies and solve just the core problem. Other constraints that work similarly are constraints on performance, usually solutions that achieve 80% or more of the result can be achieved by simplifying the objective. Such bouts of creativity only come from constraining the resources that can be marshaled at a given task. A great example of this is the SVM technique used early to recognize hand written addresses by the US post office. The computing power and training sets available were a small fraction of what we have today. Neural networks are a counter example, they seem to work better as they get larger. Perhaps that explains our own bias, or we might discover that some constraint on the networks yields vast improvements, only time will tell. My sympathies lie with occam.

Timeouts

Timeouts are a hack.
Timeouts can be worse than a hack.
Timeouts with Retry logic are almost always a recipe for disaster.
Because Timeouts are at best a hack they should be avoided.

Today, I was upgrading our micro-service backtester infra. At the core layer, there is a piece of code that is responsible for connecting services to other services. This code has clearly been a source of trouble in the past. Signs of wear and tear include: copious use of globals, lots of debugging statements, and additions of randomness in choosing values. Any micro service based architecture is always going to get stressed in one place and that is at the edges between services. Having verbose debugging on it’s own is a good thing, but almost always is an after thought. Adding globals is a good way to allow inspection of a live system. Adding randomness when choosing defaults for retries and timeouts, well that’s a whole other level.

What went wrong? At the face of it, inter process requests are simple, especially when they can rely on tcp. A tcp connection is instantiated an asynchronous request is made and a reply is expected. On a sunny day, this works! Suppose the service you are requesting is either busy or dead? How would you tell the difference. Well if it’s dead, you can’t instantiate that tcp connection. One option is on learning that the service you rely on is dead you should just die as well. But what if there was a network partition and if you had only tried one more time, you could have recovered. Enter retry logic. What if the first attempt you made to connect somehow got lost, you should probably just give up after some time and try again later. What if you connected, but subsequently the service died. Well if you don’t hear from the service within some reasonable amount of time, you should try to reach another viable service. Well between all of these many retries and timeouts, if a bunch of services all try to reach one instance in particular that instance can freeze. Then they will all try again and this cycle will persist. So adding a bit of randomness can help keep a bunch of clients from locking one service up. Essentially each client agrees to try again, but in an unspecified amount of time, so that a coordinated request is unlikely. However, this just defers the problem, because if the amount of time that it takes to service all the requests exceeds the maximum amount of randomness inserted into the configuration. The service will still feel like everyone is asking it to do something at the same time. So today, when I restarted the system, every single service piled in to request the same piece of data from one service, the service got overwhelmed, then they all tried again at random intervals, but the service was still dealing with overflow from previous requests, which it could no longer honor, since they had timed out the previous query. REPEAT.


This leads, me to my first point, timeouts are a hack. Sure, there is probably some case where it is the proper thing to do. However, if you are willing to retry the same connection 3 times for 10 seconds, you are better off trying once for 30 seconds, especially if you are on TCP, where ordering of messages is guaranteed. If you are going to try different sites/connections, you are still better off trying less times and waiting longer. If you are patient, your probability of success is the same.

Suppose during a pandemic you go to buy toilet paper. When you get there you see a line going out the door, you are not sure if there will be toilet paper left when it is your turn. You balk and then come back 15 minutes later. If you couldn’t do any useful work in those 15 minutes, you might as well have waited on line. The supply of toilet paper does not depend on how many times you balk and rejoin the line. Counting the members on the line is a much better proxy, than a timeout. Asking the other members, how much toilet paper they expect to buy is even better. You might even be able to form a consensus among your peers on whether you expect there to be toilet paper. Perhaps, they have recorded their interactions with the store and know what quality of service you can expect.

All of these alternatives have analogues in a microservice infrastructure, I will not go through them here, they deserve their own post.

We are committed to the timeout/retry gambit. As far as I know the best timeout/retry combo is implemented in tcp. It uses Exponential backoff simply put, when you encounter a failure to reach the service, you try again but each time you wait a multiple(>1) of the period you waited last time before retrying again. We still need randomness. Otherwise, we will have just scheduled collisions exponentially far apart. Instead each of the shoppers multiplies their wait period by a random number between one and two raised to the power of retries. This spreads out the requests evenly each time increasing the amount of time that the requests are spread over accounting for the evidence that repeat failure indicates there are many clients.

Fifo Allocation in KDB

One of the most common questions in trading is how to attribute profits and losses. There are many ways of doing this each have their own advantages one of the simplest is first in first out allocation. Here is the basic concept illustrated with an example, I buy 7 apples for 10 dollars, sell 4 for 11 dollars, buy 2 more for 12 and sell 3 for 11 and then sell 2 more for 12. We can think of having a list of purchases 7,2 and list of sales 4,3,2. Each sale and purchase are associated with a price. No matter the allocation strategy, we know the profit can be calculated by a weighted sum of the sales ((4*11)+(3*11)+(2*12)) minus a weighted sum of the purchases ((7*10)+(2*12)) –> 7 dollars. Allocation tries to answer the question of which trades were bad. First in first out will allocate the first 7 sales to the first 7 purchases (7 dollars profit), leaving 2 more purchased and 2 more sold and allocate those to each other (0 profit). This tells us that all the profit was made from the first 7 purchased. However, a last in first out strategy might instead see this ((4*11)-(4*10)+((2*11)-(2*12))+((1*11)-(1*10))+((2*12)-(2*10)) which suggests that the second purchase at 12 dollars was a bad trade.

We leave the discussion for more elaborate allocation methods and their respective merits for another time. In this post we will focus on the mechanics of efficiently calculating fifo allocations. Luckily, Jeff Borror covers an elegant implementation in Q for Mortals, which I will only briefly repeat here. The idea is that you can represent the buys and sells as a matrix, where each cell corresponds to the amount allocated to that purchase and sale. By convention rows will correspond to purchases and columns to sales. So in our example, we can write the allocation as

 | 4 3 2
-| -----
7| 4 3 0
2| 0 0 2

I left the corresponding order quantities in the row and columns as headers but they are actually implied. Jeff also gives us the algorithm that produces this matrix.

  • First we calculate the rolling sums of the purchase and sales
    • 7 9 for purchases
    • 4 7 9 for sales
  • We then take the cross product minimum
    • 4 7 7
      4 7 9
  • We then take the differences along the columns
    • 4 7 7
      0 0 2
  • We then take the differences along the rows
    • 4 3 0
      0 0 2
  • We are done, as a bonus here is the code in Q
    • deltas each deltas sums[buys] &\: sums[sells]

Not only is this rather clever, there is a certain logic that explains how to come to this result. The cumulative sums tells you how much max inventory you have bought or sold till this point. The minimum tells you how much you can allocate so far assuming you haven’t allocated anything. The differences down the columns subtracts the amount you have already allocated to the previous sales. The differences along the rows tells you how much you have already allocated to the previous purchases. Since you can only allocate what hasn’t yet been claimed.

Having read this far, you should feel confident you can easily do fifo allocations in KDB. I know, I did. There are even stack overflow answers that use this method. There is one problem, that occurs the moment you start dealing with a non trivial number of orders. This method uses up n^2 space. We are building a cross product of all the buys and sells. We know that the final matrix will contain mostly zeros, so we should be able to do better. We can use the traditional method for doing fifo allocation. Keep two lists of buys and sells allocated thus far and keep amending the first non zero element of each list and created a list of allocations triplets, (buy, sell, allocated). Although, this is linear implementing this in KDB is rather unKdb like. For incidental reasons, amending data structures repeatedly which is what this algorithm entails is best done by pointers, which in KDB means using globals and pass by reference semantics. It’s not long, but it’s not pretty.

codecomment
fastFifo:{[o]/takes orders signed to neg are sell
signed:abs (o unindex:group signum o);/split into buys and sells
if[any 0=count each signed 1 -1;:([]row:();col:();val:())];/exit early if there are no buys or no sells
k:(/Brains of the calculation called k (for allokate)
.[{[a;i;b;s]/args are allocations triplets(a), index(i), buys(b), sells(s)
if[0=last[get s]&last[get b];:(a;i;b;s)];/if either there are no more buys or sells to allocate, return the current allocations
l:a[;i];l[2]:$[c:b[l 0]<=s[l 1];b[l 0];s[l 1]];/l is the current allocation, if (condition c) is buys are bigger allocate the sell qty, otherwise the buy qty
(.[a;(til 3;i);:;l+(0 1 0;1 0 0)c];i+:1;@[b;l 0;-;l 2];@[s;l 1;-;l 2])}]/depedning on the c increment either buy or sell pointer
);/end definition of k
`.fifo.A set (3;count[o])#0;/create global list of allocation triplets, size of total orders is at least 1 more than max necessary size, consider a case where you have 10 sells for 1 qty, 1 buy for 10, you will have 9 allocations
`.fifo.B set signed 1;/Set the buys
`.fifo.S set signed -1;/Set the sells
res:{delete from x where 0=0^val} update row:unindex[1] row, col:unindex[-1] col, next val from flip `row`col`val!get first k over (`.fifo.A;0;`.fifo.B;`.fifo.S);/delete the cases of 0 allocation, return the original indexes of the orders after apply k until the allocations stop changing, return the allocations as a sparse representation of the matrix
delete A,B,S from `.fifo;/delete the global variables to clean up
res}/return the result

This algorithm, has sufficient performance that we could stop there. However, we could ask is there a way to get all the benefits of the elegant array solution without paying the space cost. The answer is that we can, the trick is that as we have noticed most of the matrix will actually be filled with zeros. In particular, we can see that the matrix will essentially traverse from the top left hand corner to the bottom left hand corner. If we could split the problem into small pieces and then stitch the solutions together we would have the original path of allocations.

I will now briefly sketch out how we can split this problem into pieces and then I will present an annotated version of the code.

If we just look at the first 100 buys and first 100 sells. We can simply apply Jeff’s algorithm. If we wanted to apply it to the next 100 buys and next 100 sells, we would find that we have a problem. We need to know three things, we need to know the index of the buys and sells we are currently up to and any remaining buys and sells that we have not allocated to yet in the previous iteration. Strictly speaking we can only have unallocated quantities on one side, but it is easier to simply keep track of both and letting one list be empty each time.

Here is the annotated code:

codecomment
fifoIterative2:{[o]/takes orders signed, neg are sell
signed:abs o unindex:group signum o;/split into buys and sells
iterate:{[a;b;s]/Brains of the calculation called iterate takes accumulator (a), chunks of (b), chunk of sells (s), accumulator has three slots, slot 0 is for the current count where we are up to in the buys and sells, slot 1 keeps track of the unallocated qty from the previous iteration, slot 3 keeps track of the currrent allocations
if[any 0=sum each (b:a[1;`b],:b;s:a[1;`s],:s);:a];/concatenate the unallocated qty to the respective places and check if one is an empty list, if so exit early and return the accumulator
t:sm m:deltas each deltas sums[b]&\:sums s;/apply Jeff’s algorithm and shrink matrix to sparse rep
c:(sum .[0^m;] ::) each 2 2#(::),i:last[t][`col`row];/Sum the last col, and last row allocated so that we can calculate the amount allocated to the last buy and sell respectively
a[1]:`s`b!.[i _’ (s;b);(::;0);-;c];/drop all the buys and sells that were fully allocated and subtract the amount allocated on the last buy and sell, from the first of what’s left assign this to slot 1, which holds the remainder
a[2],:update row+a[0;`b], col+a[0;`s] from t;/Append to slot 2 the current allocations, updating the indexes using the counts from slot 0
a[0]+:(`s`b!(count[s];count[b]))-count each a[1];/update the count, adding to the count the currently allocated and subtracing the remainder
a};/return a
k:max div[;200] count each signed;/calculate the max  number of rows so that each buys and sells has at most 200 elements, configurable to the maximize efficiency of the machine
res:last iterate/[(`s`b!(0;0);`s`b!(();());sm 1 1#0);(k;0N)#signed 1;(k;0N)#signed -1];/reshape the buys and sells and iterate over the chunks, select the last element which is the allocations
update row:unindex[1]row,col:unindex[-1] col from res}/reindex the buys and sells into the original indexes provided by the input o

This code actually performs much faster than the iterative version and is shorter. This code essentially has no branching except where there are no more allocations on one side and exits early.

If anyone has a better way of computing fifo allocation let me know in the comments.

Below are some timings and graphs of memory usage. The code can be found here.

Evolutes

This next article was inspired by two sources:

An Advent of Code puzzle from 2017 

And a playing with J article:

An evolute is a matrix whose cells are numbered to spiral out from the center.

17  16  15  14  13
18   5   4   3  12
19   6   1   2  11
20   7   8   9  10
21  22  23---> ...

The first question asks us to find the Manhattan distance between some positive integer cell and the 1st cell.

I see two ways to do this:

  1. Build the evolute, find the location of the integer and then calculate the distance
  2. Calculate the coordinate with respect to the center directly and take the sum of the absolute value of the coordinates. 

I will start with the second method because it is simpler. If we look at the structure of the evolute we will see that each new layer of the matrix will have an odd square in it’s bottom right corner. 1 9 25 49 ….

That means we can locate the layer by finding the nearest odd square root. Here is the code for that:

f:{j:floor sqrt x; j - not j mod 2}
q)f 10
3
q)f 25
5
q)f 24
3

Next we can find the grid coordinate or how far we are from the center for that corner. I plugged in max x so we can use the function on lists of numbers. The ‘?’ verb is being used to find index which corresponds to the steps away from the center.

g:{(1+2*til max x)?f x} 
q)g 10
1
q)g 25
2
q)g 24
1

Now all we need to do is figure out where we are on the layer.

We can be in one of 5 places:

  1. Exactly at the end of a layer
  2. Between the bottom right corner and the top right corner 
  3. Between the top right corner and the top left corner
  4. Between the top left corner and the bottom left corner
  5. Between the bottom left corner and the end of the layer

We can divide by the size of the layer to calculate this. We can also find the remainder to know how far between we are. This gives us the following code:

place:{(x-j*j) div 1+j:f x}
offset:{(x-j*j) mod 1+j:f x}
q)place 1+til 25
0 0 1 1 2 2 3 3 0 0 0 0 1 1 1 1 2 2 2 2 3 3 3 3 0
q)offset 1+til 25
0 1 0 1 0 1 0 1 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0

This allows us to write the full function to calculate the coordinate. All we need to do is to check which of those conditions are met fill in the coordinates from g, plus the appropriate offset and correct sign. 

coordinate:{[x]
x:(),x;
f:{j:floor sqrt x; j - not j mod 2};
g:{(1+2til max x)?y}; place:{(x-yy) div 1+y};
offset:{(x-yy) mod 1+y}; j:f x; brc:x=jj;
grid:g[x;j];
p:place[x;j];
o:offset[x;j];
?[brc;flip (grid;neg[grid]);
?[0=p;flip (grid+1;neg[grid]+o);
?[1=p;flip (grid+1-o;grid+1);
?[2=p;flip (neg[grid+1];1+grid-o);
flip (neg[grid+1]+o;neg[1+grid])]]]]
}
q)coordinate 1+til 9
0 0
1 1
1 1
0 1
-1 1
-1 0
-1 -1
0 -1
1 -1

At this point we can find the coordinate of any cell number. Since we can do this for all the points, we should be able to create the evolute for any dimension by generating all the coordinates and then shifting them to the appropriate column/row number and inserting them into that matrix. We can do this cleverly by sorting the row and column separately, since we want a particular orientation and sorting each of them corresponds to either a horizontal transposition or a vertical transposition. 

evol:{t:`row`col!(x div 2)+flip cordinate j:1+til x*x; (x;x)#exec val from `col xdesc `row xasc flip update val:j from t} 
q)evol 7
37 36 35 34 33 32 31
38 17 16 15 14 13 30
39 18 5 4 3 12 29
40 19 6 1 2 11 28
41 20 7 8 9 10 27
42 21 22 23 24 25 26
43 44 45 46 47 48 49

Now that we can create a evolute and we see that given an evolute we can find the coordinates, we also see that permuting the row and column indexes gives us different orientations. This leads us to see how Joey Tuttle/Eugene McDonnell created an evolute from  scratch. Let’s raze an evolute:

q)raze evol 5
17 16 15 14 13 18 5 4 3 12 19 6 1 2 11 20 7 8 9 10 21 22 23 24 25

Now let us permute it by the magnitude and then look at the differences:

q)iasc raze evol 5
12 13 8 7 6 11 16 17 18 19 14 9 4 3 2 1 0 5 10 15 20 21 22 23 24
q)deltas iasc raze evol 5
12 1 -5 -1 -1 5 5 1 1 1 -5 -5 -5 -1 -1 -1 -1 5 5 5 5 1 1 1 1
q)deltas iasc raze evol 3
4 1 -3 -1 -1 3 3 1 1

We now see a pretty clear pattern. We are repeating a pattern of 1,neg x,-1,x. 
Each time we are taking 2 elements from the pattern. We take an increasing number of them and continue until we take x elements in shot. Here is Eugene’s explanation translated into q

f:{-1+2*x}
g:{(-1;x;1;neg x)}
k:{(f x)#g x}
h:{1+til x}
j:{-1 _ raze 2#'h x}
l:{-1 rotate raze (j x)#'k x}
m:{iasc sums l x}
n:{(x;x)#m x}
q)f 5
9
q)g 5
- -1 5 1 -5
q)k 5
-1 5 1 -5 -1 5 1 -5 -1
q)h 5
1 2 3 4 5
q)j 5
1 1 2 2 3 3 4 4 5
q)l 5
-1 -1 5 1 1 -5 -5 -1 -1 -1 5 5 5 1 1 1 1 -5 -5 -5 -5 -1 -1 -1 -1
q)m 5
24 23 22 21 20 9 8 7 6 19 10 1 0 5 18 11 2 3 4 17 12 13 14 15 16
q)n 5
24 23 22 21 20
9 8 7 6 19
10 1 0 5 18
11 2 3 4 17
12 13 14 15 16

Combining this into function we get:

evolute:{(x;x)#iasc sums -1 rotate raze (-1 _ raze 2#'1+til x)#'(-1+2*x)#(1;neg x;-1;x)}

We can shorten it just a bit and make it a bit faster by seeing that grabbing parts j,k,l can make use of how kdb overloads ‘where’. In KDB ‘where’ gives you the indexes of the 1s in a boolean mask. For example:
q)where 101b
0 2
However, a little thought reveals that we are returning the index of the element the number at that index. So we return one 0, zero 1,one 2. Generalizing this, we can think that where of 1 2 3 should return one 0,two 1s and three 2s and indeed KDB does this.
q) where 1 2 3
0 1 1 2 2 2
Knowing this, we can rewrite j k and l, which make heavy use of each and remove all the razing. 
I am going to show how I built this up. 

q){((x-1)#2),1} 5
2 2 2 2 1
q){1+where ((x-1)#2),1} 5
1 1 2 2 3 3 4 4 5
q){(where 1+(where ((x-1)#2),1))} 5
0 1 2 2 3 3 4 4 4 5 5 5 6 6 6 6 7 7 7 7 8 8 8 8 8
q){(where 1+(where ((x-1)#2),1)) mod 4} 5 /so that we cycle
0 1 2 2 3 3 0 0 0 1 1 1 2 2 2 2 3 3 3 3 0 0 0 0 0
q){(1;neg x;-1;x)(where 1+(where ((x-1)#2),1)) mod 4} 5
1 -5 -1 -1 5 5 1 1 1 -5 -5 -5 -1 -1 -1 -1 5 5 5 5 1 1 1 1 1
q){-1 rotate (1;neg x;-1;x)(where 1+(where ((x-1)#2),1)) mod 4} 5
1 1 -5 -1 -1 5 5 1 1 1 -5 -5 -5 -1 -1 -1 -1 5 5 5 5 1 1 1 1

The last step recreates l. So the whole function using where looks like this:

evolute:{(x;x)#iasc sums -1 rotate (1;neg x;-1;x)(where 1+(where ((x-1)#2),1)) mod 4}

Part 2 of the advent of code is even more interesting. It describes us filling the spiral by summing the neighbors, the center square starts at 1. Here is the description from advent.

Then, in the same allocation order as shown above, they store the sum of the values in all adjacent squares, including diagonals.
So, the first few squares' values are chosen as follows:
Square 1 starts with the value 1.
Square 2 has only one adjacent filled square (with value 1), so it also stores 1.
Square 3 has both of the above squares as neighbors and stores the sum of their values, 2.
Square 4 has all three of the aforementioned squares as neighbors and stores the sum of their values, 4.
Square 5 only has the first and fourth squares as neighbors, so it gets the value 5.
Once a square is written, its value does not change. Therefore, the first few squares would receive the following values:
147  142  133  122   59
304    5    4    2   57
330   10    1    1   54
351   11   23   25   26
362  747  806--->   ...

What is the first value written that is larger than your puzzle input?

Now given our earlier coordinate function, we just need to add elements in into empty matrix of 0s and sum the neighbors at each step to determine what to add.  First, let’s write the simple pieces of finding the neighbors and summing them. Then we will conquer the composition.

/To find the neighbors according to our rule, 
/ we basically find the coordinates for (1..9)
/ then we can add to x pair so it's centered around that one
neighbors:{x+/:coordinate 1+til 9}
q)neighbors (1;0)
0 1
1 1
1 2
0 2
-1 2
-1 1
-1 0
0 0
1 0
/We need to shift our coordinate system depending on the size of the matrix
shift:{i:div[count x;2]; (i+y)}
q)shift[3 3#0;coordinate 5]
0 2
/calculate the sum of the neighbors
sumN:{sum over .[x] each neighbors y}
q)(1+evolute 5)
17 16 15 14 13
18 5 4 3 12
19 6 1 2 11
20 7 8 9 10
21 22 23 24 25
q)(1+evolute 5)[2;3]
2
/Manually sum the neghbors of 2
q)sum 4 3 12 1 2 11 8 9 10
60
q)sumN[1+evolute 5;(2;3)]
60

Okay, we are going to start with a matrix that has single 1 in the center.

center1:{.[(x;x)#0;(a;a:x div 2);:;1]} /we are ammending a 1 at the x div 2
q)center1 5
0 0 0 0 0
0 0 0 0 0
0 0 1 0 0
0 0 0 0 0
0 0 0 0 0

Now that we see how easy it is to amend an element at a particular index, we can combine this with the notion of a fold(over). The idea is that we will keep updating this matrix with new elements based on the neighbor sum:

cumEvolute:{{.[x;y;:;sumN[x;y]]}/[j;shift[j:center1[x]] coordinate 1+til (x*x)]}
q)cumEvolute 5
362 351 330 304 147
747 11 10 5 142
806 23 1 4 133
880 25 1 2 122
931 26 54 57 59
/If we want to rotate it we can simply reverse flip
q)reverse flip cumEvolute 5
147 142 133 122 59
304 5 4 2 57
330 10 1 1 54
351 11 23 25 26
362 747 806 880 931

Using cumEvolute we can find when we generate the first integer bigger than the input, by creating a stopping condition and counting the number of iterations.

stopCum:{{$[z<max over x;x;.[x;y;:;sumN[x;y]]]}/[j;shift[j:center1[x]] coordinate 1+til (x*x);y]}
q)stopCum[5;20]
0 0 0 0 0
0 11 10 5 0
0 23 1 4 0
0 0 1 2 0
0 0 0 0 0
q)max over stopCum[5;20]
23
q)max over stopCum[10;100]
122

As of Join Performance a surprising result

The Kx wiki describes in detail how to structure as of join queries for best performance.

https://code.kx.com/wiki/Reference/aj

3. There is no need to select on quote, i.e. irrespective of the number of quote records, use:

aj[`sym`time;select .. from trade where ..;quote]

instead of

aj[`sym`time;select .. from trade where ..;
             select .. from quote where ..]

The reason for this is that since the quote table is partitioned on date and grouped on sym the as of join function simply scans the sectors of the disk in linear order and grabs the first matches using binary search.

Linear access to disk is very different from random access.  So when you perform a search linearly on partitioned historical database it runs pretty fast.

However, someone asked me if you were going to query over and over, the same data whether then it made sense to do some pre-filtering on the quote table before doing multiple ajs.

At the time, I said it shouldn’t be faster. I thought that the overhead of reading from disk afresh each time was going to be much smaller than the overhead of allocating a very large amount of room in memory.

So I tested it:

taj1:{[t;d] now:.z.T;
do[10;aj[`sym`time;select from t where date=d;select from quote where date=d]];        after:.z.T;
after-now}

taj2:{[t;d]
now:.z.T;
quotecache:select from quote where date=d, sym in exec sym from t;                                do[10;aj[`sym`time;select from t where date=d;quotecache]];
after:.z.T;
after-now}

The timings, were not even close the taj2, which precaches data when run 10 times, for a trade table with only 1000 records on only 100 symbols took more than 5 minutes to run. While taj1 was taking only 5 seconds, that is approximately 2 aj per second.

When I scaled the number of symbols to 1000, the cached version didn’t comeback and I had to cancel the query after 25 minutes.  The non cached version took 10 seconds or 1 aj per second.

The intuition is correct, if you can prefilter, then not having to read in all the data twice should be faster. So I checked the meta of the cached table, it was loosing the p attribute while filtering. 

I then created version taj3:

taj3:{[t;d]
now:.z.T;
quotecache:update `p#sym from select from quote where date=d, sym in exec sym  from t;
do[10;aj[`sym`time;select from t where date=d;quotecache]];
after:.z.T;
after-now}

If you reapply the p attribute the first query for 1000 symbols costs you 3.6 seconds, but increasing the number of times you run the aj is almost free. So running 10 times only costs 4 seconds. Running 20 times the taj3 also took 4 seconds. So if you know you have a limited universe than reducing and pre-filtering is worth it if you will be running the queries again and again, JUST REMEMBER TO REAPPLY THE P ATTRIBUTE. Since the quote table is only filtered, the p attribute will be a really cheap operation, since order is preserved.

Q Idioms assemembered

This is more of an expanding list of Q idioms I have had to either assemble or remember or some combination.

  1. Cross product of two lists is faster in table form

    q)show x:til 3

    0 1 2

    q)show y:2#7

    7 7

    q)x cross y

    0 7

    0 7

    1 7

    1 7

    2 7

    2 7

    q) ([] x) cross ([] y)

    x y

    0 7

    0 7

    1 7

    1 7

    2 7

    2 7

    \t til[1000] cross 1000#5

    85

    \t flip value flip ([] til 1000) cross ([] x1:1000#5)

    33

    /and if you are happy working with it in table form

    \t ([] til 1000) cross ([] x1:1000#5)

    5

  2. Take last N observations for a column in a tableBoth of these can be used to create a list of indexes and the columns can be simply   projected on to the list of indexes

    /a) Intuitive way take last n from c

    q) N:10; C:til 100000;

    q) {[n;c]c{y-x}[til n] each til count c}[N;C] /timing 31 milli

    0

    1 0

    2 1 0

    3 2 1 0

    4 3 2 1 0

    ….

    /b) Faster way using xprev and flip

    q) \t {[n;c] flip (1+til n) xprev \\: c}[10;c]  /timing 9 millesec

  3. Create a Polynomial Function from the coefficients

    poly:{[x] (‘)[wsum[x;];xexp/:[;til count x]]};

    f:poly 0 1 2 3;

    f til 5

    0 6 34 102 228

  4. Count Non Null entries

    All in K:

    fastest:{(#x)-+/^x}

    slower:{+/~^x}

    slowest:{#*=^x}

    Q translation:

    fastest:{count[x]-sum null x}

    slower:{sum not null x}

    slowest:{count first group null x}

  5. Camel case char separated symbols:

camelCase:{[r;c]`$ssr[;r;””] each @'[h;i;:;]upper h @’ i:1+ss'[;r] h:string x}
most common case
q)c:` sv/: `a`b cross `e`fff`ggk cross `r`f
`a.e.r`a.e.f`a.fff.r`a.fff.f`a.ggk.r`a.ggk.f`b.e.r`b.e.f`b.fff.r`b.fff.f`b.ggk.r`b.ggk.f
camelCase[“.”;c]
`aER`aEF`aFffR`aFffF`aGgkR`aGgkF`bER`bEF`bFffR`bFffF`bGgkR`bGgkF

6.  Progress style bars:

p:0;do[10; p+:1; system “sleep .1”;1 “\r “,p#”#”];-1 “”;

{[p]1 “\r “,(a#”#”),((75-a:7h$75*p%100)#” “),string[p],”%”;if[p=100;-1 “”;]}

7.  AutoCorr:

autocorr:{x%first x:x{(y#x)$neg[y]#x}/:c-til c:count x-:avg x}

8.  Index of distinct elements:

idistinct:{first each group x}

Tree Tables/Parent Vector representation of Dictionaries

Introduction:

http://archive.vector.org.uk/art10500340

Vector, the Journal of the British APL Association

TreeTables and the parent vector are magical things and they are a great way to flatten deeply nested structures. In particular, I will use this representation of a dictionary in the next post to allow you to import q libraries under a different namespace, kind of like what python allows with “import x as y”.

In this post I will focus on how we can easily rewrite a dictionary into a flat table with 4 columns. Let me motivate this a bit first. Take a deeply nested dictionary:

a| `b`c!(2;`d`e!(6;`f`g!7 8))
b| `c`d!4 5
c| `d`e!(6;`f`g!7 8)

In json that might look like this:

 
  "a": 
    "b":2,
    "c": 
      "d":6,
      "e": 
        "f":7,
        "g":8
      }
    }
  },
  "b": 
    "c":4,
    "d":5
  },
  "c": 
    "d":6,
    "e": 
      "f":7,
      "g":8
    }
  }
}

As we can see each key is a letter and each value is either another dictionary or a number. In q we can apply functions into these arbitrarily nested dictionaries provided all the values conform. So for example we can add 10+d

a| `b`c!(12;`d`e!(16;`f`g!17 18))
b| `c`d!14 15
c| `d`e!(16;`f`g!17 18)

However, once your dictionaries stop conforming things can get a bit hairy. As a simple example let’s add a simple dictionary f!f to the key f in d

a| `b`c!(2;`d`e!(6;`f`g!7 8))
b| `c`d!4 5
c| `d`e!(6;`f`g!7 8)
f| (,`f)!,`f

Suddenly 10+d throws a type error. Because we can’t add 10 to the non numeric dictionary. We can always right logic to avoid those cases, but it would be simpler to simply pull out all the conforming values add 10 to them and then replace them in the original structure. A tabular tree representation aids in this. As an example the same dictionary in tabular form: (I have omitted some results in the middle for ease of reading.) Which I will call treeTable:

l p  c  d
———
0 0  :: 1
1 0  `a 1
1 0  `b 1
1 0  `c 1
1 0  `f 1
2 1  `b 1
2 1  `c 1

4 15 8  0
5 21 7  0
5 22 8  0

The column c (child) contains every key and value, column d tells you if the row is leaf node or has a dictionary below it. Pulling out all the numeric children is as simple as:

numeric:6 7 8 9h
select c from treeTable where not d, (abs type each c) in numeric
c
-
2
4
5
6
6
7
8
7
8

This table can then  be modified and assuming we can convert it back into its dictionary form we could then  work with nested data easily.

How the Magic is done:

At this point you either believe this structure is useful or you believe it isn’t but you are still interested in understanding how this transformation happens.

Let’s look at the original treeTable and unpack the meaning of the columns.

The first column  represents the level of the dictionary that the key or value is located in.
The second column is the index of the parent row in this table. The root of the table is self parenting meaning that the 0th row is at the 0th level and it’s parent is 0 (keep this in mind it will be useful in a second). All elements will have a parent.
The third column is the child and it is the value at this level of nesting. It will either be a key if there are more levels below or it will just be the value at that level.
The fourth and final column indicates whether or not this row is a dictionary. That is, whether the child should be treated as a key or a value.

To convert a dictionary into a treeTable we use breadth first search, that is the purpose of the column. We first define a primitive treeTable with only one row the root.

l p  c  d
---------
0 0  :: 1

All treeTables will have this row. If the thing we are trying to convert into a treeTable is not a dictionary but is instead SOMEKIND_OF_THING_THAT_IS_NOT_A_DICTIONARY. Then there will be only one more row in this table.

l p  c  d
---------
0 0  :: 1
1 0 SOMEKIND_OF_THING_THAT_IS_NOT_A_DICTIONARY 0

We will then know we are done because there are no more dictionaries to unpack at the last level.

If instead we get a dictionary. We will first record all the keys at that first level and return the table.

Our ability to find a value at a particular level is enabled through the use of the parent column in the treeTable. Suppose we have a simple dictionary:

 
 "d":6,
 "e": 
   "f":7,
   "g":8
  }
}

We can record this as the following two columns:

p c index
----------
0 :: 0 
0 `d 1 
0 `e 2 
2 `f 3 
2 `g 4 
1  6 5 
3  7 6 
4  8 7

The first row is the root it is self parenting. The next two rows are both top level keys. So their parent is the root. f and g both are under e so their parent is 2, but 6 is under d so it’s parent is 1. To make this easier to see, I have added the virtual index column which is always available.  Then the final two rows are under f and g respectively.

If we want the parent of particular row, and we have the parent column, which I will call p:0 0 0 2 2 1 3 4

We can index into that row to get the parent. p 6 -> 3 which is f

We can repeat this until we get to the root. p 3 -> 2; p 2-> 0 ; p 0 -> 0. To find the root in one step we can use KDB’s built-in converge function which will apply until two consecutive results are the same. This is why it was so convenient for the root to be self-parenting. So see the path to the root use scan instead of over.

p over 6 -> 0

p scan 6 -> 3 2 0. Now that we have a path, we just need to get the keys that correspond to that path, this is done by indexing the path against the child column.

c 3 2 0 -> f e ::

We now can get the unique path to any element.

The next step is using the path to index into any level of the dictionary. This is accomplished with a special object called getItems.

getItems is defined by combining the indexing at depth verb with a function that reverses the path list and checks if the path list happens to be only the root. In which case, we simply return the original item.

Using just those two ideas, we are able to construct the treeTable. The algorithm is to index one level at a time each time recording the if the level contains dictionaries or not. If it contains no dictionaries we are done and we will return the same table twice in row, which means that our function will converge. Using breadth first search we avoid any stackoverflow issues that could happen with a recursive solution, instead the function becomes tail recursive, meaning all the necessary ingredients to call the function again are returned as the output. That is why on the first call the function returns the first row of a treeTable. That way each call after simply indexes deeper into the original dictionary to return more levels of the treeTable.

The Code and An Example:

toTreeTable:{[d]
 getItems:('[;] over (.[d;];{$[x~(enlist[::]);x;1_reverse x]}));
 tTT:{[getItems;t] 
 $[98h~type t;;t:([]l:(1#0);p:0;c:(::);d:1b)];
 lev:last t[`l];
 k:exec i from t where l=lev,d;
 $[count k;;:t];
 paths:t[`c] (t[`p]\')k;
 items:getItems each paths;
 id:where bd:99h=type each items;
 p:raze (count each key each items[id])#'k[id];
 c:raze key each items[id];
 df:count[p]#1b;
 lvl:count[p]#lev+1;
 id:where not bd;
 p:p,k[id];
 c:c,items[id];
 df:df,count[id]#0b;
 lvl:lvl,count[id]#lev+1;
 t upsert flip `l`p`c`d!(lvl;p;c;df)}[getItems];
 tTT over ()}



/an example of a nested structure 


b:`c`d!4 5
e:`f`g!7 8
c:`d`e!(6;e)
a:`b`c!(2; c)
d:`a`b`c!(a;b;c)
q)toTreeTable d
l p c d
---------
0 0 :: 1
1 0 `a 1
1 0 `b 1
1 0 `c 1
2 1 `b 1
2 1 `c 1
2 2 `c 1
2 2 `d 1
2 3 `d 1
2 3 `e 1
3 5 `d 1
3 5 `e 1
3 9 `f 1
3 9 `g 1
3 4 2 0
3 6 4 0
3 7 5 0
3 8 6 0
4 11 `f 1
4 11 `g 1
4 10 6 0
4 12 7 0
4 13 8 0
5 18 7 0
5 19 8 0

 

And Back Again!

Now that we covered how to get a treeTable we can also understand how to go back to a dictionary.

We apply the opposite approach. The core function returns a dictionary. Each time we return a dictionary that is slightly deeper than the previous time. We put placeholder empty dictionaries until we build the final result. Since we know whether each row is a key or a value, we know whether the current item requires a placeholder.

toDictFromTreeTable:{[tt]
 tD:{[tt;dSoFar;lev]
 dS:exec {x!count[x]#enlist[()!()]}[c] by p from tt where l=lev, d;
 dS:dS,.[!; value exec p,c from tt where l=lev, not d];
 pR:tt[`c](-1_|:) each (tt[`p]\')[key dS];
 pC:tt[`c]key dS;
 paths:raze each {(1_x;enlist[y])}'[pR;pC];
 $[lev>1;.[;;:;]/[dSoFar;paths;value dS];first value dS]}[tt];
 tD/[()!();1+til last tt[`l]]}

Wow This is Even More General Than We Thought:

When I first built this, I tried to make sure that I covered simple dictionaries and values. So I was curious what would happen to keyed tables.  Keyed tables are special in that they are essentially dictionaries whose key and values are both dictionaries. Since a dictionary is a pair of lists and a list of dictionaries is a table. A keyed table is simply a dictionary whose key is a table and and whose value is a table. A trivial example to illustrate this point:

q)k:([]k:til 5)
k
-
0
1
2
3
4
q)v:([]v:10*til 5)
v 
--
0 
10
20
30
40
q)kv:k!v
k| v 
-| --
0| 0 
1| 10
2| 20
3| 30
4| 40
/Indexing against the key table returns the value table
q)kv[k]
v 
--
0 
10
20
30
40
/but we can also apply select only certain rows using the k column
/in this case I reverse the key table and take the first 2 rows.
q)kv[2#reverse k]
v 
--
40
30


Now what happens if we turn a key table into a treeTable:

l p c d
---------------
0 0 :: 1
1 0 (,`k)!,0 1
1 0 (,`k)!,1 1
1 0 (,`k)!,2 1
1 0 (,`k)!,3 1
1 0 (,`k)!,4 1
2 1 `v 1
2 2 `v 1
2 3 `v 1
2 4 `v 1
2 5 `v 1
3 6 0 0
3 7 10 0
3 8 20 0
3 9 30 0
3 10 40 0

It converts the key part of the table into key dictionaries that are the parents of the value dictionaries in the table. And we can turn it back:

q)toDictFromTreeTable toTreeTable kv
k| v 
-| --
0| 0 
1| 10
2| 20
3| 30
4| 40
q)kv ~toDictFromTreeTable toTreeTable kv
1b

In other words, keyedTables are treated like dictionaries, this means that if you only want to look at values, you will only see values, simply by select from the treeTable where not d. The internal dictionaries inside a keyed table are broken apart into their component dictionaries and the values are stored independently.

Tables Get Treated as singletons.

Since tables are actually lists of dictionaries, and lists are treated as values. A table is also treated as a value and placed directly into the child column.

q)t:([] til 10)
x
-
0
1
2
3
4
5
6
7
8
9
q)toTreeTable t
l p c d
---------------------------------
0 0 :: 1
1 0 +(,`x)!,0 1 2 3 4 5 6 7 8 9 0

The function from TreeTable correctly undoes the toTreeTable function
but the treeTable form is actually more nested than the original table.

q)toDictFromTreeTable toTreeTable t
x
-
0
1
2
3
4
5
6
7
8
9

We can fix this by expanding our parent vector to notate whether a current element is a dictionary, list or an atom. That way we would create a node that is the head of every list and then iterate through the indexes in the list. This is left as an exercise, or until I need this functionality.

A Short Interlude on KDB Load balancing

As you might know KDB is inherently single threaded. In other words, it can do at most one thing at a time in sequence.

This is a good thing. The benefits include:

  • Ordered and inherently sequential transactions that see a consistent view of the world.
  • Fewer Cache misses since data is processed to completion instead of being rescheduled.
  • Programs are easier to reason about.

This post isn’t really about how awesome being single threaded is. There are plenty of articles about that featuring Node.js and javascript. Linked here. Also plenty of stack overflow answers about the subject and an old paper about KDB.

OED

Synchronous: ADJECTIVE

  • Existing or occurring at the same time.

Asynchronous: ADJECTIVE

  • Not existing or occurring at the same time.
  • Computing Telecommunications
    Controlling the timing of operations by the use of pulses sent when the previous operation is completed rather than at regular intervals.

In general, in this single threaded world, you simply design components that do perform some function and they go out to other services that perform other functions. If you do everything over async IPC (Inter-Process-Communication). Then you are golden because all of your q processes are responding to the events as soon as they happen. In a simple example, if you have two q processes. One is responsible for holding realtime data and it’s derived/calculated values. And the other is responsible for calculating those derived values. The updater can keep updating its table and send async events to the model calculator to recalculate. When the model recalculates it can send an async request for the updater to include some values in it’s derived model.

However, perfect systems like this don’t get built. Instead, we have many q processes using synchronous IPC. There is a reason for this –though not a good one. It is significantly easier to reason about a synchronous request. In your mind you simply model all the requests one after another and then you do something with all of the requested values.

Unfortunately, if a request takes a long time to process you can force a service to spend all its time honoring your request. Continuing the example above, suppose someone didn’t understand the initial plan and decides to calculate all the values on the realtime service, that would cause the realtime service to spend more of it’s time processing and less of it’s time updating. We can of course turn off the ability to query a realtime service like this, but that is a bit like throwing the baby out with the bath water. Synchronous requests especially for a client are much easier to understand.

So let’s start with a quick analogy that will help explain what we want to do. Synchronous messages are like phone calls.

You call me, I answer the phone you ask a question, I answer.

Asynchronous messages are like text messages.

You text me. Sometime later I read the message. I might respond.

Ideally, what we would like is to hire a secretary who intercepts incoming calls answers the phone takes a message sends it as a text and keeps the caller on the line until the text is answered at which point the secretary can respond to the caller. Since this job is pretty easy we expect that we can have one secretary serving many callers and keeping them all on the line until their query has been resolved. Unfortunately, as of q version 3.5 we still cannot do this (though it has been promised in version 3.6).

So instead we will simulate this with another mode in KDBQ called multithreaded input mode. Multithreaded input mode allows you to answer queries concurrently (up to 1020) within one q process. However, it places many limitations on what kind of queries it can execute and you don’t want heavy calculations to actually execute concurrently for the reasons described earlier. So we settle for the next best thing. A single master forwards each query to a slave, the slave executes the query returns the result to the master and the master returns the result to the client. This allows us to scale a particular service without meddling with the client code. It turns out a bit of engineering is required to make this work and it does turn the Multithreaded input mode into a router, something e documentation says it is not built to do {YMMV}. ( I have tested it and it seems to work)

DISCLAIMER

This works on linux for version 3.x and on version 3.0-3.4 on mac, it uses UNIX sockets and mac doesn’t support abstract sockets. I have implemented something similar using the fifo system as well, but it has it’s own drawbacks and won’t be described here. I’m sure something similar with files can be done in windows systems, but I haven’t spent the time working it out at this point.

Architecture

You have one master process, this master is run with a negative port number, which will put it into multithreaded input mode. It also runs with a timer turned on so it will be updating and making sure it’s slaves are alive every so N milliseconds. I have set N a thousand so every 1 seconds (this is probably too often but it makes testing faster, since you can see the result of taking a slave offline within a second).

q master.q -p -5000 -t 1000

It runs with a master script called master.q

You also need to run some slave nodes that will listen on their own individual ports. For our example we will use 4001,4002,4003,4004.

q -p 4001

Now I will explain what is inside master.q

//the unix handles of the slave processes
handles:`:unix://4001`:unix://4002`:unix://4003`:unix://4004
h: {@[{hopen x};x;0Ni]} each handles;
han:handles!h;
han:(where han=0Ni)_han;
/ 
 we will define a function that updates all of
 open handles, so that none of them are stale. This will be executed
 in the main thread, which means that all of the slaves should 
 have returned by this point. So if they don't respond it means they
 are dead.
\
.z.ts:{ 
    alive:{@[{x (=;1b;1b)};x;0b]} each han;
    $[all[han] & count [han]=count[handles] ;;
      [dead:where not alive; hclose each han[key dead];
       inactive: handles where (key han) not in handles;
       h:{@[{hopen x};x;0Ni]} each (dead,inactive);
       hd:(where hd=0Ni)_hd:(dead,inactive)!h;
       `han upsert hd;]]};
/
Because each concurrent handle is unique we can assign a slave 
by binning the request handles over the number of active slaves.
In other words we are round robining over the slaves.
\
.z.pg:{[m] (first value (.z.w mod count han)_han) m};
.z.ps:.z.pg;

The reason we have to do all this work, is that any command that changes the global state is illegal and therefore can’t be in .z.pg, and though hopen isn’t illegal it pauses all of the concurrency for all the threads because it effectively changes global state.

As a consequence a serious limitation of this approach is that if a client opens a connection with the this multithreaded gateway while querying they can lock up all the other concurrent requests. In practice this is done by using the connection string query instead of hopen.

For example:

::5000 "2+2"

is equivalent to

h:hopen `::5000; h "2+2"

but the first case will lock up the multithreaded gateway and the second won’t.

Secondly this gateway does not support asynchronous requests with callbacks. That is you can call the async request but you won’t be able to make sure the host calls you back. There are probably some ugly workarounds, like including your host and port in the function.

This is a fully functional load balancer and allows clients to continue to write synchronous code without locking your servers. Technically each of these slaves can actually forward the requests onto other machines, since the slaves are not actually limited.

How the Mighty have Fallen, or How I learned to love GREP.

**Note This post was written over a two week period where I kept learning more about this problem. Many of the notes refer to later parts of the post, so that caveats can be explored. Also some parenthetical notes are almost incomprehensible without reading many of the linked articles on the page.  There is little I can do to summarize all of those articles, this post is already way too long.

This post is inspired by a friend of mine who was solving a particularly common “Big Data” problem. (I enjoyed the last two weeks immensely, and have a much greater appreciation for Ken Thompson, Andrew Gallant, Richard Stallman, and of course Arthur) 

[I use Big data loosely, it usually relates to a problem being too large for the machine that you have access too. A great essay can be found in Programming pearls Vol 1 Cracking the Oyster. Where a big data problem is one that uses more than a megabyte of memory. ]

Essentially, my friend had a data set that contained roughly 100 million rows and each row had a column that was a free text description. A human could easily read several records and identify which records belonged together in one category, but the sheer number of records prevented any such attempt to be made. On the other hand, a computer using only a small set of keywords would easily misclassify many of the records. It would be easier to create more categories than there were and then combine categories based on statistical analysis. My friend tried fuzzy matching using something I covered in an earlier post, however, the algorithm slows down drastically with larger datasets.

By the time he told me his problem he had already invented an on the fly hashing and categorization algorithm that I have still yet to fully grok (his algorithm might be faster but it is optimized for his specific problem).

Problem Statement: Assuming you have a set of rows with descriptions and a set of keyword/tokens, assign each row to one or more tokens.

In particular every token will have a set of rows that belong to it. Rows can be assigned to multiple tokens and tokens have many rows.   The total number of tokens was 1 million. The sample dataset would have 10 million rows [limited by the memory on one machine].

The first approach was to use pandas and use groupby methods. However, that was disastrously slow because Pandas strings are actually Python objects instead of being Numpy objects and are simply not optimized for speed.

This led me to move to pure Numpy. I started by constructing a toy problem. The toy problem would have a million rows. Each row would consist of only 5 randomly selected lowercase English letters separated by spaces. The goal would be to identify all the rows that contained each letter.

Here is that experiment:

In [1]:

import pandas as pd
import random
import string
import numpy as np
In [2]:
" ".join(random.sample(string.ascii_lowercase, 5))
Out[2]:
'i k a r m'
In [3]:
#million
I=1000000
tokens=list(string.ascii_lowercase)
df=pd.DataFrame({"d":[" ".join(random.sample(tokens, 5)) for i in range(I)], "x":range(I)})
In [4]:
df.head()
Out[4]:
d x
0 b c v l y 0
1 m l q u t 1
2 b f s i n 2
3 g p e k h 3
4 j k a b l 4
In [5]:
#convert the description column into a fixed size array. 
#being a bit conservative on the max size of this column is okay
d=df.d.as_matrix().astype('S1000')
In [6]:
#d is now just fixed strings
d
Out[6]:
array([b'b c v l y', b'm l q u t', b'b f s i n', ..., b'p f k n h',
       b'y p t o i', b'e g o a i'], 
      dtype='|S1000')
In [7]:
# We see 26 arrays with the indexes of the rows that have each char
list(map(lambda x: np.where(np.char.find(d, np.bytes_(x))>-1),tokens))
Out[7]:
[(array([     4,     24,     39, ..., 999967, 999990, 999999]),),
 (array([     0,      2,      4, ..., 999990, 999992, 999993]),),
 (array([     0,      9,     11, ..., 999978, 999980, 999991]),),
 (array([    12,     15,     18, ..., 999979, 999983, 999985]),),
 (array([     3,     13,     14, ..., 999991, 999995, 999999]),),
 (array([     2,      7,     13, ..., 999980, 999989, 999997]),),
 (array([     3,      6,      7, ..., 999992, 999994, 999999]),),
 (array([     3,      6,     10, ..., 999987, 999993, 999997]),),
 (array([     2,      8,      9, ..., 999993, 999998, 999999]),),
 (array([     4,      5,     10, ..., 999990, 999993, 999994]),),
 (array([     3,      4,      5, ..., 999988, 999996, 999997]),),
 (array([     0,      1,      4, ..., 999984, 999987, 999996]),),
 (array([     1,     12,     15, ..., 999964, 999975, 999979]),),
 (array([     2,      9,     10, ..., 999985, 999989, 999997]),),
 (array([    16,     21,     25, ..., 999996, 999998, 999999]),),
 (array([     3,      6,     12, ..., 999992, 999997, 999998]),),
 (array([     1,     15,     20, ..., 999986, 999988, 999989]),),
 (array([     5,     13,     17, ..., 999974, 999984, 999989]),),
 (array([     2,     14,     18, ..., 999990, 999994, 999996]),),
 (array([     1,      7,     19, ..., 999985, 999994, 999998]),),
 (array([     1,      7,      8, ..., 999993, 999994, 999995]),),
 (array([     0,     10,     11, ..., 999974, 999983, 999988]),),
 (array([     6,      9,     18, ..., 999987, 999990, 999991]),),
 (array([    17,     25,     28, ..., 999988, 999992, 999995]),),
 (array([     0,      5,      7, ..., 999991, 999995, 999998]),),
 (array([     8,     17,     18, ..., 999976, 999995, 999996]),)]
In [8]:
%%timeit
list(map(lambda x: np.where(np.char.find(d, np.bytes_(x))>-1),tokens))
 1 loop, best of 3: 22.5 s per loop

The toy problem took 22.5 seconds on my core i7 desktop from 2015. All the benchmarks were done on the same machine to prevent skewing of results.

I decided to abandon Numpy. 22 seconds for 26 tokens would not scale. At this rate finding the matches for a million tokens would take  (1 million / 26) * 22 seconds = 9.79344729 days. Multiplying by another 10 because my problem was only on 1 million rows would mean this problem was going to run for 90 days at least, you can’t even return most clothes after that much time.

So I turned to my favorite language optimized for big data. KDBQ.

Here is basically the same code as Numpy with comments inline:

.Q.a is just all the lowercase ascii letters

/create a function g that will generate the random string
g:{raze ” “,’5?.Q.a}
/create a table with a million rows using g and column d
t:([]d:g@’1000000#0)
/demonstrate the function on one letter “a” for the first 10 rows.
/This creates a boolean mask where 1 means that the row contains “a”
{first `boolean$ss[x;”a”]} each 10#t.d
0000000010b
/use where to save only the rows that have the letter
where {first `boolean$ss[x;”a”]} each 10#t.d
,9
/combine this function and apply it for each letter returning 26 rows with the indexes
/for the first 10 rows
{where {first `boolean$ss[y;x]}[x] each 10#t.d} each .Q.a
,9
2 4
1 3
`long$()
1 5 9
3 7 8
`long$()
,4
5 8 9
,7
1 7 9
,1
`long$()
5 8
,0
6 9
4 5 6
0 2 3
,2
`long$()
`long$()
0 2
..
/Time this function for a million rows
q)\t {where {first `boolean$ss[y;x]}[x] each t.d} each .Q.a
5335
/Try to parallelize the problem with 6 slave threads using peach
q)\t {where {first `boolean$ss[y;x]}[x] each t.d} peach .Q.a
3104
/Try to further parallelize the problem, no luck
q)\t {where {first `boolean$ss[y;x]}[x] peach t.d} peach .Q.a
3087

KDB has done a lot better than Numpy on the same problem, BUT this problem will still take: 1.33547009 days * 10 because this was only a million rows. We get 10 days, which is more reasonable. But we can do better.

I decided to take advantage of the KDBQ enumerated string AKA symbol type.

This led to the following code and timings:

q)g:{raze `$’5?.Q.a}
q)t:([]d:g@’1000000#0)
q)first t[`d]
`f`j`r`p`j
q)tokens:`$’.Q.a
q)\t {where {x in y}[x] each t.d} each tokens
3796
q)\t {where {x in y}[x] each t.d} peach tokens
2218
q)\t {where {x in y}[x] peach t.d} peach tokens
2288

Definitely an improvement, but nothing to write home about.

At this point, it was time to start looking to improve the algorithm. A quick StackOverflow search lead me to the name of the problem. MultiString/MultiPattern Search. There are two dominant algorithms: Aho-Corasick and Rabin-Karp. Aho-Corasick was developed by Aho of AWK fame and Margaret Corasick (not much is written on the internet about Margaret Corasick, except that she was clearly one of the leading female pioneers in computer science). The Algorithm is linear in the number of inputs and the number of tokens and the number of matches. The Rabin-Karp algorithm, in true Rabin style, is technically not linear but average case linear assuming a good hashing function is chosen and it involves using prime numbers so that you can roll through the string without looping back, more here.

I was going to start implementing one of these algorithms, but as I was searching for a good implementation, I found that the reference implementation was to be found in GNU Grep. So out of curiosity, I tried running my toy problem using grep.

I used KDBQ to write out my table to t.txt and deleted the first line so that I just had the rows of strings on each new line.  I then saved the another file called tokens.txt where each new line was an ascii letter. I then constructed the necessary grep to get the line numbers where a token appears.

grep -n token filename

999918: x b w a b
999923: a v u p u
999924: y v i a p
999927: c k j a g
999930: a l n i p
999935: p a b e k
999944: f a a n t
999949: b u a h p
999951: n j x u a
….

This will give you the line numbers and the matches using cut command you can return only the line numbers.

Here is what that looks like:

$ grep -n “a” t.txt | cut -d: -f1

999918 
999923
999924
….

I then decided to go read the tokens.txt file line by line and grep for the line in the table of a million rows. Each result of line numbers is written to a file in the results folder with that token name. So if we want all the rows that match a token we simply look at the row numbers in that file.

time while read LINE; do grep -n $LINE t.txt| cut -d: -f1 > results/$LINE; done <tokens.txt

real 0m0.749s
user 0m0.852s
sys 0m0.044s

Less than a second and we haven’t even parallelized the operation. Adding & so that things run in the background and waiting for the process to finish we see the toy problem for what it is, A TOY.

$ time $(while read LINE; do grep -n $LINE t.txt| cut -d: -f1 > results/$LINE & done <tokens.txt; wait)

real 0m0.216s
user 0m1.424s
sys 0m0.112s

With parallelization the toy problem takes less than a fifth of a second. Estimating the run time for the original problem we see that it will now take: (1 million/26)*.216 seconds = 2.3 hours * 10 this is now in the less than a day category.

Curious and surprised that grep was reading off the disk faster than KDBQ was reading from memory, I decided to find out why grep was so FAST.

Summarizing one of the original authors of GNU grep:

  • Like KDBQ, GREP memory maps files
  • Like KDBQ grep will often use the Boyer-Moore search algorithm that basically allows quick skipping of bytes where no match is possible.
  • Grep avoids looking at the lines and only reconstructs line breaks if there is a match

Mike Hartel’s famous quote from that article is very enlightening.

The key to making programs fast is to make them do practically nothing. ;-)

So essentially KDBQ is slower than grep because it still thinks of strings in each row as independent rows and it isn’t designed solely for ripping through strings. The tool was originally the brainchild of Ken Thompson of K&R C fame. So is it surprising that after 40 years the tool that was optimized for string searching is still the king? Not really, but we thought Arthur could do better, we eagerly await K6 which promises to be an entire OS, which will probably include something like grep. (I later realized that I could use Q in a different way to achieve the same result much faster, keep reading… )

A quick aside, for those unfamiliar with grep. Grep is a Unix command utility, that allows searching files for strings that match on a line. Many options were added as Grep developed. Grep gave you a way to search for regular expressions. That is, for things that are similar to a certain string. The extended regular expression engine was originally Aho’s invention and was called EGrep. This merged into Grep as the argument -e.  FGrep gave you a way to search for exact matches much faster which became -f. Similarly, -n allows you to print line numbers and -o allows you to print only the matching pattern.

Anyway, knowing something about the algorithms that allow grep to run so quickly, I decided that it was time to try a more realistic example. Firstly, I thought that longer words might further advantage grep since it could skip more quickly through the file, secondly, I wondered how linux would handle opening many background processes at once.

So to construct the more realistic dataset I went and found a file that contained the 20,000 most common english words.

I read the files into KDBQ and created my tokens and my table:

q)tokens:read0 `:20k.txt
/now instead always taking 5 random tokens we are randomly taking between 1 and 10 /random words without replacement, so that we get a unique set for each row.
q)g:{” ” sv (neg first 1+1?10)?tokens}
q)t:([]d:g@’1000000#0)
q)\t {where {first `boolean$ss[y;x]}[x] each t.d} peach 10#tokens
2453
q)\t {where {first `boolean$ss[y;x]}[x] peach t.d} peach 10#tokens
2498
q)\t {where {first `boolean$ss[y;x]}[x] each t.d} peach 100#tokens
12188
q)\t {where {first `boolean$ss[y;x]}[x] peach t.d} peach 100#tokens
12329
q)\t {where {first `boolean$ss[y;x]}[x] peach t.d} peach 1000#tokens
184171

What we see is that getting the list of indexes for 10 tokens takes about 2.5 seconds, with no significant difference from adding another peach. Getting the indexes for 100 tokens takes 12 seconds, which looks like an improvement, to about 1.2 seconds per 10 tokens. However, getting the indexes for 1000 tokens takes about 3 minutes for an average of about 2 seconds per 10 tokens. So KDBQ is definitely performing better on longer tokens, and tokens get longer deeper into the list so that is contributing to the speedup. However, we really are just using KDBQ to place the files on the filesystem for grep to pick them up.

q)save`:t.txt
/save only the top 100 tokens in a separate file
q)tokens:([]100#tokens)
q)save `:tokens.txt

Now we can see how grep does on this much more realistic problem.

Running grep on 100 tokens; we get the following timing result:

time $(while read LINE; do grep -n ” $LINE ” t.txt| cut -d: -f1 > results/$LINE & done <tokens.txt; wait)

real 0m4.460s
user 0m6.684s
sys 0m0.884s

So for 100 words it took 4.46 seconds, so less than 5 seconds. I decided to run it for the whole set of 20,000 words:

Here are the timings experimenting with using Fgrep and using and not using background processes:

It took 5 minutes to run on 20k tokens for a million line table.

time $(while read LINE; do grep -n ” $LINE ” t.txt| cut -d: -f1 > results/$LINE & done <20k.txt; wait)

real 4m54.079s
user 17m13.756s
sys 2m9.064s

Not running in parallel is a bit faster in terms of total CPU time, but is almost 6 times slower in wall clock time, which means you will be waiting awhile.
time while read LINE; do grep -n ” $LINE ” t.txt| cut -d: -f1 > results/$LINE; done <20k.txt

real 23m41.263s
user 14m14.508s
sys 1m57.508s

Using Fgrep was marginally faster because it searches a literal string without any pattern matching support, if you need pattern matching then you can’t use this, but the performance improvement was not all that significant. (This should have been a clue that something was wrong, see later in this post)

time $(while read LINE; do fgrep -n ” $LINE ” t.txt| cut -d: -f1 > results/$LINE & done <20k.txt; wait)

real 4m29.731s
user 17m13.004s
sys 2m2.336s

Now multiplying through for the more realistic example we find that (1 million/20k)*5 minutes=4.1 hours, multiplying by 10 we get 40 hours. So it seems that our toy problem was underestimating the difficulty of the problem.

40 hours seemed like too much time to wait. So off I went in search of something better. Luckily, I came across ripgrep. Ripgrep aims be to be a better version of Grep written in Rust instead of C. [Rust is a language that aims to be a better version of C, so this is fitting. The name rip grep, I’m sure has something to do with R.I.P grep]. The post  on Ripgrep taught me a lot about the underlying implementations of different search tools. It made me realize that some of the information I was reading from older posts was outdated. In particular, grep no longer used memory maps (–mmap), even if you asked for them explicitly. So grep was fast, but not as fast as it could be. Also, the newest grep had somehow increased it’s speed on non-ASCII characters. So some of the traditional speedups like setting up LC_ALL=C, which forces grep to use ASCII characters doesn’t actually increase the speed. Which many of these StackOverFlow answers mention here, here, and here.

Here are some timings that explicitly contradict these results:

Searching for aardvark in a 10 million row file results in no speed improvement when switching between ASCII and UTF-8

$ time LC_ALL=en_US.UTF-8 grep -F aardvark -w -o -n t.txt > results/aardvark

real 0m0.339s
user 0m0.260s
sys 0m0.076s
$ time LC_ALL=C grep -F aardvark -w -o -n t.txt > results/aardvark

real 0m0.340s
user 0m0.276s
sys 0m0.064s

Also, many posts recommended using GNU parallel, but I found it did a poor job of actually distributing work onto all the cores of my CPU and that it was probably much better if you were going to send work to multiple machines. Instead simply using background processes was going to be fine.

Once I read the article about ripgrep, it was time to test its bold claim that it was going to be faster than all the other tools. My more realistic example with 20k words became really easy for ripgrep, once I figured out that the 20k tokens can be split into  smaller 200 line files with split.  Split is a tool written by Richard Stallman that splits the original files into files with the prefix specified. Here I use ‘S’.

$split -l200 20k.txt S

The small-token-files could be used to match against the large input and output a line number and a match on each line.

For example, in our toy problem, if we want to match on letters a-z. We can break it up into chunks of 2 letters. We will get 13 files.

file 1.–> corresponds to letters a,b
1:a
3:b
4:a
4:b

file 3. –> corresponds to letter e,f
1:e
15:f
74:e
80:e

So that for each token-file we would get a corresponding file with matches. Here is that command and timing:

$ time $(for i in S*; do rg -F -f $i -n -o -w t.txt>results/groupby$i & done; wait)

real 0m36.187s
user 3m23.104s
sys 1m15.308s

Now my more realistic problem seemed like a joke. It was taking less than a minute to go through a 10 million row file with 20 thousand tokens. I was starting to become very impressed with ripgrep. However, I decided that I first needed to try and find the older versions of grep that supported mmap to make the comparison fair.

After building from source no less than 8 versions of grep, I discovered some very curious things.

  • grep lost mmap support as of 2010, so any version after 2.6 essentially was not going to work, and my current ubuntu only came with 2.25, 25 is larger than 6.
  • Old versions of grep were in general faster on large input patterns than new versions, they built the Aho-Corasick tree faster.
  • Grep versions before 2.5 did not support printing only the match (-o) option, which was necessary for us to collect all the indexes back together at the end.
  • So Grep versions 2.5.x were the only ones to support both -o and –mmap options

I ended up settling on grep 2.5.1 for a very strange reason.  I was excited that grep 2.5.1 was going really fast on some test inputs, when I discovered that when matching on multiple patterns at once and printing the line number. I would get the following output:

1:tenor
bleeding
fertilizer
colchester
2:mining
swap
everyone
miriam
persuasion
….

It seemed that grep had decided that it would be easier to just print the line number once per many matches. So I went to the next version grep 2.5.4, which printed the matches as expected.

1:tenor
1:bleeding
1:fertilizer
1:colchester
2:mining
2:swap
2:everyone
2:miriam
2:persuasion
….

However,  it was an order of magnitude slower. So I ended up attaching a small AWK program that would add the line numbers for the lines that had them implicitly. Here is that AWK program:

awk -F: ‘NF==2{line=$1} NF==1{$0=line FS $0}1’

This AWK program had almost no overhead when I tested it so it was much better than relying on grep to print all the line numbers for each match.

Once I had the right version of grep, things really started moving. I discovered that grep could match many more tokens at once than ripgrep.

In fact, my more realistic example with 20 thousand tokens ran in the same amount of time whether or not I split the token files. In fact, the overhead of starting many processes was slowing down grep2.5.1.

MOST REALISTIC PROBLEM, (yet):

This led me to go and create the most realistic example I could without generating nonsense words myself. (I might do this later, but I was feeling a bit lazy and I think this problem is close enough). I went here and got 355k words. I went and recreated my dataset. I had a file of tokens.txt which had 354950 words to be exact. I deleted about 40 words that started with a ‘(single quote) because ripgrep has a bug (bug was fixed within 14 hours of posting it, amazing!) that prevents it from properly printing out only matches that match exactly that word[-o -w]( a very low priority bug , but I figured it was worth filing and skipping those words). I created my 10 million row dataset which had a random number of words between 1 and 10 from the tokens file. I saved that into t.txt.

I then created two very similar scripts that allowed ripgrep and grep2.5.1 to finally compete.

grepper.sh:

#!/bin/bash

FILE=$1

export LC_ALL=C
pid_count=0
for F in S*
do
~/Downloads/grep-2.5.1/src/grep –mmap -F -f $F -n -w -o $FILE|awk -F: ‘NF==2{line=$1} NF==1{$0=line FS $0}1’ >results/group$F &
if (( ++pid_count > 4 )); then
wait -n
((pid_count–))
fi
done
wait

rger.sh

#!/bin/bash
FILE=$1

pid_count=0
for F in S*
do
/home/pinhus/Downloads/rg –mmap -F -f $F -n -w -o $FILE >results/group$F &
if (( ++pid_count > 4 )); then
wait -n
((pid_count–))
fi
done
wait

The two programs behave very similarly and produce identical output. So I will describe them together.

The program takes a filename argument $1 and creates a count of how many processes are running so far. For each pattern File S* it runs a literal match on each patternfile and prints the line numbers along with each match into a corresponding file inside the results folder called group{patternfile}. It runs this in the background and increments the number of running processes. It waits until there are less than 4 running processes, decrements the number of processes and adds another process for the next patternfile. When it is all done it will wait until all the pattern files are written out.

The main difference between them is that the grep version uses strictly ASCII characters whereas the ripgrep version uses UTF8 by default.

There is also a difference in how to optimally split the pattern file. After much experimentation, I found that ripgrep could comfortably handle up to about 400 lines at a time, whereas grep2.5.1 was optimal with 20k lines at a time.

Here are the respective split commands:

#grep split, assuming 355k lines token file
# -n will create N chunks and not break a lines if you add l/N
split -n l/20 tokens.txt S

#ripgrep split assuming 355k lines in token file
split -n l/888 tokens.txt S

As you can see there are only 20 chunks created for the grep version and there are 888 files for the ripgrep version.

However, the real kicker is that the grep version now runs faster, because we have taken care of the overhead of having too many process handles at once.

Here are the timings:

time bash grepper.sh t.txt

real 1m27.812s
user 2m55.840s
sys 0m2.100s

time bash rger.sh t.txt

real 3m29.019s
user 25m19.516s
sys 1m55.148s

Long Live grep!

Multiplying to the original problem of 1 million tokens we now see that grep will take about 4.5 minutes. Let’s say 5 minutes to be safe. And ripgrep will now take about 10.5 minutes.

The reason that grep now beats ripgrep is that grep uses Aho-Corasick to do multiple matching, whereas ripgrep actually uses a clever version of regex. Instead of implementing Aho-Corasick, ripgrep joins the entire pattern file with | which makes the regex engine match any of the those patterns. Ripgrep then goes one step further and tries to optimize this by searching for a common literal among all those patterns so that it can skip using the regex engine which is much slower. If it finds a common literal string among the patterns it will submit it to it’s implementation of the teddy algorithm, which takes advantage of Intel’s more recent SIMD instructions (single instruction multiple data, AKA vectorized processing). The teddy algorithm submits the literal to be checked against 16 bytes at a time and if it finds a match it will go into the slower regex engine to verify. This implementation forces our pattern files to be smaller since a pattern-file that is too big will not have a common literal or will match too often. In fact, the pattern-files have a limit that will throw an error if they are too big because they overwhelm the regex engine. Because of this limitation, ripgrep has to search the entire file many more times than grep has to. Additionally, ripgrep will perform terribly if the original patternfile is not sorted. I suspected that ripgrep was benefitting from our sorted tokens file, and was able to find common literals, which was allowing it to be much faster than a random tokens file. To test this I used round robin split,

#ripgrep split assuming 355k lines in token file
# replace the l with r. which causes the lines to be round robined across 888 files.
split -n r/888 tokens.txt S

The difference was staggering, using a 100 thousand line input file:

#pattern file is sorted and broken into sorted chunks
time bash rger.sh t100k.txt

real 0m8.205s
user 0m59.716s
sys 0m1.400s

#pattern file is deliberately round robined across chunks
time bash rger.sh t100k.txt

real 1m11.450s
user 9m8.252s
sys 0m2.092s

The difference is between 8 seconds for a 100 thousand row file and 1.2 minutes. Almost an order of magnitude difference. Sort your pattern file before using ripgrep. Also because ripgrep uses a nonnaive version of Boyer-Moore it is able to skip many more bytes. (We eagerly await this optimization from Arthur, after all SIMDs are a natural part of vectorized languages.)

On the other hand, because grep uses Aho-Corasick it can handle much larger pattern files, the chunks merely facilitate creation of a smaller tree to traverse and you don’t need to sort your pattern files. This also allows grep to read the input file many fewer times. Once again, algorithms beat brute force  (kind of, see last note where I discover that many matches were missed because of a bug in 2.5.1, most of the points stand but I think you might be better off using ripgrep).

Now to collect the outputs and create the final result, which is a token and all of it’s associated row numbers, we go back to KDBQ.

All of our output files live in the results directory and look something like this:

402:bicultural
630:bicapitate
697:bibby
861:biblicist
945:bicylindrical
951:bibliopolist
1200:bicalcarate

So we need a function that gets all the files from a directory:

tree:{ $[ x ~ k:key x; x; 11h = type k; raze (.z.s ` sv x,) each k;()]}

This is one of the utilities provided in the KDB Cookbooks/tips. It will recursively go down a directory and give you all the files. We just need 1 level, but it’s fast enough, so there is no need to write our own.

We create a table that will have two columns, an index column and a token column.

t:([]index:();token:())

We then create a small function that will read in each file split on the colon and parse the line and upsert it into the table.

{`t upsert flip 0:[(“IS”;”:”);x]}

0: will read a text file and takes arguments that specify what it expects see, as well as a delimiter. In this case a colon.
flip arranges the tuple to be aligned with the table
`t upsert will insert tuple.

By default 0: will read the two columns into rows, so we need to flip it.
q)(“IS”;”:”) 0: `:results/groupSaaa
72 165 567 698 811 838 1003 1136 12..
abiliment abience a1 4th aberuncator abbey abductors abdominohysterectomy 10..
q)flip (“IS”;”:”) 0: `:results/groupSaaa
72i `abiliment
165i `abience
567i `a1
698i `4th
811i `aberuncator
838i `abbey

Running this on each file in the results directory like so:

q)\t {`t upsert flip 0:[(“IS”;”:”);x]} each tree `:results
4361
q)\t select index by token from t
567
q)select index by token from t
token| index ..
—–| ———————————————————————-..
1080 | 1283 127164 391580 443156 543775 561432 608795 862810 889780 1029079 1..
10th | 27747 57768 238851 293795 332990 487295 564854 732020 755630 834528 83..
1st | 37915 113165 175561 203936 213171 318590 382626 411531 473224 653353 6..
2 | 26642 67272 108161 122598 192496 211019 250924 281576 300950 344445 45..
2nd | 3275 50750 93226 99089 107105 208955 212762 234366 238134 403142 48216

Reading in all the output files and aggregating by token takes less than 5 seconds.

However, once I did this aggregation, I thought there must be a better way to take advantage of a KDB only solution.

PURE KDBQ:

If you knew that in the description, the only tokens you were interested were words and that the only types of word breaks were spaces then you could do following:

/read in the raw 10 million row file
q)read0 `:t.txt
“cellulipetal outname kymograms splendiferousness preneural raploch noncontin..
“copernicus”
“overglanced”
“consists”
“epicure casern mythologization”
“angiosymphysis”
“shellfishes renniogen lactonic”
..
/define raw as this file
q)raw:read0 `:t.txt
/split the first row on spaces
q)” ” vs first raw
“cellulipetal”
“outname”
“kymograms”
“splendiferousness”
“preneural”
“raploch”
“noncontingency”
..
/cast each word to a symbol
q)`$” ” vs first raw
`cellulipetal`outname`kymograms`splendiferousness`preneural`raploch…
/apply this function to each row in raw and call it sraw
q)sraw:{`$” ” vs x} each raw
/this takes about 19 seconds
q)\t sraw:{`$” ” vs x} each raw
18744
/sraw looks like this
q)sraw
`cellulipetal`outname`kymograms`splendiferousness`preneural`raploch..
,`copernicus
,`overglanced
,`consists
`epicure`casern`mythologization
,`angiosymphysis
`shellfishes`renniogen`lactonic
..
/create an index column that tracks each row (add 1 so that it will agree with grep)
q)index:1+til count sraw
q)index
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29..
/join together each row number with every symbol in that row
q)flatrows:raze index {x,’raze y}’ sraw
q)\t flatrows:raze index {x,’raze y}’ sraw
7634
/flatrows looks like this:
q)flatrows
1 `cellulipetal
1 `outname
1 `kymograms
1 `splendiferousness
1 `preneural
1 `raploch
1 `noncontingency
..
/create a new table that has a rid and a token
q)t2:([]index:flatrows[;0]; token:flatrows[;1])
q)\t t2:([]index:flatrows[;0]; token:flatrows[;1])
4631
/t2 looks like this:
q)t2
index token
———————–
1 cellulipetal
1 outname
1 kymograms
1 splendiferousness
1 preneural
..
/aggregate and get all indexes for each token
/this takes less than 8 seconds
q)\t res2:select index by token from t2
7741
/res2 looks like
q)res2
token| index ..
—–| ———————————————————————-..
1080 | 1283 127164 391580 443156 543775 561432 608795 862810 889780 1029079 1..
10th | 27747 57768 238851 293795 332990 487295 564854 732020 755630 834528 83..
1st | 37915 113165 175561 203936 213171 318590 382626 411531 473224 653353 6..
2 | 26642 67272 108161 122598 192496 211019 250924 281576 300950 344445 45..
2nd | 3275 50750 93226 99089 107105 208955 212762 234366 238134 403142 48216
/total time without ever leaving kdb: 19+8+5+8=40 seconds!!

/If you are not interested in all the tokens you can easily subselect
/as an example asked for a 1000 random subsample of tokens from res2
q)sampleTokens:neg[1000]?(0!res2)[`token]
/sampleTokens looks like this
q)sampleTokens
`bimillionaire`abarthrosis`krone`volleyers`correct`pitarah`aeroelastic`zilch`..
/I can ask for res2 where the token is in sampleTokens
q)select from res2 where token in sampleTokens
token | index ..
————| —————————————————————..
abaciscus | 19644 81999 125546 145360 159922 213205 223650 227560 353332 41..
abactinally | 92704 265473 480112 499419 510784 539605 617031 668847 752500 8..
abarthrosis | 25638 61103 215367 281385 305352 385693 553059 585997 598293 66..
abir | 96166 214361 258595 259398 309117 387755 419601 720156 746810 8..
abstemiously| 4611 27043 61983 72530 140138 152631 227466 272363 278554 31059..
..
/this takes a millisecond for 1000 tokens
q)\t select from res2 where token in sampleTokens
1

So if you know that everything is in ASCII and you know that the words in the description are space delimited nothing beats KDBQ. Arthur you win. I did think to tokenize my inputs, but I didn’t think to wrap them with an index number, until I was collecting the results.

A disclaimer where I suggest using ripgrep until grep reintroduces mmap for large files. Old versions of grep have bugs that are not documented all that well.

On the other hand if you need UTF-8 use ripgrep and if you need regular expressions use ripgrep. The recent lack of support for –mmap in grep means that the chance of doing better with it requires you to build it yourself and deal with bugs.

(It turns out that grep2.5.1 has a bug where -w will not match word boundaries correctly and 2.5.4 fixes that bug but introduces a performance penalty by printing every line, for every match – a task much better suited to AWK. To work around this and still get the same performance, you either need to add spaces to your matches or turn off -w and deal with the overlap matches that will inevitably occur.) I ended up using the command below that mostly did away with the performance penalty, but was still slower than ripgrep.

On the other hand, the bug I found in ripgrep while writing this post was fixed within 14 hours of me reporting it. A serious benefit of a tool that is actively supported.

time cat tokens.txt | parallel –pipe –block 1000k –round-robin grep2.5.1 -Ff – –mmap -on t.txt |awk -F: ‘NF==2{line=$1} NF==1{$0=line FS $0}1’ >results/test

real 4m41.146s
user 15m35.820s
sys 0m8.096s

I didn’t get a chance to fully optimize the block size, but it’s within reach of ripgrep.

GNU Parallel was used:

O. Tange (2011): GNU Parallel – The Command-Line Power Tool,
;login: The USENIX Magazine, February 2011:42-47.