Improved factor.h/cpp code and finished corresponding unit tests
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Fri, 2 Apr 2010 15:20:06 +0000 (17:20 +0200)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Fri, 2 Apr 2010 15:20:06 +0000 (17:20 +0200)
- Fixed bug in min( const TFactor<T> &, const TFactor<T> & )
- Fixed bug in max( const TFactor<T> &, const TFactor<T> & )
- Added TFactor<T>::takeAbs()
- Added TFactor<T>::takeExp()
- Added TFactor<T>::takeLog( bool )
- Added TFactor<T>::operator-() const
- Renamed TFactor<T>::states() into TFactor<T>::nrStates()

14 files changed:
ChangeLog
include/dai/factor.h
include/dai/prob.h
src/bbp.cpp
src/bp.cpp
src/bp_dual.cpp
src/cbp.cpp
src/emalg.cpp
src/factor.cpp
src/factorgraph.cpp
src/gibbs.cpp
src/matlab/matlab.cpp
swig/example_sprinkler.m
tests/unit/factor.cpp

index 5d35e84..8d185e8 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -7,8 +7,15 @@ TODO:
 * Fixed some bugs in the MatLab interface build system
 * Fixed a bug in utils/fginfo.cpp
 * Improved factor.h/cpp:
+  - Fixed bug in min( const TFactor<T> &, const TFactor<T> & )
+  - Fixed bug in max( const TFactor<T> &, const TFactor<T> & )
+  - Added TFactor<T>::takeAbs()
+  - Added TFactor<T>::takeExp()
+  - Added TFactor<T>::takeLog( bool )
+  - Added TFactor<T>::operator-() const
   - Added TFactor<T>::sumAbs() const
   - Added TFactor<T>::operator=( const TFactor<T> &y )
+  - Renamed TFactor<T>::states() into TFactor<T>::nrStates()
   - Added get(size_t) and set(size_t,T) operators; get() is 
     equivalent to "operator[](size_t) const" and set() should
     be used instead of the non-const operator[], which has been deprecated
@@ -72,7 +79,7 @@ TODO:
 * Removed RandomDRegularGraph(), please use createGraphRegular() instead
 * Compressed Makefile
 * Added unit tests for var, smallset, varset, graph, bipgraph,
-  weightedgraph, enum, util, properties, index, prob
+  weightedgraph, enum, util, properties, index, prob, factor
 * Added unit testing framework
 * Added initialization of TRWBP weights by sampling spanning trees
 * Cleaned up MR code:
index 48e135c..9925835 100644 (file)
@@ -95,11 +95,15 @@ template <typename T> class TFactor {
 
         /// Constructs factor depending on variables in \a vars, copying the values from \a p
         TFactor( const VarSet& vars, const TProb<T> &p ) : _vs(vars), _p(p) {
-            DAI_DEBASSERT( _vs.nrStates() == _p.size() );
+            DAI_ASSERT( _vs.nrStates() == _p.size() );
         }
 
         /// Constructs factor depending on variables in \a vars, permuting the values given in \a p accordingly
         TFactor( const std::vector<Var> &vars, const std::vector<T> &p ) : _vs(vars.begin(), vars.end(), vars.size()), _p(p.size()) {
+            size_t nrStates = 1;
+            for( size_t i = 0; i < vars.size(); i++ )
+                nrStates *= vars[i].states();
+            DAI_ASSERT( nrStates == p.size() );
             Permute permindex(vars);
             for( size_t li = 0; li < p.size(); ++li )
                 _p.set( permindex.convertLinearIndex(li), p[li] );
@@ -136,6 +140,12 @@ template <typename T> class TFactor {
         /// Returns the number of possible joint states of the variables on which the factor depends, \f$\prod_{l\in L} S_l\f$
         /** \note This is equal to the length of the value vector.
          */
+        size_t nrStates() const { return _p.size(); }
+
+        /// Returns the number of possible joint states of the variables on which the factor depends, \f$\prod_{l\in L} S_l\f$
+        /** \note This is equal to the length of the value vector.
+         *  \deprecated Please use dai::TFactor::nrStates() instead.
+         */
         size_t states() const { return _p.size(); }
 
         /// Returns the Shannon entropy of \c *this, \f$-\sum_i p_i \log p_i\f$
@@ -173,11 +183,19 @@ template <typename T> class TFactor {
 
     /// \name Unary transformations
     //@{
-        /// Returns pointwise absolute value
-        TFactor<T> abs() const {
+        /// Returns negative of \c *this
+        TFactor<T> operator- () const { 
             // Note: the alternative (shorter) way of implementing this,
             //   return TFactor<T>( _vs, _p.abs() );
             // is slower because it invokes the copy constructor of TProb<T>
+            TFactor<T> x;
+            x._vs = _vs;
+            x._p = -_p;
+            return x;
+        }
+
+        /// Returns pointwise absolute value
+        TFactor<T> abs() const {
             TFactor<T> x;
             x._vs = _vs;
             x._p = _p.abs();
@@ -226,10 +244,21 @@ template <typename T> class TFactor {
     /// \name Unary operations
     //@{
         /// Draws all values i.i.d. from a uniform distribution on [0,1)
-        TFactor<T> & randomize () { _p.randomize(); return *this; }
+        TFactor<T>& randomize() { _p.randomize(); return *this; }
 
         /// Sets all values to \f$1/n\f$ where \a n is the number of states
-        TFactor<T>& setUniform () { _p.setUniform(); return *this; }
+        TFactor<T>& setUniform() { _p.setUniform(); return *this; }
+
+        /// Applies absolute value pointwise
+        TFactor<T>& takeAbs() { _p.takeAbs(); return *this; }
+
+        /// Applies exponent pointwise
+        TFactor<T>& takeExp() { _p.takeExp(); return *this; }
+
+        /// Applies logarithm pointwise
+        /** If \a zero == \c true, uses <tt>log(0)==0</tt>; otherwise, <tt>log(0)==-Inf</tt>.
+         */
+        TFactor<T>& takeLog( bool zero = false ) { _p.takeLog(zero); return *this; }
 
         /// Normalizes factor using the specified norm
         /** \throw NOT_NORMALIZABLE if the norm is zero
@@ -240,7 +269,7 @@ template <typename T> class TFactor {
     /// \name Operations with scalars
     //@{
         /// Sets all values to \a x
-        TFactor<T> & fill (T x) { _p.fill( x ); return *this; }
+        TFactor<T>& fill (T x) { _p.fill( x ); return *this; }
 
         /// Adds scalar \a x to each value
         TFactor<T>& operator+= (T x) { _p += x; return *this; }
@@ -430,7 +459,7 @@ template <typename T> class TFactor {
     //@{
         /// Returns a slice of \c *this, where the subset \a vars is in state \a varsState
         /** \pre \a vars sould be a subset of vars()
-         *  \pre \a varsState < vars.states()
+         *  \pre \a varsState < vars.nrStates()
          *
          *  The result is a factor that depends on the variables of *this except those in \a vars,
          *  obtained by setting the variables in \a vars to the joint state specified by the linear index
@@ -473,7 +502,7 @@ template<typename T> TFactor<T> TFactor<T>::slice( const VarSet& vars, size_t va
     // OPTIMIZE ME
     IndexFor i_vars (vars, _vs);
     IndexFor i_varsrem (varsrem, _vs);
-    for( size_t i = 0; i < states(); i++, ++i_vars, ++i_varsrem )
+    for( size_t i = 0; i < nrStates(); i++, ++i_vars, ++i_varsrem )
         if( (size_t)i_vars == varsState )
             result.set( i_varsrem, _p[i] );
 
@@ -548,7 +577,7 @@ template<typename T> T TFactor<T>::strength( const Var &i, const Var &j ) const
  */
 template<typename T> std::ostream& operator<< (std::ostream& os, const TFactor<T>& f) {
     os << "(" << f.vars() << ", (";
-    for( size_t i = 0; i < f.states(); i++ )
+    for( size_t i = 0; i < f.nrStates(); i++ )
         os << (i == 0 ? "" : ", ") << f[i];
     os << "))";
     return os;
@@ -574,8 +603,8 @@ template<typename T> T dist( const TFactor<T> &f, const TFactor<T> &g, typename
  *  \pre f.vars() == g.vars()
  */
 template<typename T> TFactor<T> max( const TFactor<T> &f, const TFactor<T> &g ) {
-    DAI_ASSERT( f._vs == g._vs );
-    return TFactor<T>( f._vs, max( f.p(), g.p() ) );
+    DAI_ASSERT( f.vars() == g.vars() );
+    return TFactor<T>( f.vars(), max( f.p(), g.p() ) );
 }
 
 
@@ -584,8 +613,8 @@ template<typename T> TFactor<T> max( const TFactor<T> &f, const TFactor<T> &g )
  *  \pre f.vars() == g.vars()
  */
 template<typename T> TFactor<T> min( const TFactor<T> &f, const TFactor<T> &g ) {
-    DAI_ASSERT( f._vs == g._vs );
-    return TFactor<T>( f._vs, min( f.p(), g.p() ) );
+    DAI_ASSERT( f.vars() == g.vars() );
+    return TFactor<T>( f.vars(), min( f.p(), g.p() ) );
 }
 
 
@@ -606,14 +635,14 @@ template<typename T> T MutualInfo(const TFactor<T> &f) {
 typedef TFactor<Real> Factor;
 
 
-/// Returns a binary single-variable factor \f$ \exp(hx) \f$ where \f$ x = \pm 1 \f$
+/// Returns a binary unnormalized single-variable factor \f$ \exp(hx) \f$ where \f$ x = \pm 1 \f$
 /** \param x Variable (should be binary)
  *  \param h Field strength
  */
 Factor createFactorIsing( const Var &x, Real h );
 
 
-/// Returns a binary pairwise factor \f$ \exp(J x_1 x_2) \f$ where \f$ x_1, x_2 = \pm 1 \f$
+/// Returns a binary unnormalized pairwise factor \f$ \exp(J x_1 x_2) \f$ where \f$ x_1, x_2 = \pm 1 \f$
 /** \param x1 First variable (should be binary)
  *  \param x2 Second variable (should be binary)
  *  \param J Coupling strength
index d7d6846..a44da76 100644 (file)
@@ -495,15 +495,15 @@ class TProb {
         }
 
         /// Applies absolute value pointwise
-        const TProb<T>& takeAbs() { return pwUnaryOp( fo_abs<T>() ); }
+        TProb<T>& takeAbs() { return pwUnaryOp( fo_abs<T>() ); }
 
         /// Applies exponent pointwise
-        const TProb<T>& takeExp() { return pwUnaryOp( fo_exp<T>() ); }
+        TProb<T>& takeExp() { return pwUnaryOp( fo_exp<T>() ); }
 
         /// Applies logarithm pointwise
         /** If \a zero == \c true, uses <tt>log(0)==0</tt>; otherwise, <tt>log(0)==-Inf</tt>.
          */
-        const TProb<T>& takeLog(bool zero=false) {
+        TProb<T>& takeLog(bool zero=false) {
             if( zero ) {
                 return pwUnaryOp( fo_log0<T>() );
             } else
index fd643f9..8c73591 100644 (file)
@@ -157,7 +157,7 @@ void BBP::RegenerateInds() {
         _indices[i].reserve( _fg->nbV(i).size() );
         foreach( const Neighbor &I, _fg->nbV(i) ) {
             _ind_t index;
-            index.reserve( _fg->factor(I).states() );
+            index.reserve( _fg->factor(I).nrStates() );
             for( IndexFor k(_fg->var(i), _fg->factor(I).vars()); k.valid(); ++k )
                 index.push_back( k );
             _indices[i].push_back( index );
@@ -188,7 +188,7 @@ void BBP::RegenerateU() {
     for( size_t I = 0; I < _fg->nrFactors(); I++ ) {
         _Umsg[I].resize( _fg->nbF(I).size() );
         foreach( const Neighbor &i, _fg->nbF(I) ) {
-            Prob prod( _fg->factor(I).states(), 1.0 );
+            Prob prod( _fg->factor(I).nrStates(), 1.0 );
             foreach( const Neighbor &j, _fg->nbF(I) )
                 if( i.node != j.node ) {
                     Prob n_jI( _bp_dual.msgN( j, j.dual ) );
@@ -217,7 +217,7 @@ void BBP::RegenerateS() {
                         if( k != i && k.node != j.node ) {
                             const _ind_t &ind = _index( k, k.dual );
                             Prob p( _bp_dual.msgN( k, k.dual ) );
-                            for( size_t x_I = 0; x_I < prod.states(); x_I++ )
+                            for( size_t x_I = 0; x_I < prod.nrStates(); x_I++ )
                                 prod.set( x_I, prod[x_I] * p[ind[x_I]] );
                         }
                     }
@@ -279,7 +279,7 @@ void BBP::RegeneratePsiAdjoints() {
     _adj_psi_F.reserve( _fg->nrFactors() );
     for( size_t I = 0; I < _fg->nrFactors(); I++ ) {
         Prob p( _adj_b_F_unnorm[I] );
-        DAI_ASSERT( p.size() == _fg->factor(I).states() );
+        DAI_ASSERT( p.size() == _fg->factor(I).nrStates() );
         foreach( const Neighbor &i, _fg->nbF(I) ) {
             Prob n_iI( _bp_dual.msgN( i, i.dual ) );
             const _ind_t& ind = _index( i, i.dual );
@@ -737,7 +737,7 @@ std::vector<Prob> BBP::getZeroAdjF( const FactorGraph &fg ) {
     vector<Prob> adj_2;
     adj_2.reserve( fg.nrFactors() );
     for( size_t I = 0; I < fg.nrFactors(); I++ )
-        adj_2.push_back( Prob( fg.factor(I).states(), 0.0 ) );
+        adj_2.push_back( Prob( fg.factor(I).nrStates(), 0.0 ) );
     return adj_2;
 }
 
@@ -777,7 +777,7 @@ void BBP::initCostFnAdj( const BBPCostFunction &cfn, const vector<size_t> *state
                 psi1_adj.push_back( p );
             }
             for( size_t I = 0; I < fg.nrFactors(); I++ ) {
-                size_t dim = fg.factor(I).states();
+                size_t dim = fg.factor(I).nrStates();
                 Prob p( dim, 0.0 );
                 for( size_t xI = 0; xI < dim; xI++ )
                     p.set( xI, 1 + log( _ia->beliefF(I)[xI] / fg.factor(I).p()[xI] ) );
@@ -793,7 +793,7 @@ void BBP::initCostFnAdj( const BBPCostFunction &cfn, const vector<size_t> *state
             vector<Prob> b2_adj;
             b2_adj.reserve( fg.nrFactors() );
             for( size_t I = 0; I < fg.nrFactors(); I++ ) {
-                size_t dim = fg.factor(I).states();
+                size_t dim = fg.factor(I).nrStates();
                 Prob p( dim, 0.0 );
                 for( size_t xI = 0; xI < dim; xI++ ) {
                     Real bIxI = _ia->beliefF(I)[xI];
@@ -872,7 +872,7 @@ void BBP::initCostFnAdj( const BBPCostFunction &cfn, const vector<size_t> *state
             vector<Prob> b2_adj;
             b2_adj.reserve( fg.nrVars() );
             for( size_t I = 0; I <  fg.nrFactors(); I++ ) {
-                size_t n = fg.factor(I).states();
+                size_t n = fg.factor(I).nrStates();
                 Prob delta( n, 0.0 );
 
                 size_t x_I = getFactorEntryForState( fg, I, state );
index efa45ce..4f4f35c 100644 (file)
@@ -105,7 +105,7 @@ void BP::construct() {
             newEP.newMessage = Prob( var(i).states() );
 
             if( DAI_BP_FAST ) {
-                newEP.index.reserve( factor(I).states() );
+                newEP.index.reserve( factor(I).nrStates() );
                 for( IndexFor k( var(i), factor(I).vars() ); k.valid(); ++k )
                     newEP.index.push_back( k );
             }
index bf9f714..1591492 100644 (file)
@@ -62,7 +62,7 @@ void BP_dual::regenerateBeliefs() {
     for( size_t i = 0; i < fg().nrVars(); i++ )
         _beliefs.b1.push_back( Prob( fg().var(i).states() ) );
     for( size_t I = 0; I < fg().nrFactors(); I++ )
-        _beliefs.b2.push_back( Prob( fg().factor(I).states() ) );
+        _beliefs.b2.push_back( Prob( fg().factor(I).nrStates() ) );
 }
 
 
index ef208ac..0fbbb44 100644 (file)
@@ -196,12 +196,12 @@ void CBP::runRecurse( InfAlg *bp, Real orig_logZ, vector<size_t> clamped_vars_li
             DAI_ASSERT(/*0<=xi &&*/ xi < var(i).states() );
     else
         foreach( size_t xI, xis )
-            DAI_ASSERT(/*0<=xI &&*/ xI < factor(i).states() );
+            DAI_ASSERT(/*0<=xI &&*/ xI < factor(i).nrStates() );
     // - otherwise, clamp and recurse, saving margin estimates for each
     // clamp setting. afterwards, combine estimates.
 
     // compute complement of 'xis'
-    vector<size_t> cmp_xis = complement( xis, clampingVar ? var(i).states() : factor(i).states() );
+    vector<size_t> cmp_xis = complement( xis, clampingVar ? var(i).states() : factor(i).nrStates() );
 
     /// \idea dai::CBP::runRecurse() could be implemented more efficiently with a nesting version of backupFactors/restoreFactors
     // this improvement could also be done locally: backup the clamped factor in a local variable,
@@ -317,7 +317,7 @@ bool CBP::chooseNextClampVar( InfAlg *bp, vector<size_t> &clamped_vars_list, siz
             // only pick probable values for variable
             size_t xi;
             do {
-                xi = rnd( factor(i).states() );
+                xi = rnd( factor(i).nrStates() );
                 t++;
             } while( bp->beliefF(i).p()[xi] < tiny && t < t1 );
             DAI_ASSERT( t < t1 );
@@ -435,7 +435,7 @@ bool CBP::chooseNextClampVar( InfAlg *bp, vector<size_t> &clamped_vars_list, siz
         if( clampingVar )
             DAI_ASSERT(/*0<=xi &&*/ xi < var(i).states() );
         else
-            DAI_ASSERT(/*0<=xi &&*/ xi < factor(i).states() );
+            DAI_ASSERT(/*0<=xi &&*/ xi < factor(i).nrStates() );
         DAI_IFVERB(2, "CHOOSE_BBP (num clamped = " << clamped_vars_list.size() << ") chose " << i << " state " << xi << endl);
     } else
         DAI_THROW(UNKNOWN_ENUM_VALUE);
@@ -503,7 +503,7 @@ std::pair<size_t, size_t> BBPFindClampVar( const InfAlg &in_bp, bool clampingVar
             }
         }
         DAI_ASSERT(/*0 <= argmax_var_state &&*/
-               argmax_var_state < in_bp.fg().factor(argmax_var).states() );
+               argmax_var_state < in_bp.fg().factor(argmax_var).nrStates() );
     }
     if( maxVarOut )
         *maxVarOut = max_var;
index 3732190..73ff1d2 100644 (file)
@@ -173,8 +173,8 @@ void SharedParameters::collectSufficientStatistics( InfAlg &alg ) {
         VarSet &vs = _varsets[i->first];
 
         Factor b = alg.belief(vs);
-        Prob p( b.states(), 0.0 );
-        for( size_t entry = 0; entry < b.states(); ++entry )
+        Prob p( b.nrStates(), 0.0 );
+        for( size_t entry = 0; entry < b.nrStates(); ++entry )
             p.set( entry, b[perm.convertLinearIndex(entry)] ); // apply inverse permutation
         _estimation->addSufficientStatistics( p );
     }
@@ -188,7 +188,7 @@ void SharedParameters::setParameters( FactorGraph &fg ) {
         VarSet &vs = _varsets[i->first];
 
         Factor f( vs, 0.0 );
-        for( size_t entry = 0; entry < f.states(); ++entry )
+        for( size_t entry = 0; entry < f.nrStates(); ++entry )
             f.set( perm.convertLinearIndex(entry), p[entry] );
 
         fg.setFactor( i->first, f );
index 504b506..bcc4905 100644 (file)
@@ -40,7 +40,7 @@ Factor createFactorIsing( const Var &n1, const Var &n2, Real J ) {
 
 Factor createFactorExpGauss( const VarSet &ns, Real beta ) {
     Factor fac( ns );
-    for( size_t t = 0; t < fac.states(); t++ )
+    for( size_t t = 0; t < fac.nrStates(); t++ )
         fac.set( t, std::exp(rnd_stdnormal() * beta) );
     return fac;
 }
index c4bbf26..b21d7b2 100644 (file)
@@ -86,11 +86,11 @@ std::ostream& operator<< ( std::ostream &os, const FactorGraph &fg ) {
             os << i->states() << " ";
         os << endl;
         size_t nr_nonzeros = 0;
-        for( size_t k = 0; k < fg.factor(I).states(); k++ )
+        for( size_t k = 0; k < fg.factor(I).nrStates(); k++ )
             if( fg.factor(I)[k] != (Real)0 )
                 nr_nonzeros++;
         os << nr_nonzeros << endl;
-        for( size_t k = 0; k < fg.factor(I).states(); k++ )
+        for( size_t k = 0; k < fg.factor(I).nrStates(); k++ )
             if( fg.factor(I)[k] != (Real)0 )
                 os << k << " " << setw(os.precision()+4) << fg.factor(I)[k] << endl;
     }
@@ -320,7 +320,7 @@ void FactorGraph::clampVar( size_t i, const vector<size_t> &is, bool backup ) {
 
 
 void FactorGraph::clampFactor( size_t I, const vector<size_t> &is, bool backup ) {
-    size_t st = factor(I).states();
+    size_t st = factor(I).nrStates();
     Factor newF( factor(I).vars(), (Real)0 );
 
     foreach( size_t i, is ) {
index 1e96b98..d1a0597 100644 (file)
@@ -72,7 +72,7 @@ void Gibbs::construct() {
     _factor_counts.clear();
     _factor_counts.reserve( nrFactors() );
     for( size_t I = 0; I < nrFactors(); I++ )
-        _factor_counts.push_back( _count_t( factor(I).states(), 0 ) );
+        _factor_counts.push_back( _count_t( factor(I).nrStates(), 0 ) );
 
     _sample_count = 0;
 
index 256987f..1b49b8b 100644 (file)
@@ -44,7 +44,7 @@ mxArray *Factors2mx(const vector<Factor> &Ps) {
 
         mxArray *BiP = mxCreateNumericArray(I->vars().size(), &(*(dims.begin())), mxDOUBLE_CLASS, mxREAL);
         double *BiP_data = mxGetPr(BiP);
-        for( size_t j = 0; j < I->states(); j++ )
+        for( size_t j = 0; j < I->nrStates(); j++ )
             BiP_data[j] = (*I)[j];
 
         mxSetField(Bi,0,"Member",BiMember);
index c7b83de..989e3db 100644 (file)
@@ -61,7 +61,7 @@ SprinklerNetwork = dai.FactorGraph(SprinklerFactors);
 
 % Write factorgraph to a file
 SprinklerNetwork.WriteToFile('sprinkler.fg');
-fprintf('Sprinkler network written to sprinkler.fg');
+fprintf('Sprinkler network written to sprinkler.fg\n');
 
 % Output some information about the factorgraph
 fprintf('%d variables\n', SprinklerNetwork.nrVars());
index cb0441a..3beaf4f 100644 (file)
@@ -31,18 +31,18 @@ const double tol = 1e-8;
 BOOST_AUTO_TEST_CASE( ConstructorsTest ) {
     // check constructors
     Factor x1;
-    BOOST_CHECK_EQUAL( x1.states(), 1 );
+    BOOST_CHECK_EQUAL( x1.nrStates(), 1 );
     BOOST_CHECK( x1.p() == Prob( 1, 1.0 ) );
     BOOST_CHECK( x1.vars() == VarSet() );
 
     Factor x2( 5.0 );
-    BOOST_CHECK_EQUAL( x2.states(), 1 );
+    BOOST_CHECK_EQUAL( x2.nrStates(), 1 );
     BOOST_CHECK( x2.p() == Prob( 1, 5.0 ) );
     BOOST_CHECK( x2.vars() == VarSet() );
 
     Var v1( 0, 3 );
     Factor x3( v1 );
-    BOOST_CHECK_EQUAL( x3.states(), 3 );
+    BOOST_CHECK_EQUAL( x3.nrStates(), 3 );
     BOOST_CHECK( x3.p() == Prob( 3, 1.0 / 3.0 ) );
     BOOST_CHECK( x3.vars() == VarSet( v1 ) );
     BOOST_CHECK_EQUAL( x3[0], 1.0 / 3.0 );
@@ -51,14 +51,14 @@ BOOST_AUTO_TEST_CASE( ConstructorsTest ) {
 
     Var v2( 1, 2 );
     Factor x4( VarSet( v1, v2 ) );
-    BOOST_CHECK_EQUAL( x4.states(), 6 );
+    BOOST_CHECK_EQUAL( x4.nrStates(), 6 );
     BOOST_CHECK( x4.p() == Prob( 6, 1.0 / 6.0 ) );
     BOOST_CHECK( x4.vars() == VarSet( v1, v2 ) );
     for( size_t i = 0; i < 6; i++ )
         BOOST_CHECK_EQUAL( x4[i], 1.0 / 6.0 );
 
     Factor x5( VarSet( v1, v2 ), 1.0 );
-    BOOST_CHECK_EQUAL( x5.states(), 6 );
+    BOOST_CHECK_EQUAL( x5.nrStates(), 6 );
     BOOST_CHECK( x5.p() == Prob( 6, 1.0 ) );
     BOOST_CHECK( x5.vars() == VarSet( v1, v2 ) );
     for( size_t i = 0; i < 6; i++ )
@@ -68,7 +68,7 @@ BOOST_AUTO_TEST_CASE( ConstructorsTest ) {
     for( size_t i = 0; i < 6; i++ )
         x[i] = 10.0 - i;
     Factor x6( VarSet( v1, v2 ), x );
-    BOOST_CHECK_EQUAL( x6.states(), 6 );
+    BOOST_CHECK_EQUAL( x6.nrStates(), 6 );
     BOOST_CHECK( x6.vars() == VarSet( v1, v2 ) );
     for( size_t i = 0; i < 6; i++ )
         BOOST_CHECK_EQUAL( x6[i], x[i] );
@@ -79,14 +79,14 @@ BOOST_AUTO_TEST_CASE( ConstructorsTest ) {
     x.resize( 6 );
     x[4] = 10.0 - 4; x[5] = 10.0 - 5;
     Factor x8( VarSet( v2, v1 ), &(x[0]) );
-    BOOST_CHECK_EQUAL( x8.states(), 6 );
+    BOOST_CHECK_EQUAL( x8.nrStates(), 6 );
     BOOST_CHECK( x8.vars() == VarSet( v1, v2 ) );
     for( size_t i = 0; i < 6; i++ )
         BOOST_CHECK_EQUAL( x8[i], x[i] );
 
     Prob xx( x );
     Factor x9( VarSet( v2, v1 ), xx );
-    BOOST_CHECK_EQUAL( x9.states(), 6 );
+    BOOST_CHECK_EQUAL( x9.nrStates(), 6 );
     BOOST_CHECK( x9.vars() == VarSet( v1, v2 ) );
     for( size_t i = 0; i < 6; i++ )
         BOOST_CHECK_EQUAL( x9[i], x[i] );
@@ -115,7 +115,7 @@ BOOST_AUTO_TEST_CASE( ConstructorsTest ) {
     vars.push_back( v8 );
     vars.push_back( v7 );
     Factor x11( vars, w );
-    BOOST_CHECK_EQUAL( x11.states(), 12 );
+    BOOST_CHECK_EQUAL( x11.nrStates(), 12 );
     BOOST_CHECK( x11.vars() == VarSet( vars.begin(), vars.end() ) );
     BOOST_CHECK_EQUAL( x11[0], 0.1 );
     BOOST_CHECK_EQUAL( x11[1], 3.5 );
@@ -140,7 +140,7 @@ BOOST_AUTO_TEST_CASE( ConstructorsTest ) {
 
 BOOST_AUTO_TEST_CASE( QueriesTest ) {
     Factor x( Var( 5, 5 ), 0.0 );
-    for( size_t i = 0; i < x.states(); i++ )
+    for( size_t i = 0; i < x.nrStates(); i++ )
         x.set( i, 2.0 - i );
 
     // test min, max, sum, sumAbs, maxAbs
@@ -200,67 +200,48 @@ BOOST_AUTO_TEST_CASE( QueriesTest ) {
     BOOST_CHECK( a == b );
 }
 
-/*
+
 BOOST_AUTO_TEST_CASE( UnaryTransformationsTest ) {
-    Prob x( 3 );
+    Var v( 0, 3 );
+    Factor x( v );
     x.set( 0, -2.0 );
     x.set( 1, 0.0 );
     x.set( 2, 2.0 );
 
-    Prob y = -x;
-    Prob z = x.pwUnaryTr( std::negate<Real>() );
+    Factor y = -x;
     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.set( 0, 0.0 );
-    z.set( 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.set( 0, 0.0 );
-    z.set( 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.set( 0, 2.0 );
     y = x.normalized();
@@ -282,15 +263,16 @@ BOOST_AUTO_TEST_CASE( UnaryTransformationsTest ) {
 
 
 BOOST_AUTO_TEST_CASE( UnaryOperationsTest ) {
-    Prob xorg(3);
+    Var v( 0, 3 );
+    Factor xorg( v );
     xorg.set( 0, 2.0 );
     xorg.set( 1, 0.0 );
     xorg.set( 2, 1.0 );
-    Prob y(3);
+    Factor y( v );
 
-    Prob x = xorg;
-    BOOST_CHECK( x.setUniform() == Prob(3) );
-    BOOST_CHECK( x == Prob(3) );
+    Factor x = xorg;
+    BOOST_CHECK( x.setUniform() == Factor( v ) );
+    BOOST_CHECK( x == Factor( v ) );
 
     y.set( 0, std::exp(2.0) );
     y.set( 1, 1.0 );
@@ -298,9 +280,6 @@ BOOST_AUTO_TEST_CASE( UnaryOperationsTest ) {
     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.set( 0, std::log(2.0) );
     y.set( 1, -INFINITY );
@@ -311,17 +290,11 @@ BOOST_AUTO_TEST_CASE( UnaryOperationsTest ) {
     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.set( 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.set( 0, 2.0 / 3.0 );
     y.set( 1, 0.0 / 3.0 );
@@ -351,7 +324,7 @@ BOOST_AUTO_TEST_CASE( UnaryOperationsTest ) {
 
     for( size_t repeat = 0; repeat < 10000; repeat++ ) {
         x.randomize();
-        for( size_t i = 0; i < x.size(); i++ ) {
+        for( size_t i = 0; i < x.nrStates(); i++ ) {
             BOOST_CHECK( x[i] < 1.0 );
             BOOST_CHECK( x[i] >= 0.0 );
         }
@@ -359,62 +332,21 @@ BOOST_AUTO_TEST_CASE( UnaryOperationsTest ) {
 }
 
 
-BOOST_AUTO_TEST_CASE( ScalarTransformationsTest ) {
-    Prob x(3);
-    x.set( 0, 2.0 );
-    x.set( 1, 0.0 );
-    x.set( 2, 1.0 );
-    Prob y(3);
-
-    y.set( 0, 3.0 ); y.set( 1, 1.0 ); y.set( 2, 2.0 );
-    BOOST_CHECK( (x + 1.0) == y );
-    y.set( 0, 0.0 ); y.set( 1, -2.0 ); y.set( 2, -1.0 );
-    BOOST_CHECK( (x + (-2.0)) == y );
-
-    y.set( 0, 1.0 ); y.set( 1, -1.0 ); y.set( 2, 0.0 );
-    BOOST_CHECK( (x - 1.0) == y );
-    y.set( 0, 4.0 ); y.set( 1, 2.0 ); y.set( 2, 3.0 );
-    BOOST_CHECK( (x - (-2.0)) == y );
-
-    BOOST_CHECK( (x * 1.0) == x );
-    y.set( 0, 4.0 ); y.set( 1, 0.0 ); y.set( 2, 2.0 );
-    BOOST_CHECK( (x * 2.0) == y );
-    y.set( 0, -1.0 ); y.set( 1, 0.0 ); y.set( 2, -0.5 );
-    BOOST_CHECK( (x * -0.5) == y );
-
-    BOOST_CHECK( (x / 1.0) == x );
-    y.set( 0, 1.0 ); y.set( 1, 0.0 ); y.set( 2, 0.5 );
-    BOOST_CHECK( (x / 2.0) == y );
-    y.set( 0, -4.0 ); y.set( 1, 0.0 ); y.set( 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.set( 0, 4.0 ); y.set( 1, 0.0 ); y.set( 2, 1.0 );
-    BOOST_CHECK( (x ^ 2.0) == y );
-    y.set( 0, 1.0 / std::sqrt(2.0) ); y.set( 1, INFINITY ); y.set( 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( ScalarOperationsTest ) {
-    Prob xorg(3), x(3);
+    Var v( 0, 3 );
+    Factor xorg( v ), x( v );
     xorg.set( 0, 2.0 );
     xorg.set( 1, 0.0 );
     xorg.set( 2, 1.0 );
-    Prob y(3);
+    Factor y( v );
 
     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) );
+    BOOST_CHECK( x.fill( 1.0 ) == Factor(v, 1.0) );
+    BOOST_CHECK( x == Factor(v, 1.0) );
+    BOOST_CHECK( x.fill( 2.0 ) == Factor(v, 2.0) );
+    BOOST_CHECK( x == Factor(v, 2.0) );
+    BOOST_CHECK( x.fill( 0.0 ) == Factor(v, 0.0) );
+    BOOST_CHECK( x == Factor(v, 0.0) );
 
     x = xorg;
     y.set( 0, 3.0 ); y.set( 1, 1.0 ); y.set( 2, 2.0 );
@@ -451,14 +383,14 @@ BOOST_AUTO_TEST_CASE( ScalarOperationsTest ) {
     y.set( 0, -4.0 ); y.set( 1, 0.0 ); y.set( 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) );
+    BOOST_CHECK( (x /= 0.0) == Factor(v, 0.0) );
+    BOOST_CHECK( x == Factor(v, 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) );
+    BOOST_CHECK( (x ^= 0.0) == Factor(v, 1.0) );
+    BOOST_CHECK( x == Factor(v, 1.0) );
     x = xorg;
     y.set( 0, 4.0 ); y.set( 1, 0.0 ); y.set( 2, 1.0 );
     BOOST_CHECK( (x ^= 2.0) == y );
@@ -469,13 +401,57 @@ BOOST_AUTO_TEST_CASE( ScalarOperationsTest ) {
 }
 
 
-BOOST_AUTO_TEST_CASE( VectorOperationsTest ) {
+BOOST_AUTO_TEST_CASE( ScalarTransformationsTest ) {
+    Var v( 0, 3 );
+    Factor x( v );
+    x.set( 0, 2.0 );
+    x.set( 1, 0.0 );
+    x.set( 2, 1.0 );
+    Factor y( v );
+
+    y.set( 0, 3.0 ); y.set( 1, 1.0 ); y.set( 2, 2.0 );
+    BOOST_CHECK( (x + 1.0) == y );
+    y.set( 0, 0.0 ); y.set( 1, -2.0 ); y.set( 2, -1.0 );
+    BOOST_CHECK( (x + (-2.0)) == y );
+
+    y.set( 0, 1.0 ); y.set( 1, -1.0 ); y.set( 2, 0.0 );
+    BOOST_CHECK( (x - 1.0) == y );
+    y.set( 0, 4.0 ); y.set( 1, 2.0 ); y.set( 2, 3.0 );
+    BOOST_CHECK( (x - (-2.0)) == y );
+
+    BOOST_CHECK( (x * 1.0) == x );
+    y.set( 0, 4.0 ); y.set( 1, 0.0 ); y.set( 2, 2.0 );
+    BOOST_CHECK( (x * 2.0) == y );
+    y.set( 0, -1.0 ); y.set( 1, 0.0 ); y.set( 2, -0.5 );
+    BOOST_CHECK( (x * -0.5) == y );
+
+    BOOST_CHECK( (x / 1.0) == x );
+    y.set( 0, 1.0 ); y.set( 1, 0.0 ); y.set( 2, 0.5 );
+    BOOST_CHECK( (x / 2.0) == y );
+    y.set( 0, -4.0 ); y.set( 1, 0.0 ); y.set( 2, -2.0 );
+    BOOST_CHECK( (x / -0.5) == y );
+    BOOST_CHECK( (x / 0.0) == Factor(v, 0.0) );
+
+    BOOST_CHECK( (x ^ 1.0) == x );
+    BOOST_CHECK( (x ^ 0.0) == Factor(v, 1.0) );
+    y.set( 0, 4.0 ); y.set( 1, 0.0 ); y.set( 2, 1.0 );
+    BOOST_CHECK( (x ^ 2.0) == y );
+    y.set( 0, 1.0 / std::sqrt(2.0) ); y.set( 1, INFINITY ); y.set( 2, 1.0 );
+    Factor 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( SimilarFactorOperationsTest ) {
     size_t N = 6;
-    Prob xorg(N), x(N);
+    Var v( 0, N );
+    Factor xorg( v ), x( v );
     xorg.set( 0, 2.0 ); xorg.set( 1, 0.0 ); xorg.set( 2, 1.0 ); xorg.set( 3, 0.0 ); xorg.set( 4, 2.0 ); xorg.set( 5, 3.0 );
-    Prob y(N);
+    Factor y( v );
     y.set( 0, 0.5 ); y.set( 1, -1.0 ); y.set( 2, 0.0 ); y.set( 3, 0.0 ); y.set( 4, -2.0 ); y.set( 5, 3.0 );
-    Prob z(N), r(N);
+    Factor z( v ), r( v );
 
     z.set( 0, 2.5 ); z.set( 1, -1.0 ); z.set( 2, 1.0 ); z.set( 3, 0.0 ); z.set( 4, 0.0 ); z.set( 5, 6.0 );
     x = xorg;
@@ -484,7 +460,7 @@ BOOST_AUTO_TEST_CASE( VectorOperationsTest ) {
         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.binaryOp( y, std::plus<Real>() ) == z );
     BOOST_CHECK( x == z );
 
     z.set( 0, 1.5 ); z.set( 1, 1.0 ); z.set( 2, 1.0 ); z.set( 3, 0.0 ); z.set( 4, 4.0 ); z.set( 5, 0.0 );
@@ -494,7 +470,7 @@ BOOST_AUTO_TEST_CASE( VectorOperationsTest ) {
         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.binaryOp( y, std::minus<Real>() ) == z );
     BOOST_CHECK( x == z );
 
     z.set( 0, 1.0 ); z.set( 1, 0.0 ); z.set( 2, 0.0 ); z.set( 3, 0.0 ); z.set( 4, -4.0 ); z.set( 5, 9.0 );
@@ -504,7 +480,7 @@ BOOST_AUTO_TEST_CASE( VectorOperationsTest ) {
         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.binaryOp( y, std::multiplies<Real>() ) == z );
     BOOST_CHECK( x == z );
 
     z.set( 0, 4.0 ); z.set( 1, 0.0 ); z.set( 2, 0.0 ); z.set( 3, 0.0 ); z.set( 4, -1.0 ); z.set( 5, 1.0 );
@@ -514,118 +490,314 @@ BOOST_AUTO_TEST_CASE( VectorOperationsTest ) {
         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.set( 0, 4.0 ); z.set( 1, 0.0 ); z.set( 2, INFINITY ); 
-    // z.set( 3, INFINITY );
-    z.set( 4, -1.0 ); z.set( 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.set( 3, 0.0 ); r.set( 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.set( 3, 0.0 ); r.set( 3, 0.0 );
-    BOOST_CHECK( x == r );
-
-    z.set( 0, std::sqrt(2.0) ); z.set( 1, INFINITY ); z.set( 2, 1.0 ); z.set( 3, 1.0 ); z.set( 4, 0.25 ); z.set( 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.binaryOp( y, fo_divides0<Real>() ) == z );
     BOOST_CHECK( x == z );
 }
 
 
-BOOST_AUTO_TEST_CASE( VectorTransformationsTest ) {
+BOOST_AUTO_TEST_CASE( SimilarFactorTransformationsTest ) {
     size_t N = 6;
-    Prob x(N);
+    Var v( 0, N );
+    Factor x( v );
     x.set( 0, 2.0 ); x.set( 1, 0.0 ); x.set( 2, 1.0 ); x.set( 3, 0.0 ); x.set( 4, 2.0 ); x.set( 5, 3.0 );
-    Prob y(N);
+    Factor y( v );
     y.set( 0, 0.5 ); y.set( 1, -1.0 ); y.set( 2, 0.0 ); y.set( 3, 0.0 ); y.set( 4, -2.0 ); y.set( 5, 3.0 );
-    Prob z(N), r(N);
+    Factor z( v ), r( v );
 
     z.set( 0, 2.5 ); z.set( 1, -1.0 ); z.set( 2, 1.0 ); z.set( 3, 0.0 ); z.set( 4, 0.0 ); z.set( 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>() );
+    z = x.binaryTr( y, std::plus<Real>() );
     BOOST_CHECK( r == z );
 
     z.set( 0, 1.5 ); z.set( 1, 1.0 ); z.set( 2, 1.0 ); z.set( 3, 0.0 ); z.set( 4, 4.0 ); z.set( 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>() );
+    z = x.binaryTr( y, std::minus<Real>() );
     BOOST_CHECK( r == z );
 
     z.set( 0, 1.0 ); z.set( 1, 0.0 ); z.set( 2, 0.0 ); z.set( 3, 0.0 ); z.set( 4, -4.0 ); z.set( 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>() );
+    z = x.binaryTr( y, std::multiplies<Real>() );
     BOOST_CHECK( r == z );
 
     z.set( 0, 4.0 ); z.set( 1, 0.0 ); z.set( 2, 0.0 ); z.set( 3, 0.0 ); z.set( 4, -1.0 ); z.set( 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>() );
+    z = x.binaryTr( y, fo_divides0<Real>() );
     BOOST_CHECK( r == z );
+}
 
-    z.set( 0, 4.0 ); z.set( 1, 0.0 ); z.set( 2, INFINITY ); 
-    // z.set( 3, INFINITY );
-    z.set( 4, -1.0 ); z.set( 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.set( 0, std::sqrt(2.0) ); z.set( 1, INFINITY ); z.set( 2, 1.0 ); z.set( 3, 1.0 ); z.set( 4, 0.25 ); z.set( 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_AUTO_TEST_CASE( FactorOperationsTest ) {
+    size_t N = 9;
+    Var v1( 1, 3 );
+    Var v2( 2, 3 );
+    Factor xorg( v1 ), x( v1 );
+    xorg.set( 0, 2.0 ); xorg.set( 1, 0.0 ); xorg.set( 2, -1.0 );
+    Factor y( v2 );
+    y.set( 0, 0.5 ); y.set( 1, -1.0 ); y.set( 2, 0.0 );
+    Factor r;
+
+    Factor z( VarSet( v1, v2 ) );
+    z.set( 0,  2.5 ); z.set( 1,  0.5 ); z.set( 2, -0.5 ); 
+    z.set( 3,  1.0 ); z.set( 4, -1.0 ); z.set( 5, -2.0 );
+    z.set( 6,  2.0 ); z.set( 7,  0.0 ); z.set( 8, -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.binaryOp( y, std::plus<Real>() ) == z );
+    BOOST_CHECK( x == z );
+
+    z.set( 0,  1.5 ); z.set( 1, -0.5 ); z.set( 2, -1.5 ); 
+    z.set( 3,  3.0 ); z.set( 4,  1.0 ); z.set( 5,  0.0 );
+    z.set( 6,  2.0 ); z.set( 7,  0.0 ); z.set( 8, -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.binaryOp( y, std::minus<Real>() ) == z );
+    BOOST_CHECK( x == z );
+
+    z.set( 0,  1.0 ); z.set( 1,  0.0 ); z.set( 2, -0.5 ); 
+    z.set( 3, -2.0 ); z.set( 4,  0.0 ); z.set( 5,  1.0 );
+    z.set( 6,  0.0 ); z.set( 7,  0.0 ); z.set( 8,  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.binaryOp( y, std::multiplies<Real>() ) == z );
+    BOOST_CHECK( x == z );
+
+    z.set( 0,  4.0 ); z.set( 1,  0.0 ); z.set( 2, -2.0 ); 
+    z.set( 3, -2.0 ); z.set( 4,  0.0 ); z.set( 5,  1.0 );
+    z.set( 6,  0.0 ); z.set( 7,  0.0 ); z.set( 8,  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.binaryOp( y, fo_divides0<Real>() ) == z );
+    BOOST_CHECK( x == z );
+}
+
+
+BOOST_AUTO_TEST_CASE( FactorTransformationsTest ) {
+    size_t N = 9;
+    Var v1( 1, 3 );
+    Var v2( 2, 3 );
+    Factor x( v1 );
+    x.set( 0, 2.0 ); x.set( 1, 0.0 ); x.set( 2, -1.0 );
+    Factor y( v2 );
+    y.set( 0, 0.5 ); y.set( 1, -1.0 ); y.set( 2, 0.0 );
+    Factor r;
+
+    Factor z( VarSet( v1, v2 ) );
+    z.set( 0,  2.5 ); z.set( 1,  0.5 ); z.set( 2, -0.5 ); 
+    z.set( 3,  1.0 ); z.set( 4, -1.0 ); z.set( 5, -2.0 );
+    z.set( 6,  2.0 ); z.set( 7,  0.0 ); z.set( 8, -1.0 );
+    r = x + y;
+    for( size_t i = 0; i < N; i++ )
+        BOOST_CHECK_CLOSE( r[i], z[i], tol );
+    BOOST_CHECK( r == z );
+    z = x.binaryTr( y, std::plus<Real>() );
+    BOOST_CHECK( r == z );
+
+    z.set( 0,  1.5 ); z.set( 1, -0.5 ); z.set( 2, -1.5 ); 
+    z.set( 3,  3.0 ); z.set( 4,  1.0 ); z.set( 5,  0.0 );
+    z.set( 6,  2.0 ); z.set( 7,  0.0 ); z.set( 8, -1.0 );
+    r = x - y;
+    for( size_t i = 0; i < N; i++ )
+        BOOST_CHECK_CLOSE( r[i], z[i], tol );
+    BOOST_CHECK( r == z );
+    z = x.binaryTr( y, std::minus<Real>() );
+    BOOST_CHECK( r == z );
+
+    z.set( 0,  1.0 ); z.set( 1,  0.0 ); z.set( 2, -0.5 ); 
+    z.set( 3, -2.0 ); z.set( 4,  0.0 ); z.set( 5,  1.0 );
+    z.set( 6,  0.0 ); z.set( 7,  0.0 ); z.set( 8,  0.0 );
+    r = x * y;
+    for( size_t i = 0; i < N; i++ )
+        BOOST_CHECK_CLOSE( r[i], z[i], tol );
+    BOOST_CHECK( r == z );
+    z = x.binaryTr( y, std::multiplies<Real>() );
+    BOOST_CHECK( r == z );
+
+    z.set( 0,  4.0 ); z.set( 1,  0.0 ); z.set( 2, -2.0 ); 
+    z.set( 3, -2.0 ); z.set( 4,  0.0 ); z.set( 5,  1.0 );
+    z.set( 6,  0.0 ); z.set( 7,  0.0 ); z.set( 8,  0.0 );
+    r = x / y;
+    for( size_t i = 0; i < N; i++ )
+        BOOST_CHECK_CLOSE( r[i], z[i], tol );
+    BOOST_CHECK( r == z );
+    z = x.binaryOp( y, fo_divides0<Real>() );
     BOOST_CHECK( r == z );
 }
 
 
+BOOST_AUTO_TEST_CASE( MiscOperationsTest ) {
+    Var v1(1, 2);
+    Var v2(2, 3);
+    Factor x( VarSet( v1, v2 ) );
+    x.randomize();
+
+    // slice
+    Factor y = x.slice( v1, 0 );
+    BOOST_CHECK( y.vars() == VarSet( v2 ) );
+    BOOST_CHECK_EQUAL( y.nrStates(), 3 );
+    BOOST_CHECK_EQUAL( y[0], x[0] );
+    BOOST_CHECK_EQUAL( y[1], x[2] );
+    BOOST_CHECK_EQUAL( y[2], x[4] );
+    y = x.slice( v1, 1 );
+    BOOST_CHECK( y.vars() == VarSet( v2 ) );
+    BOOST_CHECK_EQUAL( y.nrStates(), 3 );
+    BOOST_CHECK_EQUAL( y[0], x[1] );
+    BOOST_CHECK_EQUAL( y[1], x[3] );
+    BOOST_CHECK_EQUAL( y[2], x[5] );
+    y = x.slice( v2, 0 );
+    BOOST_CHECK( y.vars() == VarSet( v1 ) );
+    BOOST_CHECK_EQUAL( y.nrStates(), 2 );
+    BOOST_CHECK_EQUAL( y[0], x[0] );
+    BOOST_CHECK_EQUAL( y[1], x[1] );
+    y = x.slice( v2, 1 );
+    BOOST_CHECK( y.vars() == VarSet( v1 ) );
+    BOOST_CHECK_EQUAL( y.nrStates(), 2 );
+    BOOST_CHECK_EQUAL( y[0], x[2] );
+    BOOST_CHECK_EQUAL( y[1], x[3] );
+    y = x.slice( v2, 2 );
+    BOOST_CHECK( y.vars() == VarSet( v1 ) );
+    BOOST_CHECK_EQUAL( y.nrStates(), 2 );
+    BOOST_CHECK_EQUAL( y[0], x[4] );
+    BOOST_CHECK_EQUAL( y[1], x[5] );
+    for( size_t i = 0; i < x.nrStates(); i++ ) {
+        y = x.slice( VarSet( v1, v2 ), 0 );
+        BOOST_CHECK( y.vars() == VarSet() );
+        BOOST_CHECK_EQUAL( y.nrStates(), 1 );
+        BOOST_CHECK_EQUAL( y[0], x[0] );
+    }
+    y = x.slice( VarSet(), 0 );
+    BOOST_CHECK_EQUAL( y, x );
+
+    // embed
+    Var v3(3, 4);
+    BOOST_CHECK_THROW( x.embed( VarSet( v3 ) ), Exception );
+    BOOST_CHECK_THROW( x.embed( VarSet( v3, v2 ) ), Exception );
+    y = x.embed( VarSet( v3, v2 ) | v1 );
+    for( size_t i = 0; i < y.nrStates(); i++ )
+        BOOST_CHECK_EQUAL( y[i], x[i % 6] );
+    y = x.embed( VarSet( v1, v2 ) );
+    BOOST_CHECK_EQUAL( x, y );
+
+    // marginal
+    y = x.marginal( v1 );
+    BOOST_CHECK( y.vars() == VarSet( v1 ) );
+    BOOST_CHECK_EQUAL( y[0], (x[0] + x[2] + x[4]) / x.sum() );
+    BOOST_CHECK_EQUAL( y[1], (x[1] + x[3] + x[5]) / x.sum() );
+    y = x.marginal( v2 );
+    BOOST_CHECK( y.vars() == VarSet( v2 ) );
+    BOOST_CHECK_EQUAL( y[0], (x[0] + x[1]) / x.sum() );
+    BOOST_CHECK_EQUAL( y[1], (x[2] + x[3]) / x.sum() );
+    BOOST_CHECK_EQUAL( y[2], (x[4] + x[5]) / x.sum() );
+    y = x.marginal( VarSet() );
+    BOOST_CHECK( y.vars() == VarSet() );
+    BOOST_CHECK_EQUAL( y[0], 1.0 );
+    y = x.marginal( VarSet( v1, v2 ) );
+    BOOST_CHECK( y == x.normalized() );
+
+    y = x.marginal( v1, true );
+    BOOST_CHECK( y.vars() == VarSet( v1 ) );
+    BOOST_CHECK_EQUAL( y[0], (x[0] + x[2] + x[4]) / x.sum() );
+    BOOST_CHECK_EQUAL( y[1], (x[1] + x[3] + x[5]) / x.sum() );
+    y = x.marginal( v2, true );
+    BOOST_CHECK( y.vars() == VarSet( v2 ) );
+    BOOST_CHECK_EQUAL( y[0], (x[0] + x[1]) / x.sum() );
+    BOOST_CHECK_EQUAL( y[1], (x[2] + x[3]) / x.sum() );
+    BOOST_CHECK_EQUAL( y[2], (x[4] + x[5]) / x.sum() );
+    y = x.marginal( VarSet(), true );
+    BOOST_CHECK( y.vars() == VarSet() );
+    BOOST_CHECK_EQUAL( y[0], 1.0 );
+    y = x.marginal( VarSet( v1, v2 ), true );
+    BOOST_CHECK( y == x.normalized() );
+
+    y = x.marginal( v1, false );
+    BOOST_CHECK( y.vars() == VarSet( v1 ) );
+    BOOST_CHECK_EQUAL( y[0], x[0] + x[2] + x[4] );
+    BOOST_CHECK_EQUAL( y[1], x[1] + x[3] + x[5] );
+    y = x.marginal( v2, false );
+    BOOST_CHECK( y.vars() == VarSet( v2 ) );
+    BOOST_CHECK_EQUAL( y[0], x[0] + x[1] );
+    BOOST_CHECK_EQUAL( y[1], x[2] + x[3] );
+    BOOST_CHECK_EQUAL( y[2], x[4] + x[5] );
+    y = x.marginal( VarSet(), false );
+    BOOST_CHECK( y.vars() == VarSet() );
+    BOOST_CHECK_EQUAL( y[0], x.sum() );
+    y = x.marginal( VarSet( v1, v2 ), false );
+    BOOST_CHECK( y == x );
+
+    // maxMarginal
+    y = x.maxMarginal( v1 );
+    BOOST_CHECK( y.vars() == VarSet( v1 ) );
+    BOOST_CHECK_EQUAL( y[0], x.slice( v1, 0 ).max() / (x.slice( v1, 0 ).max() + x.slice( v1, 1 ).max()) );
+    BOOST_CHECK_EQUAL( y[1], x.slice( v1, 1 ).max() / (x.slice( v1, 0 ).max() + x.slice( v1, 1 ).max()) );
+    y = x.maxMarginal( v2 );
+    BOOST_CHECK( y.vars() == VarSet( v2 ) );
+    BOOST_CHECK_EQUAL( y[0], x.slice( v2, 0 ).max() / (x.slice( v2, 0 ).max() + x.slice( v2, 1 ).max() + x.slice( v2, 2 ).max()) );
+    BOOST_CHECK_EQUAL( y[1], x.slice( v2, 1 ).max() / (x.slice( v2, 0 ).max() + x.slice( v2, 1 ).max() + x.slice( v2, 2 ).max()) );
+    BOOST_CHECK_EQUAL( y[2], x.slice( v2, 2 ).max() / (x.slice( v2, 0 ).max() + x.slice( v2, 1 ).max() + x.slice( v2, 2 ).max()) );
+    y = x.maxMarginal( VarSet() );
+    BOOST_CHECK( y.vars() == VarSet() );
+    BOOST_CHECK_EQUAL( y[0], 1.0 );
+    y = x.maxMarginal( VarSet( v1, v2 ) );
+    BOOST_CHECK( y == x.normalized() );
+
+    y = x.maxMarginal( v1, true );
+    BOOST_CHECK( y.vars() == VarSet( v1 ) );
+    BOOST_CHECK_EQUAL( y[0], x.slice( v1, 0 ).max() / (x.slice( v1, 0 ).max() + x.slice( v1, 1 ).max()) );
+    BOOST_CHECK_EQUAL( y[1], x.slice( v1, 1 ).max() / (x.slice( v1, 0 ).max() + x.slice( v1, 1 ).max()) );
+    y = x.maxMarginal( v2, true );
+    BOOST_CHECK( y.vars() == VarSet( v2 ) );
+    BOOST_CHECK_EQUAL( y[0], x.slice( v2, 0 ).max() / (x.slice( v2, 0 ).max() + x.slice( v2, 1 ).max() + x.slice( v2, 2 ).max()) );
+    BOOST_CHECK_EQUAL( y[1], x.slice( v2, 1 ).max() / (x.slice( v2, 0 ).max() + x.slice( v2, 1 ).max() + x.slice( v2, 2 ).max()) );
+    BOOST_CHECK_EQUAL( y[2], x.slice( v2, 2 ).max() / (x.slice( v2, 0 ).max() + x.slice( v2, 1 ).max() + x.slice( v2, 2 ).max()) );
+    y = x.maxMarginal( VarSet(), true );
+    BOOST_CHECK( y.vars() == VarSet() );
+    BOOST_CHECK_EQUAL( y[0], 1.0 );
+    y = x.maxMarginal( VarSet( v1, v2 ), true );
+    BOOST_CHECK( y == x.normalized() );
+
+    y = x.maxMarginal( v1, false );
+    BOOST_CHECK( y.vars() == VarSet( v1 ) );
+    BOOST_CHECK_EQUAL( y[0], x.slice( v1, 0 ).max() );
+    BOOST_CHECK_EQUAL( y[1], x.slice( v1, 1 ).max() );
+    y = x.maxMarginal( v2, false );
+    BOOST_CHECK( y.vars() == VarSet( v2 ) );
+    BOOST_CHECK_EQUAL( y[0], x.slice( v2, 0 ).max() );
+    BOOST_CHECK_EQUAL( y[1], x.slice( v2, 1 ).max() );
+    BOOST_CHECK_EQUAL( y[2], x.slice( v2, 2 ).max() );
+    y = x.maxMarginal( VarSet(), false );
+    BOOST_CHECK( y.vars() == VarSet() );
+    BOOST_CHECK_EQUAL( y[0], x.max() );
+    y = x.maxMarginal( VarSet( v1, v2 ), false );
+    BOOST_CHECK( y == x );
+}
+
+
 BOOST_AUTO_TEST_CASE( RelatedFunctionsTest ) {
-    Prob x(3), y(3), z(3);
+    Var v( 0, 3 );
+    Factor x(v), y(v), z(v);
     x.set( 0, 0.2 );
     x.set( 1, 0.8 );
     x.set( 2, 0.0 );
@@ -646,48 +818,36 @@ BOOST_AUTO_TEST_CASE( RelatedFunctionsTest ) {
     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.set( 1, 0.7 ); x.set( 2, 0.1 );
     y.set( 0, 0.1 ); y.set( 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)") );
+    BOOST_CHECK_EQUAL( s, std::string("({x0}, (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)") );
+    BOOST_CHECK_EQUAL( s, std::string("({x0}, (0.1, 0.5, 0.4))") );
 
     z = min( x, y );
     BOOST_CHECK_EQUAL( z[0], 0.1 );
@@ -698,11 +858,38 @@ BOOST_AUTO_TEST_CASE( RelatedFunctionsTest ) {
     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 );
+    for( double J = -1.0; J <= 1.01; J += 0.1 ) {
+        Factor x = createFactorIsing( Var(0,2), Var(1,2), J ).normalized();
+        BOOST_CHECK_CLOSE( x[0], std::exp(J) / (4.0 * std::cosh(J)), tol );
+        BOOST_CHECK_CLOSE( x[1], std::exp(-J) / (4.0 * std::cosh(J)), tol );
+        BOOST_CHECK_CLOSE( x[2], std::exp(-J) / (4.0 * std::cosh(J)), tol );
+        BOOST_CHECK_CLOSE( x[3], std::exp(J) / (4.0 * std::cosh(J)), tol );
+        BOOST_CHECK_SMALL( MutualInfo( x ) - (J * std::tanh(J) - std::log(std::cosh(J))), tol );
+    }
+    Var v1( 1, 3 );
+    Var v2( 2, 4 );
+    BOOST_CHECK_SMALL( MutualInfo( (Factor(v1).randomize() * Factor(v2).randomize()).normalized() ), tol );
+    BOOST_CHECK_THROW( MutualInfo( createFactorIsing( Var(0,2), 1.0 ).normalized() ), Exception );
+    BOOST_CHECK_THROW( createFactorIsing( v1, 0.0 ), Exception );
+    BOOST_CHECK_THROW( createFactorIsing( v1, v2, 0.0 ), Exception );
+    for( double J = -1.0; J <= 1.01; J += 0.1 ) {
+        Factor x = createFactorIsing( Var(0,2), J ).normalized();
+        BOOST_CHECK_CLOSE( x[0], std::exp(-J) / (2.0 * std::cosh(J)), tol );
+        BOOST_CHECK_CLOSE( x[1], std::exp(J) / (2.0 * std::cosh(J)), tol );
+        BOOST_CHECK_SMALL( x.entropy() - (-J * std::tanh(J) + std::log(2.0 * std::cosh(J))), tol );
+    }
+    
+    x = createFactorDelta( v1, 2 );
+    BOOST_CHECK_EQUAL( x[0], 0.0 );
+    BOOST_CHECK_EQUAL( x[1], 0.0 );
+    BOOST_CHECK_EQUAL( x[2], 1.0 );
+    x = createFactorDelta( v1, 1 );
+    BOOST_CHECK_EQUAL( x[0], 0.0 );
+    BOOST_CHECK_EQUAL( x[1], 1.0 );
+    BOOST_CHECK_EQUAL( x[2], 0.0 );
+    x = createFactorDelta( v1, 0 );
+    BOOST_CHECK_EQUAL( x[0], 1.0 );
+    BOOST_CHECK_EQUAL( x[1], 0.0 );
+    BOOST_CHECK_EQUAL( x[2], 0.0 );
+    BOOST_CHECK_THROW( createFactorDelta( v1, 4 ), Exception );
 }
-*/