Counting Hamiltonian triangulations (OEIS 3122)

I'm trying to translate formulas from the 3122th sequence in OEIS database, in particular this one:

B(s, m) = sum((m! / m_1! ... m_s!) * r(1)^{m_1} ... r(s)^{m_s})
where the sum is over all partitions of s such that
s = m_1 + 2*m_2 + ... + s*m_s and m = m_1 + m_2 + ... + m_s.

My complete solution is available at GitHub.
First question. What the b(s,m) function really does isn't important, all that matters here is that I would like to store already computed values. Can one perform the memoisation in an elegant and smarter way (like Haskell here achieves)? It speeds up the computation a lot.

memoB = [[-1 for x in range(100)] for y in range(100)]

def b(s, m):
    if memoB[s][m] == -1:
        # some calculations
        memoB[s][m] = total
    return memoB[s][m]

Second (algorithmic) question.
The slowest part of my code generates all integer partitions m = m_1 + m_2 + ... + m_s such that s = m_1 + 2m_2 + ... + sm_s. For example, if s = 10, m = 5, these are:

[4, 0, 0, 0, 0, 1, 0, 0, 0, 0]
[3, 1, 0, 0, 1, 0, 0, 0, 0, 0]
[3, 0, 1, 1, 0, 0, 0, 0, 0, 0]
[2, 2, 0, 1, 0, 0, 0, 0, 0, 0]
[2, 1, 2, 0, 0, 0, 0, 0, 0, 0]
[1, 3, 1, 0, 0, 0, 0, 0, 0, 0]
[0, 5, 0, 0, 0, 0, 0, 0, 0, 0]

The current algorithm (coded below) finds the partitions that I am looking for, but is inefficient. Can it be improved? If yes, how?

```
def partitions(s, m):
stack = [] # contains partial partitions
matches = [] # stores good partitions
for i in range(1 + min(s, m)):
stack.append(([i], s - i, m - i))

while len(stack) > 0:
seq, s_sum, m_sum = stack.pop()
# seq - initial section of a partition
# s_sum - difference between s and partial sum of seq
# m_sum - difference between m and partial sum of seq

# if seq has already s summands and both partial sums agree,
# it is marked as a good partition
if len(seq) == s:
if s_sum == 0 and m_sum == 0:
matches.a

To start with, the spec is badly phrased. I know it's not your fault, but it might help you (and would almost certainly help reviewers) to clarify it. It talks about partitions of s but since m = m_1 + m_2 + ... + m_s the m_i are actually a partition of m, and it requires unpicking of the definition. I suggest the following rephrasing:


$$B(s, m) = \sum_{\lambda \vdash s,\, |\lambda| = m} \binom{m}{a_1, \ldots, a_k} r(1)^{a_1} \ldots r(k)^{a_k}$$
where the sum is over all partitions of \$s\$ into exactly \$m\$ parts, \$\lambda = 1^{a_1} 2^{a_2} \ldots k^{a_k}\$ in the frequency representation.

Memoisation support in Python

30 seconds with Google turned up https://docs.python.org/3/library/functools.html#functools.lru_cache . If you know that you're only going to want 10000 different arguments memoised then

@lru_cache(maxsize=10000)
def b(s, m):
    ...

Optimising the integer partition search

How do you form a partition of \$s\$ into \$m\$ parts? You take a part \$p\$ and add it to a partition of \$s - p\$ into \$m - 1\$ parts. To avoid duplication you enforce that the parts are decreasing.

def partitions(s, m):
    return partitions_inner(s, m, s)

def partitions_inner(total, num_parts, max_part):
    # Each of the parts must be between 1 and max_part
    if total < num_parts or num_parts * max_part < total:
        return []
    if num_parts == 1:
        return [[total]]

    results = []
    for part in range(1, max_part + 1):
        for subpartition in partitions_inner(total - part, num_parts - 1, part):
            results.append([part] + subpartition)
    return results

This should be about as efficient as your code, but has the great advantage that partitions_inner can be memoised. Converting to frequency representation is quite simple, and in fact you could refactor this to produce output directly in frequency representation.

But that's irrelevant

If you look at how the \$B(s, m)\$ are used, you see that it's only in \$A(n, s) = \sum_{m=1}^{s} \binom{n}{m} B(s, m)\$. So \$B\$ sums over all partitions with a given number of parts, and then is summed over all numbers of parts. If you cut out the middle-man you get just $$A(n, s) = \sum_{\lambda\vdash s} \binom{n}{a_1, \ldots, a_k} r(1)^{a_1} \ldots r(k)^{a_k}$$ where the sum is over all partitions of \$s\$, \$\lambda = 1^{a_1} 2^{a_2} \ldots k^{a_k}\$ in the frequency representation.

Generating all partitions is (obviously) simpler than generating all partitions with a given number of parts.

def partitions_inner(total, max_part):
    if total == 0:
        return [[]]

    results = []
    for part in range(1, min(max_part, total) + 1):
        for subpartition in partitions_inner(total - part, part):
            results.append([part] + subpartition)
    return results

Again, this can be memoised. It's also possible to trade off memory for speed by splitting up the multinomial as $$\begin{eqnarray}T(\lambda) & = & \binom{m}{a_1, \ldots, a_k} r(1)^{a_1} \ldots r(k)^{a_k} \\
A(n, s) & = & \sum_{\lambda\vdash s} \binom{n}{m} T(\lambda)\end{eqnarray}$$ where the sum is over all partitions of \$s\$, \$\lambda = 1^{a_1} 2^{a_2} \ldots k^{a_k}\$ in the frequency representation, and \$m = \sum_i a_i\$. \$T(\lambda)\$ can be memoised, but also it can be calculated incrementally by using a more complicated object to store the partitions instead of just an array of the parts:

class Partition:
    def __init__(self, maxPart, tail = None):
        self.part = maxPart
        self.tail = tail
        if tail is None:
            self.partFreq = 1
            self.numParts = 1
            tailTerm = 1
        else:
            self.partFreq = tail.partFreq + 1 if tail.part == self.part else 1
            self.numParts = tail.numParts + 1
            tailTerm = tail.term
        self.term = tailTerm * r(self.part) * self.numParts // self.partFreq

@lru_cache(maxsize=100)
def partition(n):
    accum = []
    for maxPart in range(1, n):
        for tail in partition(n - maxPart):
            if tail.part > maxPart:
                break

            accum.append(Partition(maxPart, tail))
    accum.append(Partition(n))
    return accum

(Note that I'm not a Pythonista and I haven't checked the above code against PEP8).

StackExchange Code Review Q#155933, answer score: 7