The Philosopher's Language

A long time ago, there was a philosopher. He had realized that most of the problems that the other philosophers were dealing with, was due to them not having a good enough language to express their thoughts with. A lot of their effort was put into defining new terms and concepts to add to their existing language. Still, they seemed to always quarrel about the meaning of very basic notions, whose nature should be obvious to every human. So our man set out to recreate language. Language, which has been with humans since the very beginning. Perhaps it even defines what a human is. From the ground up, he started recreating humanity's closest ally. A new Language, for those who truly seek Knowledge. A Language which would focus the powers of the mind to think only clear thoughts, and allow one to share those thoughts with others, without confusion. His ambition was greater than that of any philosopher who came before him, or any who has come after him since. His whole working life, he dedicated to this language. In the end, he had made a large book about it. When the book was being published, he had to come up with a name for it. Because he had worked so hard on the Language, he knew he scope of his ignorance. Even though he had laid down a lifetime, he could not even scratch the surface of the Language. Others would have to come after him and expand upon it. So he gave the book a name that meant that its content was Basic, but also that it was something Foundational. He named it Elements.

Turing machines

I want to dig deeper into the meaning of computation, and the limits of what programming can be. As part of that, I read through Turing's 1936 paper [1]. Turing defines a computable sequence as any sequence (possibly infinite) that is the output of a program of finite size. He goes on to describe a finite set of symbol manipulation operations that is sufficient to compute any computable sequence. This later became known as a Turing machine. The operations themselves can be assigned symbols, so a Turing machine can be designed to emulate other Turing machines, given only a description of the emulated machine as input. Such a system is called a Universal Turing machine.

In the same very productive paper, Turing also proves that deciding if an arbitrary program will output an infinite number of symbols, cannot be done in finite time. His proof is quite cognitively challenging, but Tom LeBlanc has a version that is easier to follow [2].

I also implemented one of his early examples, which prints the sequence

001011011101111...

A solution to this is described in table form:

Turing uses very stylistic letters as notation, which are ironically not very computer-friendly

My code for this:


I keep track of the tape using a dict. Both tape, head, and state are global variables.

Conclusions
Using a Turing machine is a very pure way of programming. That is to say, one works with an extremely small set of operations, and one needs to model every detail of the problem. It feels like building a house by starting collecting mud from the creek instead of buying bricks.

It is not an efficient way of programming, even if one would implement the operations directly in hardware. Modern hardware is better than being able to see a single symbol at a time. It is not efficient for the programmer either. It is tricky to debug and a hard way to model problems. The Turing machine mostly has value as a theoretical tool.

[1] On Computable Numbers (Turing 1936)
[2] Undecidability

Sampling: Computing stationary distribution

Previous:

• Idea for series sampling
• Mutating Box Trees

Sampling a series with the box tree defines a Markov chain. A valid set of states is to split the first dimension everywhere where a box ends:

Given a value on the x-axis, the sampling distribution is uniquely defined by which interval (confined by red lines) the value is in. 

Realizing this equivalence with respect to the sampling, gives us a finite number of states. We can calculate the transition probabilities from a given state by integrating the conditional distribution over the states' intervals. The conditional distribution is just a piecewise uniform distribution, so integrating over an interval is simple. This gives us the transition matrix P. 

If the stationary distribution exists, it can be found by solving:

Where we have the additional constraint that the sum of x should be 1. This can be solved by scipy.linalg.nullspace. Using this stationary distribution instead of relying on series sampling gives us a faster and also exact measure of the series entropy. 

Greedy search for transition function when using theoretical values for entropy, instead of samples. An improvement from before is that the optimization doesn't get stuck and is unable to find a better neighbour. However, it clearly gets stuck in a local valley.

This leads to a better search for interesting transition functions (low prior entropy, high posterior entropy), in the sense that it is able to keep optimizing for a long time. The shortcoming of the sampling-based measure was that sometimes the finite sample sequence was a "lucky" one that gave a good score, but perhaps was not representative for the ensemble average. Then, it was hard for it to find a better neighbour, since the current state was overestimated. Using the Markov chain model solves this. However, the search clearly gets stuck in a local valley.

Also, the measure itself should have some addition. I would like to premier transition distributions with a lot of almost-stationary states. This would make the series stay in those states for a long time before transitioning to another state. Still, we don't want it to get stuck in those states. The property of not getting stuck in some states is called ergodicity. We should look for almost non-ergodic processes.  

Zen of Python: Sparse is better than dense

In the sampling mini project, and also in the 2D interval selection problem, I had to implement a few operations on intervals one dimension and higher. Let's call intervals in an arbitrary number of dimensions boxes.

• contains: given a box and a point, checks if the point is contained in the box. 
• intersection: given two boxes, returns the overlap between them (if there is any). 
• profile: given a box and a point specified in a subset S of the box's dimensions D, return the box's expansion in D\S (the remaining dimensions), if the point is contained in the box's projection onto S. We can see this as intersect between the box and the point expanded infinitely in D\S. 
• volume: given a box, return the product of its size in all its dimensions.
• corners: given a box, return all its corners as tuples. 

A simple google search did not suffice to find an impressive python package that supports these relatively abstract operations. I decided to clean up my code to make it a bit more general. It struck me that a good data structure for boxes might be a dict. The values are 2-tuples (or empty tuples for empty intervals) and the keys, representing the dimensions, can be any hashable. Even if I until now have worked with things that are naturally contained in euclidean space, I can imagine getting use for a sparse representation of intervals whenever dealing with a product space of ordered sets. Some bonus features:

• Handle empty intervals
• Projection onto arbitrary dimension is trivial. Just a key access, like so: box[dim].

Actually, this concept can be even further generalized. The only operations above that rely on the dimension being ordered are volume and corners. The others we can get as long as we have defined intersection on each of the dimensions in use. The general concept is a product space, stored only as a collection of expansions in each dimension. 

Sampling: mutating box trees

Previous:

• Idea for series sampling

In the previous post, I sampled time series by defining transition probabilities using boxes of different sizes. We want to find the most complex time series using the least complex representation of the transitions. We need to start searching for better transition functions in the box representation.

First, we need to mutate the box representation. We can see the boxes as a non-balanced binary tree. Each node represents a box. Each box is either split or not split - i.e. each node has either 0 or 2 children. The root node is the box representing the whole space. I mutate the box tree by removing both children of a non-leaf box. I then insert as many leafs as was lost in the remove operation. So, the number of leafs are kept the same. This gives a very quickly changing box landscape:

Random mutations on the box space. All sub-boxes in a random box are deleted. Then, new leafs are inserted, keeping the number of boxes the same. 

Now we need a measure for the complexity of the box tree. I use area-wise entropy:

Where A_i is the area of the box and A is the total area. This is a measure of how much we can compress a point's location, if we know that it is in one of the boxes with equal probability. If we do a greedy search to maximize the entropy, we get:

Greedy search to maximize box entropy. Optimum is that all boxes have same area. There are 64 boxes, and the total area is 1, so ideally they should all have the area 2^(-6).

Greedy search to minimize entropy. Optimum is that one box covers the whole area and the others have 0 area. This cannot be represented by the box tree, so it produces a geometric series of boxes.

We need also a score for the resulting time series. Ideally, the series should cover a lot of the box space. I use the count-wise entropy:

Where C_i is the count of the series for box i, and C is the total number of elements in the series. With this, we can do a greedy search to maximize S_series - S_boxes:

Greedy search to maximize S_series - S_boxes. Estimating S_series by sampling is not very reliable.

The problem is that we cannot reliably find out the entropy of the series, because it relies on the stationary distribution. At this point, I realize that what we are doing actually defines a Markov chain. And for Markov chains, we can calculate the stationary distribution. So, that is a good next step in order to get a better measure of process complexity. 

Sampling: series

Previous:

• sampling points in box shaped clusters

In an earlier post, we used the following design goal for a good sampling:

"generate a set of random points with a lot of features, using minimal coding effort"

An approach based on generating a random (non-complete) binary tree of boxes gave nice results:

From previous post

We can plot the boxes for clarity:

Fractally generated boxes.

Fractally generated boxes, with random sampled points.

Now I had the idea that this could be used to generate interesting time series. We can put the current value on the x-axis, and sample the next value from the profile distribution. Illustration:

Profile / conditional distribution of y(t) when y(t-1)=0.4. Notice that box density = conditional density. 

So let's generate some time series with this. We start at y(0) = 0, and then update recursively using the conditional distribution induced by the boxes. 

Sampling directly from the conditional distribution. 

Well, that is a bit disappointing. The generated time series just looks kind of random. It is quite stateless. We can try residual sampling instead. Then we change the update equation, from:

to:

B is the box conditional distribution, and lambda is a small coefficient. This makes it more like a derivative. An other way of seeing it, in terms of the direct sampling, is that is is as if we have forced all boxes to be around the y(t)=y(t-1) line. It also gives more interesting time series:

Residual sampling - we only make small updates. Time series has more "features". However, it only explores a small part of the state space. 

This method looks more promising, but clearly we want to get better at generating conditional distributions that will make the sampling explore a larger portion of the state space. 

SAT: 2D interval selection

Previous:

• SAT: optimal interval selection

"Cut rectangular slices of pizza. The pizza is divided into cells. Each slice must contain at least L tomato cells (T) and L mushroom cells (M). Each slice must contain at most H cells [1]. Maximize the total number of cells in the cuts."

In the last post, we solved a 2D slicing problem by optimizing each row independently. We should expect to do better by actually selecting 2D intervals instead. A 2D interval is represented as:

(x_start, x_end, y_start, y_end)

To generalize the code to 2D intervals, we do not have to change much. This is thanks to the fact that we only dealt with intervals in a very abstract sense before: as elements with a size (for the score) and an overlapping relation (to forbid overlapping slices). We need only to change the way we generate the intervals, and the overlapping relation. The code:

Now we need a solution strategy. Tests show that the 2D interval selection gets tough (starts taking more than 10 seconds) when solving for larger than 20x20 matrices. However 12x12 matrices can be solved quickly (about 0.6 seconds). Solving with 12x12 sections often give a perfect score on the subproblem:

Optimization on 12x12 subproblem. Colored borders mark slices. 

Optimization on 12x12 subproblem. Colored borders mark slices.
In many cases, almost all slices are just rows or columns.

Now let's solve for the medium instance [1]. The max score is 50000. The row-wise maximization gave us 48977. The top result on google for "hashcode pizza" uses a greedy algorithm and gets 49567 on medium [2]. This time we get a score of 49980 missing only 20 cells. This clearly outperforms the other approaches. The runtime is about 10 minutes, without parallelization.

Conclusion
This result does not only say something about the power of model-based approaches, but also about the problem itself. It is quite an easy problem where we can get an almost perfect score by just solving independent subproblems.

[1] pizza.pdf
[2] sagishporer/hashcode-practice-pizza