In [24]:
from functools import lru_cache
import operator
from itertools import count, zip_longest, chain
from collections import OrderedDict
from random import randint
In [2]:
def memo_holder(optimizer):
def f(*args, **kwds):
return_memo_table = kwds.pop('memo_table', False)
pair = optimized, memo_table = optimizer(*args, **kwds)
return pair if return_memo_table else optimized
return f
Have a look at https://en.wikipedia.org/wiki/Monge_array
In [3]:
def parity_numbered_rows(matrix, parity, include_index=False):
start = 0 if parity == 'even' else 1
return [(i, r) if include_index else r
for i in range(start, len(matrix), 2)
for r in [matrix[i]]]
def argmin(iterable, only_index=True):
index, minimum = index_min = min(enumerate(iterable), key=operator.itemgetter(1))
return index if only_index else index_min
def interleaving(one, another):
for o, a in zip_longest(one, another):
yield o
if a: yield a
def is_sorted(iterable, pred=lambda l, g: l <= g):
_, *rest = iterable
return all(pred(l, g) for l, g in zip(iterable, rest))
def minima_indexes(matrix):
if len(matrix) == 1: return [argmin(matrix.pop())]
recursion = minima_indexes(parity_numbered_rows(matrix, parity='even'))
even_minima = OrderedDict((i, m) for i, m in zip(count(start=0, step=2), recursion))
odd_minima = [argmin(odd_r[start:end]) + start
for o, odd_r in parity_numbered_rows(matrix, parity='odd', include_index=True)
for start in [even_minima[o-1]]
for end in [even_minima[o+1]+1 if o+1 in even_minima else None]]
return list(interleaving(even_minima.values(), odd_minima))
def minima(matrix):
return [matrix[i][m] for i, m in enumerate(minima_indexes(matrix))]
def is_not_monge(matrix):
return any(any(matrix[r][m] > matrix[r][i] for i in range(m))
for r, m in enumerate(minima_indexes(matrix)))
The following is a Monge array:
In [70]:
matrix = [
[10, 17, 13, 28, 23],
[17, 22, 16, 29, 23],
[24, 28, 22, 34, 24],
[11, 13, 6, 17, 7],
[45, 44, 32, 37, 23],
[36, 33, 19, 21, 6],
[75, 66, 51, 53, 34],
]
minima(matrix)
Out[70]:
In [71]:
minima_indexes(matrix)
Out[71]:
In [72]:
is_not_monge(matrix)
Out[72]:
The following is not a Monge array:
In [73]:
matrix = [
[37, 23, 22, 32],
[21, 6, 7, 10],
[53, 34, 30, 31],
[32, 13, 9, 6],
[43, 21, 15, 8],
]
minima(matrix) # produces a wrong answer!!!
Out[73]:
In [74]:
minima_indexes(matrix)
Out[74]:
In [75]:
is_not_monge(matrix)
Out[75]:
In [35]:
@memo_holder
def longest_increasing_subsequence(seq):
L = []
for i, current in enumerate(seq):
"""opt, arg = max([(l, j) for (l, j) in L[:i] if l[-1] < current],
key=lambda p: len(p[0]),
default=([], tuple()))
L.append(opt + [current], (arg, i))"""
L.append(max(filter(lambda prefix: prefix[-1] < current, L[:i]), key=len, default=[]) + [current])
return max(L, key=len), L
def lis_rec(seq):
@lru_cache(maxsize=None)
def rec(i):
current = seq[i]
return max([rec(j) for j in range(i) if seq[j] < current], key=len, default=[]) + [current]
return max([rec(i) for i, _ in enumerate(seq)], key=len)
a simple test case taken from page 157:
In [36]:
seq = [5, 2, 8, 6, 3, 6, 9, 7] # see page 157
In [37]:
subseq, memo_table = longest_increasing_subsequence(seq, memo_table=True)
In [38]:
subseq
Out[38]:
memoization table shows that [2,3,6,7] is another solution:
In [39]:
memo_table
Out[39]:
In [40]:
lis_rec(seq)
Out[40]:
The following is an average case where the sequence is generated randomly:
In [41]:
length = int(5e3)
seq = [randint(0, length) for _ in range(length)]
In [42]:
%timeit longest_increasing_subsequence(seq)
In [43]:
%timeit lis_rec(seq)
worst scenario where the sequence is a sorted list, in increasing order:
In [44]:
seq = range(length)
In [45]:
%timeit longest_increasing_subsequence(seq)
In [46]:
%timeit lis_rec(seq)
In [26]:
@memo_holder
def edit_distance(xs, ys,
gap_in_xs=lambda y: 1, # cost of putting a gap in `xs` when reading `y`
gap_in_ys=lambda x: 1, # cost of putting a gap in `ys` when reading `x`
mismatch=lambda x, y: 1, # cost of mismatch (x, y) in the sense of `==`
gap = '▢',
mark=lambda s: s.swapcase(),
reduce=sum):
T = {}
T.update({ (i, 0):(xs[:i], gap * i, i) for i in range(len(xs)+1) })
T.update({ (0, j):( gap * j,ys[:j], j) for j in range(len(ys)+1) })
def combine(w, z):
a, b, c = zip(w, z)
return ''.join(a), ''.join(b), reduce(c)
for i, x in enumerate(xs, start=1):
for j, y in enumerate(ys, start=1):
T[i, j] = min(combine(T[i-1, j], (x, gap, gap_in_ys(x))),
combine(T[i, j-1], (gap, y, gap_in_xs(y))),
combine(T[i-1, j-1], (x, y, 0) if x == y else (mark(x), mark(y), mismatch(x, y))),
key=lambda t: t[2])
return T[len(xs), len(ys)], T
In [21]:
(xs, ys, cost), memo_table = edit_distance('exponential', 'polynomial', memo_table=True)
In [22]:
print('edit with cost {}:\n\n{}\n{}'.format(cost, xs, ys))
In [23]:
memo_table
Out[23]:
In [49]:
(xs, ys, cost), memo_table = edit_distance('exponential', 'polynomial', memo_table=True, mismatch=lambda x,y: 10)
In [50]:
print('edit with cost {}:\n\n{}\n{}'.format(cost, xs, ys))
In [51]:
memo_table
Out[51]:
In [49]:
@memo_holder
def matrix_product_ordering(o, w, op):
n = len(o)
T = {(i,i):(lambda: o[i], 0) for i in range(n)}
def combine(i,r,j):
t_ir, c_ir = T[i, r]
t_rj, c_rj = T[r+1, j]
return (lambda: op(t_ir(), t_rj()), c_ir + c_rj + w(i, r+1, j+1)) # w[i]*w[r+1]*w[j+1])
for d in range(1, n):
for i in range(n-d):
j = i + d
T[i, j] = min([combine(i, r, j) for r in range(i, j)], key=lambda t: t[-1])
opt, cost = T[0, n-1]
return (opt(), cost), T
def parens_str_proxy(w, **kwds):
return matrix_product_ordering(o=['▢']*(len(w)-1),
w=lambda l, c, r: w[l]*w[c]*w[r],
op=lambda a, b: '({} {})'.format(a, b),
**kwds)
In [50]:
(opt, cost), memo_table = parens_str_proxy(w={0:100, 1:20, 2:1000, 3:2, 4:50}, memo_table=True)
In [51]:
opt, cost
Out[51]:
In [52]:
{k:(thunk(), cost) for k,(thunk, cost) in memo_table.items()}
Out[52]:
In [53]:
(opt, cost), memo_table = parens_str_proxy(w={i:(i+1) for i in range(10)}, memo_table=True)
In [54]:
opt, cost
Out[54]:
In [55]:
{k:(thunk(), cost) for k,(thunk, cost) in memo_table.items()}
Out[55]:
In [56]:
from sympy import fibonacci, Matrix, init_printing
init_printing()
In [57]:
(opt, cost), memo_table = parens_str_proxy(w={i:fibonacci(i+1) for i in range(10)}, memo_table=True)
In [58]:
opt, cost
Out[58]:
In [59]:
{k:(thunk(), cost) for k,(thunk, cost) in memo_table.items()}
Out[59]:
In [60]:
def to_matrix_cost(dim, memo_table):
n, m = dim
return Matrix(n, m, lambda n,k: memo_table.get((n, k), (lambda: None, 0))[-1])
In [61]:
to_matrix_cost(dim=(9,9), memo_table=memo_table)
Out[61]:
In [211]:
@memo_holder
def longest_common_subsequence(A, B, gap_A, gap_B,
equal=lambda a, b: 1, shrink_A=lambda a: 0, shrink_B=lambda b: 0,
reduce=sum):
T = {}
T.update({(i, 0):([gap_A]*i, 0) for i in range(len(A)+1)})
T.update({(0, j):([gap_B]*j, 0) for j in range(len(B)+1)})
def combine(w, z):
alpha, beta = zip(w, z)
return list(chain.from_iterable(alpha)), reduce(beta)
for i, a in enumerate(A, start=1):
for j, b in enumerate(B, start=1):
T[i, j] = combine(T[i-1, j-1], ([a], equal(a, b))) if a == b else max(
combine(T[i, j-1], ([gap_B], -shrink_B(b))),
combine(T[i-1, j], ([gap_A], -shrink_A(a))),
key=lambda t: t[-1])
opt, cost = T[len(A), len(B)]
return (opt, cost), T
def pprint_memo_table(T, joiner, do=str):
return {k:(joiner.join(map(do, v[0])), v[1]) for k, v in T.items()}
In [224]:
(opt, cost), memo_table = longest_common_subsequence(A='ADCAAB',
B='BAABDCDCAACACBA', gap_A='▢', gap_B='○',
#shrink_B=lambda b:1,
memo_table=True)
In [225]:
print('BAABDCDCAACACBA')
print(''.join(opt))
In [226]:
pprint_memo_table(memo_table, joiner='')
Out[226]:
In [215]:
(opt, cost), memo_table = longest_common_subsequence(A=[0,1,1,2,3,5,8,13,21,34,55],
B=[1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796, 58786, 208012, 742900, 2674440,],
gap_A='▢', gap_B='○',
#shrink_B=lambda b:1,
memo_table=True)
In [216]:
print(','.join(map(str, opt)))
In [217]:
pprint_memo_table(memo_table, joiner=',')
Out[217]:
This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.