In [1]:
import sys
import numpy as np
import random
import itertools
import heapq as hq
In [2]:
sys.path.append('src/')
In [3]:
#import pyximport; pyximport.install()
from inference_lv import do_inference_list_viterbi
In [4]:
#from inference import do_inference_brute_force
In [5]:
random.seed(0)
np.random.seed(0)
In [6]:
class HeapItem: # an item in heapq (min-heap)
def __init__(self, priority, task):
self.priority = priority
self.task = task
self.string = str(priority) + ': ' + str(task)
def __lt__(self, other):
return self.priority < other.priority
def __repr__(self):
return self.string
def __str__(self):
return self.string
Brute force search.
In [7]:
def brute_force(ps, L, M, unary_params, pw_params, unary_features, pw_features, debug=True, top=5):
Cu = np.zeros(M, dtype=np.float) # unary_param[p] x unary_features[p]
Cp = np.zeros((M, M), dtype=np.float) # pw_param[pi, pj] x pw_features[pi, pj]
# a intermediate POI should NOT be the start POI, NO self-loops
for pi in range(M):
Cu[pi] = np.dot(unary_params[pi, :], unary_features[pi, :])
for pj in range(M):
Cp[pi, pj] = -np.inf if (pj == ps or pi == pj) else np.dot(pw_params[pi, pj, :], pw_features[pi, pj, :])
Q = []
poi_set = [p for p in range(M) if p != ps]
for x in itertools.product(poi_set, repeat=L-1):
y = [ps] + list(x)
score = 0
for j in range(1, L):
score += Cp[y[j - 1], y[j]] + Cu[y[j]]
if len(Q) < top:
hq.heappush(Q, HeapItem(score, np.array(y)))
else:
hq.heappushpop(Q, HeapItem(score, np.array(y))) # pop the smallest, then push
results = []
scores = []
while len(Q) > 0:
hterm = hq.heappop(Q)
results.append(hterm.task)
scores.append(hterm.priority)
# reverse the order: smallest -> largest => largest -> smallest
results.reverse()
scores.reverse()
if debug is True:
for score, y in zip(scores, results):
print(y, score)
return results
The list Viterbi algorithm described in paper Sequentially finding the N-best list in hidden Markov models (2001).
In [36]:
def list_viterbi_2001(ps, L, M, unary_params, pw_params, unary_features, pw_features, top=10):
assert(L > 1)
assert(M >= L)
assert(ps >= 0)
assert(ps < M)
assert(top > 0)
maxIter = 1e6
Cu = np.zeros(M, dtype=np.float)
Cp = np.zeros((M, M), dtype=np.float)
for pi in range(M):
Cu[pi] = np.dot(unary_params[pi, :], unary_features[pi, :])
for pj in range(M):
Cp[pi, pj] = -np.inf if (pj == ps or pi == pj) else np.dot(pw_params[pi, pj, :], pw_features[pi, pj, :])
Alpha = np.zeros((L, M), dtype=np.float)
Beta = np.zeros((L, M), dtype=np.float)
for pj in range(M): Alpha[1, pj] = Cp[ps, pj] + Cu[pj]
for t in range(2, L):
for pj in range(M):
Alpha[t, pj] = np.max([Alpha[t-1, pi] + Cp[pi, pj] + Cu[pj] for pi in range(M)])
for pi in range(M): Beta[L-1, pi] = 0
for t in range(L-1, 1, -1):
for pi in range(M):
Beta[t-1, pi] = np.max([Cp[pi, pj] + Cu[pj] + Beta[t, pj] for pj in range(M)])
Beta[0, ps] = np.max([Cp[ps, pj] + Cu[pj] + Beta[1, pj] for pj in range(M)])
Fp = np.zeros((L-1, M, M), dtype=np.float)
for t in range(L-1):
for pi in range(M):
for pj in range(M):
Fp[t, pi, pj] = Alpha[t, pi] + Cp[pi, pj] + Cu[pj] + Beta[t+1, pj]
y_best = np.ones(L, dtype=np.int) * (-1)
y_best[0] = ps
for t in range(1, L): y_best[t] = np.argmax(Fp[t-1, y_best[t-1], :])
Q = []
priority = -np.max(Alpha[L-1, :])
print(-priority)
print(Alpha[L-1, y_best[L-1]])
print(Fp[L-2, y_best[L-2], y_best[L-1]])
partition_index = -1
exclude_set = set()
hq.heappush(Q, HeapItem(priority, (y_best, partition_index, exclude_set)))
results = []
k = 0; y_last = None
while len(Q) > 0 and k < maxIter:
hitem = hq.heappop(Q)
k_priority = hitem.priority
(k_best, k_partition_index, k_exclude_set) = hitem.task
k += 1; y_last = k_best
print('OUT: %s, %.6f' % (k_best, -k_priority))
if len(set(k_best)) == L:
results.append(k_best); top -= 1
if top == 0: return results
partition_index_start = 1
if k_partition_index > 0:
assert(k_partition_index < L)
partition_index_start = k_partition_index
for parix in range(partition_index_start, L):
new_best = np.zeros(L, dtype=np.int) * (-1)
new_best[:parix] = k_best[:parix]
new_exclude_set = set({k_best[parix]})
if parix == partition_index_start: new_exclude_set = new_exclude_set | k_exclude_set
candidate_points = [p for p in range(M) if p not in new_exclude_set]
if len(candidate_points) == 0: continue
candidate_maxix = np.argmax([Fp[parix-1, k_best[parix-1], p] for p in candidate_points])
new_best[parix] = candidate_points[candidate_maxix]
for pk in range(parix+1, L):
new_best[pk] = np.argmax([Fp[pk-1, new_best[pk-1], p] for p in range(M)])
new_priority = (-k_priority)
new_priority += Fp[parix-1, k_best[parix-1], new_best[parix]] - Fp[parix-1, k_best[parix-1], k_best[parix]]
new_priority *= -1.0
hq.heappush(Q, HeapItem(new_priority, (new_best, parix, new_exclude_set)))
#print('IN : %s, %.6f' % (new_best, -new_priority))
if len(Q) == 0:
sys.stderr.write('WARN: empty queue, return the last one\n')
results.append(y_last); top -= 1
while len(Q) > 0 and top > 0:
hitem = hq.heappop(Q)
results.append(hitem.task[0]); top -= 1
return results
The list Viterbi algorithm described in paper List Viterbi decoding algorithms with applications (1994).
In [38]:
def list_viterbi_1994(ps, L, M, unary_params, pw_params, unary_features, pw_features, top=10):
assert(L > 1)
assert(M >= L)
assert(ps >= 0)
assert(ps < M)
assert(top > 0)
maxIter = 1e6
Cu = np.zeros(M, dtype=np.float)
Cp = np.zeros((M, M), dtype=np.float)
for pi in range(M):
Cu[pi] = np.dot(unary_params[pi, :], unary_features[pi, :])
for pj in range(M):
Cp[pi, pj] = -np.inf if (pj == ps or pi == pj) else np.dot(pw_params[pi, pj, :], pw_features[pi, pj, :])
Alpha = np.zeros((L, M), dtype=np.float)
Beta = np.zeros((L, M), dtype=np.float)
for pj in range(M): Alpha[1, pj] = Cp[ps, pj] + Cu[pj]
for t in range(2, L):
for pj in range(M):
Alpha[t, pj] = np.max([Alpha[t-1, pi] + Cp[pi, pj] + Cu[pj] for pi in range(M)])
for pi in range(M): Beta[L-1, pi] = 0
for t in range(L-1, 1, -1):
for pi in range(M):
Beta[t-1, pi] = np.max([Cp[pi, pj] + Cu[pj] + Beta[t, pj] for pj in range(M)])
Beta[0, ps] = np.max([Cp[ps, pj] + Cu[pj] + Beta[1, pj] for pj in range(M)])
Fp = np.zeros((L-1, M, M), dtype=np.float)
for t in range(L-1):
for pi in range(M):
for pj in range(M):
Fp[t, pi, pj] = Alpha[t, pi] + Cp[pi, pj] + Cu[pj] + Beta[t+1, pj]
y_best = np.ones(L, dtype=np.int) * (-1)
y_best[0] = ps
for t in range(1, L): y_best[t] = np.argmax(Fp[t-1, y_best[t-1], :])
Q = []
priority = -np.max(Alpha[L-1, :])
print(-priority)
print(Alpha[L-1, y_best[L-1]])
print(Fp[L-2, y_best[L-2], y_best[L-1]])
partition_index = -1
exclude_set = set()
hq.heappush(Q, HeapItem(priority, (y_best, partition_index, exclude_set)))
print(y_best)
print('-----------------------')
y1 = np.ones(L, dtype=np.int) * -1
y1[0] = y_best[0]
y1[2:] = y_best[2:]
cpoints = [p for p in range(M) if p != y_best[1]]
y1scores = [Cp[y1[0], p] + Cu[p] + Cp[p, y1[2]] for p in cpoints]
y1ix = np.argmax(y1scores)
y1[1] = cpoints[y1ix]
print(y1)
print('-----------------------')
score1 = -priority + Fp[1, y1[1], y_best[2]] - Fp[1, y_best[1], y_best[2]]
score2 = 0
for j in range(L-1):
ss = y1[j]
tt = y1[j+1]
score2 += Cp[ss, tt] + Cu[tt]
#print(-priority)
print(score1)
print(score2)
print('-----------------------')
results = []
k = 0; y_last = None
while len(Q) > 0 and k < maxIter:
hitem = hq.heappop(Q)
k_priority = hitem.priority
(k_best, k_partition_index, k_exclude_set) = hitem.task
k += 1; y_last = k_best
print('OUT: %s, %.6f' % (k_best, -k_priority))
if len(set(k_best)) == L:
results.append(k_best); top -= 1
if top == 0: return results
partition_index_start = L-1
if k_partition_index > 0:
assert(k_partition_index < L)
partition_index_start = k_partition_index
for parix in range(partition_index_start, 0, -1):
new_best = np.zeros(L, dtype=np.int) * (-1)
# new_best[parix+1:]
new_best[parix+1:] = k_best[parix+1:]
new_best[0] = ps
# new_best[parix]
new_exclude_set = set({k_best[parix]})
if parix == partition_index_start: new_exclude_set = new_exclude_set | k_exclude_set
candidate_points = [p for p in range(M) if p not in new_exclude_set]
if len(candidate_points) == 0: continue
if parix == 1:
candidate_maxix = np.argmax([Cp[ps, p] + Cu[p] + Cp[p, k_best[2]] for p in candidate_points])
new_best[parix] = candidate_points[candidate_maxix]
elif parix == L-1:
candidate_maxix = np.argmax([Alpha[L-1, p] for p in candidate_points])
new_best[parix] = candidate_points[candidate_maxix]
else:
candidate_maxix = np.argmax([Fp[parix, p, k_best[parix+1]] for p in candidate_points])
new_best[parix] = candidate_points[candidate_maxix]
# new_best[:parix]
if parix > 1:
for pk in range(parix-1, 0, -1):
if pk == 1:
new_best[pk] = np.argmax([Cp[ps, p] + Cu[p] + Cp[p, new_best[2]] for p in range(M)])
else:
new_best[pk] = np.argmax([Fp[pk, p, new_best[pk+1]] for p in range(M)])
# sequence score
new_priority = (-k_priority)
if parix == L-1:
new_priority += Alpha[L-1, new_best[L-1]] - Alpha[L-1, k_best[L-1]]
else:
new_priority += Fp[parix, new_best[parix], k_best[parix+1]] - Fp[parix, k_best[parix], k_best[parix+1]]
new_priority *= -1.0
new_score = 0
for j in range(L-1):
ss = new_best[j]
tt = new_best[j+1]
new_score += Cp[ss, tt] + Cu[tt]
#print(-new_priority); print(new_score)
assert(np.isclose(-new_priority, new_score))
new_priority = -new_score
hq.heappush(Q, HeapItem(new_priority, (new_best, parix, new_exclude_set)))
#print('IN : %s, %.6f, %d' % (new_best, -new_priority, parix))
if len(Q) == 0:
sys.stderr.write('WARN: empty queue, return the last one\n')
results.append(y_last); top -= 1
while len(Q) > 0 and top > 0:
hitem = hq.heappop(Q)
results.append(hitem.task[0]); top -= 1
return results
In [26]:
#M0 = 90
M0 = 10
n_u = 10
n_p = 5
w_u = np.random.rand(M0*n_u).reshape(M0, n_u)
f_u = np.random.rand(M0*n_u).reshape(M0, n_u)
w_p = np.random.rand(M0*M0*n_p).reshape(M0, M0, n_p)
f_p = np.random.rand(M0*M0*n_p).reshape(M0, M0, n_p)
ps0 = np.random.choice(np.arange(M0))
L0 = np.random.choice(np.arange(2, 8))
#L0 = 10
indices0 = [x for x in range(M0) if x != ps0]; np.random.shuffle(indices0)
y_true0 = [ps0] + indices0[:L0-1]
y_true_list0 = [y_true0]
for j in range(8):
np.random.shuffle(indices0); y_true_list0.append([ps0] + indices0[:L0-1])
print(ps0, L0)
In [27]:
brute_force(ps0, L0, M0, w_u, w_p, f_u, f_p, debug=True)
Out[27]:
In [28]:
do_inference_list_viterbi(ps0, L0, M0, w_u, w_p, f_u, f_p)
Out[28]:
In [37]:
list_viterbi_2001(ps0, L0, M0, w_u, w_p, f_u, f_p)
Out[37]:
In [39]:
list_viterbi_1994(ps0, L0, M0, w_u, w_p, f_u, f_p)
Out[39]:
In [ ]:
do_inference_list_viterbi(ps0, L0, M0, w_u, w_p, f_u, f_p, y_true=y_true0, y_true_list=y_true_list0) # allow sub-tours
In [ ]:
do_inference_brute_force(ps0, L0, M0, w_u, w_p, f_u, f_p, debug=True)
In [ ]:
do_inference_list_viterbi(ps0, L0, M0, w_u, w_p, f_u, f_p)
In [ ]:
a = [14, 0, 18, 27, 25]
#a = [14, 0, 27, 25, 18]
#a = [14, 0, 18, 27, 0]
priority = 0
In [ ]:
for t in range(1, L0):
ss = a[t-1]
tt = a[t]
priority += np.dot(w_p[ss, tt], f_p[ss, tt]) + np.dot(w_u[tt], f_u[tt])
In [ ]:
priority
In [ ]:
do_inference_brute_force(ps0, L0, M0, w_u, w_p, f_u, f_p, y_true=y_true0, y_true_list=y_true_list0)
In [ ]:
do_inference_list_viterbi(ps0, L0, M0, w_u, w_p, f_u, f_p, y_true=y_true0, y_true_list=y_true_list0) # allow sub-tours