Selecting distinct points from a plane using DPP

dpp_plane.py

import numpy as np
import numpy.linalg
import math

from numpy.random import random as rand

np.set_printoptions(precision=4, linewidth=200, suppress=False)

class DppProbabilityCalculator(object):
  def __init__(self, mat):
    """
    Create a DPP marginal kernel from a L matrix (mat)
    using eigenvalue scaling trick
    """
    self.matrix = mat
    (ev, evec) = np.linalg.eigh(mat)
    self.eval = ev
    self.evec = evec
    rescale = ev / (ev + 1)
    self.kernel = (evec * rescale).dot(evec.T)

  def prob_estimator(self, selected):
    return ProbEstimator(self.kernel, selected)


class ProbEstimator(object):
  def __init__(self, kernel, selection):
    self.selection = np.array(selection, dtype=np.int32)
    self.kernel = kernel

  def bruteforce(self, i):
    if i in self.selection:
      return 0.0
    mysel = np.append(self.selection, i)
    mykern = self.kernel[mysel, :][:, mysel]
    return np.linalg.det(mykern)

class Point:
  def __init__(self, x, y):
    self.x = x
    self.y = y

  def dist(self, o):
    dx = self.x - o.x
    dy = self.y - o.y
    return numpy.sqrt(dx * dx + dy * dy)


totnum = 20

points = [Point(totnum / 2 - x, totnum / 2 - y) for x in xrange(0, totnum) for y in xrange(0, totnum)]

plen = len(points)

distmatrix = np.zeros((plen, plen))

for i in xrange(0, plen):
  for j in xrange(0, plen):
    distmatrix[i, j] = points[i].dist(points[j])


print("Distance matrix is")
print(distmatrix)

maxel = np.max(distmatrix)
# this line adds no bias that items are similar
# mat = (maxel - distmatrix) / maxel

# this line adds a bias that items are (very) similar to each other
simiarity_bias = 0.9
mat = simiarity_bias + (1 - simiarity_bias) * (maxel - distmatrix) / maxel

# it is the parameter r in the paper
# chose the one you like



# this section multiplies the quality feature into L-kernel
# comment it out to get pure similarity selection

center = Point(0, 0)
cquality = np.zeros(plen)
coeff = totnum
for i in xrange(0, plen):
  cdist = points[i].dist(center)
  c = np.exp(-cdist / coeff)
  if c < 0.1:
    cquality[i] = 0.1
  else:
    cquality[i] = c


print("L-kernel before q is:")
print(mat)
cq = cquality.reshape((plen, 1))
cqx = cq.dot(cq.T)
print(cqx)

np.multiply(mat, cqx, mat)

# end of quality section

print("L-kernel is:")
print(mat)

dpp = DppProbabilityCalculator(mat)

print("DPP K-kernel is:")
print(dpp.kernel)

selected = []

probs = np.zeros((totnum, totnum))

toselect = 8

matrices = []

for trial in xrange(0, toselect):
  first = dpp.prob_estimator(selected)

  for i in xrange(0, plen):
    x = i / totnum
    y = i % totnum
    probs[x, y] = first.bruteforce(i)

  probsum = np.sum(probs)
  normProbs = probs / probsum
  matrices.append(normProbs)

  anitem = np.argmax(probs)
  print("selected item #%d" % anitem)
  selected.append(anitem)

import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import AxesGrid

fig = plt.figure(1)
grid = AxesGrid(fig, 111,
                nrows_ncols=(2, 4),
                axes_pad=0.1,
                share_all=True,
                cbar_location="top",
                cbar_mode="single",
                cbar_size="3%"
                )

xmax = np.max(matrices)
xmin = np.min(matrices)

cb = []

for i in xrange(0, toselect):
  matx = matrices[i]
  cb.append(grid[i].imshow(matx, vmax=xmax, vmin=xmin))

grid.cbar_axes[0].colorbar(cb[0])



print("Selection:")
print(selected)

plt.show()