# Question
Title: Enumerating all partial permutations of given length in lexicographic order
I need to generate all unique tuples of length *k* chosen from a series of unique, positive integers. In my case *n* choose *k* will have *n*=10, 1 \<= *k* \<= 10; and the series I am choosing from is { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 } .

For each chosen tuple, all of its permutations need to be present in the output. For example, for *k* = 2, the chosen tuple (0, 1) needs to have both (0, 1) and (1, 0) in the output. The output overall needs to be in lexicographic order, which in turn implies that lexicographic ordering of permutations within a given chosen combination is not enough.

I'm following along with

https://www.cs.utexas.edu/users/djimenez/utsa/cs3343/lecture25.html

which offers a simple recursive algorithm to generate combinations. Their recommendation to do this:

> Note that you can get all permutations of n things taken k at a time by simply calling perm (v, maxk, 0); at the base case of combinations. This generates all k! permutations of each of the n C k combinations, taking O(k! n (n C k)) = O((n+1)!/(n-k)!) time.

gets me 90% of the way there, but yields result tuples in an order that is not useful to me. For k=2 it yields

```
[0, 1]
[1, 0]
[0, 2]
[2, 0]
...

```

where I would want it to yield

```
[0, 1]
[0, 2]
[0, 3]
[0, 4]
[0, 5]
[0, 6]
[0, 7]
[0, 8]
[0, 9]
[1, 0]
[1, 2]
...

```

Is there a modification to this algorithm - or a different one - that will do this? I hope it doesn't require putting everything in memory and then applying an in-place sort.

# Answer
Here is a simple iterative solution. We maintain two arrays, an output array $L$, and a Boolean array $A$, which keeps track of the elements currently in $L$. We update $A$ as we add and remove elements from $L$.

We initialize $L$ with $0,\ldots,k-1$, which is also our first output. Now we repeatedly try to "increment" $L$. If successful, we output the incremented $L$, and try to increment $L$ again. If unsuccessful, we have gone through all combinations, and can terminate.

In order to increment $L$, we start by trying to increment position $k-1$ (this is the rightmost position). We do so by scanning the array $A$ starting at $L\[k-1\] + 1$, until we find an element $x$ which is not in $L$. If this search is successful, we replace $L\[k-1\]$ with $x$. Otherwise, we try to do the same for $L\[k-2\],\ldots,L\[0\]$, until successful; if all of these fail, we declare failure (and terminate).

Suppose that we have successfully replaced $L\[j\]$. If $j = k-1$, we're done. Otherwise, we have to update $L\[j+1\],\ldots,L\[k-1\]$, which we do by scanning $A$ starting at $0$, inserting the first $k-j-1$ elements not already in $L$.

Here is a Python implementation:

```
def gen(n, k):
    A = [True] * k + [False] * (n-k)
    L = list(range(k))
    yield tuple(L)
    while True:
        j = k - 1
        while j >= 0:
            A[L[j]] = False
            t = L[j] + 1
            while t < n and A[t]:
                t += 1
            if t < n:
                A[t] = True
                L[j] = t
                r = 0
                j += 1
                while j < k:
                    while A[r]:
                        r += 1
                    A[r] = True
                    L[j] = r
                    j += 1
                yield tuple(L)
                break
            else:
                j -= 1
        if j == -1:
            break

```

---

An even simpler recursive solution simply iterates over all possible choices of the first element; within these, iterates over all possible choices of the second element; and so on.

Here is a Python implementation:

```
def gen_rec(n, k, L=[]):
    if len(L) == k:
        yield tuple(L)
    for i in range(n):
        if i not in L:
            yield from gen_rec(n, k, L + [i])

```
> 4 votes

# Answer
The most appealing solution so far seems to be a spin on Python's own `permutations` function (source) that can be slightly simplified for this use case; thanks @Peilonrayz:

```
def combinations(r: int, n: int=10) -> Iterable[Sequence[int]]:
    pool = range(n)
    indices = list(pool)
    cycles = list(range(n, n-r, -1))
    r_range = range(r-1, -1, -1)

    yield pool[:r]

    while True:
        for i in r_range:
            cycles[i] -= 1

            if cycles[i] > 0:
                j = cycles[i]
                indices[i], indices[-j] = indices[-j], indices[i]
                yield indices[:r]
                break

            indices[i:] = indices[i+1:] + indices[i:i+1]
            cycles[i] = n - i
        else:
            return

```

As best as I can tell, the way this works:

* This is an iterative algorithm that swaps elements in a mutable list `indices` of length `n` from which `r` elements will be chosen; these elements will always be chosen from the beginning of the list
* `indices` starts with a list of elements to be chosen in their original order
* A secondary mutable list `cycles` of length `r` controls which elements are swapped
* On an inner condition based on `cycles`, only two elements are exchanged with each other and one output tuple is produced. If that condition is false, an output tuple is not produced, and instead `indices` sees a portion of its elements shifted left by one position, with the middle element being bumped to the end.

I don't think transcribing the code above into pseudocode has much value, so instead I'll output a dump that illustrates the value of these control variables over a few iterations where n=10, r=3:

```
indices=[0, 1, 3, 2, 4, 5, 6, 7, 8, 9] cycles=[10, 9, 7] i=2, indices 2 and 3 were exchanged, output [0, 1, 3]
indices=[0, 1, 4, 2, 3, 5, 6, 7, 8, 9] cycles=[10, 9, 6] i=2, indices 2 and 4 were exchanged, output [0, 1, 4]
indices=[0, 1, 5, 2, 3, 4, 6, 7, 8, 9] cycles=[10, 9, 5] i=2, indices 2 and 5 were exchanged, output [0, 1, 5]
indices=[0, 1, 6, 2, 3, 4, 5, 7, 8, 9] cycles=[10, 9, 4] i=2, indices 2 and 6 were exchanged, output [0, 1, 6]
indices=[0, 1, 7, 2, 3, 4, 5, 6, 8, 9] cycles=[10, 9, 3] i=2, indices 2 and 7 were exchanged, output [0, 1, 7]
indices=[0, 1, 8, 2, 3, 4, 5, 6, 7, 9] cycles=[10, 9, 2] i=2, indices 2 and 8 were exchanged, output [0, 1, 8]
indices=[0, 1, 9, 2, 3, 4, 5, 6, 7, 8] cycles=[10, 9, 1] i=2, indices 2 and 9 were exchanged, output [0, 1, 9]
indices=[0, 1, 9, 2, 3, 4, 5, 6, 7, 8], shifting at 2:
        [      \2, 3, 4, 5, 6, 7, 8,/9]

indices=[0, 2, 1, 3, 4, 5, 6, 7, 8, 9] cycles=[10, 8, 8] i=1, indices 1 and 2 were exchanged, output [0, 2, 1]
indices=[0, 2, 3, 1, 4, 5, 6, 7, 8, 9] cycles=[10, 8, 7] i=2, indices 2 and 3 were exchanged, output [0, 2, 3]
indices=[0, 2, 4, 1, 3, 5, 6, 7, 8, 9] cycles=[10, 8, 6] i=2, indices 2 and 4 were exchanged, output [0, 2, 4]
indices=[0, 2, 5, 1, 3, 4, 6, 7, 8, 9] cycles=[10, 8, 5] i=2, indices 2 and 5 were exchanged, output [0, 2, 5]
indices=[0, 2, 6, 1, 3, 4, 5, 7, 8, 9] cycles=[10, 8, 4] i=2, indices 2 and 6 were exchanged, output [0, 2, 6]
indices=[0, 2, 7, 1, 3, 4, 5, 6, 8, 9] cycles=[10, 8, 3] i=2, indices 2 and 7 were exchanged, output [0, 2, 7]
indices=[0, 2, 8, 1, 3, 4, 5, 6, 7, 9] cycles=[10, 8, 2] i=2, indices 2 and 8 were exchanged, output [0, 2, 8]
indices=[0, 2, 9, 1, 3, 4, 5, 6, 7, 8] cycles=[10, 8, 1] i=2, indices 2 and 9 were exchanged, output [0, 2, 9]
indices=[0, 2, 9, 1, 3, 4, 5, 6, 7, 8], shifting at 2:
        [      \1, 3, 4, 5, 6, 7, 8,/9]

indices=[0, 3, 1, 2, 4, 5, 6, 7, 8, 9] cycles=[10, 7, 8] i=1, indices 1 and 3 were exchanged, output [0, 3, 1]
indices=[0, 3, 2, 1, 4, 5, 6, 7, 8, 9] cycles=[10, 7, 7] i=2, indices 2 and 3 were exchanged, output [0, 3, 2]
indices=[0, 3, 4, 1, 2, 5, 6, 7, 8, 9] cycles=[10, 7, 6] i=2, indices 2 and 4 were exchanged, output [0, 3, 4]
indices=[0, 3, 5, 1, 2, 4, 6, 7, 8, 9] cycles=[10, 7, 5] i=2, indices 2 and 5 were exchanged, output [0, 3, 5]
indices=[0, 3, 6, 1, 2, 4, 5, 7, 8, 9] cycles=[10, 7, 4] i=2, indices 2 and 6 were exchanged, output [0, 3, 6]
indices=[0, 3, 7, 1, 2, 4, 5, 6, 8, 9] cycles=[10, 7, 3] i=2, indices 2 and 7 were exchanged, output [0, 3, 7]
indices=[0, 3, 8, 1, 2, 4, 5, 6, 7, 9] cycles=[10, 7, 2] i=2, indices 2 and 8 were exchanged, output [0, 3, 8]
indices=[0, 3, 9, 1, 2, 4, 5, 6, 7, 8] cycles=[10, 7, 1] i=2, indices 2 and 9 were exchanged, output [0, 3, 9]
indices=[0, 3, 9, 1, 2, 4, 5, 6, 7, 8], shifting at 2:
        [      \1, 2, 4, 5, 6, 7, 8,/9]

indices=[0, 4, 1, 2, 3, 5, 6, 7, 8, 9] cycles=[10, 6, 8] i=1, indices 1 and 4 were exchanged, output [0, 4, 1]
indices=[0, 4, 2, 1, 3, 5, 6, 7, 8, 9] cycles=[10, 6, 7] i=2, indices 2 and 3 were exchanged, output [0, 4, 2]
indices=[0, 4, 3, 1, 2, 5, 6, 7, 8, 9] cycles=[10, 6, 6] i=2, indices 2 and 4 were exchanged, output [0, 4, 3]
indices=[0, 4, 5, 1, 2, 3, 6, 7, 8, 9] cycles=[10, 6, 5] i=2, indices 2 and 5 were exchanged, output [0, 4, 5]
indices=[0, 4, 6, 1, 2, 3, 5, 7, 8, 9] cycles=[10, 6, 4] i=2, indices 2 and 6 were exchanged, output [0, 4, 6]
indices=[0, 4, 7, 1, 2, 3, 5, 6, 8, 9] cycles=[10, 6, 3] i=2, indices 2 and 7 were exchanged, output [0, 4, 7]
indices=[0, 4, 8, 1, 2, 3, 5, 6, 7, 9] cycles=[10, 6, 2] i=2, indices 2 and 8 were exchanged, output [0, 4, 8]
indices=[0, 4, 9, 1, 2, 3, 5, 6, 7, 8] cycles=[10, 6, 1] i=2, indices 2 and 9 were exchanged, output [0, 4, 9]
indices=[0, 4, 9, 1, 2, 3, 5, 6, 7, 8], shifting at 2:
        [      \1, 2, 3, 5, 6, 7, 8,/9]

indices=[0, 5, 1, 2, 3, 4, 6, 7, 8, 9] cycles=[10, 5, 8] i=1, indices 1 and 5 were exchanged, output [0, 5, 1]
indices=[0, 5, 2, 1, 3, 4, 6, 7, 8, 9] cycles=[10, 5, 7] i=2, indices 2 and 3 were exchanged, output [0, 5, 2]
indices=[0, 5, 3, 1, 2, 4, 6, 7, 8, 9] cycles=[10, 5, 6] i=2, indices 2 and 4 were exchanged, output [0, 5, 3]
indices=[0, 5, 4, 1, 2, 3, 6, 7, 8, 9] cycles=[10, 5, 5] i=2, indices 2 and 5 were exchanged, output [0, 5, 4]
indices=[0, 5, 6, 1, 2, 3, 4, 7, 8, 9] cycles=[10, 5, 4] i=2, indices 2 and 6 were exchanged, output [0, 5, 6]
indices=[0, 5, 7, 1, 2, 3, 4, 6, 8, 9] cycles=[10, 5, 3] i=2, indices 2 and 7 were exchanged, output [0, 5, 7]
indices=[0, 5, 8, 1, 2, 3, 4, 6, 7, 9] cycles=[10, 5, 2] i=2, indices 2 and 8 were exchanged, output [0, 5, 8]
indices=[0, 5, 9, 1, 2, 3, 4, 6, 7, 8] cycles=[10, 5, 1] i=2, indices 2 and 9 were exchanged, output [0, 5, 9]
indices=[0, 5, 9, 1, 2, 3, 4, 6, 7, 8], shifting at 2:
        [      \1, 2, 3, 4, 6, 7, 8,/9]

indices=[0, 6, 1, 2, 3, 4, 5, 7, 8, 9] cycles=[10, 4, 8] i=1, indices 1 and 6 were exchanged, output [0, 6, 1]
indices=[0, 6, 2, 1, 3, 4, 5, 7, 8, 9] cycles=[10, 4, 7] i=2, indices 2 and 3 were exchanged, output [0, 6, 2]
indices=[0, 6, 3, 1, 2, 4, 5, 7, 8, 9] cycles=[10, 4, 6] i=2, indices 2 and 4 were exchanged, output [0, 6, 3]
indices=[0, 6, 4, 1, 2, 3, 5, 7, 8, 9] cycles=[10, 4, 5] i=2, indices 2 and 5 were exchanged, output [0, 6, 4]
indices=[0, 6, 5, 1, 2, 3, 4, 7, 8, 9] cycles=[10, 4, 4] i=2, indices 2 and 6 were exchanged, output [0, 6, 5]
indices=[0, 6, 7, 1, 2, 3, 4, 5, 8, 9] cycles=[10, 4, 3] i=2, indices 2 and 7 were exchanged, output [0, 6, 7]
indices=[0, 6, 8, 1, 2, 3, 4, 5, 7, 9] cycles=[10, 4, 2] i=2, indices 2 and 8 were exchanged, output [0, 6, 8]
indices=[0, 6, 9, 1, 2, 3, 4, 5, 7, 8] cycles=[10, 4, 1] i=2, indices 2 and 9 were exchanged, output [0, 6, 9]
indices=[0, 6, 9, 1, 2, 3, 4, 5, 7, 8], shifting at 2:
        [      \1, 2, 3, 4, 5, 7, 8,/9]

indices=[0, 7, 1, 2, 3, 4, 5, 6, 8, 9] cycles=[10, 3, 8] i=1, indices 1 and 7 were exchanged, output [0, 7, 1]
indices=[0, 7, 2, 1, 3, 4, 5, 6, 8, 9] cycles=[10, 3, 7] i=2, indices 2 and 3 were exchanged, output [0, 7, 2]
indices=[0, 7, 3, 1, 2, 4, 5, 6, 8, 9] cycles=[10, 3, 6] i=2, indices 2 and 4 were exchanged, output [0, 7, 3]
indices=[0, 7, 4, 1, 2, 3, 5, 6, 8, 9] cycles=[10, 3, 5] i=2, indices 2 and 5 were exchanged, output [0, 7, 4]
indices=[0, 7, 5, 1, 2, 3, 4, 6, 8, 9] cycles=[10, 3, 4] i=2, indices 2 and 6 were exchanged, output [0, 7, 5]
indices=[0, 7, 6, 1, 2, 3, 4, 5, 8, 9] cycles=[10, 3, 3] i=2, indices 2 and 7 were exchanged, output [0, 7, 6]
indices=[0, 7, 8, 1, 2, 3, 4, 5, 6, 9] cycles=[10, 3, 2] i=2, indices 2 and 8 were exchanged, output [0, 7, 8]
indices=[0, 7, 9, 1, 2, 3, 4, 5, 6, 8] cycles=[10, 3, 1] i=2, indices 2 and 9 were exchanged, output [0, 7, 9]
indices=[0, 7, 9, 1, 2, 3, 4, 5, 6, 8], shifting at 2:
        [      \1, 2, 3, 4, 5, 6, 8,/9]

indices=[0, 8, 1, 2, 3, 4, 5, 6, 7, 9] cycles=[10, 2, 8] i=1, indices 1 and 8 were exchanged, output [0, 8, 1]
indices=[0, 8, 2, 1, 3, 4, 5, 6, 7, 9] cycles=[10, 2, 7] i=2, indices 2 and 3 were exchanged, output [0, 8, 2]
indices=[0, 8, 3, 1, 2, 4, 5, 6, 7, 9] cycles=[10, 2, 6] i=2, indices 2 and 4 were exchanged, output [0, 8, 3]
indices=[0, 8, 4, 1, 2, 3, 5, 6, 7, 9] cycles=[10, 2, 5] i=2, indices 2 and 5 were exchanged, output [0, 8, 4]
indices=[0, 8, 5, 1, 2, 3, 4, 6, 7, 9] cycles=[10, 2, 4] i=2, indices 2 and 6 were exchanged, output [0, 8, 5]
indices=[0, 8, 6, 1, 2, 3, 4, 5, 7, 9] cycles=[10, 2, 3] i=2, indices 2 and 7 were exchanged, output [0, 8, 6]
indices=[0, 8, 7, 1, 2, 3, 4, 5, 6, 9] cycles=[10, 2, 2] i=2, indices 2 and 8 were exchanged, output [0, 8, 7]
indices=[0, 8, 9, 1, 2, 3, 4, 5, 6, 7] cycles=[10, 2, 1] i=2, indices 2 and 9 were exchanged, output [0, 8, 9]
indices=[0, 8, 9, 1, 2, 3, 4, 5, 6, 7], shifting at 2:
        [      \1, 2, 3, 4, 5, 6, 7,/9]

indices=[0, 9, 1, 2, 3, 4, 5, 6, 7, 8] cycles=[10, 1, 8] i=1, indices 1 and 9 were exchanged, output [0, 9, 1]
indices=[0, 9, 2, 1, 3, 4, 5, 6, 7, 8] cycles=[10, 1, 7] i=2, indices 2 and 3 were exchanged, output [0, 9, 2]
indices=[0, 9, 3, 1, 2, 4, 5, 6, 7, 8] cycles=[10, 1, 6] i=2, indices 2 and 4 were exchanged, output [0, 9, 3]
indices=[0, 9, 4, 1, 2, 3, 5, 6, 7, 8] cycles=[10, 1, 5] i=2, indices 2 and 5 were exchanged, output [0, 9, 4]
indices=[0, 9, 5, 1, 2, 3, 4, 6, 7, 8] cycles=[10, 1, 4] i=2, indices 2 and 6 were exchanged, output [0, 9, 5]
indices=[0, 9, 6, 1, 2, 3, 4, 5, 7, 8] cycles=[10, 1, 3] i=2, indices 2 and 7 were exchanged, output [0, 9, 6]
indices=[0, 9, 7, 1, 2, 3, 4, 5, 6, 8] cycles=[10, 1, 2] i=2, indices 2 and 8 were exchanged, output [0, 9, 7]
indices=[0, 9, 8, 1, 2, 3, 4, 5, 6, 7] cycles=[10, 1, 1] i=2, indices 2 and 9 were exchanged, output [0, 9, 8]
indices=[0, 9, 8, 1, 2, 3, 4, 5, 6, 7], shifting at 2:
        [      \1, 2, 3, 4, 5, 6, 7,/8]

indices=[0, 9, 1, 2, 3, 4, 5, 6, 7, 8], shifting at 1:
        [    \1, 2, 3, 4, 5, 6, 7, 8,/9]

indices=[1, 0, 2, 3, 4, 5, 6, 7, 8, 9] cycles=[9, 9, 8] i=0, indices 0 and 1 were exchanged, output [1, 0, 2]
indices=[1, 0, 3, 2, 4, 5, 6, 7, 8, 9] cycles=[9, 9, 7] i=2, indices 2 and 3 were exchanged, output [1, 0, 3]
indices=[1, 0, 4, 2, 3, 5, 6, 7, 8, 9] cycles=[9, 9, 6] i=2, indices 2 and 4 were exchanged, output [1, 0, 4]
indices=[1, 0, 5, 2, 3, 4, 6, 7, 8, 9] cycles=[9, 9, 5] i=2, indices 2 and 5 were exchanged, output [1, 0, 5]
indices=[1, 0, 6, 2, 3, 4, 5, 7, 8, 9] cycles=[9, 9, 4] i=2, indices 2 and 6 were exchanged, output [1, 0, 6]
indices=[1, 0, 7, 2, 3, 4, 5, 6, 8, 9] cycles=[9, 9, 3] i=2, indices 2 and 7 were exchanged, output [1, 0, 7]
indices=[1, 0, 8, 2, 3, 4, 5, 6, 7, 9] cycles=[9, 9, 2] i=2, indices 2 and 8 were exchanged, output [1, 0, 8]
indices=[1, 0, 9, 2, 3, 4, 5, 6, 7, 8] cycles=[9, 9, 1] i=2, indices 2 and 9 were exchanged, output [1, 0, 9]
indices=[1, 0, 9, 2, 3, 4, 5, 6, 7, 8], shifting at 2:
        [      \2, 3, 4, 5, 6, 7, 8,/9]

indices=[1, 2, 0, 3, 4, 5, 6, 7, 8, 9] cycles=[9, 8, 8] i=1, indices 1 and 2 were exchanged, output [1, 2, 0]
indices=[1, 2, 3, 0, 4, 5, 6, 7, 8, 9] cycles=[9, 8, 7] i=2, indices 2 and 3 were exchanged, output [1, 2, 3]
indices=[1, 2, 4, 0, 3, 5, 6, 7, 8, 9] cycles=[9, 8, 6] i=2, indices 2 and 4 were exchanged, output [1, 2, 4]
indices=[1, 2, 5, 0, 3, 4, 6, 7, 8, 9] cycles=[9, 8, 5] i=2, indices 2 and 5 were exchanged, output [1, 2, 5]
indices=[1, 2, 6, 0, 3, 4, 5, 7, 8, 9] cycles=[9, 8, 4] i=2, indices 2 and 6 were exchanged, output [1, 2, 6]
indices=[1, 2, 7, 0, 3, 4, 5, 6, 8, 9] cycles=[9, 8, 3] i=2, indices 2 and 7 were exchanged, output [1, 2, 7]
indices=[1, 2, 8, 0, 3, 4, 5, 6, 7, 9] cycles=[9, 8, 2] i=2, indices 2 and 8 were exchanged, output [1, 2, 8]
indices=[1, 2, 9, 0, 3, 4, 5, 6, 7, 8] cycles=[9, 8, 1] i=2, indices 2 and 9 were exchanged, output [1, 2, 9]
indices=[1, 2, 9, 0, 3, 4, 5, 6, 7, 8], shifting at 2:
        [      \0, 3, 4, 5, 6, 7, 8,/9]

indices=[1, 3, 0, 2, 4, 5, 6, 7, 8, 9] cycles=[9, 7, 8] i=1, indices 1 and 3 were exchanged, output [1, 3, 0]
indices=[1, 3, 2, 0, 4, 5, 6, 7, 8, 9] cycles=[9, 7, 7] i=2, indices 2 and 3 were exchanged, output [1, 3, 2]
indices=[1, 3, 4, 0, 2, 5, 6, 7, 8, 9] cycles=[9, 7, 6] i=2, indices 2 and 4 were exchanged, output [1, 3, 4]
indices=[1, 3, 5, 0, 2, 4, 6, 7, 8, 9] cycles=[9, 7, 5] i=2, indices 2 and 5 were exchanged, output [1, 3, 5]
indices=[1, 3, 6, 0, 2, 4, 5, 7, 8, 9] cycles=[9, 7, 4] i=2, indices 2 and 6 were exchanged, output [1, 3, 6]
indices=[1, 3, 7, 0, 2, 4, 5, 6, 8, 9] cycles=[9, 7, 3] i=2, indices 2 and 7 were exchanged, output [1, 3, 7]
indices=[1, 3, 8, 0, 2, 4, 5, 6, 7, 9] cycles=[9, 7, 2] i=2, indices 2 and 8 were exchanged, output [1, 3, 8]
indices=[1, 3, 9, 0, 2, 4, 5, 6, 7, 8] cycles=[9, 7, 1] i=2, indices 2 and 9 were exchanged, output [1, 3, 9]
indices=[1, 3, 9, 0, 2, 4, 5, 6, 7, 8], shifting at 2:
        [      \0, 2, 4, 5, 6, 7, 8,/9]

indices=[1, 4, 0, 2, 3, 5, 6, 7, 8, 9] cycles=[9, 6, 8] i=1, indices 1 and 4 were exchanged, output [1, 4, 0]
indices=[1, 4, 2, 0, 3, 5, 6, 7, 8, 9] cycles=[9, 6, 7] i=2, indices 2 and 3 were exchanged, output [1, 4, 2]
indices=[1, 4, 3, 0, 2, 5, 6, 7, 8, 9] cycles=[9, 6, 6] i=2, indices 2 and 4 were exchanged, output [1, 4, 3]
indices=[1, 4, 5, 0, 2, 3, 6, 7, 8, 9] cycles=[9, 6, 5] i=2, indices 2 and 5 were exchanged, output [1, 4, 5]
indices=[1, 4, 6, 0, 2, 3, 5, 7, 8, 9] cycles=[9, 6, 4] i=2, indices 2 and 6 were exchanged, output [1, 4, 6]
indices=[1, 4, 7, 0, 2, 3, 5, 6, 8, 9] cycles=[9, 6, 3] i=2, indices 2 and 7 were exchanged, output [1, 4, 7]
indices=[1, 4, 8, 0, 2, 3, 5, 6, 7, 9] cycles=[9, 6, 2] i=2, indices 2 and 8 were exchanged, output [1, 4, 8]
indices=[1, 4, 9, 0, 2, 3, 5, 6, 7, 8] cycles=[9, 6, 1] i=2, indices 2 and 9 were exchanged, output [1, 4, 9]
indices=[1, 4, 9, 0, 2, 3, 5, 6, 7, 8], shifting at 2:
        [      \0, 2, 3, 5, 6, 7, 8,/9]

indices=[1, 5, 0, 2, 3, 4, 6, 7, 8, 9] cycles=[9, 5, 8] i=1, indices 1 and 5 were exchanged, output [1, 5, 0]
indices=[1, 5, 2, 0, 3, 4, 6, 7, 8, 9] cycles=[9, 5, 7] i=2, indices 2 and 3 were exchanged, output [1, 5, 2]
indices=[1, 5, 3, 0, 2, 4, 6, 7, 8, 9] cycles=[9, 5, 6] i=2, indices 2 and 4 were exchanged, output [1, 5, 3]
indices=[1, 5, 4, 0, 2, 3, 6, 7, 8, 9] cycles=[9, 5, 5] i=2, indices 2 and 5 were exchanged, output [1, 5, 4]
indices=[1, 5, 6, 0, 2, 3, 4, 7, 8, 9] cycles=[9, 5, 4] i=2, indices 2 and 6 were exchanged, output [1, 5, 6]
indices=[1, 5, 7, 0, 2, 3, 4, 6, 8, 9] cycles=[9, 5, 3] i=2, indices 2 and 7 were exchanged, output [1, 5, 7]
indices=[1, 5, 8, 0, 2, 3, 4, 6, 7, 9] cycles=[9, 5, 2] i=2, indices 2 and 8 were exchanged, output [1, 5, 8]
indices=[1, 5, 9, 0, 2, 3, 4, 6, 7, 8] cycles=[9, 5, 1] i=2, indices 2 and 9 were exchanged, output [1, 5, 9]
indices=[1, 5, 9, 0, 2, 3, 4, 6, 7, 8], shifting at 2:
        [      \0, 2, 3, 4, 6, 7, 8,/9]

indices=[1, 6, 0, 2, 3, 4, 5, 7, 8, 9] cycles=[9, 4, 8] i=1, indices 1 and 6 were exchanged, output [1, 6, 0]
indices=[1, 6, 2, 0, 3, 4, 5, 7, 8, 9] cycles=[9, 4, 7] i=2, indices 2 and 3 were exchanged, output [1, 6, 2]
indices=[1, 6, 3, 0, 2, 4, 5, 7, 8, 9] cycles=[9, 4, 6] i=2, indices 2 and 4 were exchanged, output [1, 6, 3]
indices=[1, 6, 4, 0, 2, 3, 5, 7, 8, 9] cycles=[9, 4, 5] i=2, indices 2 and 5 were exchanged, output [1, 6, 4]
indices=[1, 6, 5, 0, 2, 3, 4, 7, 8, 9] cycles=[9, 4, 4] i=2, indices 2 and 6 were exchanged, output [1, 6, 5]
indices=[1, 6, 7, 0, 2, 3, 4, 5, 8, 9] cycles=[9, 4, 3] i=2, indices 2 and 7 were exchanged, output [1, 6, 7]
indices=[1, 6, 8, 0, 2, 3, 4, 5, 7, 9] cycles=[9, 4, 2] i=2, indices 2 and 8 were exchanged, output [1, 6, 8]
indices=[1, 6, 9, 0, 2, 3, 4, 5, 7, 8] cycles=[9, 4, 1] i=2, indices 2 and 9 were exchanged, output [1, 6, 9]
indices=[1, 6, 9, 0, 2, 3, 4, 5, 7, 8], shifting at 2:
        [      \0, 2, 3, 4, 5, 7, 8,/9]

indices=[1, 7, 0, 2, 3, 4, 5, 6, 8, 9] cycles=[9, 3, 8] i=1, indices 1 and 7 were exchanged, output [1, 7, 0]
indices=[1, 7, 2, 0, 3, 4, 5, 6, 8, 9] cycles=[9, 3, 7] i=2, indices 2 and 3 were exchanged, output [1, 7, 2]
indices=[1, 7, 3, 0, 2, 4, 5, 6, 8, 9] cycles=[9, 3, 6] i=2, indices 2 and 4 were exchanged, output [1, 7, 3]
indices=[1, 7, 4, 0, 2, 3, 5, 6, 8, 9] cycles=[9, 3, 5] i=2, indices 2 and 5 were exchanged, output [1, 7, 4]
indices=[1, 7, 5, 0, 2, 3, 4, 6, 8, 9] cycles=[9, 3, 4] i=2, indices 2 and 6 were exchanged, output [1, 7, 5]
indices=[1, 7, 6, 0, 2, 3, 4, 5, 8, 9] cycles=[9, 3, 3] i=2, indices 2 and 7 were exchanged, output [1, 7, 6]
indices=[1, 7, 8, 0, 2, 3, 4, 5, 6, 9] cycles=[9, 3, 2] i=2, indices 2 and 8 were exchanged, output [1, 7, 8]
indices=[1, 7, 9, 0, 2, 3, 4, 5, 6, 8] cycles=[9, 3, 1] i=2, indices 2 and 9 were exchanged, output [1, 7, 9]
indices=[1, 7, 9, 0, 2, 3, 4, 5, 6, 8], shifting at 2:
        [      \0, 2, 3, 4, 5, 6, 8,/9]

indices=[1, 8, 0, 2, 3, 4, 5, 6, 7, 9] cycles=[9, 2, 8] i=1, indices 1 and 8 were exchanged, output [1, 8, 0]
indices=[1, 8, 2, 0, 3, 4, 5, 6, 7, 9] cycles=[9, 2, 7] i=2, indices 2 and 3 were exchanged, output [1, 8, 2]
indices=[1, 8, 3, 0, 2, 4, 5, 6, 7, 9] cycles=[9, 2, 6] i=2, indices 2 and 4 were exchanged, output [1, 8, 3]
indices=[1, 8, 4, 0, 2, 3, 5, 6, 7, 9] cycles=[9, 2, 5] i=2, indices 2 and 5 were exchanged, output [1, 8, 4]
indices=[1, 8, 5, 0, 2, 3, 4, 6, 7, 9] cycles=[9, 2, 4] i=2, indices 2 and 6 were exchanged, output [1, 8, 5]
indices=[1, 8, 6, 0, 2, 3, 4, 5, 7, 9] cycles=[9, 2, 3] i=2, indices 2 and 7 were exchanged, output [1, 8, 6]
indices=[1, 8, 7, 0, 2, 3, 4, 5, 6, 9] cycles=[9, 2, 2] i=2, indices 2 and 8 were exchanged, output [1, 8, 7]
indices=[1, 8, 9, 0, 2, 3, 4, 5, 6, 7] cycles=[9, 2, 1] i=2, indices 2 and 9 were exchanged, output [1, 8, 9]
indices=[1, 8, 9, 0, 2, 3, 4, 5, 6, 7], shifting at 2:
        [      \0, 2, 3, 4, 5, 6, 7,/9]

indices=[1, 9, 0, 2, 3, 4, 5, 6, 7, 8] cycles=[9, 1, 8] i=1, indices 1 and 9 were exchanged, output [1, 9, 0]
indices=[1, 9, 2, 0, 3, 4, 5, 6, 7, 8] cycles=[9, 1, 7] i=2, indices 2 and 3 were exchanged, output [1, 9, 2]
indices=[1, 9, 3, 0, 2, 4, 5, 6, 7, 8] cycles=[9, 1, 6] i=2, indices 2 and 4 were exchanged, output [1, 9, 3]
indices=[1, 9, 4, 0, 2, 3, 5, 6, 7, 8] cycles=[9, 1, 5] i=2, indices 2 and 5 were exchanged, output [1, 9, 4]
indices=[1, 9, 5, 0, 2, 3, 4, 6, 7, 8] cycles=[9, 1, 4] i=2, indices 2 and 6 were exchanged, output [1, 9, 5]
indices=[1, 9, 6, 0, 2, 3, 4, 5, 7, 8] cycles=[9, 1, 3] i=2, indices 2 and 7 were exchanged, output [1, 9, 6]
indices=[1, 9, 7, 0, 2, 3, 4, 5, 6, 8] cycles=[9, 1, 2] i=2, indices 2 and 8 were exchanged, output [1, 9, 7]
indices=[1, 9, 8, 0, 2, 3, 4, 5, 6, 7] cycles=[9, 1, 1] i=2, indices 2 and 9 were exchanged, output [1, 9, 8]
indices=[1, 9, 8, 0, 2, 3, 4, 5, 6, 7], shifting at 2:
        [      \0, 2, 3, 4, 5, 6, 7,/8]

indices=[1, 9, 0, 2, 3, 4, 5, 6, 7, 8], shifting at 1:
        [    \0, 2, 3, 4, 5, 6, 7, 8,/9]

...

indices=[9, 8, 7, 0, 1, 2, 3, 4, 5, 6] cycles=[1, 1, 1] i=2, indices 2 and 9 were exchanged, output [9, 8, 7]
indices=[9, 8, 7, 0, 1, 2, 3, 4, 5, 6], shifting at 2:
        [      \0, 1, 2, 3, 4, 5, 6,/7]

indices=[9, 8, 0, 1, 2, 3, 4, 5, 6, 7], shifting at 1:
        [    \0, 1, 2, 3, 4, 5, 6, 7,/8]

indices=[9, 0, 1, 2, 3, 4, 5, 6, 7, 8], shifting at 0:
        [  \0, 1, 2, 3, 4, 5, 6, 7, 8,/9]

```

As to the details on how and why exactly this works - the answer is complicated and requires theoretical understanding of *permutation cycles* that I frankly don't have.
> 4 votes

# Answer
The way Raymond Hettinger's itertools permutations works is by following the cyclical nature of permutations. And converting from a recursive to iterative function.

## Establishing a recursive pattern

We can build a recursive function to show us the underlying pattern. We can build off of Yuval Filmus's recursive solution. However we should change the code to show all the values. So we can pass `L` and a list to pick to pick from.

```
def permutations(n, k):
    def inner(l, r):
        print(l, r)
        if len(l) == n:
            yield tuple(l)
        else:
            for v in r:
                _r = r[:]
                _r.remove(v)
                yield from inner(l + [v], _r)
    return inner([], list(range(k)))

```

If we run the above code we should be able to start seeing a pattern. This pattern is the same regardless of depth. I have grouped the output by depth and highlighted some numbers to make the pattern more obvious.

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

[⓪] [①, 2, 3, 4, 5, 6]
[①] [0, ②, 3, 4, 5, 6]
[②] [0, 1, ③, 4, 5, 6]
[③] [0, 1, 2, ④, 5, 6]
[④] [0, 1, 2, 3, ⑤, 6]
[⑤] [0, 1, 2, 3, 4, ⑥]
[6] [0, 1, 2, 3, 4, 5]

[0, ①] [②, 3, 4, 5, 6]
[0, ②] [1, ③, 4, 5, 6]
[0, ③] [1, 2, ④, 5, 6]
[0, ④] [1, 2, 3, ⑤, 6]
[0, ⑤] [1, 2, 3, 4, ⑥]
[0, 6] [1, 2, 3, 4, 5]

[1, ⓪] [②, 3, 4, 5, 6]
[1, ②] [0, ③, 4, 5, 6]
[1, ③] [0, 2, ④, 5, 6]
[1, ④] [0, 2, 3, ⑤, 6]
[1, ⑤] [0, 2, 3, 4, ⑥]
[1, 6] [0, 2, 3, 4, 5]

...

[0, 1, ②] [③, 4, 5, 6]
[0, 1, ③] [2, ④, 5, 6]
[0, 1, ④] [2, 3, ⑤, 6]
[0, 1, ⑤] [2, 3, 4, ⑥]
[0, 1, 6] [2, 3, 4, 5]

[0, 2, ①] [③, 4, 5, 6]
[0, 2, ③] [1, ④, 5, 6]
[0, 2, ④] [1, 3, ⑤, 6]
[0, 2, ⑤] [1, 3, 4, ⑥]
[0, 2, 6] [1, 3, 4, 5]

...

```

## Optimizing the recursive solution

From this we can see rather than building a new list we can just swap indexes. After we have swapped to the end of the 'chunk' we have to revert the list to how it was when we entered the function so the caller doesn't get messed up.

```
def permutations(n, k):
    def inner(l, r):
        i = len(l)
        if i == n:
            yield tuple(l)
        else:
            l.append(r.pop(0))
            yield from inner(l, r)
            for j in range(len(r)):
                l[-1], r[j] = r[j], l[-1]
                yield from inner(l, r)
            r.append(l.pop())
    return inner([], list(range(k)))

```

We can join `l` and `r` together to make a single list. If we always return the first `n` values then it doesn't matter if the two are combined.

```
def permutations(n, k):
    def inner(i):
        if i == n:
            yield pool[:n]
        else:
            yield from inner(i + 1)
            for j in reversed(range(1, cycles[i])):
                pool[i], pool[-j] = pool[-j], pool[i]
                yield from inner(i + 1)
            pool[i:] = pool[i + 1:] + pool[i:i + 1]
    pool = list(range(k))
    cycles = list(range(k, k-n, -1))
    return inner(0)

```

## Converting to an iterative algorithm

First we get n=1 working. This is really simple as it's just as simple for loop. Which is the body of the recursive function anyway.

```
def permutations(n, k):
    pool = list(range(k))
    cycles = list(range(k, k-n, -1))
    i = 0
    yield pool[:n]
    for j in reversed(range(1, cycles[i])):
        pool[i], pool[-j] = pool[-j], pool[i]
        yield pool[:n]
    pool[i:] = pool[i + 1:] + pool[i:i + 1]

```

Now the biggest challenge is getting the code to work with n=1 and n=2. To do this we can use a `while True` and a `for` loop. This is somewhat complicated so I'll show you how you could write a `itertools.product` function.

```
def product(n, k):
    pool = [0] * n
    yield tuple(pool)
    while True:
        for i in reversed(range(n)):
            pool[i] += 1
            if pool[i] != k:
                yield tuple(pool)
                break
            pool[i] = 0
        else:
            return

```

We can now convert permutations to work the same way for any n.

```
def permutations(n, k):
    pool = list(range(k))
    cycles = list(range(k, k - n, -1))
    r_range = range(n - 1, -1, -1)
    yield pool[:n]
    while True:
        for i in r_range:
            cycles[i] -= 1
            if cycles[i] > 0:
                j = cycles[i]
                pool[i], pool[-j] = pool[-j], pool[i]
                yield pool[:n]
                break
            pool[i:] = pool[i + 1:] + pool[i:i + 1]
            cycles[i] = k - i
        else:
            return

```
> 3 votes

# Answer
The problem of iterating all partial permutations is closely related to the problem of generating a uniform random sequence $U$ of $k$ (u)nique integers from $0 \le U\[j\] \< n$. The usual solution to that problem is first to generate a uniform (r)andom sequence of numbers $R$ satisfying

* $|R| = k$
* $0 \le R\[j\] \< n - j$

then covert $R$ to $U$. So for example, given $R$, the $U$ would be:

$$\begin{array} {c|c} R & U \\ \hline \[0, 0, 0\] & \[0, 1, 2\] \\ \[0, 0, 1\] & \[0, 1, 3\] \\ \[1, 1, 1\] & \[1, 2, 3\] \\ \[7, 0, 7\] & \[7, 0, 8\] \\ \end{array}$$

So if $R\[a\] = x$ then $U\[a\] = y$, where $y$ is the $x'th$ smallest natural number, skipping those in $U\[0 .. a-1\]$. Generating $U$ from $R$ is typically done in $O(k(\log k)^2)$ time. $U\[j\]$ is calculated by doing a (very complicated) binary search a sorted set of the first $U\[0..j-1\]$ already calculated elements.

A partial permutation iterator can be written by enumerating $R$ in lexicographical ordering, then converting $R$ to $U$ in the usual way.

It can also be done much better though, because the iterator makes updates to $R$, $U$, and the (s)orted set of $U$ (call it $S$) in a predictable way. Each next() call of the permutation iterator changes $R$ in the following way:

* $R\[0 .. z-1\]$ doesn't change
* $R\[z\]$ is incremented by 1
* $R\[z+1 .. k\]$ is set to 0

(for some index $z$). This means that $U\[0 .. z - 1\]$ doesn't change each iteration. So each iteration we can do the following steps:

* $S$ can be shrunk to $|S| = z-1$ by removing $U\[z .. k\]$ (from the previous iteration) from it.
* $U\[z\]$ can be calculated by looking at $R\[z\]$ and the current $S$
* Insert $U\[z\]$ into $S$
* $U\[z + 1..k\]$ can be calculated by enumerating a counter $c$ from zero, and if $c$ is not in $S$ then append it to the end of $U$, stop when $U$ is full.
* $U\[z + 1 .. k\]$ can be inserted into the sorted $S$ to make $|S| = k$ (for the sake of the next iteration of the iterator) using the same technique as merge from mergesort

Each of the above steps can be done in linear $O(k)$ time.

This algorithm has a lot of overhead, so if $k$ is predictably small, it will be faster to do something like

```
for (U[0]=0; U[0] < n; U[0]++) {
    for (U[1]=0; U[1] < n; U[1]++) {
        if (U[0] == U[1]) continue;
        // do your thing here
    }
 }

```

And if $n = k$ then a total permutation iterator will also be much faster. This $O(k)$ algorithm is almost entirely theoretical, since you'd need something like $k$ to be a small fraction of $n$, but large enough to not do the above forloop technique, for it to be ideal.
> 2 votes

---
Tags: algorithms, combinatorics, permutations
---


# Question
Title: Create an optimal to express all natural number as an arithmetic expression using the alphabet $\Sigma=\{(,),1,+,\times\}$
Hope you had a fantastic christmas break :)

I am trying to find an algorithm in polynomial time that finds the shortest arithmetical expression (the one with the least amount of 1 symbols) to express any natural number, while only using the elements of $\Sigma=\\{(,),1,+,\times\\}$.

As an example there are multiple ways to express $6\in \mathbb{N}$. but we want the one with the least amount of 1 symbols.

1+1+1+1+1+1

(1+1)\*(1+1+1)

(1+1)\*1+1+1+1+1

**Algorithm:**

My current algorithm GetExpression(n) uses a recursive mechanism to get the expression for a number $n\in \mathbb{N}$.

* if $n=1$ return "1"
* trying to find a divider for n (for loop from i=2 to $\lceil \sqrt{n} \rceil$ with trying $n\mod i=0$)
* if there is a divider return "(GetExpression(divider)$\times$GetExpression($\frac{n}{divider}$))"
* if n is prime return "(1+GetExpresion(n-1))"

**Runtime:**

Due to the loop trying to find the divider going in $\mathcal{O}(n)$ and both options of recursive calls having $0\<\alpha\_1=\frac{n-1}{n}\<1$ (in the prime case) and $0\<\alpha\_2+\alpha\_3\<1$ (in the non-prime case) we can use the master theorem to see, that our agorithm runs in $\mathcal{O}(n)$ (since we only make one of both options at a time we dont need to look for $\sum\_{i=1}^3\alpha\_i\<1$ but only for $\alpha\_1\<1$ and $\sum\_{i=2}^3\alpha\_i\<1$)

(actually we need to make a special case for $divider=2$ since there couldbe $\alpha\_2+\alpha\_3=\frac{1}{2}+\frac{2}{4}=1$ if $n=4$ but we don't care right now).

**Correctness:**

My problem is how to proof that this algorithm is correct. I was thinking about a induction but i have no glue how to build it.

# Answer
> 3 votes
You can't prove that your algorithm is correct because it isn't. First counter example is n = 46. If we don't bother writing 2 = 1+1 to make it more readable, your algorithm gives

```
46 = 2 * (1 + 2 * (1 + 2 * (1 + 2 * 2)))

```

which uses 13 1's.

A better solution is

```
46 = 1 + ((1 + 2) * ((1 + 2) * (1 + 2 * 2))

```

which uses only 12 1's.

The difference is bigger for n = 235 (difference 2), n = 649 (difference 3), n = 1,081 (difference = 4), n = 7,849 (difference = 5), n = 33,793 (difference = 6), n = 221,617 (difference = 8).

---
Tags: algorithms, regular-expressions, number-theory
---


# Question
Title: Why does Lamport clocks increment on both message received and sent?
According to wikipedia, Lamport clocks need to be incremented when sending a message (`time = time + 1`) and when receiving a message (`time = max(time_stamp, time) + 1`).

In my implementation, `Increment()` is used when sending and `Witness()` is used when receiving, each of them following the previous rules.

I can't understand why not only `Increment()` actually increment the counter, that is why `Witness()` does not simply do `time = max(time_stamp, time)` instead of `time = max(time_stamp, time) + 1`.

The source of my confusion is the following scenario:

* clock start at 1
* a message is sent, clock is at 2
* this message is seen locally, clock is now at 3

In this scenario, the clock incremented twice for the same message. Ideally it would just be one. Did I break something by witnessing the clock within the same process? Did I get something wrong?

# Answer
From Leslie Lamport *Time, Clocks, and the Ordering of Events in a Distributed System*:

> The space-time diagram of Figure 1 might then yield the picture in Figure 2. Condition C 1 means that there must be a tick line between any two events on a process line, and condition C2 means that every message line must cross a tick line.

I think my confusion was on what actually is an event.

In Lamport's paper and in Wikipedia, sending and receiving the messages are the events. In that case, incrementing the counter both when sending and receiving a message makes sense, to ensure that a clock tick separates those events.

In my case, the data structure modifications are the events. The message passing between the processes is irrelevant. In that case, I believe that incrementing the counter only when the data changes makes sense.

Now, if applying Lamport's clock principle this way still qualifies as Lamport's clock and holds all the related properties, I do not know.
> 0 votes

---
Tags: clocks
---


title: Finding all vertices on negative cycles

text:
# Question
Title: Finding all vertices on negative cycles
Given a weighted digraph, I can check whether a given vertex belongs to a negative cycle in $O(|V|\cdot|E|)$ using Bellman-Ford. But what if I need to find all vertices on negative cycles? Is there a way to do it faster than Floyd-Warshall's $O(|V|^3)$?

# Answer
> 6 votes
If you don't constrain yourself to simple cycles, you can actually use the Bellman-Ford algorithm to find all the relevant vertices.

Start by running a DFS on the graph to find its strongly connected components. Let them be $G\_1,G\_2,...G\_k$. For each $G\_i$, run BF on the subgraph. If the BF algorithms detects a negative cycle (which it can after $|V\_i|+1$ iterations where $V\_i$ is the vertices in $G\_i$), add $V\_i$ to your list of vertices that belong to a negative (not necessarily simple) cycle.

If one of the strongly connected components contain a negative cycle, then because $G\_i$ is finite and so are the weights, you can loop through the negative cycle as much as you need, and then add go in a cycle through all the vertices in $G\_i$. You only need to loop an amount that will guarantee a total negative cycle.

It is easy to show that there are no cycles between any $G\_i,G\_j$ and therefore there are no negative cycles other than inside the strongly connected components.

Since you run BF on each subgraph seperately, you keep the time complexity of BF which is $O(|V|\cdot|E|)$.

# Answer
> 4 votes
I was not able to turn up any better algorithm in my research.

In practice, you could improve the running time for many graphs by first decomposing the graph into strongly connected components, then running Floyd-Warshall on each strongly connected component. This does not improve the worst-case complexity (since in the worst case the entire graph could form one large strongly connected component), but it might help on many graphs you're likely to run into in practice.

# Answer
> 1 votes
**Update**

I just realized that **my answer below may be incorrect**. The intended idea was to perform |V|-1 waves of "relaxation", as in Bellman-Ford (see CLRS). At the end, instead of returning FALSE to denote "negative weight cycle found, so Bellman-Ford failed to find a shortest path tree", we just print the cycle, then continue looking for more cycles to print.

However, it is doubtful (or at least unclear) whether this will actually find *all* negative weight cycles (or at least simple cycles). For example, if the |V|-1 waves of relaxation only recorded tree edges for *extremely* negative cycles, but not for *moderately* negative cycles, then the latter may not be picked up by our print statements, since our print statements only detect cycles that are formed by tree edges.

**Original answer**

Getting the negative weight cycle is possible with a minor modification to Bellman-Ford.

The implementation of Bellman-Ford in CLRS instantly returns FALSE whenever a negative-weight cycle is found:

```
BELLMAN-FORD(G, w, s)
  INITIALIZE-SINGLE-SOURCE(G, s)
  for i = 1 to |G.V| - 1
    for each edge(u,v) in G.E
      RELAX(u,v,w)
  for each edge(u,v) in G.E
    if v.d > u.d + w(u,v)
      return FALSE
  return TRUE

```

Instead of returning FALSE, we look at the condition `v.d > u.d + w(u,v)`, which indicates even after |V|-1 passes through all edges, `v.d` can still be improved. This can only mean that there's a path from `v` to `u` through tree edges found during the RELAX loop, for which completing the path with edge `(u,v)` would cause a negative weight cycle. Hence, to print out the cycle, just follow `v` and its "parent" pointers all the way up to `u`. That path plus `(u,v)` is your negative weight cycle.

---
Tags: algorithms, graphs, shortest-path
---


# Question
Title: How do you read the regular expression (0^∗10^+)^+?
> Give an example of a string **in** the language of $(0^\*10^+)^+$.

I've been asked to give an example of a string in this language but I'm confused on how to read this notation. I'm guessing the acceptable strings are supposed to start with 0, but that's about all I can infer.

# Answer
As I see in the comments, I think you have a bit of confusion in the syntax and semantics of regular expressions.

* $^\*$ stands for Kleene star, which means "zero or more repetitions". (the string can be empty)
* $+$ stands for Kleene plus, which means "one or more repetitions". (the string cannot be empty)

So for your example $$(0^∗10^+)^+$$ you should first look inside the brackets:

1. $0^\*$ can be empty or $0$,$00$,$000$,$\dots$
2. $10^+$ can't be empty, so it only can be of the form $10,1010,101010,10101010,\dots$

So now we can combine them and we can take for example $001010$. But this is still around brackets and has the form kleene plus ($+$), so this can be $$(001010)^+=001010,001010001010,001010001010001010,\dots$$

So we can take the string $001010001010$ as an example for this regular expression.
> 0 votes

# Answer
This regular expression describes the language of all words over $\\{0,1\\}$ which do not consist entirely of $0$-s, and in which every $1$ is followed by $0$.
> 4 votes

---
Tags: regular-expressions
---


# Question
Title: Does strong typing contribute to better performance optimization by compiler?
Do guarantees provided by the type systems enable better performance optimizations by the *modern* compilers? I found a 20+ year article suggesting that it would. Also, intuitively speaking, knowing more about the code in advance should also be useful. However, a friend claims that modern JIT compilers make such claim irrelevant. In his words:

> This is a paper from 1998. Today, JIT compilers eliminate many of the cases mentioned in the paper for dynamically typed systems. For example, JavaScript uses integer arithmetic for integer types, without knowing the type declarations in advance. I doubt if predeclared types are as significant for performance today as they were two decades ago. Native instruction set and memory management mechanisms seem to have much more influence on performance to me.

What's the situation in modern compilers?

# Answer
I think your friend somewhat presents a false dichotomy.

I will just give one example: when it first came out, the Self VM was one of the fastest dynamic language implementations. In fact, the Smalltalk VM *written in Self* that shipped as part of the Self system was one of the fastest Smalltalk VMs of its time, despite being written in a dynamic language (much more dynamic than Smalltalk, Python, or Ruby) and running on top of another VM.

Even more, at the time it was first released, Self was one of the fastest OO language implementations in general, even competing with some C++ compilers of that time.

*However*, it achieved this performance precisely *because* of its aggressively optimizing JIT compiler aided by dynamic type inference and type feedback!

It is precisely *because* JIT compilers have more and better type information than AOT compilers that JIT-compiled language implementations can be so fast.

A JIT compiler can potentially have all the same information available that an AOT compiler has. But, a JIT compiler also has *runtime* information that an AOT compiler *cannot possibly* have because of all our favorite impossibility results such as the Halting Problem, Rice's Theorem, etc.

For example, Escape Analysis is equivalent to solving the Halting Problem. So, an AOT compiler cannot know in every circumstance whether a reference will escape the local scope or not, and thus must conservatively allocate it under the assumption that it will escape. A JIT compiler, however, can simply *assume* that the reference will never escape and compile the code accordingly, and then when it *observes* the reference escaping at runtime, it recompiles the code. (This is sometimes called *Escape Detection*.)

Likewise, Class Hierarchy Analysis is equivalent to solving the Halting Problem (in languages with dynamic code loading at least). So, an AOT compiler is limited in its ability to de-virtualize and thus inline potentially overridden methods. A JIT compiler, however, can simply *look* at the class hierarchy at runtime and see whether the method is overridden or not, and thus potentially inline it. And if it observes a piece of code loaded later that overrides the method after the fact, it can re-compile the code so that the method is no longer inlined but dynamically dispatched again.

So, yes, type information is helpful, and JIT compilers *do* use type information in their optimizations, including type information empirically collected at runtime that is not available to AOT compilers.

Your friend even uses an example of using type information in their own example! Modern ECMAScript execution engines indeed use type information about numeric values collected at runtime to optimize them even though semantically ECMAScript only has IEEE Std 754-2019 *binary64* as its one and only numeric type. (The recently added `bigint` is a distinct type with distinct numeric literals.) Typically, modern ECMAScript execution engines will segregate numeric values into at least doubles and 53 bit integers, possibly even more.

Similarly, they will detect the type of arrays and represent them as machine arrays (i.e. contiguous pieces of memory) even though semantically, an array is actually just a dictionary whose keys are the strings `"0"`, `"1"`, `"2"`, etc.

The two key differences between an AOT compiler and a JIT compiler are:

* The JIT compiler has more information available: it can theoretically gather all the same static information that an AOT compiler can, but it also has runtime information available that the AOT compiler cannot possibly have.
* On the flip side, the AOT compiler has "infinite" resources available (it can take as long to compile and use as much RAM as it wants), whereas the JIT compiler is competing with the user program for the same resources *while the user is waiting for their program to start*.

However, the kinds of information and the kinds of optimizations are mostly the same for both. Both like long stretches of straight-line code, both like type information, etc.

Personally, I believe that the best possible performance will come from an approach that has not yet been heavily explored commercially: use an AOT compiler to crunch out as much static information as possible, and perform heavy and expensive optimizations, but *keep all that information around* and *keep a rich representation of the program*. Then, hand this off to a JIT compiler which can use all of the expensively calculated information and the rich representation to even further optimize at runtime.

Real world systems typically work differently. For example, C++ compilers are highly sophisticated, by they typically produce very low-level "anemic" output such as x86 machine code which has no types and loses almost all of the information and intent of the original program. On the other hand, we have typical Java implementations which have sophisticated JIT compilers but very stupid practically non-optimizing AOT compilers.
> 7 votes

# Answer
If the language allows that the same variable x could have different types on different executions of a procedure, that doesn't mean the programmer actually wants to use this feature, or does use the feature. Typically you have three cases: 1. x has always the same type (but without proof). 2. x has the same type in 99.99% of all cases. 3. It happens quite regularly that x has different types.

In the first case, which I assume you are mostly interested in, most likely you will end up with code like "if x has type T then ... else "generate and execute code for another type". The performance impact will depend on how often the test has to be done. For example if you have a loop iterating for start \<= i \<= end, there's a good chance that the code looks like "if "start" has type int and "end" has type int then "run a loop assuming that "start", "end", "i" and variables derived from them all have type int" else "generate and execute code for other cases". The performance impact will be minimal.
> 2 votes

---
Tags: compilers, type-checking, program-optimization
---


# Question
Title: Vertices reachable from negative-weight cycles in Bellman-Ford
TLDR: I want to know if there's a simple way to fill in distances for all vertices *reachable* from negative weight cycles (not just ones on the cycle itself) once Bellman-Ford has found a negative-weight cycle.

The Bellman-Ford algorithm finds single-source shortest distances $d$ for a graph of $|V|$ vertices by making $|V|-1$ passes through all $|E|$ edges with weights $w$, and then checking each edge $u \rightarrow v$, with weight $w(u,v)$, for negative-weight cycles. If $v.d \> u.d + w(u,v)$, then we know that there exists a negative weight cycle consisting of some path from $v$ to $u$, then the edge $(u,v)$. It's straightforward to find every vertex **on** that cycle and set all their $d$ to $-\infty$.

However, this also must mean that any vertices **reachable** from but not **on** this cycle must also have $d = -\infty$. One way to update all their distances is to collapse each cycle to a vertex and then perform DFS from each such "cycle vertex", setting $x.d = -\infty$ for every reachable vertex $x$.

I'm wondering if there's a faster way to update all the distances, involving perhaps a smaller modification to Bellman-Ford?

# Answer
> 1 votes
You don't need to perform DFS separately from each such vertex. Instead, add a new fake vertex $s\_0$, with an edge from $s\_0$ to each vertex in a negative cycle; use DFS to find all vertices reachable from $s\_0$; and set their distances to $-\infty$. This needs only a single DFS, rather than many DFS's. It is simple and fast; there is no hope for something that is asymptotically faster, as it takes linear time, and in the worst case, you might need to visit all vertices, taking linear time in the worst case.

---
Tags: graphs, shortest-path, weighted-graphs, dag
---


# Question
Title: How to interpret parametric formulation of information bottleneck?
I'm reading this paper on latent representations with the information bottleneck https://arxiv.org/pdf/1804.06216.pdf and in section three, the authors write that the parametric formulation of the information bottleneck is

$\max\limits\_{\phi,\theta} - I\_\phi(t;x) + \lambda I\_{\phi,\theta}(t;y).$

My question is not about the information bottleneck entirely, but how to understand the first term $\max\limits\_{\phi,\theta}$ where $\phi$ and $\theta$ are parameters of a model. This notation doesn't make sense to me as I usually expect something to follow the $\max$. Does it just mean the largest value over all the parameters in the model?

# Answer
The notation stands for the maximum value of $-I\_{\phi}(t;x) + \lambda I\_{\phi,\theta}(t;y)$, where the maximum is taken over $\phi,\theta$. The result will be a function of $t,\lambda,x,y$.

Often in such a notation, the possible values of the variables being maximized over is restricted, and sometimes this forms part of the notation (say, $\max\_{0 \leq \phi,\theta \leq 1}$). In other times, the possible values are restricted, but this is not reflected in the notation, and should be gleaned from context.
> 1 votes

---
Tags: machine-learning, information-theory
---


title: Linear separability

text:
# Question
Title: Linear separability
I have an assignment in linear programming and this question was asked: given an input of two data sets $P\_1,P\_2$ containing points $(x,y)$, create a linear program that finds an equation of a line the separates between the two data sets.

Now I saw in Wikipedia the article about linear separability and I didn't understand what are the variables of $w\_1,w\_2,\ldots,w\_n$ stand for, and why does their algorithm work. If anyone could clarify this for me it would be great.

# Answer
The general equation of a line in the plane is $$ ax + by = c. $$ (The equation $y = ax + b$ cannot handle vertical lines.) This line separates the plane into two parts: those satisfying $ax+by \> c$, and those satisfying $ax+by \< c$. Your two datasets are linearly separable if there exist $a,b,c$ such that all points in $P\_1$ satisfy $ax+by \> c$, and all points in $P\_2$ satisfy $ax+by \< c$. (If all points in $P\_1$ satisfy $ax+by\<c$ and all points in $P\_2$ satisfy $ax+by\>c$, then we can negate $a,b,c$ to get to the situation above.)

In order to find $a,b,c$ using linear programming, we need to replace the strict inequalities $\>,\<$ with non-strict inequalities $\ge,\le$. The idea is the for small enough $\epsilon\>0$, the points in $P\_1$ would actually satisfy $ax+by \geq c+\epsilon$, and the points in $P\_2$ would satisfy $ax+by \leq c-\epsilon$. While we don't know $\epsilon$, if we multiply the two inequalities by $1/\epsilon$, we effectively make $\epsilon$ equal to $1$. Therefore we look for $a,b,c$ such that $ax+by \geq c+1$ for all points in $P\_1$, and $ax+by \leq c-1$ for all points in $P\_2$. This is a linear program, which can be solved using standard algorithms.
> 2 votes

---
Tags: algorithms, linear-programming
---


# Question
Title: number of ways evaluation of expression such that value not changed
one example:

> How many ways we can do possible value-preserving parenthesis the following expression in such a way that value not changed after parenthesis with one constraint that parenthesis among one variable is not correct (i.e: $(i)+j+k$ is not acceptable) ?

as an example for expression $i+j+k$, parenthesis like $i+(j+k)$, $((i+j)+k)$ are correct parenthesis but $(i)+j+k$ is not. I take a complete short example as follows.

# Answer
1. What's wrong with (i) + j + k ?
2. I live in a world where floating-point arithmetic has rounding errors, and where integer arithmetic has overflow checks, so (i + j) + k and i + (j + k) are most definitely not the same, for example if i is a large positive integer, and j, k are two large negative integers then the second one can produce an overflow (crash) while the first one doesn't.
> 1 votes

---
Tags: algorithms, data-structures, algorithm-analysis, combinatorics, permutations
---


# Question
Title: Trapdoor functions that depends only on NP-hard?
As far as I can tell, there are no examples of trapdoor functions which hardness assumptions is merely hardness of some NP-hard problem in worst case. That is, a trapdoor function in which, it is proved that, if there is an PPT algorithm that can invert it with non-negligible probability, then there is a PT algorithm which can solve the worst case of some NP-hard problem.

If I am wrong on the above, can you show me an example of such function?

If I am not, why is it so hard to make such a trapdoor function? Is there a concrete reason why it's hard to develop such function (beside "we haven't figure it out yet")?

This supposedly "duplicate" question (Why hasn't there been an encryption algorithm that is based on the known NP-Hard problems?) doesn't actually answer this question. That question isn't why worst case NP-hard isn't enough to make an encryption out of it. This question is essentially about why there had been no successful attempt at reducing average-case hardness down to worst-case NP-hard, for the purpose of making trapdoor function. There had been problems where average-case is just as hard as worst-case, so the fact that not all problems has worst-case=average-case doesn't mean you can't use different problems where that is true.

# Answer
> 1 votes
The crux of the problem lies in worst to average case reductions. The standard construction of a function which is average case hard is taking a function which is worst case hard and applying all kinds of fancy error correcting codes on its truth table. See Ta-Shma's notes sketching the STV construction. As far as I know, we don't know how to do that for NP (STV starts with a PSPACE hard function).

Even worse, we believe this can't be done for NP. See Bogdanov Trevisan, where it is shown that if an NP complete problem has a worst to average case reduction then the hierarchy collapses. Your problem is not exactly worst to average case reduction, but reducing an NP hard problem to inverting a one way function. However, these are related and Trevisan's paper also rules out such non-adaptive reductions.

More directly related to your question is AGGM, which also rules out adaptive reductions for the case where the size of the inverse image is easy to compute. You might also want to check out this paper. These results are 10-15 years old (STV is around 20), I have no idea if anything stronger was proven since.

---
Tags: complexity-theory, np-hard
---


# Question
Title: Difficulties proving when a language is in P or is NP-complete
I have some difficulties in understanding how to prove when a language is in P or is NP-Complete. Specifically, consider the following decision problems w.r.t undirected graphs $G = (V, E)$:

1. $L\_1 = \\{ \langle G = (V, E), k\rangle \mid \text{$G$ consists of connected} \\ \text{components of size at most $\log|V|$ each, and has a vertex cover of size $k$}\\}$
2. $L\_2 = \\{\langle G = (V, E), k \rangle \mid \text{$G$ consists of two connected} \\ \text{ components of size $\frac{|V|}{2}$ each, and has a vertex cover of size $k$} \\}$

I think that $L\_1$ is in $\text{P}$, and $L\_2$ is in $\text{NP-complete}$. Am I right? How can I prove that?

# Answer
> 2 votes
The language $L\_1$ is in P since you can brute force an optimal vertex cover in each connected component.

The language $L\_2$ is NP-hard by reduction from vertex cover. Roughly, given an arbitrary graph, we have to make it connected and to add another connected component, controlling the change of the minimum vertex cover throughout the process.

---
Tags: algorithms, complexity-theory, np-complete, np
---


# Question
Title: Maximum Independent Set of special Directed Graph
I was given this *special type* of Directed Graph and was asked to find it's Maximum Independent Set.

Graph Properties :

* Graph has $N$ vertices and $N$ edges
* There can be **no edge from a vertex $v$ to itself**
* Each vertex $v$ **has only 1 outgoing edge and 0 or more* incoming edges*

I framed this problem as a linked list with potentially some cycles. So far I've tried a DFS approach, two Greedy ones and right now I am thinking of a DP one trying to convert the linked list to a Binary Tree by removing the cycles somehow.

All of the approaches below respect the fact that **there can be no edges between vertices of the under construction independent set**.So when a vertex is under consideration to be a member of it I check if it points to a vertex of the under construction independent set.

First Greedy approach :

Sort the vertices based on the number of incoming edges and keep taking the ones with the least incoming edges let's call them $v\_1,v\_2...$ as well as the ones who point where $v\_1,v\_2...$ do as they are basically ***free*** and so on.

Second Greedy approach :

Take the ones with no incoming edges let's call them $z\_1,z\_2...$ as well as the ones who point where $z\_1,z\_2...$ do as they are basically ***free*** as well and from there see which vertices have the most incoming edges and take the ones who point at them.

Both approaches seem to fail by 1 in some test cases.

As for the DFS one it was something like take 2 a time where in $v\_1 \rightarrow v\_2 \rightarrow v\_3 \rightarrow v\_4 \rightarrow v\_5$ I would take $v\_1,v\_3,v\_5..$ and so on.

I also tried prioritizing the ones with 0 incoming edges in the DFS approach but no luck and I think I am in the right direction but I am missing a small detail.

Any thoughts or insight would be really helpful.

# Answer
> 3 votes
The type of graph you are describing is identical to the one encountered in Pollard's rho algorithm. Connected components of your graph are cycles with directed trees feeding into them (possibly none). It follows that it suffices to solve your problem for such graphs.

The first step is to compute, for each directed tree, the maximum size of an independent set, depending on whether the root is included or not. Denote the maximum size of an independent set in the first case by $\alpha\_+(T)$, and in the second case by $\alpha\_-(T)$, and let $\alpha(T) = \max(\alpha\_+(T),\alpha\_-(T))$. If $T$ is a single vertex then $\alpha\_+(T) = 1$ and $\alpha\_-(T) = 0$. If $T$ is composed of a root and subtrees $T\_1,\ldots,T\_d$, then $\alpha\_+(T) = \alpha\_-(T\_1) + \cdots + \alpha\_-(T\_d) + 1$, and $\alpha\_-(T) = \alpha(T\_1) + \cdots + \alpha(T\_d)$.

Now suppose that we have a cycle of length $\ell$, and let $T\_1,\ldots,T\_\ell$ be the directed trees whose roots are the vertices of the cycle (the trees could be single vertices). Let $\alpha\_{\pm\_1,\pm\_2,i}$ be the maximum size of an independent set in $T\_1,\ldots,T\_i$ which contains the root of $T\_1$ according to the first $\pm$, the root of $T\_i$ according to the second $\pm$, and which doesn't contain two adjacent roots. First, $\alpha\_{\pm,\pm,1} = \alpha\_\pm(T\_1)$, and $\alpha\_{\pm,\mp,1} = 0$. Second, $\alpha\_{\pm,+,i+1} = \alpha\_{\pm,-,i} + \alpha\_+(T\_{i+1})$ and $\alpha\_{\pm,-,i+1} = \max(\alpha\_{\pm,+,i},\alpha\_{\pm,-,i}) + \alpha\_-(T\_{i+1})$. Finally, the maximum size of an independent set of the entire component is $\max(\alpha\_{+,-,\ell},\alpha\_{-,+,\ell},\alpha\_{-,-,\ell})$.

---
Tags: graphs, greedy-algorithms, depth-first-search, minmax-algorithm
---


# Question
Title: Determine if a language is Decidable or semi decidable
Consider the language $L = \\{\langle M \rangle: \text{ $M$ accepts at most two single-letter words}\\}$, where $\langle M\rangle$ is the encoding of Turing machine $M$.

We need to determine, without using Rice's Theorem, whether:

1. $L\in \text{R}$, i.e., $L$ is decidable.
2. $L \in \text{RE}\setminus \text{R}$, i.e., $L$ is semi-decidable but not decidable
3. $L \notin \text{RE}$, i.e., $L$ is not semi-decidable.

Any hints or some steps on how to approach this problem would be appreciated.

Also, if $L\notin \text{R}$, can we prove that using Rice's Theorem?

# Answer
The language $L$ is not in $\text{RE}$. Consider the complement of the halting problem $\overline{Halt\_{TM}} = \\{ \langle M, w\rangle: \text{ $M$ does not halt on $w$}\\}$ which is not in $\text{RE}$. We show that $\overline{Halt\_{TM}}$ reduces to $L$ (and therefore, $L \notin \text{RE}$). From this point, I am assuming w.l.o.g that every machine has an input alphabet that contains at least three letters $a, b, $ and $c$.

A mapping reduction from $\overline{Halt\_{TM}}$ to $L$ operates as follows. On input $\langle M, w\rangle$ of the reduction, the reduction outputs $\langle K\rangle$, where $K$ is a machine that operates as follows. On input $x$ of $K$, if $x = a$ or $x = b$, $K$ accepts. Otherwise, if $x\notin \\{ a, b\\}$, then $K$ simulates the run of $M$ on $w$. If $M$ halts on $w$, then $K$ accepts $x$. The reduction is correct. Indeed, if $M$ does not halt on $w$, then $K$ accepts the inputs $x = a$, and $x = b$, and does not halt on any other input. In particular, $K$ accepts at most two single-letter words. Conversely, if $M$ Halts on $w$, then $K$ accepts all inputs $x$, in particular, $K$ accepts at least three single letter words $a, b,$ and $c$.

Since $\text{R}\subseteq \text{RE}$, the above also proves that $L$ is not decidable. Yet, if you want to prove that $L$ is not decidable using Rice's theorem, then you note that the property $\text{P} = \\{ M: \text{ $M$ accepts at most two single-letter words} \\}$ is a non-trivial semantic property:

1. $\text{P}$ is nontrivial: it is easy to show that there are machines that accept at most two single-letter words, and there are machines that accept at least three single-letter words. I leave the details to you.
2. $\text{P}$ is semantic: let $M\_1$ and $M\_2$ be two machines with $L(M\_1) = L(M\_2)$. We need to show that both $M\_1$ and $M\_2$ are in $\text{P}$, or both are not in $\text{P}$. This is also easy as $M\_1$ and $M\_1$ recognise the same language, in particular, they accept the same single-letter words.
> 0 votes

---
Tags: turing-machines, undecidability, semi-decidability, rice-theorem
---


# Question
Title: How to derive the worst case time complexity of Heapify algorithm?
I would like to know how to derive the time complexity for the *Heapify Algorithm for Heap Data Structure*.

I am asking this question in the light of the book "Fundamentals of Computer Algorithms" by Ellis Horowitz et al. I am adding some screenshots of the algorithm as well as the derivation given in the book.

Algorithm for Heapify():

**procedure** $HEAPIFY(A,n)$
//Readjust the elements in A(1:n) to form a heap//
**integer** $n,i$
**for** $i\leftarrow\lfloor n/2 \rfloor$ **to** $1$ **by** $-1$ **do**
**call** $ADJUST(A, i, n)$
**repeat**
**end** $HEAPIFY$

Derivation for worst case complexity:

I understood the first part and last part of this calculation but I cannot figure out how $2^{i-1}(k-i)$ changed into $i 2^{k-i-1}$.

All the derivations I can find in the internet takes a different approach by considering the height of the tree differently. I know that approach also leads to the same answer but I would like to know about this approach.

You might need the following information:

$2^k-1 = n$ or approximately $2^k = n$, where $k$ is the number of levels, starting from the root node and the level of root is 1 (not 0) and $n$ is the number of nodes.

Also the worst case time complexity of the Adjust() function is proportional to the height of the sub-tree it is called, that is $O(log n)$, where $n$ is the total number of elements in the sub-tree.

# Answer
> 2 votes
It's a variable substitution.

First, realize that in the leftmost side of the equation, the last term of the sum is zero (because when $i = k$, $k-i = 0$). So, the range of the first summation can be written as $1 \le i \le k-1$. Now, substitute $i$ with $k-i$. $i$ iterates over the set ${1, 2, ... , k-1}$ and $k-i$ iterates over the set ${k-1, ... 2, 1}$, they are the same set, so, we can do this.

---
Tags: algorithm-analysis, heaps, mathematical-analysis, heap-sort, derivation
---


# Question
Title: Mixed Integer Linear Programming application problem
I have a Mixed-Integer Linear Programming question.

There's this game called Islanders. Pretty graphics, it's about maximizing your score building a small city in a confined island.

Here's a quick rundown of the rules:

1. There are various building types
2. Every building type has a different size
3. Every building type has a different range (explained later)
4. Every building may give a non-zero base score when placed
5. Every building may have give a positive, zero or negative score based on what buildings are within its range when placed
6. Buildings already placed are not affected by newer buildings (their score contribution is not updated)

I thought it would be a fun exercise to formulate these rules into an MILP and see what comes out of it.

I have written the attached Python 3 program, using the PuLP library (link). I have used the default solver (CBC), provided with the installation of PuLP.

There are currently 4 building types encoded, with their base score declared at line 61, the 4x4 cross-score matrix at line 64 and their dimension and range in lines 72 and 75.

In line 82 there is currently stated that there is 1 building of the first type and 3 buildings of the 2nd type to be placed optimally (N=5)

There are several decision variables:

1. `px` and `py` are each building coordinates
2. `t` is the build order (integer)
3. `placed_before` is a NxN array saying if building i is to be placed before building j (binary)
4. `east_of`, `north_of` have similar logic (binary)
5. `covers` states whether building j is in range of building i (binary)
6. `x_overlaps` and `y_overlaps` check whether two buildings are built on top of each other, which is restricted (binary)
7. `gives_score` is a NxN array saying if building i contributes to the score of building j (is whthin range is is placed before it, binary)

The constraints are as follows:

1. Lines 140-147 state the building order constraints, which form variable `t`
2. 149-161 decide if buildings are east-west and north-east, used for later
3. 163-183 establish whether buildings overlaps each other, in the east-west (x) direction and the north-south (y) direction
4. 185-219 establish if building i covers building j, filling in the `covers` variable
5. 222-237 enforces the separation of two buildings, if they overlap in both x and y directions
6. 239-250 finds out if a building i contributes to the score of building j

Currently, the cost function under optimization is `gives_score`, to maximize the number of buildings which give score to another.

The issue:

* Currently with 1 of the first building and 4 of the second, the optimization does not finish in reasonable time. 380 constraints are generated.
* If instead you go for 3 of the second building type, the problem becomes tractable again, finishing in approx. 10s. 228 constraints are generated.
* Also, if you disable the constraints of lines 222-237, the optimization manages to run, but the buildings overlap.

Any suggestions on the problem formulation or explanations on why this is such a hard problem?

Thank you in advance.

**Bonus:** Here is an example optimal result for 1 of type-1 building and 3 of type-2:

# Answer
> 1 votes
While I can't answer if your model is okay or if it has some subtle flaw, I think I can confirm that general approach works and the issue is with CBC solver.

I reimplemented the game on my own: https://gist.github.com/afish/56d710cf440bf652fa2ea51f356f7378

CPLEX finds solution for city center and 4 houses in 3 seconds. CBC can't do it (I killed it after couple minutes). I imagine the same situation might be with your code.

I don't know what solvers you can use but seems like CBC is not a good choice here. Just for the record, I also tried MIPCL, SCIP, and GLPK - none of them could solve the problem in a reasonable (~seconds) time.

---
Tags: integer-programming, python
---


# Question
Title: Algorithm to find diameter of a tree using BFS/DFS. Why does it work?
This link provides an algorithm for finding the diameter of an undirected tree **using BFS/DFS**. Summarizing:

> Run BFS on any node s in the graph, remembering the node u discovered last. Run BFS from u remembering the node v discovered last. d(u,v) is the diameter of the tree.

Why does it work ?

Page 2 of this provides a reasoning, but it is confusing. I am quoting the initial portion of the proof:

> Run BFS on any node s in the graph, remembering the node u discovered last. Run BFS from u remembering the node v discovered last. d(u,v) is the diameter of the tree.
>
> Correctness: Let a and b be any two nodes such that d(a,b) is the diameter of the tree. There is a unique path from a to b. Let t be the first node on that path discovered by BFS. If the paths $p\_1$ from s to u and $p\_2$ from a to b do not share edges, then the path from t to u includes s. So
>
> $d(t,u) \ge d(s,u)$
>
> $d(t,u) \ge d(s,a)$
>
> ....(more inequalities follow ..)

The inequalities do not make sense to me.

# Answer
All parts of proving the claim hinge on 2 crucial properties of trees with undirected edges:

* 1-connectedness (ie. between any 2 nodes in a tree there is exactly one path)
* any node can serve as the root of the tree.

Choose an arbitrary tree node $s$. Assume $u, v \in V(G)$ are nodes with $d(u,v) = diam(G)$. Assume further that the algorithm finds a node $x$ starting at $s$ first, some node $y$ starting at $x$ next. wlog $d(s,u) \geq d(s,v)$. note that $d(s,x) \geq d(s,y)$ must hold, unless the algorithm's first stage wouldn't end up at $x$. We will see that $d(x,y) = d(u,v)$.

The most general configuration of all nodes involved can be seen in the following pseudo-graphics ( possibly $s = z\_{uv}$ or $s = z\_{xy}$ or both ):

```
(u)                                            (x)
  \                                            /
   \                                          /
    \                                        /
     ( z_uv )---------( s )----------( z_xy )
    /                                        \
   /                                          \
  /                                            \
(v)                                            (y)

```

we know that:

1. $d(z\_{uv},y) \leq d(z\_{uv},v)$. otherwise $d(u,v) \< diam(G)$ contradicting the assumption.
2. $d(z\_{uv},x) \leq d(z\_{uv},u)$. otherwise $d(u,v) \< diam(G)$ contradicting the assumption.
3. $d(s,z\_{xy}) + d(z\_{xy},x) \geq d(s,z\_{uv}) + d(z\_{uv},u)$, otherwise stage 1 of the algorithm wouldn't have stopped at $x$.
4. $d(z\_{xy},y) \geq d(v,z\_{uv}) + d(z\_{uv},z\_{xy})$, otherwise stage 2 of the algorithm wouldn't have stopped at $y$.

1) and 2) imply $\\, \\ d(u,v) = d(z\_{uv},v) + d(z\_{uv},u) \\ \qquad\geq d(z\_{uv},x) + d(z\_{uv},y) = d(x,y) + 2\\, d(z\_{uv}, z\_{xy}) \\ \qquad\qquad\geq d(x,y)$.

3) and 4) imply $\\, \\ d(z\_{xy},y) + d(s,z\_{xy}) + d(z\_{xy},x) \\ \qquad\geq d(s,z\_{uv}) + d(z\_{uv},u) + d(v,z\_{uv}) + d(z\_{uv},z\_{xy}) \qquad\qquad\qquad\qquad \\ \\, $ equivalent to $\\, \\ d(x,y) = d(z\_{xy},y) + d(z\_{xy},x) \\ \qquad\geq 2\*\\,d(s,z\_{uv}) + d(v,z\_{uv}) + d(u,z\_{uv}) \\ \qquad\qquad\geq d(u,v)$.

therefore $d(u,v) = d(x,y)$.

analogue proofs hold for the alternative configurations

```
                 (u)                          (x)
                   \                          /
                    \                        /
                     \                      /
     ( s )---------( z_uv )----------( z_xy )
                     /                      \
                    /                        \
                   /                          \
                 (v)                          (y)

```

and

```
                          (x)        (u)
                          /            \
                         /              \
                        /                \
     ( s )---------( z_xy )----------( z_uv )
                        \                /
                         \              /
                          \            /
                          (y)        (v)

```

these are all possible configurations. in particular, $x \not\in path(s,u), x \not\in path(s,v)$ due to the result of stage 1 of the algorithm and $y \not\in path(x,u), y \not\in path(x,v)$ due to stage 2.
> 14 votes

# Answer
The intuition behind is very easy to understand. Suppose I have to find longest path that exists between any two nodes in the given tree.

After drawing some diagrams we can observe that the longest path will always occur between two leaf nodes( nodes having only one edge linked). This can also be proved by contradiction that if longest path is between two nodes and either or both of two nodes is not a leaf node then we can extend the path to get a longer path.

So one way is to first check what nodes are leaf nodes, then start BFS from one of the leaf node to get the node farthest from it.

Instead of first finding which nodes are leaf nodes , we start BFS from a random node and then see which node is farthest from it. Let the node farthest be x. It is clear that x is a leaf node. Now if we start BFS from x and check farthest node from it, we will get our answer.

But what is the guarantee that x will be a end point of a maximum path?

Let's see by an example :-

```
       1
    / /\ \
   6 2  4 8
         \ \
          5 9
           \
            7

```

Suppose I started my BFS from 6. The node at maximum distance from 6 is node 7. Using BFS we can get this node. Now we star BFS from node 7 to get node 9 at maximum distance. Path from node 7 to node 9 is clearly the longest path.

What if BFS that started from node 6 gave 2 as the node at maximum distance. Then when we will BFS from 2 we will get 7 as node at maximum distance and longest path will be then 2-\>1-\>4-\>5-\>7 with length 4. But the actual longest path length is 5. This cannot happen because BFS from node 6 will never give node 2 as node at maximum distance.

Hope that helps.
> 14 votes

# Answer
**Update 3 and corrected answer**

There's an error in the linked solution set (see update 2 below), but it can be easily corrected with @Yuval Filmus's suggestion in the question's comment, which further allows us to rule out one of the possibilities mentioned in the solution set.

This proof hinges on two facts: one, that there is exactly one path between any two nodes in a binary tree; and two, $d(x,y)=d(y,x)$.

Let $s$ be the root of the first BFS and let $u$ be the final vertex discovered (which means that $u$ must be the farthest point from $s$). Also, let $d(a,b)$ be a diameter. Clearly, $a$, $b$, and $u$ must all be leaves, since otherwise $d(a,b)$ wouldn't be a diameter, and $u$ wouldn't be discovered last (if any of the nodes had children, they would no longer be endpoints of paths of "maximal length", as you could just extend the path by exploring their children). Let $t$ be the lowest common ancestor of $a$ and $b$.

We claim that $t$ must be an ancestor of $u$. Suppose that it were not, then let $r$ be the lowest common ancestor of $u$ and $t$.

$d(u,s)=d(u,r)+d(r,s)$, so we have $d(u,s) \geq d(u,r)$.

Since $u$ was discovered last, $d(u,s) \geq d(a,s)$. Since $d(a,s)=d(a,r)+d(r,s)$, and $d(u,s)=d(u,r)+d(r,s)$, therefore, $d(u,r) \geq d(a,r)$.

Since $r$ is an ancestor of $t$, and $r \neq t$ (since $t$ is not an ancestor of $u$, yet $r$ is), $d(t,r) \geq 1$. Since $d(a,r)=d(a,t)+d(t,r)$, $d(a,r) \> d(a,t)$.

Thus, since $d(u,r) \> d(a,t)$, $d(u,r) + d(r,t) + d(t,b) \> d(a,t) + d(t,b)$ $\to d(u,b) \> d(a,b)$, contradicting the fact that $d(a,b)$ is a diameter (since diameters are supposed to be paths of maximal length).

Hence, $t$ is an ancestor of $u$.

Since $u$ was discovered last, $d(s,u) \geq d(s,a)$. Since $d(s,u)=d(s,t)+d(t,u)$, and $d(s,a)=d(s,t)+d(t,a)$, $d(t,u) \geq d(t,a)$. Since $d(t,u) \geq d(t,a)$, $d(b,t)+d(t,u) \geq d(b,t)+d(t,a) \to d(b,u) \geq d(b,a)$.

Since $d(a,b)$ is a diameter, $d(b,a) \geq d(b,u)$.

Therefore, since $d(b,u) \geq d(b,a)$ and $d(b,u) \leq d(b,a)$, $d(u,b) = d(a,b)$.

Therefore, $u$ is the endpoint of *some* diameter, and thus the second BFS works (the longest path from an endpoint of a diameter must be a diameter).

**Update 2**

Note that the linked solution set contains an error:

> If the paths $p\_1$ from $s$ to $u$ and $p\_2$ from $a$ to $b$ do not share edges, then the path from $t$ to $u$ includes $s$.

Consider the following counterexample: suppose $s$ is the root of the first BFS tree, and $t$ is the lowest common ancestor of leaves $a$, $u$, and $b$. No edges are shared; only the vertex $t$ is shared. Let $d(a,t) = d(b,t) = d(u,t) = 2$, and let $d(s,t)=1$. Finally, suppose $u$ is discovered last in the BFS. Then $d(a,b)$ is a diameter, $t$ is the first node discovered on that path, and the path from $t$ to $u$ does *not* include $s$.

To fix this, it's probably necessary to change "do not share edges" to "do not share vertices, as suggested by @Yuval Filmus.

**Update**

As @j\_random\_hacker points out in the comments, Lemma 1 below is not sufficient to show that $u$ is an endpoint of some diameter. Hence, the below proof is incomplete.

**Original incorrect answer**

Suppose we have $d(a, b)$ being a diameter, but we don't know which vertices $a$ and $b$ are, so we start a BFS from some vertex $s$.

If $s = a$, then the first BFS would yield either $b$ or some node $b'$ equally distant from $a$, and the second would either go back to $a$ or to some equally distant node $a'$, and the scheme obviously works. Similarly, if $s = b$, then two BFS would work for the same reason.

Otherwise, we have $s \neq a, b$.

**Lemma 0**: Both $a$ and $b$ are leaf nodes in the tree rooted at $s$.

**Proof 0**: If they weren't leaf nodes, we could increase $d(a,b)$ by extending the endpoints to leaf nodes, contradicting $d(a, b)$ being a diameter.

**Lemma 1**: At least one of $d(s,a)$ and $d(s, b)$ is the largest possible value of $d(s,u)$ for all $u$.

**Proof 1**: Suppose we have a $u$ that violates the lemma. $u$ cannot be a descendent of $a$ or $b$ since they are leaves from Lemma 0, and $u$ cannot be an ancestor, since that would make it *closer* to $s$, and it would no longer be able to violate the lemma. Let $t$ be the lowest common ancestor of $a$ and $b$. If $t = s$, then both $d(a, u)$ and $d(b, u)$ would be greater than $d(a, b)$, a contradiction. Otherwise, $t \neq s$ and either $u$ is a descendant of $t$, or not.

If $u$ is a descendant of $t$, then since $d(s,u)$ is greater than both $d(s, a)$ and $d(s, b)$, we have $d(t,u)$ is greater than both $d(t, a)$ and $d(t, b)$ and thus $d(a, u)$ and $d(b,u)$ are both greater than $d(a,b)$, a contradiction since $d(a,b)$ is a diameter.

If $u$ is not a descendant of $t$, then let $w$ be the unique ancestor of $u$ such that $d(s,t) = d(s,w)$. We know $w$ exists because $d(s,t) \< d(s,u)$. Then, we have $d(t, b) \< d(w,u) \< d(t,w) + d(w,u) = d(t,u)$ and thus $d(a, u) = d(a, t) + d(t, u) \> d(a, t) + d(t, b) = d(a, b)$, again contradicting $d(a, b)$ being a diameter.

**Main result**: Starting from any root $s$ and performing a BFS will result in some $a$ being discovered last. Using that $a$ as the root of a second $BFS$ will result in $b$ being discovered last, with $d(a,b)$ guaranteed to be a diameter.

**proof**: From lemma 1, the first BFS will find an $a$ that is furthest from $s$, which is guaranteed to be one of the endpoints of a diameter. The second BFS will find that diameter, since in a tree there's exactly one simple path between any two nodes.
> 7 votes

# Answer
It works because-

1. For any randomly chosen vertex $u$, the farthest vertex from it lies on one end of the diameter.
2. If vertex $a$ is one end of the diameter, the farthest vertex from it is the other end.

The first BFS takes you to one end of the diameter and the second BFS from that end takes you to the other end.

**Proof**

**#2** is easy to see why it's true. We will see why **#1** is true by contradiction. Let's suppose, the path from vertex $a$ to $b$ is the diameter of the tree and $v$, which is not on the diameter, is the farthest vertex from $u$. Also, let's assume paths $(a,b)$ and $(u,v)$ cross over at $t$ ($t$ can more generally be a path with multiple vertices). Because we don't have any cycles, all paths between vertices $a$, $b$, $u$ and $v$ pass through $t$.

1. By our assumptions $d(a,b)$ is the diameter and $d(u,v)$ is larger than both $d(u,a)$ and $d(u,b)$.
2. Path $(u,t)$ is the common piece in all paths originating from $u$ and hence we can deduce that $d(v, t) \> d(a, t)$ ($v$ is farther away from $t$ as compared to $a$).
3. This implies that $v$ is also farther away from $b$ as compared to $a$ and $d(v,b) \> d(a,b)$. *(add $d(t,b)$ to both sides in previous step)*
4. This contradicts our initial assumption that $(a,b)$ is a diameter because we have found a path $(v,b)$ which is longer. Symmetrically, we can also see why $d(v,a) \> d(a,b)$.
5. Hence, the farthest vertex from $u$ must either be $a$ or $b$ if $(a,b)$ is indeed the diameter.
> 3 votes

# Answer
The claim is equivalent to claiming that every deepest leaf in a tree rooted at $s$ is an endpoint to a diameter.

Suppose $p$ is an arbitrary path with end points $(a,b)$, and that $d$ is some leaf that is one of the deepest. Let $x$ be the deepest common ancestor of $a$ and $b$.

Here is how to construct a path ending in $d$ that it at least as long as path $p$, assuming WLOG $\text{depth}(a) \ge \text{depth}(b)$:

```
   x
  / \       (b,d) is at least as long
 / \ \
a   d b

   x
  / \       (a,d) is at least as long
 / / \
a d   b

   /\
  x  \      (a,d) is at least as long
 / \  \
a   b  d

```

Therefore a diameter can be constructed for every vertex $d$ from any other diameter.

(Edited: the previous proof only established that every diameter ended in a deepest vertex, not that every deepest vertex was part of a diameter, embarrassing mistake.)
> 2 votes

# Answer
By the definition of BFS, the distance (from the starting node) of each node explored is either equal to the distance of the previous node explored or greater by 1. Thus, the last node explored by BFS will be among those farthest from the starting node.

Thus, the algorithm of using BFS twice amounts to "Pick an arbitrary node $x$. Find the node $a$ farthest from $x$ (last node found by BFS starting from $x$). Find the node $b$ farthest from $a$ (last node found by BFS starting from $a$).", which thus finds two nodes of maximum distance from eachother.
> 0 votes

# Answer
One key thing to know is that a tree is always planar, which means it can be laid out on a plane, so often ordinary 2-dimensional thinking works. In this case, the algorithm says starting anywhere, go as far away as you can. The distance from that point to as far as you can get away from that point is the longest distance in the tree, and therefore the diameter.

This method would also work to find the diameter of a real, physical island if we defined that as the diameter of the smallest circle that would fully enclose the island.
> 0 votes

# Answer
@op, the way the cases are defined in the PDF may be a bit off.

I think that the two cases should be:

1. $p\_1$ does not intersect with $p\_2$, i.e. there are no common vertices between paths $p\_1$ and $p\_2$. In this case, define $t$ as the first node on $p\_2$ discovered by the first BFS starting at $s$.
2. $p\_1$ and $p\_2$ have at least one common vertex. In this case, define $t$ to be the first node on $p\_2$ discovered by the first BFS that is also on $p\_1$.

The rest of the proof in the PDF should follow.

With this definition, the figure shown by OP falls into Case 2.
> 0 votes

# Answer
First run a DFS from a random node then the diameter of a tree is the path between the deepest leaves of a node in its DFS subtree:
> -1 votes

---
Tags: algorithms, graphs, search-algorithms, trees, graph-traversal
---


# Question
Title: λ -terms that correspond to Haskell functions
Evening All

I already have a "grasp" of haskell (not terrible, about 6 month experience) and am trying to learn the fundamentals that sit behind it, thus am now turning my attention to trying to learn Lambda Calculus. I thought, perhaps mistakenly, that it might be a good idea to try to work out the λ -terms that corresponds to some functions in Haskell. I am trying the foldr and map functions.

soooo foldr should have the property

```
              foldr f u [N1, . . . , Nk] →f N1 (f N2 (. . . (f Nk u)))

```

and map would have (should imagine...)

```
              map f [N1, . . . , Nk] → [f N1, f N2, . . . , f Nk]

```

where → represents reduction by any number of β -steps....

But am struggling to flesh this out to be honest. Any help would be greatly appreciated...

Kind Regards.

# Answer
You are looking for Church encodings of datastructures in $\lambda$-calculus, and in particular of lists.
> 1 votes

---
Tags: lambda-calculus, haskell
---


# Question
Title: Is there a formalization for the trade-off between the expressiveness and the speed of a programming language?
I do not have a strong background in theoretical computer science, so the terminology I use may not be the most ideal.

Basically the "trend" across programming languages I have noticed is that there appears to be an trade-off in the speed of a programming language (how long it takes to execute a program) and its expressiveness (saying you want something done vs saying you *how* you want something done, more expressive languages focus less on the "how", from this definition).

For example, in "higher-level" languages like Python or Java - it is very easy to do something like finding the index of an item in a list with a single line like `list.indexOf(item)`. In comparison, doing the same action in a "lower-level" language like C can result in code that is much faster, but at the cost of having write the code yourself which may be several times as long.

Again, this is really just a trend I have noticed and not a hard and fast rule. Languages like C++ that are both very fast yet can also be pretty expressive at times blur the lines.

So I am curious, is there any sort of formalization in programming language theory - or even just rough rules of thumb backed by empirical studies - explicitly listing out the trade off between the expressiveness and execution speed of a programming language?

Or perhaps I am just overthinking it, and having slower code is just a unavoidable byproduct of having a more expressive language as more things are taken care of "behind the scenes" / abstracted away from during compilation, interpretation, etc.?

# Answer
> 1 votes
This is a *really* interesting problem and I don't think there are definitive answers. There are two important, and distinct dimensions to this problem.

* What does expressivity actually mean? In the light of the Church-Turing thesis, which states that all programming languages are the same (can be compiled into each other). Nevertheless it seems pretty clear that e.g. assembler programming is quite different from programming in Python. There is a research field that studies this. It was pioneered by (1), but much ink has been spilled since.
* In computational complexity theory speedup theorems show that for any algorithm solving a problem, there is a faster algorithm solving the same problem. Examples are *Blum's speedup theorem* and the *linear speedup theorem*. Basically, what those speedup theorems tell you is that you can always speed up computers (Turing machines) a bit by adding more powerful commands. In processors we went from 4 bit CPUs to 8 bit, then 16 bit, the 32 bit and nowadays 64 bit.

Neither quite answers your question, I'm afraid!

---

(1) M. Felleisen, *On the Expressive Power of Programming Languages.*

---
Tags: programming-languages
---


# Question
Title: NP-hardness for one-dimensional facility location problem with entrance fee for each customer
We have $n$ customers, $(x\_1, \dots, x\_n)$, sorted on the read line. For convenience, we also use $x\_i$ to denote its coordinate on the line. We need to locate $m$ facilities on the real line. **We note that these facilities can be located anywhere on the real line**. Each facility $j\in \[m\]$ is associated with an entrance fee $f\_j$, such as the ticket money for a swimming pool. Given a location profile $(y\_1, \dots, y\_m)\in \mathcal{R}^m$ for the facilities($(y\_1, \dots, y\_m)$ is not necessarily ordered), the cost for customer $i$ at facility $j$ is $c\_{ij} = |x\_i- y\_j|+f\_j$, which can be understood as the aggregate money of the travelling cost and the entrance fee when you taxi to a swimming pool and purchase a ticket to swim inside. And the customer will always choose the facility so as to minimize her cost.

Our goal is to find a location profile $(y\_1, \dots, y\_m)$ for the given $m$ facilities such that the total minimum cost

$$\sum\_{i\in \[n\]}\min\_{j\in \[m\]}c\_{ij}$$ is minimized.

**An easy but critical obserbation is that there is an optimal solution where each facility is located in the median agent of the continuous region of agents it serves**. If we have identical entrance fees, we only need to find the optimal $m$ continuous partitions, which can be solved by dynamic programming in $O(mn^2)$. For the general case, similar DP algortihm runs in $O(2^mn^2)$. The algorithm, however, is exponential in $m$ and only makes sense after we've proved the problem is NP-hard. Another observation which may be helpful is that for if we have known the optimal $m$ continuous partitions, we just assign the facility of the $k$-th smallest entrance fee to the partition of the $k$-th most customers.

So is this problem NP-hard? Or is there an algorithm running in polynomial time?

# Answer
> 0 votes
$\require{cancel}$

The answer shows a reduction to the rainbow path problem and thus doesn't prove anything about the complexity of the given problem. See this paper which contains hardness results for related problems (theorem 2.3 states that rainbow s-t path is hard).

The problem is $\cancel{\text{still polynomial}}$, as it can still be reduced to some version of shortest path. Observe that there always exists an optimal solution where the facilities are a subset of $\\{x\_1,...,x\_n\\}$ (given an optimal solution and a facility $y\_j$ in it, let $A\_j$ denote the set of clients attending $y\_j$, then replacing $y\_j$ with the median of $A\_j$ does not increase the value of the solution). Now construct a directed graph whose vertices are $V=\\{s,t, v\_{il}| i\in\[n\], l\in\[m\] \\}$ and whose edges are given by:

$$E=\\{(v\_{il},v\_{jk})| i\<j\in\[n\], \text{and } l,k\in\[m\]\\}\cup\\{(s,v\_{il})|i\in\[n\],l\in\[m\]\\}\\\cup\\{(v\_{n,l},t) |l\in\[m\]\\}$$

Interpret the edge $(v\_{il},v\_{jk})$ as choosing a solution where the clients $x\_{i+1},...,x\_j$ are members of the facility $y\_k$. With this in mind, we set the weights of the edges to be:

$$ w(v\_{il},v\_{jk})=\sum\limits\_{q=i+1}^j |x\_q-m\_{ij}|+(j-i)f\_k\\ w(s,v\_{il})=\sum\limits\_{q=1}^i|x\_q-m\_{0,i}|+if\_l\\ w(v\_{n,l},t)=0, \text{ for all $l\in\[m\]$} $$

Where $m\_{ij}$ is the median of $x\_{i+1},...,x\_j$. An optimal solution for your problem is the weight of the lightest path from $s$ to $t$ of length $\le m$ which uses at most one edge from the set $E\_i=\left\\{(v,v\_{l,i})|v\in V, l\in \[n\]\right\\}$ for every $i\in\[m\]$ (this implies that we do not use the same facility twice). This version of shortest (lightest) path can be solved in polynomial time (think of it as having $m$ different color edges, and we want a path that does not use the same color twice, I believe that a sort of modified Dijkstra should work).

# Answer
> 0 votes
This is not intended to be a full answer but rather to give a possible direction.

Assume $f\_j$ are natural numbers. Then let $K = \sum\_j f\_j$.

Solve the following clustering problem for points $p\_k$: $$ \arg\min\_k \sum\_i \min\_k |x\_i - p\_k| $$

Then we need to assign each cluster a facility, such that facility $j$ is associated with $f\_j$ clusters. Clearly, in order to minimize our objective these clusters need to be adjacent to each other.

Then, to determine facility $j$, we need to select the median point of clusters $c\_j^1, c\_j^2, \dots, c\_j^{f\_j}$.

Threats to validity: while $x\_i$ may have belonged to a cluster, it is not guaranteed that it will be associated with the facility that that cluster is assigned to.

---
Tags: algorithms, complexity-theory, np-hard, polynomial-time-reductions
---


# Question
Title: Flow graph with non zero lower bound or 0 capacity
I am afraid the question title might not be sufficiently accurate but I could not come up with something more accurate

Here is the problem

Given 'n' machines

* Each machine has a set of capabilities
* Each machine has max availability (A(m))

Given 'm' tasks

* Each task requires a set of capabilities
* Each task takes a certain time (D(t))
* A task has to be performed on one machine only

The problem is to determine whether all tasks can be completed.

I get stuck with the 'one machine only' requirement. The only flow graphs I can come up with do not guarantee a task is not linked to more than one machine.

It's sort of a bipartite matching problem but with capacities \> 1

I also ran into XOR-like behavior in flow networks which is similar but has the 'xor' requirement on the 'source' end where I would need it on the target end.

Would anyone have any tips? Is it at all possible to model this as a flow graph?

Tx!

Peter

PS Trying to clarify requirements with more concrete example

Assume digital print systems and print jobs

* Each digital press can run for a number of hours
* Each press has some finishing possibilities: e.g. 'sheet cutter', 'laminate', 'laser cutter', 'page folding', ....
* Each print job requires a number of hours
* Each print job needs one or more of the finishing possibilities

Given a set of machines, the availability for each and the finishing possibilities and also a set of print jobs (duration, finishing options needed) determine whether all print jobs can finish

So e.g.

* Printer p1 is available for a duration of 10 hours and has features f1 and f2
* Printer p2 is available for a duration of 10 hours and has features f2 and f3
* Job1, requiring features f1 and f2 takes 8 hours to run
* Job2, requiring features f2 and f3 takes 8 hours to run
* Job3, requiring feature f2 requires 4 hours to run

A printer that is available for e.g. 10 hours can run 10 x 1 hour jobs or 5 x 2 hour jobs, or 1 x 8 hr job and 1 x 2 hr job, ... Jobs always have to run on a single printer

The flow diagrams I could come up with always result in cases where

8 hours of p1 is assigned to job1 (leaving 2 hours for printer p1) 8 hours of p2 is assigned to job2 (leaving 2 hours for printer p2)

(so far so good)

but then

The 2 hours of p1 and p2 left are used to flow to job3 and the max flow seems to indicate the three jobs can be run (which is not ok)

# Answer
Your problem is NP-hard. In the special case where no job requires any particular features, and all printers have the same availability, this becomes just the bin packing problem, which is (strongly) NP-complete.

You could try adapting standard algorithms for bin packing to your situation. For instance, one approach would be to use integer linear programming and hope the ILP solver can handle the resulting problem instance.
> 0 votes

---
Tags: network-flow, bipartite-matching, ford-fulkerson
---