diff --git a/cpp/phoc.cc b/cpp/phoc.cc index 530df53..24a2d4a 100644 --- a/cpp/phoc.cc +++ b/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 +static inline T logm1exp1m(T x) { + const T aux = -std::expm1(x); + return aux > 0 ? std::log(aux) : -std::numeric_limits::infinity(); +} + +class ComputePairIndependece { public: template inline T operator()(Int n, const T* a, const T* b) { @@ -37,19 +43,24 @@ class ComputePair { } }; -class ComputePairMin { +class ComputePairUpperBound { public: template inline T operator()(Int n, const T* a, const T* b) { - T result = std::numeric_limits::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::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& X, ConstTensor tX(X); \ ConstTensor tY(Y); \ MutableTensor tR(R); \ - return cphoc(tX, tY, &tR, ComputePair()); \ + return cphoc(tX, tY, &tR, ComputePairIndependece()); \ } \ \ int pphoc_##STYPE(const TTYPE* X, TTYPE* R) { \ ConstTensor tX(X); \ MutableTensor 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 tX(X); \ ConstTensor tY(Y); \ MutableTensor 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 tX(X); \ MutableTensor tR(R); \ - return pphoc(tX, &tR, ComputePairMin()); \ + return pphoc(tX, &tR, ComputePairUpperBound()); \ } DEFINE_WRAPPER(f32, THFloatTensor) diff --git a/cpp/phoc.h b/cpp/phoc.h index 0ad771f..90a6124 100644 --- a/cpp/phoc.h +++ b/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); diff --git a/prob_phoc/default_impl.py b/prob_phoc/default_impl.py new file mode 100644 index 0000000..fef248b --- /dev/null +++ b/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 diff --git a/prob_phoc/phoc.py b/prob_phoc/phoc.py index 66c96ed..5d64677 100644 --- a/prob_phoc/phoc.py +++ b/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 diff --git a/prob_phoc/test.py b/prob_phoc/test.py index e53b5af..55fddff 100644 --- a/prob_phoc/test.py +++ b/prob_phoc/test.py @@ -6,46 +6,113 @@ import torch from prob_phoc.phoc import cphoc, pphoc +from prob_phoc.default_impl import compute_upper_bound -class ProbPhocTest(unittest.TestCase): - def test_cphoc_arbitrary(self): +def cphoc_upper_bound(xx, yy): + assert xx.size(1) == yy.size(1) + out = xx.new(xx.size(0), yy.size(0)).zero_() + for i, x in enumerate(xx): + for j, y in enumerate(yy): + out[i, j] = compute_upper_bound(x, y) + return out + + +def pphoc_upper_bound(xx): + N = xx.size(0) + out = cphoc_upper_bound(xx, xx) + out2 = [] + for i in range(N): + for j in range(i + 1, N): + out2.append(out[i, j]) + if isinstance(xx, torch.FloatTensor): + return torch.FloatTensor(out2) + else: + return torch.DoubleTensor(out2) + + +class ProbPHOCTest(unittest.TestCase): + @staticmethod + def test_cphoc_arbitrary_independence(): + x = torch.DoubleTensor([[0.7, 0.4], [0.9, 0.6]]).log_() + y = torch.DoubleTensor([[0.5, 0.9], [0.8, 0.9], [0.1, 0.2]]).log_() + expected_z = np.log( + np.asarray( + [ + [ + (.7 * .5 + .3 * .5) * (.4 * .9 + .6 * .1), + (.7 * .8 + .3 * .2) * (.4 * .9 + .6 * .1), + (.7 * .1 + .3 * .9) * (.4 * .2 + .6 * .8), + ], + [ + (.9 * .5 + .1 * .5) * (.6 * .9 + .4 * .1), + (.9 * .8 + .1 * .2) * (.6 * .9 + .4 * .1), + (.9 * .1 + .1 * .9) * (.6 * .2 + .4 * .8), + ], + ], + dtype=np.float32, + ) + ) + # Test float64 + z = cphoc(x, y, method="independence").numpy() + np.testing.assert_almost_equal(z, expected_z) + # Test float32 + z = cphoc( + x.type("torch.FloatTensor"), + y.type("torch.FloatTensor"), + method="independence", + ) + z = z.numpy() + np.testing.assert_almost_equal(z, expected_z.astype(np.float32)) + + @staticmethod + def test_cphoc_arbitrary_upper_bound(): x = torch.DoubleTensor([[0.7, 0.4], [0.9, 0.6]]).log_() y = torch.DoubleTensor([[0.5, 0.9], [0.8, 0.9], [0.1, 0.2]]).log_() - expected_z = np.log(np.asarray([ - [ - (.7 * .5 + .3 * .5) * (.4 * .9 + .6 * .1), - (.7 * .8 + .3 * .2) * (.4 * .9 + .6 * .1), - (.7 * .1 + .3 * .9) * (.4 * .2 + .6 * .8) - ], - [ - (.9 * .5 + .1 * .5) * (.6 * .9 + .4 * .1), - (.9 * .8 + .1 * .2) * (.6 * .9 + .4 * .1), - (.9 * .1 + .1 * .9) * (.6 * .2 + .4 * .8), - ] - ], dtype=np.float32)) + expected_z = cphoc_upper_bound(x, y).numpy() # Test float64 - z = cphoc(x, y).numpy() + z = cphoc(x, y, method="upper_bound").numpy() np.testing.assert_almost_equal(z, expected_z) # Test float32 - z = cphoc(x.type('torch.FloatTensor'), y.type('torch.FloatTensor')) + z = cphoc( + x.type("torch.FloatTensor"), + y.type("torch.FloatTensor"), + method="upper_bound", + ) z = z.numpy() np.testing.assert_almost_equal(z, expected_z.astype(np.float32)) - def test_phoc_arbitrary(self): + @staticmethod + def test_pphoc_independence(): x = torch.DoubleTensor([[.7, .4], [.9, .4], [.5, .2]]).log_() - expected_y = np.log(np.asarray([ - (.7 * .9 + .3 * .1) * (.4 * .4 + .6 * .6), - (.7 * .5 + .3 * .5) * (.4 * .2 + .6 * .8), - (.9 * .5 + .1 * .5) * (.4 * .2 + .6 * .8), - ], dtype=np.float64)) + expected_y = np.log( + np.asarray( + [ + (.7 * .9 + .3 * .1) * (.4 * .4 + .6 * .6), + (.7 * .5 + .3 * .5) * (.4 * .2 + .6 * .8), + (.9 * .5 + .1 * .5) * (.4 * .2 + .6 * .8), + ], + dtype=np.float64, + ) + ) # Test float64 (double) - y = pphoc(x).numpy() + y = pphoc(x, method="independence").numpy() np.testing.assert_almost_equal(y, expected_y) # Test float32 (float) - y = pphoc(x.type('torch.FloatTensor')).numpy() + y = pphoc(x.type("torch.FloatTensor"), method="independence").numpy() np.testing.assert_almost_equal(y, expected_y.astype(np.float32)) + @staticmethod + def test_pphoc_arbitrary_upper_bound(): + x = torch.DoubleTensor([[.7, .4], [.9, .4], [.5, .2]]).log_() + expected_z = pphoc_upper_bound(x).numpy() + # Test float64 + z = pphoc(x, method="upper_bound").numpy() + np.testing.assert_almost_equal(z, expected_z) + # Test float32 + z = pphoc(x.type("torch.FloatTensor"), method="upper_bound").numpy() + np.testing.assert_almost_equal(z, expected_z.astype(np.float32)) + -if __name__ == '__main__': +if __name__ == "__main__": unittest.main()