Improved prob.h/cpp code and unit tests
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Mon, 29 Mar 2010 15:02:32 +0000 (17:02 +0200)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Mon, 29 Mar 2010 15:02:32 +0000 (17:02 +0200)
  - TProb<T>::accumulate() now also applies op2 to init
  - Fixed bug by renaming TProb<T>::operator<=() to TProb<T>::operator<()
  - TProb<T>& operator/= (T x) now yields 0 when dividing by 0
  - Changed format of TProb<T> when streamed to an ostream

ChangeLog
include/dai/prob.h
tests/unit/prob.cpp

index 1ef55ba..3c6095a 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -2,7 +2,11 @@ git HEAD
 --------
 
 * Improved prob.h/cpp:
-  - Added operator==( const TProb<T> &q )
+  - Added TProb<T>::operator==( const TProb<T> &q )
+  - TProb<T>::accumulate() now also applies op2 to init
+  - Fixed bug by renaming TProb<T>::operator<=() to TProb<T>::operator<()
+  - TProb<T>& operator/= (T x) now yields 0 when dividing by 0
+  - Changed format of TProb<T> when streamed to an ostream
 * Improved index.h/cpp:
   - Added multifor::reset()
 * Improved properties.h/cpp:
index df133ad..2036e17 100644 (file)
@@ -305,8 +305,16 @@ template <typename T> class TProb {
         size_t size() const { return _p.size(); }
 
         /// Accumulate over all values, similar to std::accumulate
+        /** The following calculation is done:
+         *  \code
+         *  T t = op2(init);
+         *  for( const_iterator it = begin(); it != end(); it++ )
+         *      t = op1( t, op2(*it) );
+         *  return t;
+         *  \endcode
+         */
         template<typename binOp, typename unOp> T accumulate( T init, binOp op1, unOp op2 ) const {
-            T t = init;
+            T t = op2(init);
             for( const_iterator it = begin(); it != end(); it++ )
                 t = op1( t, op2(*it) );
             return t;
@@ -374,7 +382,7 @@ template <typename T> class TProb {
         /// Lexicographical comparison
         /** \pre <tt>this->size() == q.size()</tt>
          */
-        bool operator<=( const TProb<T>& q ) const {
+        bool operator<( const TProb<T>& q ) const {
             DAI_DEBASSERT( size() == q.size() );
             return lexicographical_compare( begin(), end(), q.begin(), q.end() );
         }
@@ -528,11 +536,10 @@ template <typename T> class TProb {
                 return *this;
         }
 
-        /// Divides each entry by scalar \a x
+        /// Divides each entry by scalar \a x, where division by 0 yields 0
         TProb<T>& operator/= (T x) {
-            DAI_DEBASSERT( x != 0 );
             if( x != 1 )
-                return pwUnaryOp( std::bind2nd( std::divides<T>(), x ) );
+                return pwUnaryOp( std::bind2nd( fo_divides0<T>(), x ) );
             else
                 return *this;
         }
@@ -695,9 +702,10 @@ template<typename T> T dist( const TProb<T> &p, const TProb<T> &q, typename TPro
 /** \relates TProb
  */
 template<typename T> std::ostream& operator<< (std::ostream& os, const TProb<T>& p) {
-    os << "[";
-    std::copy( p.p().begin(), p.p().end(), std::ostream_iterator<T>(os, " ") );
-    os << "]";
+    os << "(";
+    for( typename std::vector<T>::const_iterator it = p.begin(); it != p.end(); it++ )
+        os << (it != p.begin() ? ", " : "") << *it;
+    os << ")";
     return os;
 }
 
index 00045ee..e038866 100644 (file)
 using namespace dai;
 
 
+const double tol = 1e-8;
+
+
 #define BOOST_TEST_MODULE ProbTest
 
 
 #include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
 
 
 BOOST_AUTO_TEST_CASE( ConstructorsTest ) {
@@ -115,28 +119,631 @@ BOOST_AUTO_TEST_CASE( IteratorTest ) {
 
 
 BOOST_AUTO_TEST_CASE( QueriesTest ) {
+    Prob x( 5, 0.0 );
+    for( size_t i = 0; i < x.size(); i++ )
+        x[i] = 2.0 - i;
+
+    // test accumulate, min, max, sum, sumAbs, maxAbs
+    BOOST_CHECK_EQUAL( x.sum(), 0.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( 0.0, std::plus<Real>(), fo_id<Real>() ), 0.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( 1.0, std::plus<Real>(), fo_id<Real>() ), 1.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( -1.0, std::plus<Real>(), fo_id<Real>() ), -1.0 );
+    BOOST_CHECK_EQUAL( x.max(), 2.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( -INFINITY, fo_max<Real>(), fo_id<Real>() ), 2.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( 3.0, fo_max<Real>(), fo_id<Real>() ), 3.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( -5.0, fo_max<Real>(), fo_id<Real>() ), 2.0 );
+    BOOST_CHECK_EQUAL( x.min(), -2.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( INFINITY, fo_min<Real>(), fo_id<Real>() ), -2.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( -3.0, fo_min<Real>(), fo_id<Real>() ), -3.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( 5.0, fo_min<Real>(), fo_id<Real>() ), -2.0 );
+    BOOST_CHECK_EQUAL( x.sumAbs(), 6.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( 0.0, std::plus<Real>(), fo_abs<Real>() ), 6.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( 1.0, std::plus<Real>(), fo_abs<Real>() ), 7.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( -1.0, std::plus<Real>(), fo_abs<Real>() ), 7.0 );
+    BOOST_CHECK_EQUAL( x.maxAbs(), 2.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( 0.0, fo_max<Real>(), fo_abs<Real>() ), 2.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( 1.0, fo_max<Real>(), fo_abs<Real>() ), 2.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( -1.0, fo_max<Real>(), fo_abs<Real>() ), 2.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( 3.0, fo_max<Real>(), fo_abs<Real>() ), 3.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( -3.0, fo_max<Real>(), fo_abs<Real>() ), 3.0 );
+    x[1] = 1.0;
+    BOOST_CHECK_EQUAL( x.maxAbs(), 2.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( 0.0, fo_max<Real>(), fo_abs<Real>() ), 2.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( 1.0, fo_max<Real>(), fo_abs<Real>() ), 2.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( -1.0, fo_max<Real>(), fo_abs<Real>() ), 2.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( 3.0, fo_max<Real>(), fo_abs<Real>() ), 3.0 );
+    BOOST_CHECK_EQUAL( x.accumulate( -3.0, fo_max<Real>(), fo_abs<Real>() ), 3.0 );
+    for( size_t i = 0; i < x.size(); i++ )
+        x[i] = i ? (1.0 / i) : 0.0;
+    BOOST_CHECK_EQUAL( x.accumulate( 0.0, std::plus<Real>(), fo_inv0<Real>() ), 10.0 );
+    x /= x.sum();
+
+    // test entropy
+    BOOST_CHECK( x.entropy() < Prob(5).entropy() );
+    for( size_t i = 1; i < 100; i++ )
+        BOOST_CHECK_CLOSE( Prob(i).entropy(), std::log(i), tol );
+
+    // test hasNaNs and hasNegatives
+    BOOST_CHECK( !Prob( 3, 0.0 ).hasNaNs() );
+    Real c = 0.0;
+    BOOST_CHECK( Prob( 3, c / c ).hasNaNs() );
+    BOOST_CHECK( !Prob( 3, 0.0 ).hasNegatives() );
+    BOOST_CHECK( !Prob( 3, 1.0 ).hasNegatives() );
+    BOOST_CHECK( Prob( 3, -1.0 ).hasNegatives() );
+    x[0] = 0.0; x[1] = 0.0; x[2] = -1.0; x[3] = 1.0; x[4] = 100.0;
+    BOOST_CHECK( x.hasNegatives() );
+    x[2] = -INFINITY;
+    BOOST_CHECK( x.hasNegatives() );
+    x[2] = INFINITY;
+    BOOST_CHECK( !x.hasNegatives() );
+    x[2] = -1.0;
+
+    // test argmax
+    BOOST_CHECK( x.argmax() == std::make_pair( (size_t)4, (Real)100.0 ) );
+    x[4] = 0.5;
+    BOOST_CHECK( x.argmax() == std::make_pair( (size_t)3, (Real)1.0 ) );
+    x[3] = -2.0;
+    BOOST_CHECK( x.argmax() == std::make_pair( (size_t)4, (Real)0.5 ) );
+    x[4] = -1.0;
+    BOOST_CHECK( x.argmax() == std::make_pair( (size_t)0, (Real)0.0 ) );
+    x[0] = -2.0;
+    BOOST_CHECK( x.argmax() == std::make_pair( (size_t)1, (Real)0.0 ) );
+    x[1] = -3.0;
+    BOOST_CHECK( x.argmax() == std::make_pair( (size_t)2, (Real)-1.0 ) );
+    x[2] = -2.0;
+    BOOST_CHECK( x.argmax() == std::make_pair( (size_t)4, (Real)-1.0 ) );
+
+    // test draw
+    for( size_t i = 0; i < x.size(); i++ )
+        x[i] = i ? (1.0 / i) : 0.0;
+    for( size_t repeat = 0; repeat < 10000; repeat++ ) {
+        BOOST_CHECK( x.draw() < x.size() );
+        BOOST_CHECK( x.draw() != 0 );
+    }
+    x[2] = 0.0;
+    for( size_t repeat = 0; repeat < 10000; repeat++ ) {
+        BOOST_CHECK( x.draw() < x.size() );
+        BOOST_CHECK( x.draw() != 0 );
+        BOOST_CHECK( x.draw() != 2 );
+    }
+    x[4] = 0.0;
+    for( size_t repeat = 0; repeat < 10000; repeat++ ) {
+        BOOST_CHECK( x.draw() < x.size() );
+        BOOST_CHECK( x.draw() != 0 );
+        BOOST_CHECK( x.draw() != 2 );
+        BOOST_CHECK( x.draw() != 4 );
+    }
+    x[1] = 0.0;
+    for( size_t repeat = 0; repeat < 10000; repeat++ )
+        BOOST_CHECK( x.draw() == 3 );
+
+    // test <, ==
+    Prob a(3, 1.0), b(3, 1.0);
+    BOOST_CHECK( !(a < b) );
+    BOOST_CHECK( !(b < a) );
+    BOOST_CHECK( a == b );
+    a[0] = 0.0;
+    BOOST_CHECK( a < b );
+    BOOST_CHECK( !(b < a) );
+    BOOST_CHECK( !(a == b) );
+    b[2] = 0.0;
+    BOOST_CHECK( a < b );
+    BOOST_CHECK( !(b < a) );
+    BOOST_CHECK( !(a == b) );
+    b[0] = 0.0;
+    BOOST_CHECK( !(a < b) );
+    BOOST_CHECK( b < a );
+    BOOST_CHECK( !(a == b) );
+    a[1] = 0.0;
+    BOOST_CHECK( a < b );
+    BOOST_CHECK( !(b < a) );
+    BOOST_CHECK( !(a == b) );
+    b[1] = 0.0;
+    BOOST_CHECK( !(a < b) );
+    BOOST_CHECK( b < a );
+    BOOST_CHECK( !(a == b) );
+    a[2] = 0.0;
+    BOOST_CHECK( !(a < b) );
+    BOOST_CHECK( !(b < a) );
+    BOOST_CHECK( a == b );
 }
 
 
 BOOST_AUTO_TEST_CASE( UnaryTransformationsTest ) {
+    Prob x( 3 );
+    x[0] = -2.0;
+    x[1] = 0.0;
+    x[2] = 2.0;
+
+    Prob y = -x;
+    Prob z = x.pwUnaryTr( std::negate<Real>() );
+    BOOST_CHECK_EQUAL( y[0], 2.0 );
+    BOOST_CHECK_EQUAL( y[1], 0.0 );
+    BOOST_CHECK_EQUAL( y[2], -2.0 );
+    BOOST_CHECK( y == z );
+
+    y = x.abs();
+    z = x.pwUnaryTr( fo_abs<Real>() );
+    BOOST_CHECK_EQUAL( y[0], 2.0 );
+    BOOST_CHECK_EQUAL( y[1], 0.0 );
+    BOOST_CHECK_EQUAL( y[2], 2.0 );
+    BOOST_CHECK( y == z );
+
+    y = x.exp();
+    z = x.pwUnaryTr( fo_exp<Real>() );
+    BOOST_CHECK_CLOSE( y[0], std::exp(-2.0), tol );
+    BOOST_CHECK_EQUAL( y[1], 1.0 );
+    BOOST_CHECK_CLOSE( y[2], 1.0 / y[0], tol );
+    BOOST_CHECK( y == z );
+
+    y = x.log(false);
+    z = x.pwUnaryTr( fo_log<Real>() );
+    BOOST_CHECK( isnan( y[0] ) );
+    BOOST_CHECK_EQUAL( y[1], -INFINITY );
+    BOOST_CHECK_CLOSE( y[2], std::log(2.0), tol );
+    BOOST_CHECK( !(y == z) );
+    y[0] = z[0] = 0.0;
+    BOOST_CHECK( y == z );
+
+    y = x.log(true);
+    z = x.pwUnaryTr( fo_log0<Real>() );
+    BOOST_CHECK( isnan( y[0] ) );
+    BOOST_CHECK_EQUAL( y[1], 0.0 );
+    BOOST_CHECK_EQUAL( y[2], std::log(2.0) );
+    BOOST_CHECK( !(y == z) );
+    y[0] = z[0] = 0.0;
+    BOOST_CHECK( y == z );
+
+    y = x.inverse(false);
+    z = x.pwUnaryTr( fo_inv<Real>() );
+    BOOST_CHECK_EQUAL( y[0], -0.5 );
+    BOOST_CHECK_EQUAL( y[1], INFINITY );
+    BOOST_CHECK_EQUAL( y[2], 0.5 );
+    BOOST_CHECK( y == z );
+
+    y = x.inverse(true);
+    z = x.pwUnaryTr( fo_inv0<Real>() );
+    BOOST_CHECK_EQUAL( y[0], -0.5 );
+    BOOST_CHECK_EQUAL( y[1], 0.0 );
+    BOOST_CHECK_EQUAL( y[2], 0.5 );
+    BOOST_CHECK( y == z );
+
+    x[0] = 2.0;
+    y = x.normalized();
+    BOOST_CHECK_EQUAL( y[0], 0.5 );
+    BOOST_CHECK_EQUAL( y[1], 0.0 );
+    BOOST_CHECK_EQUAL( y[2], 0.5 );
+
+    y = x.normalized( Prob::NORMPROB );
+    BOOST_CHECK_EQUAL( y[0], 0.5 );
+    BOOST_CHECK_EQUAL( y[1], 0.0 );
+    BOOST_CHECK_EQUAL( y[2], 0.5 );
+
+    x[0] = -2.0;
+    y = x.normalized( Prob::NORMLINF );
+    BOOST_CHECK_EQUAL( y[0], -1.0 );
+    BOOST_CHECK_EQUAL( y[1], 0.0 );
+    BOOST_CHECK_EQUAL( y[2], 1.0 );
 }
 
 
 BOOST_AUTO_TEST_CASE( UnaryOperationsTest ) {
+    Prob xorg(3);
+    xorg[0] = 2.0;
+    xorg[1] = 0.0;
+    xorg[2] = 1.0;
+    Prob y(3);
+
+    Prob x = xorg;
+    BOOST_CHECK( x.setUniform() == Prob(3) );
+    BOOST_CHECK( x == Prob(3) );
+
+    y[0] = std::exp(2.0);
+    y[1] = 1.0;
+    y[2] = std::exp(1.0);
+    x = xorg;
+    BOOST_CHECK( x.takeExp() == y );
+    BOOST_CHECK( x == y );
+    x = xorg;
+    BOOST_CHECK( x.pwUnaryOp( fo_exp<Real>() ) == y );
+    BOOST_CHECK( x == y );
+
+    y[0] = std::log(2.0);
+    y[1] = -INFINITY;
+    y[2] = 0.0;
+    x = xorg;
+    BOOST_CHECK( x.takeLog() == y );
+    BOOST_CHECK( x == y );
+    x = xorg;
+    BOOST_CHECK( x.takeLog(false) == y );
+    BOOST_CHECK( x == y );
+    x = xorg;
+    BOOST_CHECK( x.pwUnaryOp( fo_log<Real>() ) == y );
+    BOOST_CHECK( x == y );
+
+    y[1] = 0.0;
+    x = xorg;
+    BOOST_CHECK( x.takeLog(true) == y );
+    BOOST_CHECK( x == y );
+    x = xorg;
+    BOOST_CHECK( x.pwUnaryOp( fo_log0<Real>() ) == y );
+    BOOST_CHECK( x == y );
+
+    y[0] = 2.0 / 3.0;
+    y[1] = 0.0 / 3.0;
+    y[2] = 1.0 / 3.0;
+    x = xorg;
+    BOOST_CHECK_EQUAL( x.normalize(), 3.0 );
+    BOOST_CHECK( x == y );
+
+    x = xorg;
+    BOOST_CHECK_EQUAL( x.normalize( Prob::NORMPROB ), 3.0 );
+    BOOST_CHECK( x == y );
+
+    y[0] = 2.0 / 2.0;
+    y[1] = 0.0 / 2.0;
+    y[2] = 1.0 / 2.0;
+    x = xorg;
+    BOOST_CHECK_EQUAL( x.normalize( Prob::NORMLINF ), 2.0 );
+    BOOST_CHECK( x == y );
+
+    xorg[0] = -2.0;
+    y[0] = 2.0;
+    y[1] = 0.0;
+    y[2] = 1.0;
+    x = xorg;
+    BOOST_CHECK( x.takeAbs() == y );
+    BOOST_CHECK( x == y );
+
+    for( size_t repeat = 0; repeat < 10000; repeat++ ) {
+        x.randomize();
+        for( size_t i = 0; i < x.size(); i++ ) {
+            BOOST_CHECK( x[i] < 1.0 );
+            BOOST_CHECK( x[i] >= 0.0 );
+        }
+    }
 }
 
 
-BOOST_AUTO_TEST_CASE( ScalarOperationsTest ) {
+BOOST_AUTO_TEST_CASE( ScalarTransformationsTest ) {
+    Prob x(3);
+    x[0] = 2.0;
+    x[1] = 0.0;
+    x[2] = 1.0;
+    Prob y(3);
+
+    y[0] = 3.0; y[1] = 1.0; y[2] = 2.0;
+    BOOST_CHECK( (x + 1.0) == y );
+    y[0] = 0.0; y[1] = -2.0; y[2] = -1.0;
+    BOOST_CHECK( (x + (-2.0)) == y );
+
+    y[0] = 1.0; y[1] = -1.0; y[2] = 0.0;
+    BOOST_CHECK( (x - 1.0) == y );
+    y[0] = 4.0; y[1] = 2.0; y[2] = 3.0;
+    BOOST_CHECK( (x - (-2.0)) == y );
+
+    BOOST_CHECK( (x * 1.0) == x );
+    y[0] = 4.0; y[1] = 0.0; y[2] = 2.0;
+    BOOST_CHECK( (x * 2.0) == y );
+    y[0] = -1.0; y[1] = 0.0; y[2] = -0.5;
+    BOOST_CHECK( (x * -0.5) == y );
+
+    BOOST_CHECK( (x / 1.0) == x );
+    y[0] = 1.0; y[1] = 0.0; y[2] = 0.5;
+    BOOST_CHECK( (x / 2.0) == y );
+    y[0] = -4.0; y[1] = 0.0; y[2] = -2.0;
+    BOOST_CHECK( (x / -0.5) == y );
+    BOOST_CHECK( (x / 0.0) == Prob(3, 0.0) );
+
+    BOOST_CHECK( (x ^ 1.0) == x );
+    BOOST_CHECK( (x ^ 0.0) == Prob(3, 1.0) );
+    y[0] = 4.0; y[1] = 0.0; y[2] = 1.0;
+    BOOST_CHECK( (x ^ 2.0) == y );
+    y[0] = 1.0 / std::sqrt(2.0); y[1] = INFINITY; y[2] = 1.0;
+    Prob z = (x ^ -0.5);
+    BOOST_CHECK_CLOSE( z[0], y[0], tol );
+    BOOST_CHECK_EQUAL( z[1], y[1] );
+    BOOST_CHECK_CLOSE( z[2], y[2], tol );
 }
 
 
-BOOST_AUTO_TEST_CASE( VectorTransformationsTest ) {
+BOOST_AUTO_TEST_CASE( ScalarOperationsTest ) {
+    Prob xorg(3), x(3);
+    xorg[0] = 2.0;
+    xorg[1] = 0.0;
+    xorg[2] = 1.0;
+    Prob y(3);
+
+    x = xorg;
+    BOOST_CHECK( x.fill( 1.0 ) == Prob(3, 1.0) );
+    BOOST_CHECK( x == Prob(3, 1.0) );
+    BOOST_CHECK( x.fill( 2.0 ) == Prob(3, 2.0) );
+    BOOST_CHECK( x == Prob(3, 2.0) );
+    BOOST_CHECK( x.fill( 0.0 ) == Prob(3, 0.0) );
+    BOOST_CHECK( x == Prob(3, 0.0) );
+
+    x = xorg;
+    y[0] = 3.0; y[1] = 1.0; y[2] = 2.0;
+    BOOST_CHECK( (x += 1.0) == y );
+    BOOST_CHECK( x == y );
+    y[0] = 1.0; y[1] = -1.0; y[2] = 0.0;
+    BOOST_CHECK( (x += -2.0) == y );
+    BOOST_CHECK( x == y );
+
+    x = xorg;
+    y[0] = 1.0; y[1] = -1.0; y[2] = 0.0;
+    BOOST_CHECK( (x -= 1.0) == y );
+    BOOST_CHECK( x == y );
+    y[0] = 3.0; y[1] = 1.0; y[2] = 2.0;
+    BOOST_CHECK( (x -= -2.0) == y );
+    BOOST_CHECK( x == y );
+
+    x = xorg;
+    BOOST_CHECK( (x *= 1.0) == x );
+    BOOST_CHECK( x == x );
+    y[0] = 4.0; y[1] = 0.0; y[2] = 2.0;
+    BOOST_CHECK( (x *= 2.0) == y );
+    BOOST_CHECK( x == y );
+    y[0] = -1.0; y[1] = 0.0; y[2] = -0.5;
+    BOOST_CHECK( (x *= -0.25) == y );
+    BOOST_CHECK( x == y );
+
+    x = xorg;
+    BOOST_CHECK( (x /= 1.0) == x );
+    BOOST_CHECK( x == x );
+    y[0] = 1.0; y[1] = 0.0; y[2] = 0.5;
+    BOOST_CHECK( (x /= 2.0) == y );
+    BOOST_CHECK( x == y );
+    y[0] = -4.0; y[1] = 0.0; y[2] = -2.0;
+    BOOST_CHECK( (x /= -0.25) == y );
+    BOOST_CHECK( x == y );
+    BOOST_CHECK( (x /= 0.0) == Prob(3, 0.0) );
+    BOOST_CHECK( x == Prob(3, 0.0) );
+
+    x = xorg;
+    BOOST_CHECK( (x ^= 1.0) == x );
+    BOOST_CHECK( x == x );
+    BOOST_CHECK( (x ^= 0.0) == Prob(3, 1.0) );
+    BOOST_CHECK( x == Prob(3, 1.0) );
+    x = xorg;
+    y[0] = 4.0; y[1] = 0.0; y[2] = 1.0;
+    BOOST_CHECK( (x ^= 2.0) == y );
+    BOOST_CHECK( x == y );
+    y[0] = 0.5; y[1] = INFINITY; y[2] = 1.0;
+    BOOST_CHECK( (x ^= -0.5) == y );
+    BOOST_CHECK( x == y );
 }
 
 
 BOOST_AUTO_TEST_CASE( VectorOperationsTest ) {
+    size_t N = 6;
+    Prob xorg(N), x(N);
+    xorg[0] = 2.0; xorg[1] = 0.0; xorg[2] = 1.0; xorg[3] = 0.0; xorg[4] = 2.0; xorg[5] = 3.0;
+    Prob y(N);
+    y[0] = 0.5; y[1] = -1.0; y[2] = 0.0; y[3] = 0.0; y[4] = -2.0; y[5] = 3.0;
+    Prob z(N), r(N);
+
+    z[0] = 2.5; z[1] = -1.0; z[2] = 1.0; z[3] = 0.0; z[4] = 0.0; z[5] = 6.0;
+    x = xorg;
+    r = (x += y);
+    for( size_t i = 0; i < N; i++ )
+        BOOST_CHECK_CLOSE( r[i], z[i], tol );
+    BOOST_CHECK( x == z );
+    x = xorg;
+    BOOST_CHECK( x.pwBinaryOp( y, std::plus<Real>() ) == z );
+    BOOST_CHECK( x == z );
+
+    z[0] = 1.5; z[1] = 1.0; z[2] = 1.0; z[3] = 0.0; z[4] = 4.0; z[5] = 0.0;
+    x = xorg;
+    r = (x -= y);
+    for( size_t i = 0; i < N; i++ )
+        BOOST_CHECK_CLOSE( r[i], z[i], tol );
+    BOOST_CHECK( x == z );
+    x = xorg;
+    BOOST_CHECK( x.pwBinaryOp( y, std::minus<Real>() ) == z );
+    BOOST_CHECK( x == z );
+
+    z[0] = 1.0; z[1] = 0.0; z[2] = 0.0; z[3] = 0.0; z[4] = -4.0; z[5] = 9.0;
+    x = xorg;
+    r = (x *= y);
+    for( size_t i = 0; i < N; i++ )
+        BOOST_CHECK_CLOSE( r[i], z[i], tol );
+    BOOST_CHECK( x == z );
+    x = xorg;
+    BOOST_CHECK( x.pwBinaryOp( y, std::multiplies<Real>() ) == z );
+    BOOST_CHECK( x == z );
+
+    z[0] = 4.0; z[1] = 0.0; z[2] = 0.0; z[3] = 0.0; z[4] = -1.0; z[5] = 1.0;
+    x = xorg;
+    r = (x /= y);
+    for( size_t i = 0; i < N; i++ )
+        BOOST_CHECK_CLOSE( r[i], z[i], tol );
+    BOOST_CHECK( x == z );
+    x = xorg;
+    BOOST_CHECK( x.pwBinaryOp( y, fo_divides0<Real>() ) == z );
+    BOOST_CHECK( x == z );
+
+    z[0] = 4.0; z[1] = 0.0; z[2] = INFINITY; /*z[3] = INFINITY;*/ z[4] = -1.0; z[5] = 1.0;
+    x = xorg;
+    r = (x.divide( y ));
+    BOOST_CHECK_CLOSE( r[0], z[0], tol );
+    BOOST_CHECK_CLOSE( r[1], z[1], tol );
+    BOOST_CHECK_EQUAL( r[2], z[2] );
+    BOOST_CHECK( isnan(r[3]) );
+    BOOST_CHECK_CLOSE( r[4], z[4], tol );
+    BOOST_CHECK_CLOSE( r[5], z[5], tol );
+    x[3] = 0.0; r[3] = 0.0;
+    BOOST_CHECK( x == r );
+    x = xorg;
+    r = x.pwBinaryOp( y, std::divides<Real>() );
+    BOOST_CHECK_CLOSE( r[0], z[0], tol );
+    BOOST_CHECK_CLOSE( r[1], z[1], tol );
+    BOOST_CHECK_EQUAL( r[2], z[2] );
+    BOOST_CHECK( isnan(r[3]) );
+    BOOST_CHECK_CLOSE( r[4], z[4], tol );
+    BOOST_CHECK_CLOSE( r[5], z[5], tol );
+    x[3] = 0.0; r[3] = 0.0;
+    BOOST_CHECK( x == r );
+
+    z[0] = std::sqrt(2.0); z[1] = INFINITY; z[2] = 1.0; z[3] = 1.0; z[4] = 0.25; z[5] = 27.0;
+    x = xorg;
+    r = (x ^= y);
+    BOOST_CHECK_CLOSE( r[0], z[0], tol );
+    BOOST_CHECK_EQUAL( r[1], z[1] );
+    BOOST_CHECK_CLOSE( r[2], z[2], tol );
+    BOOST_CHECK_CLOSE( r[3], z[3], tol );
+    BOOST_CHECK_CLOSE( r[4], z[4], tol );
+    BOOST_CHECK_CLOSE( r[5], z[5], tol );
+    BOOST_CHECK( x == z );
+    x = xorg;
+    BOOST_CHECK( x.pwBinaryOp( y, fo_pow<Real>() ) == z );
+    BOOST_CHECK( x == z );
+}
+
+
+BOOST_AUTO_TEST_CASE( VectorTransformationsTest ) {
+    size_t N = 6;
+    Prob x(N);
+    x[0] = 2.0; x[1] = 0.0; x[2] = 1.0; x[3] = 0.0; x[4] = 2.0; x[5] = 3.0;
+    Prob y(N);
+    y[0] = 0.5; y[1] = -1.0; y[2] = 0.0; y[3] = 0.0; y[4] = -2.0; y[5] = 3.0;
+    Prob z(N), r(N);
+
+    z[0] = 2.5; z[1] = -1.0; z[2] = 1.0; z[3] = 0.0; z[4] = 0.0; z[5] = 6.0;
+    r = x + y;
+    for( size_t i = 0; i < N; i++ )
+        BOOST_CHECK_CLOSE( r[i], z[i], tol );
+    z = x.pwBinaryTr( y, std::plus<Real>() );
+    BOOST_CHECK( r == z );
+
+    z[0] = 1.5; z[1] = 1.0; z[2] = 1.0; z[3] = 0.0; z[4] = 4.0; z[5] = 0.0;
+    r = x - y;
+    for( size_t i = 0; i < N; i++ )
+        BOOST_CHECK_CLOSE( r[i], z[i], tol );
+    z = x.pwBinaryTr( y, std::minus<Real>() );
+    BOOST_CHECK( r == z );
+
+    z[0] = 1.0; z[1] = 0.0; z[2] = 0.0; z[3] = 0.0; z[4] = -4.0; z[5] = 9.0;
+    r = x * y;
+    for( size_t i = 0; i < N; i++ )
+        BOOST_CHECK_CLOSE( r[i], z[i], tol );
+    z = x.pwBinaryTr( y, std::multiplies<Real>() );
+    BOOST_CHECK( r == z );
+
+    z[0] = 4.0; z[1] = 0.0; z[2] = 0.0; z[3] = 0.0; z[4] = -1.0; z[5] = 1.0;
+    r = x / y;
+    for( size_t i = 0; i < N; i++ )
+        BOOST_CHECK_CLOSE( r[i], z[i], tol );
+    z = x.pwBinaryTr( y, fo_divides0<Real>() );
+    BOOST_CHECK( r == z );
+
+    z[0] = 4.0; z[1] = 0.0; z[2] = INFINITY; /*z[3] = INFINITY;*/ z[4] = -1.0; z[5] = 1.0;
+    r = x.divided_by( y );
+    BOOST_CHECK_CLOSE( r[0], z[0], tol );
+    BOOST_CHECK_CLOSE( r[1], z[1], tol );
+    BOOST_CHECK_EQUAL( r[2], z[2] );
+    BOOST_CHECK( isnan(r[3]) );
+    BOOST_CHECK_CLOSE( r[4], z[4], tol );
+    BOOST_CHECK_CLOSE( r[5], z[5], tol );
+    z = x.pwBinaryTr( y, std::divides<Real>() );
+    BOOST_CHECK_CLOSE( r[0], z[0], tol );
+    BOOST_CHECK_CLOSE( r[1], z[1], tol );
+    BOOST_CHECK_EQUAL( r[2], z[2] );
+    BOOST_CHECK( isnan(r[3]) );
+    BOOST_CHECK_CLOSE( r[4], z[4], tol );
+    BOOST_CHECK_CLOSE( r[5], z[5], tol );
+
+    z[0] = std::sqrt(2.0); z[1] = INFINITY; z[2] = 1.0; z[3] = 1.0; z[4] = 0.25; z[5] = 27.0;
+    r = x ^ y;
+    BOOST_CHECK_CLOSE( r[0], z[0], tol );
+    BOOST_CHECK_EQUAL( r[1], z[1] );
+    BOOST_CHECK_CLOSE( r[2], z[2], tol );
+    BOOST_CHECK_CLOSE( r[3], z[3], tol );
+    BOOST_CHECK_CLOSE( r[4], z[4], tol );
+    BOOST_CHECK_CLOSE( r[5], z[5], tol );
+    z = x.pwBinaryTr( y, fo_pow<Real>() );
+    BOOST_CHECK( r == z );
 }
 
 
 BOOST_AUTO_TEST_CASE( RelatedFunctionsTest ) {
+    Prob x(3), y(3), z(3);
+    x[0] = 0.2;
+    x[1] = 0.8;
+    x[2] = 0.0;
+    y[0] = 0.0;
+    y[1] = 0.6;
+    y[2] = 0.4;
+
+    z = min( x, y );
+    BOOST_CHECK_EQUAL( z[0], 0.0 );
+    BOOST_CHECK_EQUAL( z[1], 0.6 );
+    BOOST_CHECK_EQUAL( z[2], 0.0 );
+    z = max( x, y );
+    BOOST_CHECK_EQUAL( z[0], 0.2 );
+    BOOST_CHECK_EQUAL( z[1], 0.8 );
+    BOOST_CHECK_EQUAL( z[2], 0.4 );
+
+    BOOST_CHECK_EQUAL( dist( x, x, Prob::DISTL1 ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( y, y, Prob::DISTL1 ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTL1 ), 0.2 + 0.2 + 0.4 );
+    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTL1 ), 0.2 + 0.2 + 0.4 );
+    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTL1 ), x.innerProduct( y, 0.0, std::plus<Real>(), fo_absdiff<Real>() ) );
+    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTL1 ), y.innerProduct( x, 0.0, std::plus<Real>(), fo_absdiff<Real>() ) );
+    BOOST_CHECK_EQUAL( dist( x, x, Prob::DISTLINF ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( y, y, Prob::DISTLINF ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTLINF ), 0.4 );
+    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTLINF ), 0.4 );
+    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTLINF ), x.innerProduct( y, 0.0, fo_max<Real>(), fo_absdiff<Real>() ) );
+    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTLINF ), y.innerProduct( x, 0.0, fo_max<Real>(), fo_absdiff<Real>() ) );
+    BOOST_CHECK_EQUAL( dist( x, x, Prob::DISTTV ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( y, y, Prob::DISTTV ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTTV ), 0.5 * (0.2 + 0.2 + 0.4) );
+    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTTV ), 0.5 * (0.2 + 0.2 + 0.4) );
+    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTTV ), x.innerProduct( y, 0.0, std::plus<Real>(), fo_absdiff<Real>() ) / 2.0 );
+    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTTV ), y.innerProduct( x, 0.0, std::plus<Real>(), fo_absdiff<Real>() ) / 2.0 );
+    BOOST_CHECK_EQUAL( dist( x, x, Prob::DISTKL ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( y, y, Prob::DISTKL ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTKL ), INFINITY );
+    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTKL ), INFINITY );
+    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTKL ), x.innerProduct( y, 0.0, std::plus<Real>(), fo_KL<Real>() ) );
+    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTKL ), y.innerProduct( x, 0.0, std::plus<Real>(), fo_KL<Real>() ) );
+    BOOST_CHECK_EQUAL( dist( x, x, Prob::DISTHEL ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( y, y, Prob::DISTHEL ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTHEL ), 0.5 * (0.2 + std::pow(std::sqrt(0.8) - std::sqrt(0.6), 2.0) + 0.4) );
+    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTHEL ), 0.5 * (0.2 + std::pow(std::sqrt(0.8) - std::sqrt(0.6), 2.0) + 0.4) );
+    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTHEL ), x.innerProduct( y, 0.0, std::plus<Real>(), fo_Hellinger<Real>() ) / 2.0 );
+    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTHEL ), y.innerProduct( x, 0.0, std::plus<Real>(), fo_Hellinger<Real>() ) / 2.0 );
+    x[1] = 0.7; x[2] = 0.1;
+    y[0] = 0.1; y[1] = 0.5;
+    BOOST_CHECK_CLOSE( dist( x, y, Prob::DISTKL ), 0.2 * std::log(0.2 / 0.1) + 0.7 * std::log(0.7 / 0.5) + 0.1 * std::log(0.1 / 0.4), tol );
+    BOOST_CHECK_CLOSE( dist( y, x, Prob::DISTKL ), 0.1 * std::log(0.1 / 0.2) + 0.5 * std::log(0.5 / 0.7) + 0.4 * std::log(0.4 / 0.1), tol );
+    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTKL ), x.innerProduct( y, 0.0, std::plus<Real>(), fo_KL<Real>() ) );
+    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTKL ), y.innerProduct( x, 0.0, std::plus<Real>(), fo_KL<Real>() ) );
+
+    std::stringstream ss;
+    ss << x;
+    std::string s;
+    std::getline( ss, s );
+    BOOST_CHECK_EQUAL( s, std::string("(0.2, 0.7, 0.1)") );
+    std::stringstream ss2;
+    ss2 << y;
+    std::getline( ss2, s );
+    BOOST_CHECK_EQUAL( s, std::string("(0.1, 0.5, 0.4)") );
+
+    z = min( x, y );
+    BOOST_CHECK_EQUAL( z[0], 0.1 );
+    BOOST_CHECK_EQUAL( z[1], 0.5 );
+    BOOST_CHECK_EQUAL( z[2], 0.1 );
+    z = max( x, y );
+    BOOST_CHECK_EQUAL( z[0], 0.2 );
+    BOOST_CHECK_EQUAL( z[1], 0.7 );
+    BOOST_CHECK_EQUAL( z[2], 0.4 );
+
+    BOOST_CHECK_CLOSE( x.innerProduct( y, 0.0, std::plus<Real>(), std::multiplies<Real>() ), 0.2*0.1 + 0.7*0.5 + 0.1*0.4, tol );
+    BOOST_CHECK_CLOSE( y.innerProduct( x, 0.0, std::plus<Real>(), std::multiplies<Real>() ), 0.2*0.1 + 0.7*0.5 + 0.1*0.4, tol );
+    BOOST_CHECK_CLOSE( x.innerProduct( y, 1.0, std::plus<Real>(), std::multiplies<Real>() ), 1.0 + 0.2*0.1 + 0.7*0.5 + 0.1*0.4, tol );
+    BOOST_CHECK_CLOSE( y.innerProduct( x, 1.0, std::plus<Real>(), std::multiplies<Real>() ), 1.0 + 0.2*0.1 + 0.7*0.5 + 0.1*0.4, tol );
+    BOOST_CHECK_CLOSE( x.innerProduct( y, -1.0, std::plus<Real>(), std::multiplies<Real>() ), -1.0 + 0.2*0.1 + 0.7*0.5 + 0.1*0.4, tol );
+    BOOST_CHECK_CLOSE( y.innerProduct( x, -1.0, std::plus<Real>(), std::multiplies<Real>() ), -1.0 + 0.2*0.1 + 0.7*0.5 + 0.1*0.4, tol );
 }