Some small documentation updates
[libdai.git] / src / bbp.cpp
index eb03614..b4bd072 100644 (file)
@@ -1,21 +1,11 @@
-/*  Copyright (C) 2009  Frederik Eaton [frederik at ofb dot net]
-
-    This file is part of libDAI.
-
-    libDAI is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    libDAI is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with libDAI; if not, write to the Free Software
-    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
-*/
+/*  This file is part of libDAI - http://www.libdai.org/
+ *
+ *  libDAI is licensed under the terms of the GNU General Public License version
+ *  2, or (at your option) any later version. libDAI is distributed without any
+ *  warranty. See the file COPYING for more details.
+ *
+ *  Copyright (C) 2009  Frederik Eaton [frederik at ofb dot net]
+ */
 
 
 #include <dai/bp.h>
@@ -35,29 +25,6 @@ using namespace std;
 typedef BipartiteGraph::Neighbor Neighbor;
 
 
-Prob unnormAdjoint( const Prob &w, Real Z_w, const Prob &adj_w ) {
-    assert( w.size() == adj_w.size() );
-    Prob adj_w_unnorm( w.size(), 0.0 );
-    Real s = 0.0;
-    for( size_t i = 0; i < w.size(); i++ )
-        s += w[i] * adj_w[i];
-    for( size_t i = 0; i < w.size(); i++ )
-        adj_w_unnorm[i] = (adj_w[i] - s) / Z_w;
-    return adj_w_unnorm;
-//  THIS WOULD BE ABOUT 50% SLOWER:  return (adj_w - (w * adj_w).sum()) / Z_w;
-}
-
-
-std::vector<size_t> getGibbsState( const InfAlg &ia, size_t iters ) {
-    PropertySet gibbsProps;
-    gibbsProps.Set("iters", iters);
-    gibbsProps.Set("verbose", size_t(0));
-    Gibbs gibbs( ia.fg(), gibbsProps );
-    gibbs.run();
-    return gibbs.state();
-}
-
-
 /// Returns the entry of the I'th factor corresponding to a global state
 size_t getFactorEntryForState( const FactorGraph &fg, size_t I, const vector<size_t> &state ) {
     size_t f_entry = 0;
@@ -72,6 +39,91 @@ size_t getFactorEntryForState( const FactorGraph &fg, size_t I, const vector<siz
 }
 
 
+bool BBPCostFunction::needGibbsState() const {
+    switch( (size_t)(*this) ) {
+        case CFN_GIBBS_B:
+        case CFN_GIBBS_B2:
+        case CFN_GIBBS_EXP:
+        case CFN_GIBBS_B_FACTOR:
+        case CFN_GIBBS_B2_FACTOR:
+        case CFN_GIBBS_EXP_FACTOR:
+            return true;
+        default:
+            return false;
+    }
+}
+
+
+Real BBPCostFunction::evaluate( const InfAlg &ia, const vector<size_t> *stateP ) const {
+    Real cf = 0.0;
+    const FactorGraph &fg = ia.fg();
+
+    switch( (size_t)(*this) ) {
+        case CFN_BETHE_ENT: // ignores state
+            cf = -ia.logZ();
+            break;
+        case CFN_VAR_ENT: // ignores state
+            for( size_t i = 0; i < fg.nrVars(); i++ )
+                cf += -ia.beliefV(i).entropy();
+            break;
+        case CFN_FACTOR_ENT: // ignores state
+            for( size_t I = 0; I < fg.nrFactors(); I++ )
+                cf += -ia.beliefF(I).entropy();
+            break;
+        case CFN_GIBBS_B:
+        case CFN_GIBBS_B2:
+        case CFN_GIBBS_EXP: {
+            DAI_ASSERT( stateP != NULL );
+            vector<size_t> state = *stateP;
+            DAI_ASSERT( state.size() == fg.nrVars() );
+            for( size_t i = 0; i < fg.nrVars(); i++ ) {
+                Real b = ia.beliefV(i)[state[i]];
+                switch( (size_t)(*this) ) {
+                    case CFN_GIBBS_B:
+                        cf += b;
+                        break;
+                    case CFN_GIBBS_B2:
+                        cf += b * b / 2.0;
+                        break;
+                    case CFN_GIBBS_EXP:
+                        cf += exp( b );
+                        break;
+                    default:
+                        DAI_THROW(UNKNOWN_ENUM_VALUE);
+                }
+            }
+            break;
+        } case CFN_GIBBS_B_FACTOR:
+          case CFN_GIBBS_B2_FACTOR:
+          case CFN_GIBBS_EXP_FACTOR: {
+            DAI_ASSERT( stateP != NULL );
+            vector<size_t> state = *stateP;
+            DAI_ASSERT( state.size() == fg.nrVars() );
+            for( size_t I = 0; I < fg.nrFactors(); I++ ) {
+                size_t x_I = getFactorEntryForState( fg, I, state );
+                Real b = ia.beliefF(I)[x_I];
+                switch( (size_t)(*this) ) {
+                    case CFN_GIBBS_B_FACTOR:
+                        cf += b;
+                        break;
+                    case CFN_GIBBS_B2_FACTOR:
+                        cf += b * b / 2.0;
+                        break;
+                    case CFN_GIBBS_EXP_FACTOR:
+                        cf += exp( b );
+                        break;
+                    default:
+                        DAI_THROW(UNKNOWN_ENUM_VALUE);
+                }
+            }
+            break;
+        } default:
+            DAI_THROWE(UNKNOWN_ENUM_VALUE, "Unknown cost function " + std::string(*this));
+    }
+    return cf;
+}
+
+
 #define LOOP_ij(body) {                             \
     size_t i_states = _fg->var(i).states();         \
     size_t j_states = _fg->var(j).states();         \
@@ -105,8 +157,8 @@ 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() );
-            for( IndexFor k(_fg->var(i), _fg->factor(I).vars()); k >= 0; ++k )
+            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 );
         }
@@ -115,49 +167,49 @@ void BBP::RegenerateInds() {
 
 
 void BBP::RegenerateT() {
-    // _T[i][_I]
-    _T.resize( _fg->nrVars() );
+    // _Tmsg[i][_I]
+    _Tmsg.resize( _fg->nrVars() );
     for( size_t i = 0; i < _fg->nrVars(); i++ ) {
-        _T[i].resize( _fg->nbV(i).size() );
+        _Tmsg[i].resize( _fg->nbV(i).size() );
         foreach( const Neighbor &I, _fg->nbV(i) ) {
             Prob prod( _fg->var(i).states(), 1.0 );
             foreach( const Neighbor &J, _fg->nbV(i) )
                 if( J.node != I.node )
                     prod *= _bp_dual.msgM( i, J.iter );
-            _T[i][I.iter] = prod;
+            _Tmsg[i][I.iter] = prod;
         }
     }
 }
 
 
 void BBP::RegenerateU() {
-    // _U[I][_i]
-    _U.resize( _fg->nrFactors() );
+    // _Umsg[I][_i]
+    _Umsg.resize( _fg->nrFactors() );
     for( size_t I = 0; I < _fg->nrFactors(); I++ ) {
-        _U[I].resize( _fg->nbF(I).size() );
+        _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 ) );
                     const _ind_t &ind = _index( j, j.dual );
                     // multiply prod by n_jI
                     for( size_t x_I = 0; x_I < prod.size(); x_I++ )
-                        prod[x_I] *= n_jI[ind[x_I]];
+                        prod.set( x_I, prod[x_I] * n_jI[ind[x_I]] );
                 }
-            _U[I][i.iter] = prod;
+            _Umsg[I][i.iter] = prod;
         }
     }
 }
 
 
 void BBP::RegenerateS() {
-    // _S[i][_I][_j]
-    _S.resize( _fg->nrVars() );
+    // _Smsg[i][_I][_j]
+    _Smsg.resize( _fg->nrVars() );
     for( size_t i = 0; i < _fg->nrVars(); i++ ) {
-        _S[i].resize( _fg->nbV(i).size() );
+        _Smsg[i].resize( _fg->nbV(i).size() );
         foreach( const Neighbor &I, _fg->nbV(i) ) {
-            _S[i][I.iter].resize( _fg->nbF(I).size() );
+            _Smsg[i][I.iter].resize( _fg->nbF(I).size() );
             foreach( const Neighbor &j, _fg->nbF(I) )
                 if( i != j ) {
                     Factor prod( _fg->factor(I) );
@@ -165,14 +217,14 @@ 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++ )
-                                prod.p()[x_I] *= p[ind[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]] );
                         }
                     }
                     // "Marginalize" onto i|j (unnormalized)
                     Prob marg;
                     marg = prod.marginal( VarSet(_fg->var(i), _fg->var(j)), false ).p();
-                    _S[i][I.iter][j.iter] = marg;
+                    _Smsg[i][I.iter][j.iter] = marg;
                 }
         }
     }
@@ -180,19 +232,19 @@ void BBP::RegenerateS() {
 
 
 void BBP::RegenerateR() {
-    // _R[I][_i][_J]
-    _R.resize( _fg->nrFactors() );
+    // _Rmsg[I][_i][_J]
+    _Rmsg.resize( _fg->nrFactors() );
     for( size_t I = 0; I < _fg->nrFactors(); I++ ) {
-        _R[I].resize( _fg->nbF(I).size() );
+        _Rmsg[I].resize( _fg->nbF(I).size() );
         foreach( const Neighbor &i, _fg->nbF(I) ) {
-            _R[I][i.iter].resize( _fg->nbV(i).size() );
+            _Rmsg[I][i.iter].resize( _fg->nbV(i).size() );
             foreach( const Neighbor &J, _fg->nbV(i) ) {
                 if( I != J ) {
                     Prob prod( _fg->var(i).states(), 1.0 );
                     foreach( const Neighbor &K, _fg->nbV(i) )
                         if( K.node != I && K.node != J.node )
                             prod *= _bp_dual.msgM( i, K.iter );
-                    _R[I][i.iter][J.iter] = prod;
+                    _Rmsg[I][i.iter][J.iter] = prod;
                 }
             }
         }
@@ -217,7 +269,7 @@ void BBP::RegeneratePsiAdjoints() {
     _adj_psi_V.reserve( _fg->nrVars() );
     for( size_t i = 0; i < _fg->nrVars(); i++ ) {
         Prob p( _adj_b_V_unnorm[i] );
-        assert( p.size() == _fg->var(i).states() );
+        DAI_ASSERT( p.size() == _fg->var(i).states() );
         foreach( const Neighbor &I, _fg->nbV(i) )
             p *= _bp_dual.msgM( i, I.iter );
         p += _init_adj_psi_V[i];
@@ -227,13 +279,13 @@ 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] );
-        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 );
             // multiply prod with n_jI
             for( size_t x_I = 0; x_I < p.size(); x_I++ )
-                p[x_I] *= n_iI[ind[x_I]];
+                p.set( x_I, p[x_I] * n_iI[ind[x_I]] );
         }
         p += _init_adj_psi_F[I];
         _adj_psi_F.push_back( p );
@@ -267,19 +319,19 @@ void BBP::RegenerateParMessageAdjoints() {
                         const _ind_t &ind = _index( j, j.dual );
                         // multiply prod with n_jI
                         for( size_t x_I = 0; x_I < prod.size(); x_I++ )
-                            prod[x_I] *= n_jI[ind[x_I]];
+                            prod.set( x_I, prod[x_I] * n_jI[ind[x_I]] );
                     }
                 Prob marg( _fg->var(i).states(), 0.0 );
                 const _ind_t &ind = _index( i, I.iter );
                 for( size_t r = 0; r < prod.size(); r++ )
-                    marg[ind[r]] += prod[r];
+                    marg.set( ind[r], marg[ind[r]] + prod[r] );
                 _new_adj_n[i][I.iter] = marg;
                 upMsgN( i, I.iter );
             }
 
             { // calculate adj_m
                 Prob prod( _adj_b_V_unnorm[i] );
-                assert( prod.size() == _fg->var(i).states() );
+                DAI_ASSERT( prod.size() == _fg->var(i).states() );
                 foreach( const Neighbor &J, _fg->nbV(i) )
                     if( J.node != I.node )
                         prod *= _bp_dual.msgM(i,J.iter);
@@ -304,7 +356,7 @@ void BBP::RegenerateSeqMessageAdjoints() {
         foreach( const Neighbor &I, _fg->nbV(i) ) {
             // calculate adj_m
             Prob prod( _adj_b_V_unnorm[i] );
-            assert( prod.size() == _fg->var(i).states() );
+            DAI_ASSERT( prod.size() == _fg->var(i).states() );
             foreach( const Neighbor &J, _fg->nbV(i) )
                 if( J.node != I.node )
                     prod *= _bp_dual.msgM( i, J.iter );
@@ -324,18 +376,34 @@ void BBP::RegenerateSeqMessageAdjoints() {
                     const _ind_t& ind = _index( j, j.dual );
                     // multiply prod with n_jI
                     for( size_t x_I = 0; x_I < prod.size(); x_I++ )
-                        prod[x_I] *= n_jI[ind[x_I]];
+                        prod.set( x_I, prod[x_I] * n_jI[ind[x_I]] );
                 }
             Prob marg( _fg->var(i).states(), 0.0 );
             const _ind_t &ind = _index( i, I.iter );
             for( size_t r = 0; r < prod.size(); r++ )
-                marg[ind[r]] += prod[r];
+                marg.set( ind[r], marg[ind[r]] + prod[r] );
             sendSeqMsgN( i, I.iter,marg );
         }
     }
 }
 
 
+void BBP::Regenerate() {
+    RegenerateInds();
+    RegenerateT();
+    RegenerateU();
+    RegenerateS();
+    RegenerateR();
+    RegenerateInputs();
+    RegeneratePsiAdjoints();
+    if( props.updates == Properties::UpdateType::PAR )
+        RegenerateParMessageAdjoints();
+    else
+        RegenerateSeqMessageAdjoints();
+    _iters = 0;
+}
+
+
 void BBP::calcNewN( size_t i, size_t _I ) {
     _adj_psi_V[i] += T(i,_I) * _adj_n_unnorm[i][_I];
     Prob &new_adj_n_iI = _new_adj_n[i][_I];
@@ -343,10 +411,10 @@ void BBP::calcNewN( size_t i, size_t _I ) {
     size_t I = _fg->nbV(i)[_I];
     foreach( const Neighbor &j, _fg->nbF(I) )
         if( j != i ) {
-            const Prob &p = _S[i][_I][j.iter];
+            const Prob &p = _Smsg[i][_I][j.iter];
             const Prob &_adj_m_unnorm_jI = _adj_m_unnorm[j][j.dual];
             LOOP_ij(
-                new_adj_n_iI[xi] += p[xij] * _adj_m_unnorm_jI[xj];
+                new_adj_n_iI.set( xi, new_adj_n_iI[xi] + p[xij] * _adj_m_unnorm_jI[xj] );
             );
             /* THE FOLLOWING WOULD BE ABOUT TWICE AS SLOW:
             Var vi = _fg->var(i);
@@ -363,7 +431,7 @@ void BBP::calcNewM( size_t i, size_t _I ) {
     const Prob &adj = _adj_m_unnorm[i][_I];
     const _ind_t &ind = _index(i,_I);
     for( size_t x_I = 0; x_I < p.size(); x_I++ )
-        p[x_I] *= adj[ind[x_I]];
+        p.set( x_I, p[x_I] * adj[ind[x_I]] );
     _adj_psi_F[I] += p;
     /* THE FOLLOWING WOULD BE SLIGHTLY SLOWER:
     _adj_psi_F[I] += (Factor( _fg->factor(I).vars(), U(I, I.dual) ) * Factor( _fg->var(i), _adj_m_unnorm[i][_I] )).p();
@@ -372,7 +440,7 @@ void BBP::calcNewM( size_t i, size_t _I ) {
     _new_adj_m[i][_I] = Prob( _fg->var(i).states(), 0.0 );
     foreach( const Neighbor &J, _fg->nbV(i) )
         if( J != I )
-            _new_adj_m[i][_I] += _R[I][I.dual][J.iter] * _adj_n_unnorm[i][J.iter];
+            _new_adj_m[i][_I] += _Rmsg[I][I.dual][J.iter] * _adj_n_unnorm[i][J.iter];
 }
 
 
@@ -457,7 +525,7 @@ void BBP::setSeqMsgM( size_t i, size_t _I, const Prob &p ) {
 void BBP::sendSeqMsgN( size_t i, size_t _I, const Prob &f ) {
     Prob f_unnorm = unnormAdjoint( _bp_dual.msgN(i,_I), _bp_dual.zN(i,_I), f );
     const Neighbor &I = _fg->nbV(i)[_I];
-    assert( I.iter == _I );
+    DAI_ASSERT( I.iter == _I );
     _adj_psi_V[i] += f_unnorm * T( i, _I );
 #if 0
     if(f_unnorm.sumAbs() > pv_thresh) {
@@ -479,8 +547,8 @@ void BBP::sendSeqMsgN( size_t i, size_t _I, const Prob &f ) {
                 DAI_DMSG("in sendSeqMsgN loop");
                 DAI_PV(J);
                 DAI_PV(f_unnorm);
-                DAI_PV(_R[J][J.dual][_I]);
-                DAI_PV(f_unnorm * _R[J][J.dual][_I]);
+                DAI_PV(_Rmsg[J][J.dual][_I]);
+                DAI_PV(f_unnorm * _Rmsg[J][J.dual][_I]);
             }
 #endif
             incrSeqMsgM( i, J.iter, f_unnorm * R(J, J.dual, _I) );
@@ -503,7 +571,7 @@ void BBP::sendSeqMsgM( size_t j, size_t _I ) {
     Prob um( U(I, _j) );
     const _ind_t &ind = _index(j, _I);
     for( size_t x_I = 0; x_I < um.size(); x_I++ )
-        um[x_I] *= _adj_m_unnorm_jI[ind[x_I]];
+        um.set( x_I, um[x_I] * _adj_m_unnorm_jI[ind[x_I]] );
     um *= 1 - props.damping;
     _adj_psi_F[I] += um;
 
@@ -518,10 +586,10 @@ void BBP::sendSeqMsgM( size_t j, size_t _I ) {
 //     DAI_PV(_fg->nbF(I).size());
     foreach( const Neighbor &i, _fg->nbF(I) ) {
         if( i.node != j ) {
-            const Prob &S = _S[i][i.dual][_j];
+            const Prob &S = _Smsg[i][i.dual][_j];
             Prob msg( _fg->var(i).states(), 0.0 );
             LOOP_ij(
-                msg[xi] += S[xij] * _adj_m_unnorm_jI[xj];
+                msg.set( xi, msg[xi] + S[xij] * _adj_m_unnorm_jI[xj] );
             );
             msg *= 1.0 - props.damping;
             /* THE FOLLOWING WOULD BE ABOUT TWICE AS SLOW:
@@ -538,7 +606,7 @@ void BBP::sendSeqMsgM( size_t j, size_t _I ) {
                 DAI_PV(_I);
                 DAI_PV(_fg->nbF(I).size());
                 DAI_PV(_fg->factor(I).p());
-                DAI_PV(_S[i][i.dual][_j]);
+                DAI_PV(_Smsg[i][i.dual][_j]);
 
                 DAI_PV(i);
                 DAI_PV(i.dual);
@@ -546,7 +614,7 @@ void BBP::sendSeqMsgM( size_t j, size_t _I ) {
                 DAI_PV(_fg->nbV(i).size());
             }
 #endif
-            assert( _fg->nbV(i)[i.dual].node == I );
+            DAI_ASSERT( _fg->nbV(i)[i.dual].node == I );
             sendSeqMsgN( i, i.dual, msg );
         }
     }
@@ -554,6 +622,19 @@ void BBP::sendSeqMsgM( size_t j, size_t _I ) {
 }
 
 
+Prob BBP::unnormAdjoint( const Prob &w, Real Z_w, const Prob &adj_w ) {
+    DAI_ASSERT( w.size() == adj_w.size() );
+    Prob adj_w_unnorm( w.size(), 0.0 );
+    Real s = 0.0;
+    for( size_t i = 0; i < w.size(); i++ )
+        s += w[i] * adj_w[i];
+    for( size_t i = 0; i < w.size(); i++ )
+        adj_w_unnorm.set( i, (adj_w[i] - s) / Z_w );
+    return adj_w_unnorm;
+//  THIS WOULD BE ABOUT 50% SLOWER:  return (adj_w - (w * adj_w).sum()) / Z_w;
+}
+
+
 Real BBP::getUnMsgMag() {
     Real s = 0.0;
     size_t e = 0;
@@ -595,7 +676,7 @@ void BBP::getMsgMags( Real &s, Real &new_s ) {
 //             argmax_var_state = argmax_state.first;
 //         }
 //     }
-//     assert(/*0 <= argmax_var_state &&*/
+//     DAI_ASSERT(/*0 <= argmax_var_state &&*/
 //            argmax_var_state < _fg->var(argmax_var).states());
 //     return tuple<size_t,size_t,Real>(argmax_var,argmax_var_state,max_var);
 // }
@@ -613,7 +694,7 @@ void BBP::getArgmaxMsgM( size_t &out_i, size_t &out__I, Real &mag ) {
                 out__I = I.iter;
             }
         }
-    assert( found );
+    DAI_ASSERT( found );
 }
 
 
@@ -652,27 +733,11 @@ Real BBP::getTotalMsgN() {
 }
 
 
-void BBP::Regenerate() {
-    RegenerateInds();
-    RegenerateT();
-    RegenerateU();
-    RegenerateS();
-    RegenerateR();
-    RegenerateInputs();
-    RegeneratePsiAdjoints();
-    if( props.updates == Properties::UpdateType::PAR )
-        RegenerateParMessageAdjoints();
-    else
-        RegenerateSeqMessageAdjoints();
-    _iters = 0;
-}
-
-
 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;
 }
 
@@ -686,9 +751,160 @@ std::vector<Prob> BBP::getZeroAdjV( const FactorGraph &fg ) {
 }
 
 
+void BBP::initCostFnAdj( const BBPCostFunction &cfn, const vector<size_t> *stateP ) {
+    const FactorGraph &fg = _ia->fg();
+
+    switch( (size_t)cfn ) {
+        case BBPCostFunction::CFN_BETHE_ENT: {
+            vector<Prob> b1_adj;
+            vector<Prob> b2_adj;
+            vector<Prob> psi1_adj;
+            vector<Prob> psi2_adj;
+            b1_adj.reserve( fg.nrVars() );
+            psi1_adj.reserve( fg.nrVars() );
+            b2_adj.reserve( fg.nrFactors() );
+            psi2_adj.reserve( fg.nrFactors() );
+            for( size_t i = 0; i < fg.nrVars(); i++ ) {
+                size_t dim = fg.var(i).states();
+                int c = fg.nbV(i).size();
+                Prob p(dim,0.0);
+                for( size_t xi = 0; xi < dim; xi++ )
+                    p.set( xi, (1 - c) * (1 + log( _ia->beliefV(i)[xi] )) );
+                b1_adj.push_back( p );
+
+                for( size_t xi = 0; xi < dim; xi++ )
+                    p.set( xi, -_ia->beliefV(i)[xi] );
+                psi1_adj.push_back( p );
+            }
+            for( size_t I = 0; I < fg.nrFactors(); I++ ) {
+                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] ) );
+                b2_adj.push_back( p );
+
+                for( size_t xI = 0; xI < dim; xI++ )
+                    p.set( xI, -_ia->beliefF(I)[xI] / fg.factor(I).p()[xI] );
+                psi2_adj.push_back( p );
+            }
+            init( b1_adj, b2_adj, psi1_adj, psi2_adj );
+            break;
+        } case BBPCostFunction::CFN_FACTOR_ENT: {
+            vector<Prob> b2_adj;
+            b2_adj.reserve( fg.nrFactors() );
+            for( size_t I = 0; I < fg.nrFactors(); I++ ) {
+                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];
+                    if( bIxI < 1.0e-15 )
+                        p.set( xI, -1.0e10 );
+                    else
+                        p.set( xI, 1 + log( bIxI ) );
+                }
+                b2_adj.push_back(p);
+            }
+            init_F( b2_adj );
+            break;
+        } case BBPCostFunction::CFN_VAR_ENT: {
+            vector<Prob> b1_adj;
+            b1_adj.reserve( fg.nrVars() );
+            for( size_t i = 0; i < fg.nrVars(); i++ ) {
+                size_t dim = fg.var(i).states();
+                Prob p( dim, 0.0 );
+                for( size_t xi = 0; xi < fg.var(i).states(); xi++ ) {
+                    Real bixi = _ia->beliefV(i)[xi];
+                    if( bixi < 1.0e-15 )
+                        p.set( xi, -1.0e10 );
+                    else
+                        p.set( xi, 1 + log( bixi ) );
+                }
+                b1_adj.push_back( p );
+            }
+            init_V( b1_adj );
+            break;
+        } case BBPCostFunction::CFN_GIBBS_B:
+          case BBPCostFunction::CFN_GIBBS_B2:
+          case BBPCostFunction::CFN_GIBBS_EXP: {
+            // cost functions that use Gibbs sample, summing over variable marginals
+            vector<size_t> state;
+            if( stateP == NULL )
+                state = getGibbsState( _ia->fg(), 2*_ia->Iterations() );
+            else
+                state = *stateP;
+            DAI_ASSERT( state.size() == fg.nrVars() );
+
+            vector<Prob> b1_adj;
+            b1_adj.reserve(fg.nrVars());
+            for( size_t i = 0; i < state.size(); i++ ) {
+                size_t n = fg.var(i).states();
+                Prob delta( n, 0.0 );
+                DAI_ASSERT(/*0<=state[i] &&*/ state[i] < n);
+                Real b = _ia->beliefV(i)[state[i]];
+                switch( (size_t)cfn ) {
+                    case BBPCostFunction::CFN_GIBBS_B:
+                        delta.set( state[i], 1.0 );
+                        break;
+                    case BBPCostFunction::CFN_GIBBS_B2:
+                        delta.set( state[i], b );
+                        break;
+                    case BBPCostFunction::CFN_GIBBS_EXP:
+                        delta.set( state[i], exp(b) );
+                        break;
+                    default:
+                        DAI_THROW(UNKNOWN_ENUM_VALUE);
+                }
+                b1_adj.push_back( delta );
+            }
+            init_V( b1_adj );
+            break;
+        } case BBPCostFunction::CFN_GIBBS_B_FACTOR:
+          case BBPCostFunction::CFN_GIBBS_B2_FACTOR:
+          case BBPCostFunction::CFN_GIBBS_EXP_FACTOR: {
+            // cost functions that use Gibbs sample, summing over factor marginals
+            vector<size_t> state;
+            if( stateP == NULL )
+                state = getGibbsState( _ia->fg(), 2*_ia->Iterations() );
+            else
+                state = *stateP;
+            DAI_ASSERT( state.size() == fg.nrVars() );
+
+            vector<Prob> b2_adj;
+            b2_adj.reserve( fg.nrVars() );
+            for( size_t I = 0; I <  fg.nrFactors(); I++ ) {
+                size_t n = fg.factor(I).nrStates();
+                Prob delta( n, 0.0 );
+
+                size_t x_I = getFactorEntryForState( fg, I, state );
+                DAI_ASSERT(/*0<=x_I &&*/ x_I < n);
+
+                Real b = _ia->beliefF(I)[x_I];
+                switch( (size_t)cfn ) {
+                    case BBPCostFunction::CFN_GIBBS_B_FACTOR:
+                        delta.set( x_I, 1.0 );
+                        break;
+                    case BBPCostFunction::CFN_GIBBS_B2_FACTOR:
+                        delta.set( x_I, b );
+                        break;
+                    case BBPCostFunction::CFN_GIBBS_EXP_FACTOR:
+                        delta.set( x_I, exp( b ) );
+                        break;
+                    default:
+                        DAI_THROW(UNKNOWN_ENUM_VALUE);
+                }
+                b2_adj.push_back( delta );
+            }
+            init_F( b2_adj );
+            break;
+        } default:
+            DAI_THROW(UNKNOWN_ENUM_VALUE);
+    }
+}
+
+
 void BBP::run() {
     typedef BBP::Properties::UpdateType UT;
-    Real &tol = props.tol;
+    Real tol = props.tol;
     UT &updates = props.updates;
 
     Real tic = toc();
@@ -766,17 +982,16 @@ void BBP::run() {
         }
     }
     if( props.verbose >= 3 )
-        cerr << "BBP::run() took " << toc()-tic << " seconds " << doneIters() << " iterations" << endl;
+        cerr << "BBP::run() took " << toc()-tic << " seconds " << Iterations() << " iterations" << endl;
 }
 
 
-double numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const PropertySet &bbp_props, bbp_cfn_t cfn, double h ) {
+Real numericBBPTest( const InfAlg &bp, const std::vector<size_t> *state, const PropertySet &bbp_props, const BBPCostFunction &cfn, Real h ) {
+    BBP bbp( &bp, bbp_props );
     // calculate the value of the unperturbed cost function
-    Real cf0 = getCostFn( bp, cfn, state );
-
+    Real cf0 = cfn.evaluate( bp, state );
     // run BBP to estimate adjoints
-    BBP bbp( &bp, bbp_props );
-    initBBPCostFnAdj( bbp, bp, cfn, state );
+    bbp.initCostFnAdj( cfn, state );
     bbp.run();
 
     Real d = 0;
@@ -787,7 +1002,7 @@ double numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const Prop
 
         // for each variable i
         for( size_t i = 0; i < fg.nrVars(); i++ ) {
-            vector<double> adj_est;
+            vector<Real> adj_est;
             // for each value xi
             for( size_t xi = 0; xi < fg.var(i).states(); xi++ ) {
                 // Clone 'bp' (which may be any InfAlg)
@@ -796,10 +1011,11 @@ double numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const Prop
                 // perturb it
                 size_t n = bp_prb->fg().var(i).states();
                 Prob psi_1_prb( n, 1.0 );
-                psi_1_prb[xi] += h;
+                psi_1_prb.set( xi, psi_1_prb[xi] + h );
 //                 psi_1_prb.normalize();
                 size_t I = bp_prb->fg().nbV(i)[0]; // use first factor in list of neighbors of i
-                bp_prb->fg().factor(I) *= Factor( bp_prb->fg().var(i), psi_1_prb );
+                Factor tmp = bp_prb->fg().factor(I) * Factor( bp_prb->fg().var(i), psi_1_prb );
+                bp_prb->fg().setFactor( I, tmp );
 
                 // call 'init' on the perturbed variables
                 bp_prb->init( bp_prb->fg().var(i) );
@@ -808,7 +1024,7 @@ double numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const Prop
                 bp_prb->run();
 
                 // calculate new value of cost function
-                Real cf_prb = getCostFn( *bp_prb, cfn, state );
+                Real cf_prb = cfn.evaluate( *bp_prb, state );
 
                 // use to estimate adjoint for i
                 adj_est.push_back( (cf_prb - cf0) / h );
@@ -816,12 +1032,12 @@ double numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const Prop
                 // free cloned InfAlg
                 delete bp_prb;
             }
-            Prob p_adj_est( adj_est.begin(), adj_est.end() );
+            Prob p_adj_est( adj_est );
             // compare this numerical estimate to the BBP estimate; sum the distances
             cout << "i: " << i
                  << ", p_adj_est: " << p_adj_est
                  << ", bbp.adj_psi_V(i): " << bbp.adj_psi_V(i) << endl;
-            d += dist( p_adj_est, bbp.adj_psi_V(i), Prob::DISTL1 );
+            d += dist( p_adj_est, bbp.adj_psi_V(i), DISTL1 );
         }
     }
     /*    if(1) {
@@ -838,7 +1054,7 @@ double numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const Prop
         for(size_t i=0; i<bp_dual.nrVars(); i++) {
             // for each factor I ~ i
             foreach(size_t I, bp_dual.nbV(i)) {
-                vector<double> adj_n_est;
+                vector<Real> adj_n_est;
                 // for each value xi
                 for(size_t xi=0; xi<bp_dual.var(i).states(); xi++) {
                     BP_dual bp_dual_prb(bp_dual);
@@ -852,7 +1068,7 @@ double numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const Prop
                     adj_n_est.push_back((cf_prb-cf0)/h);
                 }
 
-                vector<double> adj_m_est;
+                vector<Real> adj_m_est;
                 // for each value xi
                 for(size_t xi=0; xi<bp_dual.var(i).states(); xi++) {
                     BP_dual bp_dual_prb(bp_dual);
@@ -866,19 +1082,19 @@ double numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const Prop
                     adj_m_est.push_back((cf_prb-cf0)/h);
                 }
 
-                Prob p_adj_n_est(adj_n_est.begin(), adj_n_est.end());
+                Prob p_adj_n_est( adj_n_est );
                 // compare this numerical estimate to the BBP estimate; sum the distances
                 cerr << "i: " << i << ", I: " << I
                      << ", adj_n_est: " << p_adj_n_est
                      << ", bbp.adj_n(i,I): " << bbp.adj_n(i,I) << endl;
-                d += dist(p_adj_n_est, bbp.adj_n(i,I), Prob::DISTL1);
+                d += dist(p_adj_n_est, bbp.adj_n(i,I), DISTL1);
 
-                Prob p_adj_m_est(adj_m_est.begin(), adj_m_est.end());
+                Prob p_adj_m_est( adj_m_est );
                 // compare this numerical estimate to the BBP estimate; sum the distances
                 cerr << "i: " << i << ", I: " << I
                      << ", adj_m_est: " << p_adj_m_est
                      << ", bbp.adj_m(I,i): " << bbp.adj_m(I,i) << endl;
-                d += dist(p_adj_m_est, bbp.adj_m(I,i), Prob::DISTL1);
+                d += dist(p_adj_m_est, bbp.adj_m(I,i), DISTL1);
             }
         }
     }
@@ -886,7 +1102,7 @@ double numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const Prop
     /*    if(0) {
         // verify bbp.adj_b_V
         for(size_t i=0; i<bp_dual.nrVars(); i++) {
-            vector<double> adj_b_V_est;
+            vector<Real> adj_b_V_est;
             // for each value xi
             for(size_t xi=0; xi<bp_dual.var(i).states(); xi++) {
                 BP_dual bp_dual_prb(bp_dual);
@@ -900,12 +1116,12 @@ double numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const Prop
                 // add it to list of adjoints
                 adj_b_V_est.push_back((cf_prb-cf0)/h);
             }
-            Prob p_adj_b_V_est(adj_b_V_est.begin(), adj_b_V_est.end());
+            Prob p_adj_b_V_est( adj_b_V_est );
             // compare this numerical estimate to the BBP estimate; sum the distances
             cerr << "i: " << i
                  << ", adj_b_V_est: " << p_adj_b_V_est
                  << ", bbp.adj_b_V(i): " << bbp.adj_b_V(i) << endl;
-            d += dist(p_adj_b_V_est, bbp.adj_b_V(i), Prob::DISTL1);
+            d += dist(p_adj_b_V_est, bbp.adj_b_V(i), DISTL1);
         }
     }
     */
@@ -915,242 +1131,6 @@ double numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const Prop
 }
 
 
-bool needGibbsState( bbp_cfn_t cfn ) {
-    switch( (size_t)cfn ) {
-        case bbp_cfn_t::CFN_GIBBS_B:
-        case bbp_cfn_t::CFN_GIBBS_B2:
-        case bbp_cfn_t::CFN_GIBBS_EXP:
-        case bbp_cfn_t::CFN_GIBBS_B_FACTOR:
-        case bbp_cfn_t::CFN_GIBBS_B2_FACTOR:
-        case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR:
-            return true;
-        default:
-            return false;
-    }
-}
-
-
-void initBBPCostFnAdj( BBP &bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const vector<size_t> *stateP ) {
-    const FactorGraph &fg = ia.fg();
-
-    switch( (size_t)cfn_type ) {
-        case bbp_cfn_t::CFN_BETHE_ENT: {
-            vector<Prob> b1_adj;
-            vector<Prob> b2_adj;
-            vector<Prob> psi1_adj;
-            vector<Prob> psi2_adj;
-            b1_adj.reserve( fg.nrVars() );
-            psi1_adj.reserve( fg.nrVars() );
-            b2_adj.reserve( fg.nrFactors() );
-            psi2_adj.reserve( fg.nrFactors() );
-            for( size_t i = 0; i < fg.nrVars(); i++ ) {
-                size_t dim = fg.var(i).states();
-                int c = fg.nbV(i).size();
-                Prob p(dim,0.0);
-                for( size_t xi = 0; xi < dim; xi++ )
-                    p[xi] = (1 - c) * (1 + log( ia.beliefV(i)[xi] ));
-                b1_adj.push_back( p );
-
-                for( size_t xi = 0; xi < dim; xi++ )
-                    p[xi] = -ia.beliefV(i)[xi];
-                psi1_adj.push_back( p );
-            }
-            for( size_t I = 0; I < fg.nrFactors(); I++ ) {
-                size_t dim = fg.factor(I).states();
-                Prob p( dim, 0.0 );
-                for( size_t xI = 0; xI < dim; xI++ )
-                    p[xI] = 1 + log( ia.beliefF(I)[xI] / fg.factor(I).p()[xI] );
-                b2_adj.push_back( p );
-
-                for( size_t xI = 0; xI < dim; xI++ )
-                    p[xI] = -ia.beliefF(I)[xI] / fg.factor(I).p()[xI];
-                psi2_adj.push_back( p );
-            }
-            bbp.init( b1_adj, b2_adj, psi1_adj, psi2_adj );
-            break;
-        } case bbp_cfn_t::CFN_FACTOR_ENT: {
-            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();
-                Prob p( dim, 0.0 );
-                for( size_t xI = 0; xI < dim; xI++ ) {
-                    double bIxI = ia.beliefF(I)[xI];
-                    if( bIxI < 1.0e-15 )
-                        p[xI] = -1.0e10;
-                    else
-                        p[xI] = 1 + log( bIxI );
-                }
-                b2_adj.push_back(p);
-            }
-            bbp.init( bbp.getZeroAdjV(fg), b2_adj );
-            break;
-        } case bbp_cfn_t::CFN_VAR_ENT: {
-            vector<Prob> b1_adj;
-            b1_adj.reserve( fg.nrVars() );
-            for( size_t i = 0; i < fg.nrVars(); i++ ) {
-                size_t dim = fg.var(i).states();
-                Prob p( dim, 0.0 );
-                for( size_t xi = 0; xi < fg.var(i).states(); xi++ ) {
-                    double bixi = ia.beliefV(i)[xi];
-                    if( bixi < 1.0e-15 )
-                        p[xi] = -1.0e10;
-                    else
-                        p[xi] = 1 + log( bixi );
-                }
-                b1_adj.push_back( p );
-            }
-            bbp.init( b1_adj );
-            break;
-        } case bbp_cfn_t::CFN_GIBBS_B:
-          case bbp_cfn_t::CFN_GIBBS_B2:
-          case bbp_cfn_t::CFN_GIBBS_EXP: {
-            // cost functions that use Gibbs sample, summing over variable marginals
-            vector<size_t> state;
-            if( stateP == NULL )
-                state = getGibbsState( ia, 2*ia.Iterations() );
-            else
-                state = *stateP;
-            assert( state.size() == fg.nrVars() );
-
-            vector<Prob> b1_adj;
-            b1_adj.reserve(fg.nrVars());
-            for( size_t i = 0; i < state.size(); i++ ) {
-                size_t n = fg.var(i).states();
-                Prob delta( n, 0.0 );
-                assert(/*0<=state[i] &&*/ state[i] < n);
-                double b = ia.beliefV(i)[state[i]];
-                switch( (size_t)cfn_type ) {
-                    case bbp_cfn_t::CFN_GIBBS_B:
-                        delta[state[i]] = 1.0;
-                        break;
-                    case bbp_cfn_t::CFN_GIBBS_B2:
-                        delta[state[i]] = b;
-                        break;
-                    case bbp_cfn_t::CFN_GIBBS_EXP:
-                        delta[state[i]] = exp(b);
-                        break;
-                    default:
-                        DAI_THROW(UNKNOWN_ENUM_VALUE);
-                }
-                b1_adj.push_back( delta );
-            }
-            bbp.init( b1_adj );
-            break;
-        } case bbp_cfn_t::CFN_GIBBS_B_FACTOR:
-          case bbp_cfn_t::CFN_GIBBS_B2_FACTOR:
-          case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR: {
-            // cost functions that use Gibbs sample, summing over factor marginals
-            vector<size_t> state;
-            if( stateP == NULL )
-                state = getGibbsState( ia, 2*ia.Iterations() );
-            else
-                state = *stateP;
-            assert( state.size() == fg.nrVars() );
-
-            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();
-                Prob delta( n, 0.0 );
-
-                size_t x_I = getFactorEntryForState( fg, I, state );
-                assert(/*0<=x_I &&*/ x_I < n);
-
-                double b = ia.beliefF(I)[x_I];
-                switch( (size_t)cfn_type ) {
-                    case bbp_cfn_t::CFN_GIBBS_B_FACTOR:
-                        delta[x_I] = 1.0;
-                        break;
-                    case bbp_cfn_t::CFN_GIBBS_B2_FACTOR:
-                        delta[x_I] = b;
-                        break;
-                    case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR:
-                        delta[x_I] = exp( b );
-                        break;
-                    default:
-                        DAI_THROW(UNKNOWN_ENUM_VALUE);
-                }
-                b2_adj.push_back( delta );
-            }
-            bbp.init( bbp.getZeroAdjV(fg), b2_adj );
-            break;
-        } default:
-            DAI_THROW(UNKNOWN_ENUM_VALUE);
-    }
-}
-
-
-Real getCostFn( const InfAlg &ia, bbp_cfn_t cfn_type, const vector<size_t> *stateP ) {
-    double cf = 0.0;
-    const FactorGraph &fg = ia.fg();
-
-    switch( (size_t)cfn_type ) {
-        case bbp_cfn_t::CFN_BETHE_ENT: // ignores state
-            cf = -ia.logZ();
-            break;
-        case bbp_cfn_t::CFN_VAR_ENT: // ignores state
-            for( size_t i = 0; i < fg.nrVars(); i++ )
-                cf += -ia.beliefV(i).entropy();
-            break;
-        case bbp_cfn_t::CFN_FACTOR_ENT: // ignores state
-            for( size_t I = 0; I < fg.nrFactors(); I++ )
-                cf += -ia.beliefF(I).entropy();
-            break;
-        case bbp_cfn_t::CFN_GIBBS_B:
-        case bbp_cfn_t::CFN_GIBBS_B2:
-        case bbp_cfn_t::CFN_GIBBS_EXP: {
-            assert( stateP != NULL );
-            vector<size_t> state = *stateP;
-            assert( state.size() == fg.nrVars() );
-            for( size_t i = 0; i < fg.nrVars(); i++ ) {
-                double b = ia.beliefV(i)[state[i]];
-                switch( (size_t)cfn_type ) {
-                    case bbp_cfn_t::CFN_GIBBS_B:
-                        cf += b;
-                        break;
-                    case bbp_cfn_t::CFN_GIBBS_B2:
-                        cf += b * b / 2.0;
-                        break;
-                    case bbp_cfn_t::CFN_GIBBS_EXP:
-                        cf += exp( b );
-                        break;
-                    default:
-                        DAI_THROW(UNKNOWN_ENUM_VALUE);
-                }
-            }
-            break;
-        } case bbp_cfn_t::CFN_GIBBS_B_FACTOR:
-          case bbp_cfn_t::CFN_GIBBS_B2_FACTOR:
-          case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR: {
-            assert( stateP != NULL );
-            vector<size_t> state = *stateP;
-            assert( state.size() == fg.nrVars() );
-            for( size_t I = 0; I < fg.nrFactors(); I++ ) {
-                size_t x_I = getFactorEntryForState( fg, I, state );
-                double b = ia.beliefF(I)[x_I];
-                switch( (size_t)cfn_type ) {
-                    case bbp_cfn_t::CFN_GIBBS_B_FACTOR:
-                        cf += b;
-                        break;
-                    case bbp_cfn_t::CFN_GIBBS_B2_FACTOR:
-                        cf += b * b / 2.0;
-                        break;
-                    case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR:
-                        cf += exp( b );
-                        break;
-                    default:
-                        DAI_THROW(UNKNOWN_ENUM_VALUE);
-                }
-            }
-            break;
-        } default:
-            DAI_THROWE(UNKNOWN_ENUM_VALUE, "Unknown cost function " + std::string(cfn_type));
-    }
-    return cf;
-}
-
-
 } // end of namespace dai
 
 
@@ -1161,39 +1141,45 @@ namespace dai {
 
 void BBP::Properties::set(const PropertySet &opts)
 {
-    const std::set<PropertyKey> &keys = opts.allKeys();
-    std::set<PropertyKey>::const_iterator i;
-    for(i=keys.begin(); i!=keys.end(); i++) {
-        if(*i == "verbose") continue;
-        if(*i == "maxiter") continue;
-        if(*i == "tol") continue;
-        if(*i == "damping") continue;
-        if(*i == "updates") continue;
-        DAI_THROWE(UNKNOWN_PROPERTY_TYPE, "BBP: Unknown property " + *i);
+    const std::set<PropertyKey> &keys = opts.keys();
+    std::string errormsg;
+    for( std::set<PropertyKey>::const_iterator i = keys.begin(); i != keys.end(); i++ ) {
+        if( *i == "verbose" ) continue;
+        if( *i == "maxiter" ) continue;
+        if( *i == "tol" ) continue;
+        if( *i == "damping" ) continue;
+        if( *i == "updates" ) continue;
+        errormsg = errormsg + "BBP: Unknown property " + *i + "\n";
+    }
+    if( !errormsg.empty() )
+        DAI_THROWE(UNKNOWN_PROPERTY, errormsg);
+    if( !opts.hasKey("maxiter") )
+        errormsg = errormsg + "BBP: Missing property \"maxiter\" for method \"BBP\"\n";
+    if( !opts.hasKey("tol") )
+        errormsg = errormsg + "BBP: Missing property \"tol\" for method \"BBP\"\n";
+    if( !opts.hasKey("damping") )
+        errormsg = errormsg + "BBP: Missing property \"damping\" for method \"BBP\"\n";
+    if( !opts.hasKey("updates") )
+        errormsg = errormsg + "BBP: Missing property \"updates\" for method \"BBP\"\n";
+    if( !errormsg.empty() )
+        DAI_THROWE(NOT_ALL_PROPERTIES_SPECIFIED,errormsg);
+    if( opts.hasKey("verbose") ) {
+        verbose = opts.getStringAs<size_t>("verbose");
+    } else {
+        verbose = 0;
     }
-    if(!opts.hasKey("verbose"))
-        DAI_THROWE(NOT_ALL_PROPERTIES_SPECIFIED,"BBP: Missing property \"verbose\" for method \"BBP\"");
-    if(!opts.hasKey("maxiter"))
-        DAI_THROWE(NOT_ALL_PROPERTIES_SPECIFIED,"BBP: Missing property \"maxiter\" for method \"BBP\"");
-    if(!opts.hasKey("tol"))
-        DAI_THROWE(NOT_ALL_PROPERTIES_SPECIFIED,"BBP: Missing property \"tol\" for method \"BBP\"");
-    if(!opts.hasKey("damping"))
-        DAI_THROWE(NOT_ALL_PROPERTIES_SPECIFIED,"BBP: Missing property \"damping\" for method \"BBP\"");
-    if(!opts.hasKey("updates"))
-        DAI_THROWE(NOT_ALL_PROPERTIES_SPECIFIED,"BBP: Missing property \"updates\" for method \"BBP\"");
-    verbose = opts.getStringAs<size_t>("verbose");
     maxiter = opts.getStringAs<size_t>("maxiter");
-    tol = opts.getStringAs<double>("tol");
-    damping = opts.getStringAs<double>("damping");
+    tol = opts.getStringAs<Real>("tol");
+    damping = opts.getStringAs<Real>("damping");
     updates = opts.getStringAs<UpdateType>("updates");
 }
 PropertySet BBP::Properties::get() const {
     PropertySet opts;
-    opts.Set("verbose", verbose);
-    opts.Set("maxiter", maxiter);
-    opts.Set("tol", tol);
-    opts.Set("damping", damping);
-    opts.Set("updates", updates);
+    opts.set("verbose", verbose);
+    opts.set("maxiter", maxiter);
+    opts.set("tol", tol);
+    opts.set("damping", damping);
+    opts.set("updates", updates);
     return opts;
 }
 string BBP::Properties::toString() const {