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

SAT: optimal interval selection

We are tasked with the following:

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. Example pizza:

TTTTT
TMMMT
TTTTT

L = 1, H = 6

Optimal solution:

TT | T | TT
TM | M | MT
TT | T | TT

I want to find (almost) optimal solutions using SAT. In this first version, I reduce it to a 1D problem, assuming that we only cut row by row. In other words, no slice cover more than one row.

Some problem analysis:

• Without involving the SAT solver, we can calculate which slices are allowed in the solution.
• Also without invoking in SAT, we can calculate which slices cannot overlap. 
• Determining which cells are covered by which slices can also be done easily without modelling. 
• With all this data known beforehand, we can reduce the slicing to something very abstract. We can see it as a selection problem: select some elements of a set with some mutual constraints and some objective. 

The key here is that the mutual constraints and the objective can be easily expressed using clauses. First we enumerate all feasible slices. Then, for each pair of slices (s_0, s_1), we add the following clause if they overlap:


For each cell, c, we find the slices s_0, s_1,...s_k that overlap it. Then for each cell we add an indicator:

Now we want to find the largest S such that adding

is satisfiable. For this, we can use pysat's ITotalizer [2]. It is based on technology that allows for some optimizations when trying different values of S [3]. I applied it to a binary search. The code:

For the instances medium.in and big.in in the input data [1], this algorithm can solve row-wise optimality in reasonable time. Sample output medium.in (max score is 250):

Row time: 0.510
Row score: 250

Row time: 0.484
Row score: 242

Row time: 0.519
Row score: 247

Sample output big.in (max score is 1000):

Row time: 4.096
Row score: 880

Row time: 3.108
Row score: 904

Row time: 3.336
Row score: 881

For medium I get a total score of 48977, which is not particularly good, because it can be beaten by a greedy algorithm (49567) [4]. The max score is 50000 (250x200). Still, it is interesting that we get so close to the real optimum with such a silly constraint as row-wise maximization. 

[1] pizza.pdf
[2] PySAT: SAT technology in Python
[3] Martins, Ruben, et al. "Incremental cardinality constraints for MaxSAT." International Conference on Principles and Practice of Constraint Programming. Springer, Cham, 2014.
[4] sagishporer/hashcode-practice-pizza

SAT: parity board

We get the following challenge:

I take "clear" to mean that all tiles should be the same color, either black or white. One "move" is to click a tile, which changes the color of it, and also its 4-connected neighbors. Illustration:

Clicking a tile changes the color of it, and its 4-connected neighbors.

This can be solved with SAT. Some analysis of the problem:

• Clicking a tile twice is the same as not clicking it at all. The only thing that matters is parity of clicks. If we assume that we are looking for an optimal solution, we can assume that each tile is clicked 0 or 1 times. So, they can be modeled as boolean variables. 
• Parity is also the only thing that matters for clearing the board. Suppose we want to make all tiles white. Then, black tiles should switch color an odd number of times. White tiles should switch color an even number of times. 
• If we can count the number of color switches each tile has done, we can add cardinality constraints on the click variables that cover each tile. 
• We also have to add a cardinality constraint on the whole set of click variables. 

For small numbers like this, we can do parity as a disjunction of equality clauses. 

This produces the following solution:

0  0  1  0  0

0  0  1  0  0

1  0  0  1  1

0  1  0  1  1

0  0  0  1  0

(1 means click this tile, 0 means no click)

SAT: disjunction operator

Previous:

1. getting started
2. langford pairs
3. how not to "negate"

I want to create an operator that turns the disjunction of a set of CNFs, into one CNF. That is to say, for c ... d given, find g such that:

I got no relevant hits on "disjoint" or "disjunction" in the pysat documentation. The implementation of the operator is not very complicated, however.

The idea is to do a Tseytin transformation. For each CNF in the input, we create a variable p that is equivalent to the satisfaction of the CNF. Call this variable an indicator variable. We have:



This relationship can be turned into a CNF (standard AND-constraint):



So the method is: create indicators and indicator clauses for all input CNFs. We add all indicator clauses together into a joint CNF. Lastly, we need to add a clause with all indicators, since we need to make sure that at least one indicator is true. The code:

Now we can test it. Let's see if we can make a CNF that is true if and only if the sum of the variables is even. The code:

I modified enumeration_test so that if outputs which variable sums are present in the satisfying assignments, and in the unsatisfying assignments, respectively. The output is as expected:

Satisfying assignments:
{0, 2, 4, 6, 8, 10}

Unsatisfying assignments:
{1, 3, 5, 7, 9}

In conclusion, we could do disjunction as expected. Note that auxiliary variables were present here (implicit from the "equals" clauses, which use sequential counters), as in the negation experiment. But in this case, it worked fine.