What I (wish I) did at Geneity

This page is documents some of the really cool things that I wish I had the time to pursue during my day-to-day, but instead did during my spare time. Some of them went on to be implemented, others I still have hope for, and then there’s the ideas that I’ve since come to consider sub-optimal. The overarching theme is a desire to replace code with data and algorithms on that data.

I’d read the introduction first and then feel free to jump around, most of the sections should stand alone. Personally I recommend market DSL, my first and most academic idea.

Teach Yourself Algorithms in Five Minutes

If I was to describe the purpose of the algorithms we develop in a single word it would be “derivatives” (with “related contingencies” in a distant second). Okay, let’s back up a bit and talk about sports betting in general, here’s a crudely drawn representation of a very small offering for a soccer match:

Soccer match.

Selections are an outcome that can be bet on, for example “the home team will win”. They have a probability of occurring which is overrounded to produce a price which represents how much the bet could win. Markets are groups of related selections, if exactly one selection can win and the probability sums to one I call it a proper market. As developers we’re interested in probabilities and proper markets because prices and markets are just transformations with fewer invariants.

Another way of grouping markets is driving and derivative, where driving markets are the markets that contain enough information to calculate the derivations. The set of driving markets varies per sport and model but it’s typically a supremacy market such as “who will score the most points?” and an expectancy market such as “how many points will be scored in total?”. Algorithms take the burden of ensuring that all prices are consistent off the traders — the domain experts who set the prices for an event — by reducing their workload to just pricing the driving markets, allowing them to trade more events. Even this work can be outsourced by receiving the driving markets from an odds feed.

We call the process of fitting the model to the prices of the driving markets parameter derivation:

Parameter derivation.

Once we have parameters we can use the to build the model. The model provides data structures from which we can extract probabilities of outcomes, a process we call calculation:

Calculation.

It’s the relationship between structures and extractors that give rise to proper markets, structures contain all the possibilities at a moment in time, so an extractor that assigns each probability within it to a single selection is guaranteed to calculate a proper market.

Efficient Extractors

Some say that Python is a terrible choice of language for high-performance numerical computation. You’ll get no argument from me (or anyone on my team), but we’re stuck with it for now, so I try to push the calculations into NumPy. Sean Parent’s advice of “no raw loops” is something to live by, and this time it manifests as a rare example of optimization:

>>> import numpy

>>> # cs is a correct score grid
>>> # cs[x, y] is the probability that the match finishes x-y.
>>> def who_will_win(cs):
...     ps = [0, 0, 0] # home, draw, away.
...     for (x, y), p in numpy.ndenumerate(cs):
...         if x > y:   ps[0] += p
...         elif y > x: ps[2] += p
...         else:       ps[1] += p
...     return ps

>>> # `.T` transposes the array so we can write home on the x-axis.
>>> cs = numpy.array([
...     [0.10, 0.15, 0.05, 0.05],
...     [0.10, 0.15, 0.10, 0.05],
...     [0.05, 0.05, 0.05, 0.05],
...     [0.00, 0.00, 0.05, 0.00],
... ]).T

>>> who_will_win(cs)
[0.45, 0.30, 0.25]

>>> import timeit
>>> timeit.timeit('who_will_win(cs)',
...               'from __main__ import cs, who_will_win')
11.737718105316162
“Who will win?” for-based extractor.

“Who will win?" is a proper market, that means we can calculate it assigning each cell a bin and summing the cells in each bin, aka numpy.bincount.

>>> # TODO: Memoize, and use in-place operations.
>>> def who_will_win_bins(sz):
...     bins = numpy.indices(sz)
...     bins[0] *= -1
...     bins = bins.sum(axis=0)
...     return numpy.sign(bins) + 1

>>> who_will_win_bins(cs.shape)
array([[1, 2, 2, 2],
       [0, 1, 2, 2],
       [0, 0, 1, 2],
       [0, 0, 0, 1]])

>>> bins = who_will_win_bins(cs.shape)
>>> numpy.bincount(bins.flat, weights=cs.flat, minlength=3)
array([ 0.45,  0.3 ,  0.25])

>>> timeit.timeit('numpy.bincount(bins.flat, weights=cs.flat, minlength=3)',
...               'from __main__ import numpy, bins, cs')
6.1542229652404785
“Who will win?” bincount-based extractor.

Well this is a little faster, but my experiments in C led me to expect a 100× speed-up. Let’s try with a more realistic 21×21 grid:

>>> cs = numpy.random.uniform(0, 1, (21, 21))
>>> cs /= cs.sum()
>>> bins = who_will_win_bins(cs.shape)

>>> timeit.timeit('who_will_win(cs)',
...               'from __main__ import cs, who_will_win')
357.38719296455383

>>> timeit.timeit('numpy.bincount(bins.flat, weights=cs.flat, minlength=3)',
...               'from __main__ import numpy, bins, cs')
7.571918964385986

>>> who_will_win(cs) == numpy.bincount(bins.flat, weights=cs.flat, minlength=3)
array([ True,  True,  True], dtype=bool)
Extracting from a realistic grid.

For a few minutes I was worried that the for-based one wasn’t going to terminate, but we got there in the end! Of course if we had been able to write in C (or C++, Fortran, etc) we could get even more speed:

#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

struct Extent {
	size_t w, h;
};

void extract(struct Extent e, uint8_t *bs, float *restrict ps, float *restrict out)
{
	uint8_t *be = bs + e.w * e.h;
	while (bs != be) out[*bs++] += *ps++;
}

int main(void)
{
	struct Extent e = { 21, 21 };
	uint8_t www[e.w * e.h];
	uint8_t *bs = www;
	for (size_t j = 0; j < e.h; ++j) {
		memset(bs, 2, j);
		bs[j] = 1;
		memset(bs + j + 1, 0, e.w - j);
		bs += e.w;
	}

	float cs[e.w * e.h];
	float sum = 0;
	for (size_t n = 0; n < e.w * e.h; ++n) {
		sum += cs[n] = ((float)rand()) / RAND_MAX;
	}
	for (size_t n = 0; n < e.w * e.h; ++n) cs[n] /= sum;

	float out[3];

	struct timespec begin, end;
	clock_gettime(CLOCK_REALTIME, &begin);

	for (size_t i = 0; i < 1000000; ++i) {
		memset(out, 0, sizeof out);
		extract(e, www, cs, out);
	}

	clock_gettime(CLOCK_REALTIME, &end);
	long s = end.tv_sec - begin.tv_sec;
	long ns = end.tv_nsec - begin.tv_nsec;
	if (ns < 0) { s -= 1; ns += 1000000000; }
	printf("%ld.%09ld\n", s, ns);
}

// gcc -O3 -o extract extract.c
// ./extract
// 0.830870301
“Who will win?” C implementation.

Calculating a million markets in a second is overkill as even our biggest events only have two hundred or so. I wonder what we could put all that extra CPU to use for, maybe something involving an arbitrary number of markets and selections…

Market DSL

This was my first major breakthrough and forged my belief that (at least) my whole team could stop writing code if we could just find the right abstractions. Now I am more experienced this idea is no longer what I would propose, but it’s close. The original prototype is on GitHub.

Resulting is the process that decides which selections have won, which have lost, and which are not yet decided. You could use the algorithms to do this, taking zero as loss, one as win, and everything else as undecided. Unfortunately there are markets that we offer that the algorithms don’t model, and so instead we have two implementations that both understand what each market means. This is obviously wasteful, but worse it also leads to subtle errors when the definitions are slightly different. This led to me design a language that you could express the meaning of a market in:

>>> # "How many goals are there?"
>>> # len(collection) -> len(collection)
>>> total_goals = len(goals)

>>> # "Which team is associated with the first goal?"
>>> # nth(collection, n) -> collection[n]
>>> # attr(element, k) -> getattr(element, k)
>>> first_goal = attr(nth(goals, 0), 'team')

>>> # "Which team has the maximum total goals?"
>>> # partition(collection, k) -> groupby(attrgetter(k), collection)
>>> # map(mapping, fn) -> {k: fn(v) for k, v in mapping.items()}
>>> # max(mapping) -> {k for k, v in mapping.items() if v == max(mapping.values())}
>>> most_goals = max(map(partition(goals, 'team'), len))
Defining markets.

In this language resulting is eval and “which selections are in this market?” is members·typeof:

>>> typeof(total_goals)
>>> typeof(len)
int
>>> # soccer (incorrectly) assumes that int means total goals.
>>> # 4+ means "four or more", soccer games don't tend to go that high.
>>> members(int, soccer)
{'0', '1', '2', '3', '4+'}

>>> typeof(first_goal)
>>> typeof(goals.team)
team
>>> members(team, soccer)
{'home', 'away'}

>>> typeof(most_goals)
>>> typeof(partition(goals, 'team'))
team -> [goal]
>>> typeof(map(team -> [goal], len))
team -> int
>>> typeof(max(team -> int))
{team}
>>> members({team}, soccer)
{('home',), ('away',), ('home', 'away')}

>>> prematch = State(goals=[])
>>> # goals is actually a list of who scored and when, etc.
>>> inplay2_1 = State(goals=['home', 'home', 'away'])

>>> eval(total_goals, prematch)
0

>>> eval(first_goal, prematch)
>>> eval(nth(goals, 0), prematch)
Unknown
>>> eval(attr(Unknown, 'team'))
Unknown

>>> eval(most_goals, inplay2_1)
>>> eval(partition(goals, 'team'), inplay2_1)
{'home': ['home', 'home'], 'away': ['away']}
>>> eval(map({'home': ['home', 'home'], 'away': ['away']}, len))
{'home': 2, 'away': 1}
>>> eval(max({'home': 2, 'away': 1}))
{'home'}
typeof, members and eval.

Another particularly interesting walker estimates the relative magnitudes of the probabilities of selections in a given state, and that can be used to estimate the direction each probability will move given a state change:

>>> import functools32
>>> import itertools
>>> import operator

>>> @functools32.singledispatch
... def magnitude(expr, env):
...     raise NotImplementedError

>>> @magnitude.register(types.len)
... def magnitude_len(expr, env):
...     current = eval(expr, env)
...
...     magnitudes = {}
...     for target in members(typeof(expr), env):
...         if isinstance(target, str) and target.endswith('+'):
...             magnitudes[target] = int(target[:-1]) - current
...         elif target < current:
...             magnitudes[target] = float('inf')
...         else:
...             magnitudes[target] = target - current
...
...     return magnitudes

>>> # WARNING: Assumes scores increment by one.
>>> @magnitude.register(types.max)
... def magnitude_max(expr, env):
...     current = eval(expr.coll, env)
...     kmax, vmax = max(current.iteritems(), key=operator.itemgetter(1))
...     uniquemax = sum(1 for v in current.itervalues() if v == vmax) == 1
...
...     magnitudes = {}
...     for ks in members(typeof(expr), env):
...         # WARNING: Doesn't work if kmax=A, ks={B, C}.
...         if len(ks) > 1 or (ks == (kmax,) and uniquemax):
...             target = vmax
...         else:
...             target = vmax + 1
...
...         magnitudes[ks] = sum(target - current[k] for k in ks)
...
...     return magnitudes


>>> def trend(expr, env0, env1):
...     magnitude0 = magnitude(expr, env0)
...     magnitude1 = magnitude(expr, env1)
...     return {k: cmp(magnitude0[k], magnitude1[k]) for k in magnitude0}

>>> inplay1_0 = State(goals=['home'])

>>> trend(total_goals, prematch, inplay1_0)
{0: -1, 1: 1, 2: 1, 3: 1, '4+': 1}

>>> trend(most_goals, prematch, inplay1_0)
{('home',): 1, ('away',): -1, ('home', 'away'): -1}
Estimating relative probabilities (magnitude) and directions (trend).

We can use this as an unsophisticated way to verify that our algorithms behave sanely for all states. This function is also particularly applicable for solving a problem we face with our odds feeds. We receive odds from several sources and aggregate them before fitting the algorithm to them, the problem is that not all feeds update at the same time so we need to discard prices that predate the latest change in state but many of the feeds don’t provide the state the prices apply to. We can compare the directions of the deltas of the prices to the predicted directions as a quick and dirty way to detect this situation. It is also somewhat possible to infer a change of state from deltas too, but because many factors can influence prices we can only suspend betting until we receive the accompanying state change; this is still useful to protect our clients from bad bets.

It is even possible to use this model to extract markets from the algorithms, but the difficulty I had in making this efficient convinced me that this couldn’t be the unifying theory I hoped it would be, it’s just too powerful. However tt did lead me most of the way to the current design, the market definitions fall into one of a small number of groups, so I have taken those as basis functions.

Sports as FSMs

This section is based on a draft email that walks the reader through my thought process in developing a replacement of our event state management system. Open questions include how to model time; what to do about the explosion of states if we model several scores, for example goals, corners and cards in soccer; and what to do about cyclic machines as in tennis.

We have thousands of lines of code that answer the questions “what can happen?” and “what has happened?” throughout the system, obviously algorithms focuses on the first, but we have to interact with plenty of code that answers the second. Of course the problem with code is that it’s easy to let it get out of sync. Nowhere is this more clear than tennis where the different formats are supported to varying degrees, for the longest time the serving rules in tiebreaks were subtly wrong but only if a trader was manually controlling the state.

We model the state of an event as an ordered list of incidents such as “the home team scored at minute X” or “the serving player scored”. Given a history the two fundamental operations are:

get-state(history: [Incident]) -> State
Returns a summary of history.
set-state(history: [Incident], target: State): [Incident]
Adjusts history such that get-state returns target.

Incidents remind me of transitions in a finite-state machine, consider a three-set tennis match, we can model it as an FSM as follows:

State machine for three-set tennis.

There are many ways to encode this into a program, I have chosen to write a function because we’d need an infinitely large table to represent advantage sets once we start modeling games:

>>> def tennis_set((s1, s2)):
...     if s1 == 2 or s2 == 2:
...         return {}
...     else:
...         return {1: (s1+1, s2), 2: (s1, s2+1)}
>>> tennis_set.initial = 0, 0
Codified tennis FSM.

We can trivially build get-state and set-state from this function:

>>> def get_state(fsm, history):
...     state = fsm.initial
...     for transition in history:
...         state = fsm(state)[transition]
...     return state

>>> get_state(tennis_set, [])
(0, 0)
>>> get_state(tennis_set, [1, 2])
(1, 1)

>>> def set_state(fsm, history, target):
...     open = {(get_state(fsm, history), ())}
...     closed = set()
...     while open:
...         current, path = open.pop()
...         if current == target:
...             return history + list(path)
...         else:
...             closed.add(current)
...             for transition, next in fsm(current).iteritems():
...                 if next not in closed:
...                     open.add((next, path + (transition,)))
...      raise ValueError("unreachable")

>>> set_state(tennis_set, [1], (2, 0))
[1, 1]

>>> set_state(tennis_set, [], (2, 0))
[1, 1]

>>> set_state(tennis_set, [], (3, 0))
Traceback (most recent call last):
  ...
ValueError: unreachable

>>> set_state(tennis_set, [], (1, 1))
[2, 1]
get-state and set-state.

The weakness of using functions to define transitions is apparent, we should be able to prune the search tree, there’s no need to search past (1, 0) if the target is (0, 2) but that knowledge isn’t reified anywhere.

Even worse, set-state has a problem, when there are multiple paths to a state it picks an arbitrary one, this is actually how the existing system behaves, but if one day the algorithms start caring about the order of events we’d need to solve it properly. One approach is to record an ordered list of states, and to store the transitions between each state only if known.

Another approach — proposed in my email — is to stored all possible transitions, but this is complicated by cyclic graphs such as tennis games with deuces, and the benefits of this are smaller; knowing that there multiple possible transitions affects correctness whereas knowing which transitions they are only affects completeness. The upshot is that resulting does a slightly better job because it can eliminate some cases that didn’t happen.

Failure Reporting

The migration to microservices required a rewrite of the part of our code responsible for interfacing with the ORM. This code was tied up with the calculation process as a whole, and so I decided to take the chance to improve our error reporting. The previous system was implemented as for loops with trys to isolate failures and logging and continues to skip the markets that didn’t need calculating. This meant that to answer a simple question such as “why isn’t this market priced?” you’d have to consult the logs on the production boxes, even for user errors such as “the algorithm can’t price that market” or “the algorithm parameters are not sufficient”.

def calculate(state, model, markets):
    calculated = {}
    for market in markets:
        if state.result_known(market):
            continue

        if not model.has_structure(market):
            continue

        selections = model.calculate(market)

        if not is_valid(selections):
            continue

        calculated[market] = selections

    return calculated
Abridged legacy calculation process.

We could have addressed this by introducing a list of reasons why markets were skipped and making sure to append to it whenever we make a decision, but that could only be enforced with code reviews and I much prefer having code police itself. When you think about it, the bodies of our loops are essentially functions on a single market that returns either a processed market or a reason why that market could not be processed, and we want to apply those functions in series to all the non-skipped markets, a.k.a. the either and list monads.

> import Control.Applicative
> import Control.Monad

> data Team = Home | Away
> data State = Soccer { goals :: [Team] }
> data Market = Market { name :: String, structure :: String }
              | Calculated { name :: String, selections :: [Float] } deriving (Show)

> let resultKnown s Market { name="winner", structure=_ } = False ;
      resultKnown s Market { name="first", structure="goal" } = (length $ goals s) >= 1 ;
      resultKnown s m = False

> let hasStructure mod Market { name=_, structure="goal" } = True ;
      hasStructure mod m = False

> let isValid Calculated { name=_, selections=s } = sum s == 1

> let calculate mod m = Right Calculated { name=name m, selections=[0.25, 0.5, 0.25] }

> let markets = [ Market { name="winner", structure="goal" }
                , Market { name="winner", structure="corner" }
                , Market { name="first", structure="goal" }
                ]

> let state = Soccer { goals=[Home] }
> let mod = undefined

> -- In Python we derive `msg` from the names of the functions, each
> -- resultKnown function has a specific human-readable name.
> let require p msg m = if (p m) then (Right m) else (Left (m, msg))

> fmap (pure >=>
        require (not . resultKnown state) "result known" >=>
        require (hasStructure mod) "no structure" >=>
        calculate mod >=>
        require isValid "invalid sum"
       ) markets

[ Right (Calculated {name = "winner", selections = [0.25,0.5,0.25]})
, Left (Market {name = "winner", structure = "corner"},"no structure")
, Left (Market {name = "first", structure = "goal"},"result known")
]
Monadic calculation process.

The Left instances contain human-readable strings that are displayed to the traders as tooltips on the buttons that control whether a market is offered or not. Of course in Python we’re using exceptions to early-exit with reasons and implementing all the control flow ourselves, but apart from being wordy there’s no difference.

Testing DSLs

The algorithms codebase is particularly uniform, there are a handful of functions to test and a large dataset to test them on. Inspired by Jay Fields’ Expectations we have a family of test-suite generating functions that provide the “when” when given data structures representing “given” and “then”:

# results(state, market_results)
ht_1_0 = {'home': 1, 'away': 0, 'time': 45}

results(ht_1_0, {
    who_will_win: {'home': Unknown, 'draw': Unknown, 'away': Unknown},
    who_will_score_goal(1): {'home': Win, 'away': Lose},
})


# extracts(structure, market_selections)
# The structure constructors return symbolic data.
grid2x2 = Grid(2, 2)

extracts(grid2x2, {
    who_will_win: {'home': {(1,0)}, 'draw': {(0,0), (1,1)}, 'away': {(0,1)}},
    gt_x_goals(1): {'yes': {(1,1)}, 'no': {(0,0), (1,0), (0,1)}},
})


# derives(algorithm, state, drivers, parameters)
prematch = {'home': 0, 'away': 0, 'time': 0}

derives(soccer, prematch, {
   who_will_win: {'home': 0.300, 'draw': 0.500, 'away': 0.200},
   gt_x_goals(2): {'yes': 0.500, 'no': 0.500},
}, {
    'supremacy': 0.1,
    'expectancy': 2.5,
})
Test DSLs.

The idea is to remove all the boilerplate so that there’s no excuse for not writing tests when a new market is added. The terseness also reduces the noise that a code reviewer could be distracted by, which in conjunction with well-named “given”s make it easier to detect erroneous tests.

The format of the expected values is designed for convenience, for example the function extracts tests returns probabilities not indices, so the test runner generates a handful of grids with random probabilities and checks that the extracted probabilities match the sum of the cells.