3
\$\begingroup\$

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 + 2*m_2 + ... + s*m_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.append(seq)
            continue

        # if one of the partial sums drops below zero, 
        # the seq can't be continued
        if min(s_sum, m_sum) < 0 and max(s_sum, m_sum) > 0:
            continue

        # we extend the stack by all "one step continuations" of seq
        factor = len(seq) + 1
        bound = max(1, 1 + min(s_sum // factor, m_sum))
        stack.extend([(seq + [i], s_sum - i * factor, m_sum - i) for i in range(bound)])
    return matches

I know that one shouldn't name functions with single letters, but these are taken from the OEIS description directly.

\$\endgroup\$
  • 2
    \$\begingroup\$ Can you quote a meaningful extract of what this code should accomplish out of the link you gave? \$\endgroup\$ – 409_Conflict Feb 21 '17 at 14:18
  • \$\begingroup\$ @MathiasEttinger, I have shortened the code - I am mainly having problems with efficient implementation of the B(s, m) function from oeis.org/A003122 \$\endgroup\$ – Santiago Feb 21 '17 at 14:53
  • 2
    \$\begingroup\$ This is not really what I asked for, it's quite the opposite. Here at Code Review, we prefer to see all the ins and outs of the code, so the first version was fine as it allowed us to see the big picture. What I was asking was: "can you add a better (textual) description" of what the code is supposed to do? \$\endgroup\$ – 409_Conflict Feb 21 '17 at 15:00
  • \$\begingroup\$ @MathiasEttinger I'm not sure whether I have accomplished what you had asked for, but I tried to make my questions more clear. \$\endgroup\$ – Santiago Feb 22 '17 at 8:58
7
\$\begingroup\$

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).

\$\endgroup\$
  • \$\begingroup\$ Thank you very much for detailed answer. However I think that there is something wrong with your very first definition of $A(n, s)$: for $n = 3$, $s = 2$ it gives ${3 \choose 2} \cdot 1^2$ (from $2 = 1 + 1$) plus ${3 \choose 0,1} 1^0 \cdot 4^1$, e.g. $3 + 24 = 27$ instead of $39$. Have I missed something? \$\endgroup\$ – Santiago Feb 24 '17 at 11:17
  • \$\begingroup\$ @Santiago, I don't know how you get 27 or 39. I get 15 by both approaches. \$\endgroup\$ – Peter Taylor Feb 24 '17 at 12:19
  • \$\begingroup\$ here are my calculations: \$A(3, 2) = \frac{3!}{1!} \sum_{\lambda \vdash 2} \prod_i \frac{r(i)^{a_i}}{a_i!} = 6 \left(\underbrace{\frac{1^2}{2!}}_{\lambda = 1^2} + \underbrace{\frac{1^0}{0!}\frac{4^1}{1!}}_{\lambda = 1^0 2^1}\right) = 6 \left(4 + \frac 12\right) = 27\$ @PeterTaylor \$\endgroup\$ – Santiago Feb 24 '17 at 18:35
  • \$\begingroup\$ @Santiago, ah! That's not my very first definition of \$A(n, s)\$, which gives 15. You're right: when I tried to factor out common terms I got mixed up between m and s. \$\endgroup\$ – Peter Taylor Feb 24 '17 at 23:22
  • \$\begingroup\$ unfortunately it doesn't give the correct answer for all \$n \ge 5\$ because \$a(3, 5) = 7191\$ (in my original code) and not \$7195\$ (as your formula suggests). \$\endgroup\$ – Santiago Feb 26 '17 at 16:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.