Added different approximation to probabilistic PhocNet

cpp/phoc.cc

@@ -19,7 +19,13 @@ static inline T logsumexp(T a, T b) {
    return a + std::log1p(std::exp(b - a));
 }
 
-class ComputePair {
+template <typename T>
+static inline T logm1exp1m(T x) {
+  const T aux = -std::expm1(x);
+  return aux > 0 ? std::log(aux) : -std::numeric_limits<T>::infinity();
+}
+
+class ComputePairIndependece {
  public:
   template <typename Int, typename T>
   inline T operator()(Int n, const T* a, const T* b) {
@@ -37,19 +43,24 @@ class ComputePair {
   }
 };
 
-class ComputePairMin {
+class ComputePairUpperBound {
  public:
   template <typename Int, typename T>
   inline T operator()(Int n, const T* a, const T* b) {
-    T result = std::numeric_limits<T>::max();
-    for (Int i = 0; i < n; ++i) {
-      const T pa0 = -std::expm1(a[i]);
-      const T pb0 = -std::expm1(b[i]);
-      const T lh1 = a[i] + b[i];
-      const T lh0 = (pa0 > 0 && pb0 > 0)
-          ? std::log(pa0) + std::log(pb0)
-          : -std::numeric_limits<T>::infinity();
-      result = std::min(result, logsumexp(lh0, lh1));
+    if (n <= 0) { return 0; }
+    T ma = std::max(a[0], -std::expm1(a[0]));
+    T mb = std::max(b[0], -std::expm1(b[0]));
+    T result = std::max(a[0] + b[0], -(std::expm1(a[0]) + std::expm1(b[0])));
+    for (Int i = 1; i < n; ++i) {
+        const T a1 = a[i];
+        const T b1 = b[i];
+        const T a0 = logm1exp1m(a1);
+        const T b0 = logm1exp1m(b1);
+        const T min0 = std::min({a0 + b0, a0 + mb, b0 + ma, result});
+        const T min1 = std::min({a1 + b1, a1 + mb, b1 + ma, result});
+        result = std::max(min0, min1);
+        ma = std::min(ma, std::max(a0, a1));
+        mb = std::min(mb, std::max(b0, b1));
     }
     return result;
   }
@@ -137,26 +148,26 @@ static inline int pphoc(const ConstTensor<TT>& X,
     ConstTensor<TTYPE> tX(X);                                           \
     ConstTensor<TTYPE> tY(Y);                                           \
     MutableTensor<TTYPE> tR(R);                                         \
-    return cphoc(tX, tY, &tR, ComputePair());                           \
+    return cphoc(tX, tY, &tR, ComputePairIndependece());                \
   }                                                                     \
                                                                         \
   int pphoc_##STYPE(const TTYPE* X, TTYPE* R) {                         \
     ConstTensor<TTYPE> tX(X);                                           \
     MutableTensor<TTYPE> tR(R);                                         \
-    return pphoc(tX, &tR, ComputePair());                               \
+    return pphoc(tX, &tR, ComputePairIndependece());                    \
   }                                                                     \
                                                                         \
-  int cphoc_min_##STYPE(const TTYPE* X, const TTYPE* Y, TTYPE* R) {     \
+  int cphoc_max_##STYPE(const TTYPE* X, const TTYPE* Y, TTYPE* R) {     \
     ConstTensor<TTYPE> tX(X);                                           \
     ConstTensor<TTYPE> tY(Y);                                           \
     MutableTensor<TTYPE> tR(R);                                         \
-    return cphoc(tX, tY, &tR, ComputePairMin());                        \
+    return cphoc(tX, tY, &tR, ComputePairUpperBound());                 \
   }                                                                     \
                                                                         \
-  int pphoc_min_##STYPE(const TTYPE* X, TTYPE* R) {                     \
+  int pphoc_max_##STYPE(const TTYPE* X, TTYPE* R) {                     \
     ConstTensor<TTYPE> tX(X);                                           \
     MutableTensor<TTYPE> tR(R);                                         \
-    return pphoc(tX, &tR, ComputePairMin());                            \
+    return pphoc(tX, &tR, ComputePairUpperBound());                     \
   }
 
 DEFINE_WRAPPER(f32, THFloatTensor)

cpp/phoc.h

@@ -7,10 +7,10 @@ int pphoc_f32(const THFloatTensor* X, THFloatTensor* R);
 int pphoc_f64(const THDoubleTensor* X, THDoubleTensor* R);
 
 
-int cphoc_min_f32(const THFloatTensor* X, const THFloatTensor* Y,
+int cphoc_max_f32(const THFloatTensor* X, const THFloatTensor* Y,
                   THFloatTensor* R);
-int cphoc_min_f64(const THDoubleTensor* X, const THDoubleTensor* Y,
+int cphoc_max_f64(const THDoubleTensor* X, const THDoubleTensor* Y,
                   THDoubleTensor* R);
 
-int pphoc_min_f32(const THFloatTensor* X, THFloatTensor* R);
-int pphoc_min_f64(const THDoubleTensor* X, THDoubleTensor* R);
+int pphoc_max_f32(const THFloatTensor* X, THFloatTensor* R);
+int pphoc_max_f64(const THDoubleTensor* X, THDoubleTensor* R);

prob_phoc/default_impl.py

@@ -0,0 +1,52 @@
+from __future__ import absolute_import
+
+import math
+from scipy.misc import logsumexp
+import numpy as np
+import torch
+
+
+def _logmexpm1(x):
+    try:
+        return math.log(-math.expm1(x))
+    except ValueError:
+        return -np.inf
+
+
+def compute_independent(a, b):
+    assert torch.is_tensor(a)
+    assert torch.is_tensor(b)
+    assert len(a) == len(b)
+    a, b = a.view(-1), b.view(-1)
+
+    result = 0.0
+    for a1, b1 in zip(a, b):
+        h1 = a1 + b1
+        try:
+            h0 = math.log(-math.expm1(a1)) + math.log(-math.expm1(b1))
+        except ValueError:
+            h0 = -np.inf
+        result += logsumexp([h0, h1])
+
+    return result
+
+
+def compute_upper_bound(a, b):
+    assert torch.is_tensor(a)
+    assert torch.is_tensor(b)
+    assert len(a) == len(b)
+    a, b = a.view(-1), b.view(-1)
+    if len(a) == 0:
+        return 0.0
+
+    ma = max(a[0], -math.expm1(a[0]))
+    mb = max(b[0], -math.expm1(b[0]))
+    result = max(a[0] + b[0], -(math.expm1(a[0]) + math.expm1(b[0])))
+    for a1, b1 in zip(a[1:], b[1:]):
+        a0, b0 = _logmexpm1(a1), _logmexpm1(b1)
+        aux0 = min(a0 + b0, a0 + mb, b0 + ma, result)
+        aux1 = min(a1 + b1, a1 + mb, b1 + ma, result)
+        result = max(aux0, aux1)
+        ma = min(ma, max(a0, a1))
+        mb = min(mb, max(b0, b1))
+    return result

prob_phoc/phoc.py

@@ -6,81 +6,56 @@
 import warnings
 
 try:
-    from ._ext import cphoc_f32, cphoc_f64, cphoc_min_f32, cphoc_min_f64
-    from ._ext import pphoc_f32, pphoc_f64, pphoc_min_f32, pphoc_min_f64
+    from ._ext import cphoc_f32, cphoc_f64, cphoc_max_f32, cphoc_max_f64
+    from ._ext import pphoc_f32, pphoc_f64, pphoc_max_f32, pphoc_max_f64
 except ImportError as ex:
-    warnings.warn('The C++ implementation of prob_phoc could not be imported '
-                  '(%s). Python implementation will be used. ' % str(ex))
-
-    import math
-    from scipy.misc import logsumexp
-
-    def _compute(a, b):
-        result = 0.0
-        for i in range(a.size(0)):
-            h1 = a[i] + b[i]
-            try:
-                h0 = math.log(-math.expm1(a[i])) + math.log(-math.expm1(b[i]))
-            except ValueError:
-                h0 = -np.inf
-            result += logsumexp([h0, h1])
-
-        return result
-
-    def _compute_min(a, b):
-        result = float('inf')
-        for i in range(a.size(0)):
-            h1 = a[i] + b[i]
-            try:
-                h0 = math.log(-math.expm1(a[i])) + math.log(-math.expm1(b[i]))
-            except ValueError:
-                h0 = -np.inf
-            result = min(result, logsumexp([h0, h1]))
-
-        return result
+    warnings.warn(
+        "The C++ implementation of prob_phoc could not be imported "
+        "(%s). Python implementation will be used. " % str(ex)
+    )
 
+    from .default_impl import compute_independent, compute_upper_bound
 
     def _cphoc(x, y, out):
         for i in range(x.size(0)):
             for j in range(y.size(0)):
-                out[i, j] = _compute(x[i], y[j])
-
+                out[i, j] = compute_independent(x[i], y[j])
 
     def _pphoc(x, out):
         k = 0
         for i in range(x.size(0)):
             for j in range(i + 1, x.size(0)):
-                out[k] = _compute(x[i], x[j])
+                out[k] = compute_independent(x[i], x[j])
                 k += 1
 
-    def _cphoc_min(x, y, out):
+    def _cphoc_max(x, y, out):
         for i in range(x.size(0)):
             for j in range(y.size(0)):
-                out[i, j] = _compute_min(x[i], y[j])
-
+                out[i, j] = compute_upper_bound(x[i], y[j])
 
-    def _pphoc_min(x, out):
+    def _pphoc_max(x, out):
         k = 0
         for i in range(x.size(0)):
             for j in range(i + 1, x.size(0)):
-                out[k] = _compute_min(x[i], x[j])
+                out[k] = compute_upper_bound(x[i], x[j])
                 k += 1
 
     cphoc_f32 = cphoc_f64 = _cphoc
     pphoc_f32 = pphoc_f64 = _pphoc
-    cphoc_min_f32 = cphoc_min_f64 = _cphoc_min
-    pphoc_min_f32 = pphoc_min_f64 = _pphoc_min
+    cphoc_max_f32 = cphoc_max_f64 = _cphoc_max
+    pphoc_max_f32 = pphoc_max_f64 = _pphoc_max
 
 
-def cphoc(x, y, out=None, product=True):
+def cphoc(x, y, out=None, method="independence"):
     """Computes probabilistic PHOC relevance scores between each pair of inputs.
     """
+    assert method in ["independence", "upper_bound"]
     if isinstance(x, np.ndarray):
         x = torch.from_numpy(x)
     if isinstance(y, np.ndarray):
         y = torch.from_numpy(y)
     assert torch.is_tensor(x) and torch.is_tensor(y)
-    assert x.type() in ['torch.FloatTensor', 'torch.DoubleTensor']
+    assert x.type() in ["torch.FloatTensor", "torch.DoubleTensor"]
     assert x.type() == y.type()
     assert x.dim() == 2 and y.dim() == 2
     assert x.size(1) == y.size(1)
@@ -89,41 +64,42 @@ def cphoc(x, y, out=None, product=True):
     if out is None:
         out = x.new(x.size(0), y.size(0)).zero_()
 
-    if x.type() == 'torch.FloatTensor':
-        if product:
+    if x.type() == "torch.FloatTensor":
+        if method == "independence":
             cphoc_f32(x, y, out)
         else:
-            cphoc_min_f32(x, y, out)
+            cphoc_max_f32(x, y, out)
     else:
-        if product:
+        if method == "independence":
             cphoc_f64(x, y, out)
         else:
-            cphoc_min_f64(x, y, out)
+            cphoc_max_f64(x, y, out)
 
     return out
 
 
-def pphoc(x, out=None, product=True):
+def pphoc(x, out=None, method="independence"):
     """Pairwise probabilistic PHOC relevance scores."""
+    assert method in ["independence", "upper_bound"]
     if isinstance(x, np.ndarray):
         x = torch.from_numpy(x)
     assert torch.is_tensor(x)
-    assert x.type() in ['torch.FloatTensor', 'torch.DoubleTensor']
+    assert x.type() in ["torch.FloatTensor", "torch.DoubleTensor"]
     assert x.dim() == 2
     x = x.cpu()
 
     if out is None:
-        out = x.new(x.size(0) * (x.size(0) - 1) // 2,).zero_()
+        out = x.new(x.size(0) * (x.size(0) - 1) // 2).zero_()
 
-    if x.type() == 'torch.FloatTensor':
-        if product:
+    if x.type() == "torch.FloatTensor":
+        if method == "independence":
             pphoc_f32(x, out)
         else:
-            pphoc_min_f32(x, out)
+            pphoc_max_f32(x, out)
     else:
-        if product:
+        if method == "independence":
             pphoc_f64(x, out)
         else:
-            pphoc_min_f64(x, out)
+            pphoc_max_f64(x, out)
 
     return out