Merge branch 'master' of git@git.tuebingen.mpg.de:libdai
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Thu, 1 Apr 2010 12:38:52 +0000 (14:38 +0200)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Thu, 1 Apr 2010 12:38:52 +0000 (14:38 +0200)
33 files changed:
ChangeLog
Makefile
Makefile.CYGWIN
Makefile.LINUX
Makefile.MACOSX
Makefile.WINDOWS
README
examples/example_sprinkler.cpp
include/dai/bipgraph.h
include/dai/doc.h
include/dai/factor.h
include/dai/prob.h
include/dai/properties.h
include/dai/smallset.h
include/dai/weightedgraph.h
src/bbp.cpp
src/bp.cpp
src/bp_dual.cpp
src/daialg.cpp
src/emalg.cpp
src/factor.cpp
src/factorgraph.cpp
src/fbp.cpp
src/gibbs.cpp
src/jtree.cpp
src/mr.cpp
src/treeep.cpp
src/trwbp.cpp
src/weightedgraph.cpp
swig/dai.i
tests/testdai.cpp
tests/unit/prob.cpp
utils/fginfo.cpp

index 3c6095a..06ff5d4 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,12 +1,21 @@
 git HEAD
 --------
 
+* Fixed some bugs in the MatLab interface build system
+* Fixed a bug in utils/fginfo.cpp
+* Improved factor.h/cpp:
+  - 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
 * Improved prob.h/cpp:
   - Added TProb<T>::operator==( const TProb<T> &q )
   - TProb<T>::accumulate() now also applies op2 to init
   - Fixed bug by renaming TProb<T>::operator<=() to TProb<T>::operator<()
   - TProb<T>& operator/= (T x) now yields 0 when dividing by 0
   - Changed format of TProb<T> when streamed to an ostream
+  - 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
 * Improved index.h/cpp:
   - Added multifor::reset()
 * Improved properties.h/cpp:
@@ -29,6 +38,7 @@ git HEAD
   - More error checking in RootedTree constructor
   - The vertices in DEdge and UEdge are now named "first" and "second"
     for compatibility with a std::pair<size_t,size_t>
+  - Added GraphEL( const GraphAL& G ) constructor
 * Improved bipgraph.h/cpp:
   - Fixed bug in BipartiteGraph::eraseNode1()
   - Fixed bug in BipartiteGraph::eraseNode2()
@@ -53,7 +63,7 @@ git HEAD
     no longer has a default value in order to avoid confusion with the
      SmallSet::SmallSet( const T &t1, const T &t2 )
     constructor.
-* Removed RandomDRegularGraph()
+* Removed RandomDRegularGraph(), please use createGraphRegular() instead
 * Compressed Makefile
 * Added unit tests for var, smallset, varset, graph, bipgraph,
   weightedgraph, enum, util, properties, index, prob
index 36a5788..40e1c9b 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -41,11 +41,6 @@ endif
 ifdef WITH_MATLAB
   TARGETS:=$(TARGETS) matlabs
   # Specify the same C++ compiler and flags to mex
-  ifneq ($(OS),WINDOWS)
-    MEXFLAGS=CXX\#$(CC) CXXFLAGS\#'$(CCFLAGS)'
-  else
-    MEXFLAGS=CXX\#$(CC) CXXFLAGS\#"$(CCFLAGS)"
-  endif
   ifdef NEW_MATLAB
     MEXFLAGS:=$(MEXFLAGS) -largeArrayDims
   else
@@ -239,16 +234,16 @@ tests/testbbp$(EE) : tests/testbbp.cpp $(HEADERS) $(LIB)/libdai$(LE)
 ###################
 
 matlab/dai$(ME) : $(SRC)/matlab/dai.cpp $(HEADERS) matlab$(OE) $(LIB)/libdai$(LE)
-       $(MEX) -o$@ $< matlab$(OE) $(LIB)/libdai$(LE)
+       $(MEX) -output $@ $< matlab$(OE) $(LIB)/libdai$(LE)
 
 matlab/dai_readfg$(ME) : $(SRC)/matlab/dai_readfg.cpp $(HEADERS) factorgraph$(OE) matlab$(OE) exceptions$(OE) bipgraph$(OE)
-       $(MEX) -o$@ $< factorgraph$(OE) matlab$(OE) exceptions$(OE) bipgraph$(OE)
+       $(MEX) -output $@ $< factorgraph$(OE) matlab$(OE) exceptions$(OE) bipgraph$(OE)
 
 matlab/dai_writefg$(ME) : $(SRC)/matlab/dai_writefg.cpp $(HEADERS) factorgraph$(OE) matlab$(OE) exceptions$(OE) bipgraph$(OE)
-       $(MEX) -o$@ $< factorgraph$(OE) matlab$(OE) exceptions$(OE) bipgraph$(OE)
+       $(MEX) -output $@ $< factorgraph$(OE) matlab$(OE) exceptions$(OE) bipgraph$(OE)
 
 matlab/dai_potstrength$(ME) : $(SRC)/matlab/dai_potstrength.cpp $(HEADERS) matlab$(OE) exceptions$(OE)
-       $(MEX) -o$@ $< matlab$(OE) exceptions$(OE)
+       $(MEX) -output $@ $< matlab$(OE) exceptions$(OE)
 
 matlab$(OE) : $(SRC)/matlab/matlab.cpp $(INC)/matlab/matlab.h $(HEADERS)
        $(MEX) -c $<
index 0cbc567..287283c 100644 (file)
@@ -57,7 +57,7 @@ MATLABDIR=/agbs/share/sw/matlab
 # The following should resolve to the MatLab mex compile command
 MEX=$(MATLABDIR)/bin/mex
 # Specify the same C++ compiler and flags to mex
-MEXFLAGS=CXX\#$(CC) CXXFLAGS\#'$(CCFLAGS)'
+MEXFLAGS:=CXX\#$(CC) CXXFLAGS\#'$(CCFLAGS)'
 
 # SWIG PYTHON INTERFACE
 # The following should resolve to the SWIG command
index 70145b3..2631e21 100644 (file)
@@ -30,7 +30,13 @@ ME=.mexglx
 
 # COMPILER
 # Compile using GNU C++ Compiler
-CC=g++
+ifdef WITH_MATLAB
+# MatLab R2008b only works with older versions of g++
+  CC=g++-4.1
+else
+# If we don't build the MatLab interface, we prefer a newer version of g++
+  CC=g++
+endif
 # Output filename option of the compiler
 CCO=-o
 # Flags for the C++ compiler
@@ -57,7 +63,7 @@ MATLABDIR=/agbs/share/sw/matlab
 # The following should resolve to the MatLab mex compile command
 MEX=$(MATLABDIR)/bin/mex
 # Specify the same C++ compiler and flags to mex
-MEXFLAGS=CXX\#$(CC) CXXFLAGS\#'$(CCFLAGS)'
+MEXFLAGS:=CXX\#$(CC) CXXFLAGS\#'$(CCFLAGS)'
 
 # SWIG PYTHON INTERFACE
 # The following should resolve to the SWIG command
index baf5abd..e6cbb97 100644 (file)
@@ -56,7 +56,7 @@ MATLABDIR=/agbs/share/sw/matlab
 # The following should resolve to the MatLab mex compile command
 MEX=$(MATLABDIR)/bin/mex
 # Specify the same C++ compiler and flags to mex
-MEXFLAGS=CXX\#$(CC) CXXFLAGS\#'$(CCFLAGS)'
+MEXFLAGS:=CXX\#$(CC) CXXFLAGS\#'$(CCFLAGS)'
 
 # SWIG PYTHON INTERFACE
 # The following should resolve to the SWIG command
index f404676..90a371e 100644 (file)
@@ -58,7 +58,7 @@ MATLABDIR=c:\matlab
 # The following should resolve to the MatLab mex compile command
 MEX=$(MATLABDIR)\bin\mex
 # Specify the same C++ compiler and flags to mex
-MEXFLAGS=CXX\#$(CC) CXXFLAGS\#"$(CCFLAGS)"
+MEXFLAGS:=CXX\#$(CC) CXXFLAGS\#'$(CCFLAGS)'
 
 # SWIG PYTHON INTERFACE
 # The following should resolve to the SWIG command
diff --git a/README b/README
index 65b534e..9dd8b2e 100644 (file)
--- a/README
+++ b/README
@@ -280,7 +280,7 @@ at the MatLab prompt), which performs exact inference by the junction tree
 algorithm and approximate inference by belief propagation on the ALARM network:
 
   cd path_to_libdai/matlab
-  [psi] = dai_readfg ('../examples/alarm.fg');
+  [psi] = dai_readfg ('../tests/alarm.fg');
   [logZ,q,md,qv,qf] = dai (psi, 'JTREE', '[updates=HUGIN,verbose=0]')
   [logZ,q,md,qv,qf] = dai (psi, 'BP', '[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0]')
 
index 74a7441..0b4d74a 100644 (file)
@@ -27,33 +27,33 @@ int main() {
 
     // Define probability distribution for C
     Factor P_C( C );
-    P_C[0] = 0.5;   // C = 0
-    P_C[1] = 0.5;   // C = 1
+    P_C.set(0, 0.5);   // C = 0
+    P_C.set(1, 0.5);   // C = 1
 
     // Define conditional probability of S given C
     Factor P_S_given_C( VarSet( S, C ) );
-    P_S_given_C[0] = 0.5;   // C = 0, S = 0
-    P_S_given_C[1] = 0.9;   // C = 1, S = 0
-    P_S_given_C[2] = 0.5;   // C = 0, S = 1
-    P_S_given_C[3] = 0.1;   // C = 1, S = 1
+    P_S_given_C.set(0, 0.5);   // C = 0, S = 0
+    P_S_given_C.set(1, 0.9);   // C = 1, S = 0
+    P_S_given_C.set(2, 0.5);   // C = 0, S = 1
+    P_S_given_C.set(3, 0.1);   // C = 1, S = 1
 
     // Define conditional probability of R given C
     Factor P_R_given_C( VarSet( R, C ) );
-    P_R_given_C[0] = 0.8;   // C = 0, R = 0
-    P_R_given_C[1] = 0.2;   // C = 1, R = 0
-    P_R_given_C[2] = 0.2;   // C = 0, R = 1
-    P_R_given_C[3] = 0.8;   // C = 1, R = 1
+    P_R_given_C.set(0, 0.8);   // C = 0, R = 0
+    P_R_given_C.set(1, 0.2);   // C = 1, R = 0
+    P_R_given_C.set(2, 0.2);   // C = 0, R = 1
+    P_R_given_C.set(3, 0.8);   // C = 1, R = 1
 
     // Define conditional probability of W given S and R
     Factor P_W_given_S_R( VarSet( S, R ) | W );
-    P_W_given_S_R[0] = 1.0;  // S = 0, R = 0, W = 0
-    P_W_given_S_R[1] = 0.1;  // S = 1, R = 0, W = 0
-    P_W_given_S_R[2] = 0.1;  // S = 0, R = 1, W = 0
-    P_W_given_S_R[3] = 0.01; // S = 1, R = 1, W = 0
-    P_W_given_S_R[4] = 0.0;  // S = 0, R = 0, W = 1
-    P_W_given_S_R[5] = 0.9;  // S = 1, R = 0, W = 1
-    P_W_given_S_R[6] = 0.9;  // S = 0, R = 1, W = 1
-    P_W_given_S_R[7] = 0.99; // S = 1, R = 1, W = 1
+    P_W_given_S_R.set(0, 1.0);  // S = 0, R = 0, W = 0
+    P_W_given_S_R.set(1, 0.1);  // S = 1, R = 0, W = 0
+    P_W_given_S_R.set(2, 0.1);  // S = 0, R = 1, W = 0
+    P_W_given_S_R.set(3, 0.01); // S = 1, R = 1, W = 0
+    P_W_given_S_R.set(4, 0.0);  // S = 0, R = 0, W = 1
+    P_W_given_S_R.set(5, 0.9);  // S = 1, R = 0, W = 1
+    P_W_given_S_R.set(6, 0.9);  // S = 0, R = 1, W = 1
+    P_W_given_S_R.set(7, 0.99); // S = 1, R = 1, W = 1
 
     // Build factor graph consisting of those four factors
     vector<Factor> SprinklerFactors;
index f51c0fd..8ab63e5 100644 (file)
@@ -350,11 +350,13 @@ class BipartiteGraph {
 
         /// Calculates second-order neighbors (i.e., neighbors of neighbors) of node \a n1 of type 1.
         /** If \a include == \c true, includes \a n1 itself, otherwise excludes \a n1.
+         *  \note In libDAI versions 0.2.4 and earlier, this function used to return a std::vector<size_t>
          */
         SmallSet<size_t> delta1( size_t n1, bool include = false ) const;
 
         /// Calculates second-order neighbors (i.e., neighbors of neighbors) of node \a n2 of type 2.
         /** If \a include == \c true, includes \a n2 itself, otherwise excludes \a n2.
+         *  \note In libDAI versions 0.2.4 and earlier, this function used to return a std::vector<size_t>
          */
         SmallSet<size_t> delta2( size_t n2, bool include = false ) const;
 
index cfa6382..d6f5ed2 100644 (file)
  *  at the MatLab prompt), which performs exact inference by the junction tree
  *  algorithm and approximate inference by belief propagation on the ALARM network:
  *  <pre>  cd path_to_libdai/matlab
- *  [psi] = dai_readfg ('../examples/alarm.fg');
+ *  [psi] = dai_readfg ('../tests/alarm.fg');
  *  [logZ,q,md,qv,qf] = dai (psi, 'JTREE', '[updates=HUGIN,verbose=0]')
  *  [logZ,q,md,qv,qf] = dai (psi, 'BP', '[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0]')</pre>
  *  where "path_to_libdai" has to be replaced with the directory in which libDAI
index 7c30825..cfb616c 100644 (file)
@@ -102,7 +102,7 @@ template <typename T> class TFactor {
         TFactor( const std::vector<Var> &vars, const std::vector<T> &p ) : _vs(vars.begin(), vars.end(), vars.size()), _p(p.size()) {
             Permute permindex(vars);
             for( size_t li = 0; li < p.size(); ++li )
-                _p[permindex.convertLinearIndex(li)] = p[li];
+                _p.set( permindex.convertLinearIndex(li), p[li] );
         }
     //@}
 
@@ -118,8 +118,15 @@ template <typename T> class TFactor {
         T operator[] (size_t i) const { return _p[i]; }
 
         /// Returns a reference to the \a i 'th entry of the value vector
+        /// \deprecated Please use dai::TFactor::set() instead
         T& operator[] (size_t i) { return _p[i]; }
 
+        /// Gets \a i 'th entry of the value vector
+        T get( size_t i ) const { return _p[i]; }
+
+        /// Sets \a i 'th entry of the value vector to \a val
+        void set( size_t i, T val ) { _p.set( i, val ); }
+
         /// Returns constant reference to variable set (i.e., the variables on which the factor depends)
         const VarSet& vars() const { return _vs; }
 
@@ -460,7 +467,7 @@ template<typename T> TFactor<T> TFactor<T>::slice( const VarSet& vars, size_t va
     IndexFor i_varsrem (varsrem, _vs);
     for( size_t i = 0; i < states(); i++, ++i_vars, ++i_varsrem )
         if( (size_t)i_vars == varsState )
-            result._p[i_varsrem] = _p[i];
+            result.set( i_varsrem, _p[i] );
 
     return result;
 }
@@ -473,7 +480,7 @@ template<typename T> TFactor<T> TFactor<T>::marginal(const VarSet &vars, bool no
 
     IndexFor i_res( res_vars, _vs );
     for( size_t i = 0; i < _p.size(); i++, ++i_res )
-        res._p[i_res] += _p[i];
+        res.set( i_res, res[i_res] + _p[i] );
 
     if( normed )
         res.normalize( TProb<T>::NORMPROB );
@@ -490,7 +497,7 @@ template<typename T> TFactor<T> TFactor<T>::maxMarginal(const VarSet &vars, bool
     IndexFor i_res( res_vars, _vs );
     for( size_t i = 0; i < _p.size(); i++, ++i_res )
         if( _p[i] > res._p[i_res] )
-            res._p[i_res] = _p[i];
+            res.set( i_res, _p[i] );
 
     if( normed )
         res.normalize( TProb<T>::NORMPROB );
index 2036e17..81ca241 100644 (file)
@@ -193,10 +193,15 @@ template<typename T> struct fo_absdiff : public std::binary_function<T, T, T> {
  *
  *  \tparam T Should be a scalar that is castable from and to dai::Real and should support elementary arithmetic operations.
  */
-template <typename T> class TProb {
+template <typename T>
+class TProb {
+    public:
+        /// Type of data structure used for storing the values
+        typedef std::vector<T> container_type;
+
     private:
-        /// The vector
-        std::vector<T> _p;
+        /// The data structure that stores the values
+        container_type _p;
 
     public:
         /// Enumerates different ways of normalizing a probability measure.
@@ -221,10 +226,10 @@ template <typename T> class TProb {
         TProb() : _p() {}
 
         /// Construct uniform probability distribution over \a n outcomes (i.e., a vector of length \a n with each entry set to \f$1/n\f$)
-        explicit TProb( size_t n ) : _p(std::vector<T>(n, (T)1 / n)) {}
+        explicit TProb( size_t n ) : _p( n, (T)1 / n ) {}
 
         /// Construct vector of length \a n with each entry set to \a p
-        explicit TProb( size_t n, T p ) : _p(n, p) {}
+        explicit TProb( size_t n, T p ) : _p( n, p ) {}
 
         /// Construct vector from a range
         /** \tparam TIterator Iterates over instances that can be cast to \a T
@@ -240,7 +245,7 @@ template <typename T> class TProb {
 
         /// Construct vector from another vector
         /** \tparam S type of elements in \a v (should be castable to type \a T)
-         *  \param v vector used for initialization
+         *  \param v vector used for initialization.
          */
         template <typename S>
         TProb( const std::vector<S> &v ) : _p() {
@@ -250,13 +255,13 @@ template <typename T> class TProb {
     //@}
 
         /// Constant iterator over the elements
-        typedef typename std::vector<T>::const_iterator const_iterator;
+        typedef typename container_type::const_iterator const_iterator;
         /// Iterator over the elements
-        typedef typename std::vector<T>::iterator iterator;
+        typedef typename container_type::iterator iterator;
         /// Constant reverse iterator over the elements
-        typedef typename std::vector<T>::const_reverse_iterator const_reverse_iterator;
+        typedef typename container_type::const_reverse_iterator const_reverse_iterator;
         /// Reverse iterator over the elements
-        typedef typename std::vector<T>::reverse_iterator reverse_iterator;
+        typedef typename container_type::reverse_iterator reverse_iterator;
 
     /// \name Iterator interface
     //@{
@@ -283,14 +288,8 @@ template <typename T> class TProb {
 
     /// \name Queries
     //@{
-        /// Returns a const reference to the wrapped vector
-        const std::vector<T> & p() const { return _p; }
-
-        /// Returns a reference to the wrapped vector
-        std::vector<T> & p() { return _p; }
-
-        /// Returns a copy of the \a i 'th entry
-        T operator[]( size_t i ) const {
+        /// Gets \a i 'th entry
+        T get( size_t i ) const { 
 #ifdef DAI_DEBUG
             return _p.at(i);
 #else
@@ -298,7 +297,24 @@ template <typename T> class TProb {
 #endif
         }
 
+        /// Sets \a i 'th entry to \a val
+        void set( size_t i, T val ) {
+            DAI_DEBASSERT( i < _p.size() );
+            _p[i] = val;
+        }
+
+        /// Returns a const reference to the wrapped vector
+        const container_type& p() const { return _p; }
+
+        /// Returns a reference to the wrapped vector
+        container_type& p() { return _p; }
+
+        /// Returns a copy of the \a i 'th entry
+        T operator[]( size_t i ) const { return get(i); }
+
         /// Returns reference to the \a i 'th entry
+        /** \deprecated Please use dai::TProb::set() instead
+         */
         T& operator[]( size_t i ) { return _p[i]; }
 
         /// Returns length of the vector (i.e., the number of entries)
@@ -341,7 +357,7 @@ template <typename T> class TProb {
         /// Returns \c true if one or more entries are NaN
         bool hasNaNs() const {
             bool foundnan = false;
-            for( typename std::vector<T>::const_iterator x = _p.begin(); x != _p.end(); x++ )
+            for( const_iterator x = _p.begin(); x != _p.end(); x++ )
                 if( isnan( *x ) ) {
                     foundnan = true;
                     break;
@@ -364,7 +380,7 @@ template <typename T> class TProb {
                 arg = i;
               }
             }
-            return std::make_pair(arg,max);
+            return std::make_pair( arg, max );
         }
 
         /// Returns a random index, according to the (normalized) distribution described by *this
@@ -372,7 +388,7 @@ template <typename T> class TProb {
             Real x = rnd_uniform() * sum();
             T s = 0;
             for( size_t i = 0; i < size(); i++ ) {
-                s += _p[i];
+                s += get(i);
                 if( s > x )
                     return i;
             }
@@ -423,7 +439,7 @@ template <typename T> class TProb {
             else
                 return pwUnaryTr( fo_log<T>() );
         }
-        
+
         /// Returns pointwise inverse
         /** If \a zero == \c true, uses <tt>1/0==0</tt>; otherwise, <tt>1/0==Inf</tt>.
          */
@@ -507,7 +523,7 @@ template <typename T> class TProb {
     /// \name Operations with scalars
     //@{
         /// Sets all entries to \a x
-        TProb<T> & fill(T x) {
+        TProb<T> & fill( T x ) {
             std::fill( _p.begin(), _p.end(), x );
             return *this;
         }
@@ -703,7 +719,7 @@ template<typename T> T dist( const TProb<T> &p, const TProb<T> &q, typename TPro
  */
 template<typename T> std::ostream& operator<< (std::ostream& os, const TProb<T>& p) {
     os << "(";
-    for( typename std::vector<T>::const_iterator it = p.begin(); it != p.end(); it++ )
+    for( typename TProb<T>::const_iterator it = p.begin(); it != p.end(); it++ )
         os << (it != p.begin() ? ", " : "") << *it;
     os << ")";
     return os;
index 289b7a3..d75f80d 100644 (file)
@@ -98,6 +98,13 @@ class PropertySet : private std::map<PropertyKey, PropertyValue> {
             return *this; 
         }
 
+        /// Sets a property (a key \a key with a corresponding value \a val)
+        /** \deprecated Please use dai::PropertySet::set(const PropertyKey&, const PropertyValue&) instead
+         */
+        PropertySet& Set( const PropertyKey& key, const PropertyValue& val ) { 
+            return set( key, val );
+        }
+
         /// Set properties according to \a newProps, overriding properties that already exist with new values
         PropertySet& set( const PropertySet& newProps ) {
             const std::map<PropertyKey, PropertyValue> *m = &newProps;
@@ -106,6 +113,13 @@ class PropertySet : private std::map<PropertyKey, PropertyValue> {
             return *this;
         }
 
+        /// Set properties according to \a newProps, overriding properties that already exist with new values
+        /** \deprecated Please use dai::PropertySet::set(const PropertySet&) instead
+         */
+        PropertySet& Set( const PropertySet& newProps ) {
+            return set( newProps );
+        }
+
         /// Shorthand for (temporarily) adding properties
         /** \par Example:
             \code
@@ -148,6 +162,17 @@ PropertySet p()("method","BP")("verbose",1)("tol",1e-9)
                 }
             }
         }
+
+        /// Converts the type of the property value corresponding with \a key from string to \a ValueType (if necessary)
+        /** The implementation makes use of boost::lexical_cast
+         *  \tparam ValueType Type to which the value should be cast
+         *  \throw IMPOSSIBLE_TYPECAST if the type cast cannot be done
+         *  \deprecated Please use dai::PropertySet::convertTo() instead
+         */
+        template<typename ValueType>
+        void ConvertTo( const PropertyKey& key ) { 
+            convertTo<ValueType>( key );
+        }
     //@}
 
     //@}
@@ -194,6 +219,14 @@ PropertySet p()("method","BP")("verbose",1)("tol",1e-9)
             return x->second;
         }
 
+        /// Gets the value corresponding to \a key
+        /** \throw OBJECT_NOT_FOUND if the key cannot be found in \c *this
+         *  \deprecated Please use dai::PropertySet::get() instead
+         */
+        const PropertyValue& Get( const PropertyKey& key ) const {
+            return get( key );
+        }
+
         /// Gets the value corresponding to \a key, cast to \a ValueType
         /** \tparam ValueType Type to which the value should be cast
          *  \throw OBJECT_NOT_FOUND if the key cannot be found in \c *this
@@ -209,6 +242,17 @@ PropertySet p()("method","BP")("verbose",1)("tol",1e-9)
             }
         }
 
+        /// Gets the value corresponding to \a key, cast to \a ValueType
+        /** \tparam ValueType Type to which the value should be cast
+         *  \throw OBJECT_NOT_FOUND if the key cannot be found in \c *this
+         *  \throw IMPOSSIBLE_TYPECAST if the type cast cannot be done
+         *  \deprecated Please use dai::PropertySet::getAs() instead
+         */
+        template<typename ValueType>
+        ValueType GetAs( const PropertyKey& key ) const {
+            return getAs<ValueType>( key );
+        }
+
         /// Gets the value corresponding to \a key, cast to \a ValueType, converting from a string if necessary
         /** If the type of the value is already equal to \a ValueType, no conversion is done.
          *  Otherwise, the type of the value should be a std::string, in which case boost::lexical_cast is
index dc49a90..d00ec2e 100644 (file)
@@ -65,6 +65,7 @@ class SmallSet {
          *  \param begin Points to first element to be added.
          *  \param end Points just beyond last element to be added.
          *  \param sizeHint For efficiency, the number of elements can be speficied by \a sizeHint.
+         *  \note The \a sizeHint parameter used to be optional in libDAI versions 0.2.4 and earlier.
          */
         template <typename TIterator>
         SmallSet( TIterator begin, TIterator end, size_t sizeHint ) {
index 78362ea..fd58358 100644 (file)
@@ -27,6 +27,7 @@
 #include <climits>   // Work-around for bug in boost graph library
 #include <dai/util.h>
 #include <dai/exceptions.h>
+#include <dai/graph.h>
 
 #include <boost/graph/adjacency_list.hpp>
 #include <boost/graph/prim_minimum_spanning_tree.hpp>
@@ -40,10 +41,18 @@ namespace dai {
 class DEdge {
     public:
         /// First node index (source of edge)
-        size_t first;
+        union {
+            size_t first;
+            /// \deprecated Please use member dai::DEdge::first instead
+            size_t n1;
+        };
 
         /// Second node index (target of edge)
-        size_t second;
+        union {
+            size_t second;
+            /// \deprecated Please use member dai::DEdge::second instead
+            size_t n2;
+        };
 
         /// Default constructor
         DEdge() : first(0), second(0) {}
@@ -71,10 +80,18 @@ class DEdge {
 class UEdge {
     public:
         /// First node index
-        size_t first;
+        union {
+            size_t first;
+            /// \deprecated Please use member dai::UEdge::first instead
+            size_t n1;
+        };
 
         /// Second node index
-        size_t second;
+        union {
+            size_t second;
+            /// \deprecated Please use member dai::UEdge::second instead
+            size_t n2;
+        };
 
         /// Default constructor
         UEdge() : first(0), second(0) {}
@@ -121,6 +138,14 @@ class GraphEL : public std::set<UEdge> {
         GraphEL( InputIterator begin, InputIterator end ) {
             insert( begin, end );
         }
+
+        /// Construct from GraphAL
+        GraphEL( const GraphAL& G ) {
+            for( size_t n1 = 0; n1 < G.nrNodes(); n1++ )
+                foreach( const GraphAL::Neighbor n2, G.nb(n1) )
+                    if( n1 < n2 )
+                        insert( UEdge( n1, n2 ) );
+        }
 };
 
 
@@ -235,6 +260,39 @@ template<typename T> RootedTree MaxSpanningTree( const WeightedGraph<T> &G, bool
 }
 
 
+/// Constructs a minimum spanning tree from the (non-negatively) weighted graph \a G using Prim's algorithm.
+/** \param G Weighted graph that should have non-negative weights.
+ *  \note Uses implementation from Boost Graph Library.
+ *  \note The vertices of \a G must be in the range [0,N) where N is the number of vertices of \a G.
+ *  \deprecated Please use dai::MinSpanningTree(const WeightedGraph&, bool) instead
+ */
+template<typename T> RootedTree MinSpanningTree( const WeightedGraph<T> &G ) {
+    return MinSpanningTree( G, true );
+}
+
+
+/// Constructs a minimum spanning tree from the (non-negatively) weighted graph \a G using Prim's algorithm.
+/** \param G Weighted graph that should have non-negative weights.
+ *  \note Uses implementation from Boost Graph Library.
+ *  \note The vertices of \a G must be in the range [0,N) where N is the number of vertices of \a G.
+ *  \deprecated Please use dai::MinSpanningTree(const WeightedGraph&, bool) instead
+ */
+template<typename T> RootedTree MaxSpanningTree( const WeightedGraph<T> &G ) {
+    return MaxSpanningTree( G, true );
+}
+
+
+/// Constructs a random undirected graph of \a N nodes, where each node has connectivity \a d
+/** Algorithm 1 in [\ref StW99].
+ *  Draws a random graph of size \a N and uniform degree \a d
+ *  from an almost uniform probability distribution over these graphs
+ *  (which becomes uniform in the limit that \a d is small and \a N goes
+ *  to infinity).
+ *  \deprecated Please use dai::createGraphRegular(size_t, size_t) instead
+ */
+GraphEL RandomDRegularGraph( size_t N, size_t d );
+
+
 } // end of namespace dai
 
 
index b0e0139..fd643f9 100644 (file)
@@ -195,7 +195,7 @@ void BBP::RegenerateU() {
                     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]] );
                 }
             _Umsg[I][i.iter] = prod;
         }
@@ -218,7 +218,7 @@ void BBP::RegenerateS() {
                             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]];
+                                prod.set( x_I, prod[x_I] * p[ind[x_I]] );
                         }
                     }
                     // "Marginalize" onto i|j (unnormalized)
@@ -285,7 +285,7 @@ void BBP::RegeneratePsiAdjoints() {
             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 );
@@ -319,12 +319,12 @@ 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 );
             }
@@ -376,12 +376,12 @@ 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 );
         }
     }
@@ -414,7 +414,7 @@ void BBP::calcNewN( size_t i, size_t _I ) {
             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);
@@ -431,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();
@@ -571,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;
 
@@ -589,7 +589,7 @@ void BBP::sendSeqMsgM( size_t j, size_t _I ) {
             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:
@@ -629,7 +629,7 @@ Prob BBP::unnormAdjoint( const Prob &w, Real Z_w, const Prob &adj_w ) {
     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;
+        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;
 }
@@ -769,22 +769,22 @@ void BBP::initCostFnAdj( const BBPCostFunction &cfn, const vector<size_t> *state
                 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] ));
+                    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[xi] = -_ia->beliefV(i)[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).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] );
+                    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[xI] = -_ia->beliefF(I)[xI] / fg.factor(I).p()[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 );
@@ -798,9 +798,9 @@ void BBP::initCostFnAdj( const BBPCostFunction &cfn, const vector<size_t> *state
                 for( size_t xI = 0; xI < dim; xI++ ) {
                     Real bIxI = _ia->beliefF(I)[xI];
                     if( bIxI < 1.0e-15 )
-                        p[xI] = -1.0e10;
+                        p.set( xI, -1.0e10 );
                     else
-                        p[xI] = 1 + log( bIxI );
+                        p.set( xI, 1 + log( bIxI ) );
                 }
                 b2_adj.push_back(p);
             }
@@ -815,9 +815,9 @@ void BBP::initCostFnAdj( const BBPCostFunction &cfn, const vector<size_t> *state
                 for( size_t xi = 0; xi < fg.var(i).states(); xi++ ) {
                     Real bixi = _ia->beliefV(i)[xi];
                     if( bixi < 1.0e-15 )
-                        p[xi] = -1.0e10;
+                        p.set( xi, -1.0e10 );
                     else
-                        p[xi] = 1 + log( bixi );
+                        p.set( xi, 1 + log( bixi ) );
                 }
                 b1_adj.push_back( p );
             }
@@ -843,13 +843,13 @@ void BBP::initCostFnAdj( const BBPCostFunction &cfn, const vector<size_t> *state
                 Real b = _ia->beliefV(i)[state[i]];
                 switch( (size_t)cfn ) {
                     case BBPCostFunction::CFN_GIBBS_B:
-                        delta[state[i]] = 1.0;
+                        delta.set( state[i], 1.0 );
                         break;
                     case BBPCostFunction::CFN_GIBBS_B2:
-                        delta[state[i]] = b;
+                        delta.set( state[i], b );
                         break;
                     case BBPCostFunction::CFN_GIBBS_EXP:
-                        delta[state[i]] = exp(b);
+                        delta.set( state[i], exp(b) );
                         break;
                     default:
                         DAI_THROW(UNKNOWN_ENUM_VALUE);
@@ -881,13 +881,13 @@ void BBP::initCostFnAdj( const BBPCostFunction &cfn, const vector<size_t> *state
                 Real b = _ia->beliefF(I)[x_I];
                 switch( (size_t)cfn ) {
                     case BBPCostFunction::CFN_GIBBS_B_FACTOR:
-                        delta[x_I] = 1.0;
+                        delta.set( x_I, 1.0 );
                         break;
                     case BBPCostFunction::CFN_GIBBS_B2_FACTOR:
-                        delta[x_I] = b;
+                        delta.set( x_I, b );
                         break;
                     case BBPCostFunction::CFN_GIBBS_EXP_FACTOR:
-                        delta[x_I] = exp( b );
+                        delta.set( x_I, exp( b ) );
                         break;
                     default:
                         DAI_THROW(UNKNOWN_ENUM_VALUE);
@@ -1011,7 +1011,7 @@ Real numericBBPTest( const InfAlg &bp, const std::vector<size_t> *state, const P
                 // 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 );
index 91951a2..efa45ce 100644 (file)
@@ -175,9 +175,9 @@ Prob BP::calcIncomingMessageProduct( size_t I, bool without_i, size_t i ) const
 
                 for( size_t r = 0; r < prod.size(); ++r )
                     if( props.logdomain )
-                        prod[r] += prod_j[ind[r]];
+                        prod.set( r, prod[r] + prod_j[ind[r]] );
                     else
-                        prod[r] *= prod_j[ind[r]];
+                        prod.set( r, prod[r] * prod_j[ind[r]] );
             }
     }
     return prod;
@@ -203,23 +203,23 @@ void BP::calcNewMessage( size_t i, size_t _I ) {
 
         // Marginalize onto i
         if( !DAI_BP_FAST ) {
-            /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
+            // UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION
             if( props.inference == Properties::InfType::SUMPROD )
                 marg = Fprod.marginal( var(i) ).p();
             else
                 marg = Fprod.maxMarginal( var(i) ).p();
         } else {
-            /* OPTIMIZED VERSION */
+            // OPTIMIZED VERSION 
             marg = Prob( var(i).states(), 0.0 );
             // ind is the precalculated IndexFor(i,I) i.e. to x_I == k corresponds x_i == ind[k]
             const ind_t ind = index(i,_I);
             if( props.inference == Properties::InfType::SUMPROD )
                 for( size_t r = 0; r < prod.size(); ++r )
-                    marg[ind[r]] += prod[r];
+                    marg.set( ind[r], marg[ind[r]] + prod[r] );
             else
                 for( size_t r = 0; r < prod.size(); ++r )
                     if( prod[r] > marg[ind[r]] )
-                        marg[ind[r]] = prod[r];
+                        marg.set( ind[r], prod[r] );
             marg.normalize();
         }
     }
index d282d2d..bf9f714 100644 (file)
@@ -93,14 +93,14 @@ void BP_dual::calcNewM( size_t i, size_t _I ) {
             Prob &n = msgN(j,j.dual);
             IndexFor ind( fg().var(j), fg().factor(I).vars() );
             for( size_t x = 0; ind.valid(); x++, ++ind )
-                prod[x] *= n[ind];
+                prod.set( x, prod[x] * n[ind] );
         }
     // Marginalize onto i
     Prob marg( fg().var(i).states(), 0.0 );
     // ind is the precalculated Index(i,I) i.e. to x_I == k corresponds x_i == ind[k]
     IndexFor ind( fg().var(i), fg().factor(I).vars() );
     for( size_t x = 0; ind.valid(); x++, ++ind )
-        marg[ind] += prod[x];
+        marg.set( ind, marg[ind] + prod[x] );
 
     _msgs.Zm[i][_I] = marg.normalize();
     _msgs.m[i][_I] = marg;
@@ -142,7 +142,7 @@ void BP_dual::calcBeliefF( size_t I ) {
         IndexFor ind( fg().var(j), fg().factor(I).vars() );
         Prob n( msgN(j,j.dual) );
         for( size_t x = 0; ind.valid(); x++, ++ind )
-            prod[x] *= n[ind];
+            prod.set( x, prod[x] * n[ind] );
     }
     _beliefs.Zb2[I] = prod.normalize();
     _beliefs.b2[I] = prod;
index 02993c6..0e6b218 100644 (file)
@@ -61,9 +61,9 @@ Factor calcMarginal( const InfAlg &obj, const VarSet &vs, bool reInit ) {
                 logZ0 = logZ;
 
         if( logZ == -INFINITY )
-            Pvs[s] = 0;
+            Pvs.set( s, 0 );
         else
-            Pvs[s] = exp(logZ - logZ0); // subtract logZ0 to avoid very large numbers
+            Pvs.set( s, exp(logZ - logZ0) ); // subtract logZ0 to avoid very large numbers
 
         // restore clamped factors
         clamped->restoreFactors( vs );
@@ -132,7 +132,7 @@ vector<Factor> calcPairBeliefs( const InfAlg & obj, const VarSet& vs, bool reIni
 
                         // we assume that j.label() < k.label()
                         // i.e. we make an assumption here about the indexing
-                        pairbelief[j_val + (k_val * nj->states())] = Z_xj;
+                        pairbelief.set( j_val + (k_val * nj->states()), Z_xj );
 
                         // restore clamped factors
                         clamped->restoreFactors( vs );
@@ -190,9 +190,9 @@ vector<Factor> calcPairBeliefs( const InfAlg & obj, const VarSet& vs, bool reIni
                         Factor b_k = clamped->belief(vvs[k]);
                         for( size_t k_val = 0; k_val < vvs[k].states(); k_val++ )
                             if( vvs[j].label() < vvs[k].label() )
-                                pairbeliefs[j * N + k][j_val + (k_val * vvs[j].states())] = Z_xj * b_k[k_val];
+                                pairbeliefs[j * N + k].set( j_val + (k_val * vvs[j].states()), Z_xj * b_k[k_val] );
                             else
-                                pairbeliefs[j * N + k][k_val + (j_val * vvs[k].states())] = Z_xj * b_k[k_val];
+                                pairbeliefs[j * N + k].set( k_val + (j_val * vvs[k].states()), Z_xj * b_k[k_val] );
                     }
 
                 // restore clamped factors
index 3e50ec2..3732190 100644 (file)
@@ -64,12 +64,14 @@ Prob CondProbEstimation::estimate() {
     for( size_t parent = 0; parent < _stats.size(); parent += _target_dim ) {
         // calculate norm
         size_t top = parent + _target_dim;
-        Real norm = std::accumulate( &(_stats[parent]), &(_stats[parent]) + _target_dim, 0.0 );
+        Real norm = 0.0;
+        for( size_t i = parent; i < top; ++i )
+            norm += _stats[i];
         if( norm != 0.0 )
             norm = 1.0 / norm;
         // normalize
         for( size_t i = parent; i < top; ++i )
-            _stats[i] *= norm;
+            _stats.set( i, _stats[i] * norm );
     }
     // reset _stats to _initial_stats
     Prob result = _stats;
@@ -173,7 +175,7 @@ void SharedParameters::collectSufficientStatistics( InfAlg &alg ) {
         Factor b = alg.belief(vs);
         Prob p( b.states(), 0.0 );
         for( size_t entry = 0; entry < b.states(); ++entry )
-            p[entry] = b[perm.convertLinearIndex(entry)]; // apply inverse permutation
+            p.set( entry, b[perm.convertLinearIndex(entry)] ); // apply inverse permutation
         _estimation->addSufficientStatistics( p );
     }
 }
@@ -187,7 +189,7 @@ void SharedParameters::setParameters( FactorGraph &fg ) {
 
         Factor f( vs, 0.0 );
         for( size_t entry = 0; entry < f.states(); ++entry )
-            f[perm.convertLinearIndex(entry)] = p[entry];
+            f.set( perm.convertLinearIndex(entry), p[entry] );
 
         fg.setFactor( i->first, f );
     }
index 46c73fa..504b506 100644 (file)
@@ -41,7 +41,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++ )
-        fac[t] = std::exp(rnd_stdnormal() * beta);
+        fac.set( t, std::exp(rnd_stdnormal() * beta) );
     return fac;
 }
 
@@ -50,7 +50,7 @@ Factor createFactorPotts( const Var &n1, const Var &n2, Real J ) {
     Factor fac( VarSet( n1, n2 ), 1.0 );
     DAI_ASSERT( n1.states() == n2.states() );
     for( size_t s = 0; s < n1.states(); s++ )
-        fac[ s * (n1.states() + 1) ] = std::exp(J);
+        fac.set( s * (n1.states() + 1), std::exp(J) );
     return fac;
 }
 
@@ -58,7 +58,7 @@ Factor createFactorPotts( const Var &n1, const Var &n2, Real J ) {
 Factor createFactorDelta( const Var &v, size_t state ) {
     Factor fac( v, 0.0 );
     DAI_ASSERT( state < v.states() );
-    fac[state] = 1.0;
+    fac.set( state, 1.0 );
     return fac;
 }
 
index 76c1c18..c4bbf26 100644 (file)
@@ -188,7 +188,7 @@ std::istream& operator>> ( std::istream& is, FactorGraph &fg ) {
             is >> val;
 
             // store value, but permute indices first according to internal representation
-            facs.back()[permindex.convertLinearIndex( li )] = val;
+            facs.back().set( permindex.convertLinearIndex( li ), val );
         }
     }
 
@@ -292,7 +292,7 @@ vector<VarSet> FactorGraph::Cliques() const {
 void FactorGraph::clamp( size_t i, size_t x, bool backup ) {
     DAI_ASSERT( x <= var(i).states() );
     Factor mask( var(i), (Real)0 );
-    mask[x] = (Real)1;
+    mask.set( x, (Real)1 );
 
     map<size_t, Factor> newFacs;
     foreach( const Neighbor &I, nbV(i) )
@@ -309,7 +309,7 @@ void FactorGraph::clampVar( size_t i, const vector<size_t> &is, bool backup ) {
 
     foreach( size_t i, is ) {
         DAI_ASSERT( i <= n.states() );
-        mask_n[i] = (Real)1;
+        mask_n.set( i, (Real)1 );
     }
 
     map<size_t, Factor> newFacs;
@@ -325,7 +325,7 @@ void FactorGraph::clampFactor( size_t I, const vector<size_t> &is, bool backup )
 
     foreach( size_t i, is ) {
         DAI_ASSERT( i <= st );
-        newF[i] = factor(I)[i];
+        newF.set( i, factor(I)[i] );
     }
 
     setFactor( I, newF, backup );
index 0135899..bf140b6 100644 (file)
@@ -94,9 +94,9 @@ Prob FBP::calcIncomingMessageProduct( size_t I, bool without_i, size_t i ) const
 
                 for( size_t r = 0; r < prod.size(); ++r )
                     if( props.logdomain )
-                        prod[r] += prod_j[ind[r]];
+                        prod.set( r, prod[r] + prod_j[ind[r]] );
                     else
-                        prod[r] *= prod_j[ind[r]];
+                        prod.set( r, prod[r] * prod_j[ind[r]] );
             }
     }
     return prod;
@@ -134,11 +134,11 @@ void FBP::calcNewMessage( size_t i, size_t _I ) {
         const ind_t ind = index(i,_I);
         if( props.inference == Properties::InfType::SUMPROD )
             for( size_t r = 0; r < prod.size(); ++r )
-                marg[ind[r]] += prod[r];
+                marg.set( ind[r], marg[ind[r]] + prod[r] );
         else
             for( size_t r = 0; r < prod.size(); ++r )
                 if( prod[r] > marg[ind[r]] )
-                    marg[ind[r]] = prod[r];
+                    marg.set( ind[r], prod[r] );
         marg.normalize();
     }
 
index cc18e14..1e96b98 100644 (file)
@@ -131,7 +131,7 @@ Prob Gibbs::getVarDist( size_t i ) {
         size_t I_skip = getFactorEntryDiff( I, i );
         size_t I_entry = getFactorEntry(I) - (_state[i] * I_skip);
         for( size_t st_i = 0; st_i < i_states; st_i++ ) {
-            i_given_MB[st_i] *= f_I[I_entry];
+            i_given_MB.set( st_i, i_given_MB[st_i] * f_I[I_entry] );
             I_entry += I_skip;
         }
     }
index c4d13fb..5c034f8 100644 (file)
@@ -492,7 +492,7 @@ Factor JTree::calcMarginal( const VarSet& vs ) {
                     for( VarSet::const_iterator n = vsrem.begin(); n != vsrem.end(); n++ )
                         if( Qa[T[i].second].vars() >> *n ) {
                             Factor piet( *n, 0.0 );
-                            piet[s(*n)] = 1.0;
+                            piet.set( s(*n), 1.0 );
                             Qa[T[i].second] *= piet;
                         }
 
@@ -504,7 +504,7 @@ Factor JTree::calcMarginal( const VarSet& vs ) {
                 logZ += log(Qa[T[0].first].normalize());
 
                 Factor piet( vsrem, 0.0 );
-                piet[s] = exp(logZ);
+                piet.set( s, exp(logZ) );
                 Pvs += piet * Qa[T[0].first].marginal( vs / vsrem, false );      // OPTIMIZE ME
 
                 // Restore clamped beliefs
index 763ee66..d175f6b 100644 (file)
@@ -299,8 +299,8 @@ Real MR::calcCavityCorrelations() {
             foreach( const Neighbor &j, G.nb(i) ) {
                 // Create weights for magnetization of some spin
                 Prob p( 2, 0.0 );
-                p[0] = -1.0;
-                p[1] = 1.0;
+                p.set( 0, -1.0 );
+                p.set( 1, 1.0 );
 
                 // BBP cost function would be the magnetization of spin j
                 vector<Prob> b1_adj;
@@ -378,7 +378,7 @@ Real MR::run() {
 
 Factor MR::beliefV( size_t i ) const {
     if( supported ) {
-        Prob x(2);
+        Real x[2];
         x[0] = 0.5 - Mag[i] / 2.0;
         x[1] = 0.5 + Mag[i] / 2.0;
 
index b9f17df..9030b24 100644 (file)
@@ -319,11 +319,8 @@ void TreeEP::TreeEPSubTree::HUGIN_with_I( std::vector<Factor> &Qa, std::vector<F
         for( size_t i = _RTree.size(); (i--) != 0; ) {
             // clamp variables in nsrem
             for( VarSet::const_iterator n = _nsrem.begin(); n != _nsrem.end(); n++ )
-                if( _Qa[_RTree[i].second].vars() >> *n ) {
-                    Factor delta( *n, 0.0 );
-                    delta[s(*n)] = 1.0;
-                    _Qa[_RTree[i].second] *= delta;
-                }
+                if( _Qa[_RTree[i].second].vars() >> *n )
+                    _Qa[_RTree[i].second] *= createFactorDelta( *n, s(*n) );
             Factor new_Qb = _Qa[_RTree[i].second].marginal( _Qb[i].vars(), false );
             _Qa[_RTree[i].first] *= new_Qb / _Qb[i];
             _Qb[i] = new_Qb;
index 34c2943..17041a4 100644 (file)
@@ -119,12 +119,11 @@ Prob TRWBP::calcIncomingMessageProduct( size_t I, bool without_i, size_t i ) con
                 // ind is the precalculated IndexFor(j,I) i.e. to x_I == k corresponds x_j == ind[k]
                 const ind_t &ind = index(j, _I);
 
-                for( size_t r = 0; r < prod.size(); ++r ) {
+                for( size_t r = 0; r < prod.size(); ++r )
                     if( props.logdomain )
-                        prod[r] += prod_j[ind[r]];
+                        prod.set( r, prod[r] + prod_j[ind[r]] );
                     else
-                        prod[r] *= prod_j[ind[r]];
-                }
+                        prod.set( r, prod[r] * prod_j[ind[r]] );
             }
         }
     
index 3674f5a..d69d6ad 100644 (file)
@@ -75,4 +75,9 @@ RootedTree::RootedTree( const GraphEL &T, size_t Root ) {
 }
 
 
+GraphEL RandomDRegularGraph( size_t N, size_t d ) {
+    return GraphEL( createGraphRegular( N, d ) );
+}
+
+
 } // end of namespace dai
index 84e1ffe..c276d0c 100644 (file)
 %include "../include/dai/prob.h"
 %template(Prob) dai::TProb<dai::Real>;
 %extend dai::TProb<dai::Real> {
-        inline dai::Real __getitem__(int i) const {return (*self)[i];} /* for python */
-        inline void __setitem__(int i,dai::Real d) {(*self)[i] = d;}   /* for python */
-        inline dai::Real __paren(int i) const {return (*self)[i];}     /* for octave */
-        inline void __paren_asgn(int i,dai::Real d) {(*self)[i] = d;}  /* for octave */
+        inline dai::Real __getitem__(int i) const {return (*self).get(i);} /* for python */
+        inline void __setitem__(int i,dai::Real d) {(*self).set(i,d);}   /* for python */
+        inline dai::Real __paren(int i) const {return (*self).get(i);}     /* for octave */
+        inline void __paren_asgn(int i,dai::Real d) {(*self).set(i,d);}  /* for octave */
 };
 %include "../include/dai/factor.h"
 %extend dai::TFactor<dai::Real> {
-        inline dai::Real __getitem__(int i) const {return (*self)[i];} /* for python */
-        inline void __setitem__(int i,dai::Real d) {(*self)[i] = d;}   /* for python */
-        inline dai::Real __paren__(int i) const {return (*self)[i];}     /* for octave */
-        inline void __paren_asgn__(int i,dai::Real d) {(*self)[i] = d;}  /* for octave */
+        inline dai::Real __getitem__(int i) const {return (*self).get(i);} /* for python */
+        inline void __setitem__(int i,dai::Real d) {(*self).set(i,d);}   /* for python */
+        inline dai::Real __paren__(int i) const {return (*self).get(i);}     /* for octave */
+        inline void __paren_asgn__(int i,dai::Real d) {(*self).set(i,d);}  /* for octave */
 };
 
 %template(Factor) dai::TFactor<dai::Real>;
index 4cfa78b..f47ef4d 100644 (file)
@@ -62,8 +62,7 @@ class TestDAI {
 
             if( name == "LDPC" ) {
                 // special case: simulating a Low Density Parity Check code
-                Prob zero(2,0.0);
-                zero[0] = 1.0;
+                Real zero[2] = {1.0, 0.0};
                 for( size_t i = 0; i < fg.nrVars(); i++ )
                     varMarginals.push_back( Factor(fg.var(i), zero) );
                 allMarginals = varMarginals;
index e038866..8dcf61a 100644 (file)
@@ -48,10 +48,10 @@ BOOST_AUTO_TEST_CASE( ConstructorsTest ) {
     BOOST_CHECK_EQUAL( x3[1], 1.0 );
     BOOST_CHECK_EQUAL( x3[2], 1.0 );
     BOOST_CHECK_EQUAL( x3[3], 1.0 );
-    x3[0] = 0.5;
-    x3[1] = 1.0;
-    x3[2] = 2.0;
-    x3[3] = 4.0;
+    x3.set( 0, 0.5 );
+    x3.set( 1, 1.0 );
+    x3.set( 2, 2.0 );
+    x3.set( 3, 4.0 );
 
     Prob x4( x3.begin(), x3.end() );
     BOOST_CHECK_EQUAL( x4.size(), 4 );
@@ -90,7 +90,7 @@ BOOST_AUTO_TEST_CASE( IteratorTest ) {
     Prob x( 5, 0.0 );
     size_t i;
     for( i = 0; i < x.size(); i++ )
-        x[i] = i;
+        x.set( i, i );
 
     i = 0;
     for( Prob::const_iterator cit = x.begin(); cit != x.end(); cit++, i++ )
@@ -121,7 +121,7 @@ BOOST_AUTO_TEST_CASE( IteratorTest ) {
 BOOST_AUTO_TEST_CASE( QueriesTest ) {
     Prob x( 5, 0.0 );
     for( size_t i = 0; i < x.size(); i++ )
-        x[i] = 2.0 - i;
+        x.set( i, 2.0 - i );
 
     // test accumulate, min, max, sum, sumAbs, maxAbs
     BOOST_CHECK_EQUAL( x.sum(), 0.0 );
@@ -146,7 +146,7 @@ BOOST_AUTO_TEST_CASE( QueriesTest ) {
     BOOST_CHECK_EQUAL( x.accumulate( -1.0, fo_max<Real>(), fo_abs<Real>() ), 2.0 );
     BOOST_CHECK_EQUAL( x.accumulate( 3.0, fo_max<Real>(), fo_abs<Real>() ), 3.0 );
     BOOST_CHECK_EQUAL( x.accumulate( -3.0, fo_max<Real>(), fo_abs<Real>() ), 3.0 );
-    x[1] = 1.0;
+    x.set( 1, 1.0 );
     BOOST_CHECK_EQUAL( x.maxAbs(), 2.0 );
     BOOST_CHECK_EQUAL( x.accumulate( 0.0, fo_max<Real>(), fo_abs<Real>() ), 2.0 );
     BOOST_CHECK_EQUAL( x.accumulate( 1.0, fo_max<Real>(), fo_abs<Real>() ), 2.0 );
@@ -154,7 +154,7 @@ BOOST_AUTO_TEST_CASE( QueriesTest ) {
     BOOST_CHECK_EQUAL( x.accumulate( 3.0, fo_max<Real>(), fo_abs<Real>() ), 3.0 );
     BOOST_CHECK_EQUAL( x.accumulate( -3.0, fo_max<Real>(), fo_abs<Real>() ), 3.0 );
     for( size_t i = 0; i < x.size(); i++ )
-        x[i] = i ? (1.0 / i) : 0.0;
+        x.set( i, i ? (1.0 / i) : 0.0 );
     BOOST_CHECK_EQUAL( x.accumulate( 0.0, std::plus<Real>(), fo_inv0<Real>() ), 10.0 );
     x /= x.sum();
 
@@ -170,50 +170,50 @@ BOOST_AUTO_TEST_CASE( QueriesTest ) {
     BOOST_CHECK( !Prob( 3, 0.0 ).hasNegatives() );
     BOOST_CHECK( !Prob( 3, 1.0 ).hasNegatives() );
     BOOST_CHECK( Prob( 3, -1.0 ).hasNegatives() );
-    x[0] = 0.0; x[1] = 0.0; x[2] = -1.0; x[3] = 1.0; x[4] = 100.0;
+    x.set( 0, 0.0 ); x.set( 1, 0.0 ); x.set( 2, -1.0 ); x.set( 3, 1.0 ); x.set( 4, 100.0 );
     BOOST_CHECK( x.hasNegatives() );
-    x[2] = -INFINITY;
+    x.set( 2, -INFINITY );
     BOOST_CHECK( x.hasNegatives() );
-    x[2] = INFINITY;
+    x.set( 2, INFINITY );
     BOOST_CHECK( !x.hasNegatives() );
-    x[2] = -1.0;
+    x.set( 2, -1.0 );
 
     // test argmax
     BOOST_CHECK( x.argmax() == std::make_pair( (size_t)4, (Real)100.0 ) );
-    x[4] = 0.5;
+    x.set( 4, 0.5 );
     BOOST_CHECK( x.argmax() == std::make_pair( (size_t)3, (Real)1.0 ) );
-    x[3] = -2.0;
+    x.set( 3, -2.0 );
     BOOST_CHECK( x.argmax() == std::make_pair( (size_t)4, (Real)0.5 ) );
-    x[4] = -1.0;
+    x.set( 4, -1.0 );
     BOOST_CHECK( x.argmax() == std::make_pair( (size_t)0, (Real)0.0 ) );
-    x[0] = -2.0;
+    x.set( 0, -2.0 );
     BOOST_CHECK( x.argmax() == std::make_pair( (size_t)1, (Real)0.0 ) );
-    x[1] = -3.0;
+    x.set( 1, -3.0 );
     BOOST_CHECK( x.argmax() == std::make_pair( (size_t)2, (Real)-1.0 ) );
-    x[2] = -2.0;
+    x.set( 2, -2.0 );
     BOOST_CHECK( x.argmax() == std::make_pair( (size_t)4, (Real)-1.0 ) );
 
     // test draw
     for( size_t i = 0; i < x.size(); i++ )
-        x[i] = i ? (1.0 / i) : 0.0;
+        x.set( i, i ? (1.0 / i) : 0.0 );
     for( size_t repeat = 0; repeat < 10000; repeat++ ) {
         BOOST_CHECK( x.draw() < x.size() );
         BOOST_CHECK( x.draw() != 0 );
     }
-    x[2] = 0.0;
+    x.set( 2, 0.0 );
     for( size_t repeat = 0; repeat < 10000; repeat++ ) {
         BOOST_CHECK( x.draw() < x.size() );
         BOOST_CHECK( x.draw() != 0 );
         BOOST_CHECK( x.draw() != 2 );
     }
-    x[4] = 0.0;
+    x.set( 4, 0.0 );
     for( size_t repeat = 0; repeat < 10000; repeat++ ) {
         BOOST_CHECK( x.draw() < x.size() );
         BOOST_CHECK( x.draw() != 0 );
         BOOST_CHECK( x.draw() != 2 );
         BOOST_CHECK( x.draw() != 4 );
     }
-    x[1] = 0.0;
+    x.set( 1, 0.0 );
     for( size_t repeat = 0; repeat < 10000; repeat++ )
         BOOST_CHECK( x.draw() == 3 );
 
@@ -222,27 +222,27 @@ BOOST_AUTO_TEST_CASE( QueriesTest ) {
     BOOST_CHECK( !(a < b) );
     BOOST_CHECK( !(b < a) );
     BOOST_CHECK( a == b );
-    a[0] = 0.0;
+    a.set( 0, 0.0 );
     BOOST_CHECK( a < b );
     BOOST_CHECK( !(b < a) );
     BOOST_CHECK( !(a == b) );
-    b[2] = 0.0;
+    b.set( 2, 0.0 );
     BOOST_CHECK( a < b );
     BOOST_CHECK( !(b < a) );
     BOOST_CHECK( !(a == b) );
-    b[0] = 0.0;
+    b.set( 0, 0.0 );
     BOOST_CHECK( !(a < b) );
     BOOST_CHECK( b < a );
     BOOST_CHECK( !(a == b) );
-    a[1] = 0.0;
+    a.set( 1, 0.0 );
     BOOST_CHECK( a < b );
     BOOST_CHECK( !(b < a) );
     BOOST_CHECK( !(a == b) );
-    b[1] = 0.0;
+    b.set( 1, 0.0 );
     BOOST_CHECK( !(a < b) );
     BOOST_CHECK( b < a );
     BOOST_CHECK( !(a == b) );
-    a[2] = 0.0;
+    a.set( 2, 0.0 );
     BOOST_CHECK( !(a < b) );
     BOOST_CHECK( !(b < a) );
     BOOST_CHECK( a == b );
@@ -251,9 +251,9 @@ BOOST_AUTO_TEST_CASE( QueriesTest ) {
 
 BOOST_AUTO_TEST_CASE( UnaryTransformationsTest ) {
     Prob x( 3 );
-    x[0] = -2.0;
-    x[1] = 0.0;
-    x[2] = 2.0;
+    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>() );
@@ -282,7 +282,8 @@ BOOST_AUTO_TEST_CASE( UnaryTransformationsTest ) {
     BOOST_CHECK_EQUAL( y[1], -INFINITY );
     BOOST_CHECK_CLOSE( y[2], std::log(2.0), tol );
     BOOST_CHECK( !(y == z) );
-    y[0] = z[0] = 0.0;
+    y.set( 0, 0.0 );
+    z.set( 0, 0.0 );
     BOOST_CHECK( y == z );
 
     y = x.log(true);
@@ -291,7 +292,8 @@ BOOST_AUTO_TEST_CASE( UnaryTransformationsTest ) {
     BOOST_CHECK_EQUAL( y[1], 0.0 );
     BOOST_CHECK_EQUAL( y[2], std::log(2.0) );
     BOOST_CHECK( !(y == z) );
-    y[0] = z[0] = 0.0;
+    y.set( 0, 0.0 );
+    z.set( 0, 0.0 );
     BOOST_CHECK( y == z );
 
     y = x.inverse(false);
@@ -308,7 +310,7 @@ BOOST_AUTO_TEST_CASE( UnaryTransformationsTest ) {
     BOOST_CHECK_EQUAL( y[2], 0.5 );
     BOOST_CHECK( y == z );
 
-    x[0] = 2.0;
+    x.set( 0, 2.0 );
     y = x.normalized();
     BOOST_CHECK_EQUAL( y[0], 0.5 );
     BOOST_CHECK_EQUAL( y[1], 0.0 );
@@ -319,7 +321,7 @@ BOOST_AUTO_TEST_CASE( UnaryTransformationsTest ) {
     BOOST_CHECK_EQUAL( y[1], 0.0 );
     BOOST_CHECK_EQUAL( y[2], 0.5 );
 
-    x[0] = -2.0;
+    x.set( 0, -2.0 );
     y = x.normalized( Prob::NORMLINF );
     BOOST_CHECK_EQUAL( y[0], -1.0 );
     BOOST_CHECK_EQUAL( y[1], 0.0 );
@@ -329,18 +331,18 @@ BOOST_AUTO_TEST_CASE( UnaryTransformationsTest ) {
 
 BOOST_AUTO_TEST_CASE( UnaryOperationsTest ) {
     Prob xorg(3);
-    xorg[0] = 2.0;
-    xorg[1] = 0.0;
-    xorg[2] = 1.0;
+    xorg.set( 0, 2.0 );
+    xorg.set( 1, 0.0 );
+    xorg.set( 2, 1.0 );
     Prob y(3);
 
     Prob x = xorg;
     BOOST_CHECK( x.setUniform() == Prob(3) );
     BOOST_CHECK( x == Prob(3) );
 
-    y[0] = std::exp(2.0);
-    y[1] = 1.0;
-    y[2] = std::exp(1.0);
+    y.set( 0, std::exp(2.0) );
+    y.set( 1, 1.0 );
+    y.set( 2, std::exp(1.0) );
     x = xorg;
     BOOST_CHECK( x.takeExp() == y );
     BOOST_CHECK( x == y );
@@ -348,9 +350,9 @@ BOOST_AUTO_TEST_CASE( UnaryOperationsTest ) {
     BOOST_CHECK( x.pwUnaryOp( fo_exp<Real>() ) == y );
     BOOST_CHECK( x == y );
 
-    y[0] = std::log(2.0);
-    y[1] = -INFINITY;
-    y[2] = 0.0;
+    y.set( 0, std::log(2.0) );
+    y.set( 1, -INFINITY );
+    y.set( 2, 0.0 );
     x = xorg;
     BOOST_CHECK( x.takeLog() == y );
     BOOST_CHECK( x == y );
@@ -361,7 +363,7 @@ BOOST_AUTO_TEST_CASE( UnaryOperationsTest ) {
     BOOST_CHECK( x.pwUnaryOp( fo_log<Real>() ) == y );
     BOOST_CHECK( x == y );
 
-    y[1] = 0.0;
+    y.set( 1, 0.0 );
     x = xorg;
     BOOST_CHECK( x.takeLog(true) == y );
     BOOST_CHECK( x == y );
@@ -369,9 +371,9 @@ BOOST_AUTO_TEST_CASE( UnaryOperationsTest ) {
     BOOST_CHECK( x.pwUnaryOp( fo_log0<Real>() ) == y );
     BOOST_CHECK( x == y );
 
-    y[0] = 2.0 / 3.0;
-    y[1] = 0.0 / 3.0;
-    y[2] = 1.0 / 3.0;
+    y.set( 0, 2.0 / 3.0 );
+    y.set( 1, 0.0 / 3.0 );
+    y.set( 2, 1.0 / 3.0 );
     x = xorg;
     BOOST_CHECK_EQUAL( x.normalize(), 3.0 );
     BOOST_CHECK( x == y );
@@ -380,17 +382,17 @@ BOOST_AUTO_TEST_CASE( UnaryOperationsTest ) {
     BOOST_CHECK_EQUAL( x.normalize( Prob::NORMPROB ), 3.0 );
     BOOST_CHECK( x == y );
 
-    y[0] = 2.0 / 2.0;
-    y[1] = 0.0 / 2.0;
-    y[2] = 1.0 / 2.0;
+    y.set( 0, 2.0 / 2.0 );
+    y.set( 1, 0.0 / 2.0 );
+    y.set( 2, 1.0 / 2.0 );
     x = xorg;
     BOOST_CHECK_EQUAL( x.normalize( Prob::NORMLINF ), 2.0 );
     BOOST_CHECK( x == y );
 
-    xorg[0] = -2.0;
-    y[0] = 2.0;
-    y[1] = 0.0;
-    y[2] = 1.0;
+    xorg.set( 0, -2.0 );
+    y.set( 0, 2.0 );
+    y.set( 1, 0.0 );
+    y.set( 2, 1.0 );
     x = xorg;
     BOOST_CHECK( x.takeAbs() == y );
     BOOST_CHECK( x == y );
@@ -407,39 +409,39 @@ BOOST_AUTO_TEST_CASE( UnaryOperationsTest ) {
 
 BOOST_AUTO_TEST_CASE( ScalarTransformationsTest ) {
     Prob x(3);
-    x[0] = 2.0;
-    x[1] = 0.0;
-    x[2] = 1.0;
+    x.set( 0, 2.0 );
+    x.set( 1, 0.0 );
+    x.set( 2, 1.0 );
     Prob y(3);
 
-    y[0] = 3.0; y[1] = 1.0; y[2] = 2.0;
+    y.set( 0, 3.0 ); y.set( 1, 1.0 ); y.set( 2, 2.0 );
     BOOST_CHECK( (x + 1.0) == y );
-    y[0] = 0.0; y[1] = -2.0; y[2] = -1.0;
+    y.set( 0, 0.0 ); y.set( 1, -2.0 ); y.set( 2, -1.0 );
     BOOST_CHECK( (x + (-2.0)) == y );
 
-    y[0] = 1.0; y[1] = -1.0; y[2] = 0.0;
+    y.set( 0, 1.0 ); y.set( 1, -1.0 ); y.set( 2, 0.0 );
     BOOST_CHECK( (x - 1.0) == y );
-    y[0] = 4.0; y[1] = 2.0; y[2] = 3.0;
+    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[0] = 4.0; y[1] = 0.0; y[2] = 2.0;
+    y.set( 0, 4.0 ); y.set( 1, 0.0 ); y.set( 2, 2.0 );
     BOOST_CHECK( (x * 2.0) == y );
-    y[0] = -1.0; y[1] = 0.0; y[2] = -0.5;
+    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[0] = 1.0; y[1] = 0.0; y[2] = 0.5;
+    y.set( 0, 1.0 ); y.set( 1, 0.0 ); y.set( 2, 0.5 );
     BOOST_CHECK( (x / 2.0) == y );
-    y[0] = -4.0; y[1] = 0.0; y[2] = -2.0;
+    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[0] = 4.0; y[1] = 0.0; y[2] = 1.0;
+    y.set( 0, 4.0 ); y.set( 1, 0.0 ); y.set( 2, 1.0 );
     BOOST_CHECK( (x ^ 2.0) == y );
-    y[0] = 1.0 / std::sqrt(2.0); y[1] = INFINITY; y[2] = 1.0;
+    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] );
@@ -449,9 +451,9 @@ BOOST_AUTO_TEST_CASE( ScalarTransformationsTest ) {
 
 BOOST_AUTO_TEST_CASE( ScalarOperationsTest ) {
     Prob xorg(3), x(3);
-    xorg[0] = 2.0;
-    xorg[1] = 0.0;
-    xorg[2] = 1.0;
+    xorg.set( 0, 2.0 );
+    xorg.set( 1, 0.0 );
+    xorg.set( 2, 1.0 );
     Prob y(3);
 
     x = xorg;
@@ -463,38 +465,38 @@ BOOST_AUTO_TEST_CASE( ScalarOperationsTest ) {
     BOOST_CHECK( x == Prob(3, 0.0) );
 
     x = xorg;
-    y[0] = 3.0; y[1] = 1.0; y[2] = 2.0;
+    y.set( 0, 3.0 ); y.set( 1, 1.0 ); y.set( 2, 2.0 );
     BOOST_CHECK( (x += 1.0) == y );
     BOOST_CHECK( x == y );
-    y[0] = 1.0; y[1] = -1.0; y[2] = 0.0;
+    y.set( 0, 1.0 ); y.set( 1, -1.0 ); y.set( 2, 0.0 );
     BOOST_CHECK( (x += -2.0) == y );
     BOOST_CHECK( x == y );
 
     x = xorg;
-    y[0] = 1.0; y[1] = -1.0; y[2] = 0.0;
+    y.set( 0, 1.0 ); y.set( 1, -1.0 ); y.set( 2, 0.0 );
     BOOST_CHECK( (x -= 1.0) == y );
     BOOST_CHECK( x == y );
-    y[0] = 3.0; y[1] = 1.0; y[2] = 2.0;
+    y.set( 0, 3.0 ); y.set( 1, 1.0 ); y.set( 2, 2.0 );
     BOOST_CHECK( (x -= -2.0) == y );
     BOOST_CHECK( x == y );
 
     x = xorg;
     BOOST_CHECK( (x *= 1.0) == x );
     BOOST_CHECK( x == x );
-    y[0] = 4.0; y[1] = 0.0; y[2] = 2.0;
+    y.set( 0, 4.0 ); y.set( 1, 0.0 ); y.set( 2, 2.0 );
     BOOST_CHECK( (x *= 2.0) == y );
     BOOST_CHECK( x == y );
-    y[0] = -1.0; y[1] = 0.0; y[2] = -0.5;
+    y.set( 0, -1.0 ); y.set( 1, 0.0 ); y.set( 2, -0.5 );
     BOOST_CHECK( (x *= -0.25) == y );
     BOOST_CHECK( x == y );
 
     x = xorg;
     BOOST_CHECK( (x /= 1.0) == x );
     BOOST_CHECK( x == x );
-    y[0] = 1.0; y[1] = 0.0; y[2] = 0.5;
+    y.set( 0, 1.0 ); y.set( 1, 0.0 ); y.set( 2, 0.5 );
     BOOST_CHECK( (x /= 2.0) == y );
     BOOST_CHECK( x == y );
-    y[0] = -4.0; y[1] = 0.0; y[2] = -2.0;
+    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) );
@@ -506,10 +508,10 @@ BOOST_AUTO_TEST_CASE( ScalarOperationsTest ) {
     BOOST_CHECK( (x ^= 0.0) == Prob(3, 1.0) );
     BOOST_CHECK( x == Prob(3, 1.0) );
     x = xorg;
-    y[0] = 4.0; y[1] = 0.0; y[2] = 1.0;
+    y.set( 0, 4.0 ); y.set( 1, 0.0 ); y.set( 2, 1.0 );
     BOOST_CHECK( (x ^= 2.0) == y );
     BOOST_CHECK( x == y );
-    y[0] = 0.5; y[1] = INFINITY; y[2] = 1.0;
+    y.set( 0, 0.5 ); y.set( 1, INFINITY ); y.set( 2, 1.0 );
     BOOST_CHECK( (x ^= -0.5) == y );
     BOOST_CHECK( x == y );
 }
@@ -518,12 +520,12 @@ BOOST_AUTO_TEST_CASE( ScalarOperationsTest ) {
 BOOST_AUTO_TEST_CASE( VectorOperationsTest ) {
     size_t N = 6;
     Prob xorg(N), x(N);
-    xorg[0] = 2.0; xorg[1] = 0.0; xorg[2] = 1.0; xorg[3] = 0.0; xorg[4] = 2.0; xorg[5] = 3.0;
+    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);
-    y[0] = 0.5; y[1] = -1.0; y[2] = 0.0; y[3] = 0.0; y[4] = -2.0; y[5] = 3.0;
+    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);
 
-    z[0] = 2.5; z[1] = -1.0; z[2] = 1.0; z[3] = 0.0; z[4] = 0.0; z[5] = 6.0;
+    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;
     r = (x += y);
     for( size_t i = 0; i < N; i++ )
@@ -533,7 +535,7 @@ BOOST_AUTO_TEST_CASE( VectorOperationsTest ) {
     BOOST_CHECK( x.pwBinaryOp( y, std::plus<Real>() ) == z );
     BOOST_CHECK( x == z );
 
-    z[0] = 1.5; z[1] = 1.0; z[2] = 1.0; z[3] = 0.0; z[4] = 4.0; z[5] = 0.0;
+    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 );
     x = xorg;
     r = (x -= y);
     for( size_t i = 0; i < N; i++ )
@@ -543,7 +545,7 @@ BOOST_AUTO_TEST_CASE( VectorOperationsTest ) {
     BOOST_CHECK( x.pwBinaryOp( y, std::minus<Real>() ) == z );
     BOOST_CHECK( x == z );
 
-    z[0] = 1.0; z[1] = 0.0; z[2] = 0.0; z[3] = 0.0; z[4] = -4.0; z[5] = 9.0;
+    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 );
     x = xorg;
     r = (x *= y);
     for( size_t i = 0; i < N; i++ )
@@ -553,7 +555,7 @@ BOOST_AUTO_TEST_CASE( VectorOperationsTest ) {
     BOOST_CHECK( x.pwBinaryOp( y, std::multiplies<Real>() ) == z );
     BOOST_CHECK( x == z );
 
-    z[0] = 4.0; z[1] = 0.0; z[2] = 0.0; z[3] = 0.0; z[4] = -1.0; z[5] = 1.0;
+    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 );
     x = xorg;
     r = (x /= y);
     for( size_t i = 0; i < N; i++ )
@@ -563,7 +565,7 @@ BOOST_AUTO_TEST_CASE( VectorOperationsTest ) {
     BOOST_CHECK( x.pwBinaryOp( y, fo_divides0<Real>() ) == z );
     BOOST_CHECK( x == z );
 
-    z[0] = 4.0; z[1] = 0.0; z[2] = INFINITY; /*z[3] = INFINITY;*/ z[4] = -1.0; z[5] = 1.0;
+    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 );
@@ -572,7 +574,7 @@ BOOST_AUTO_TEST_CASE( VectorOperationsTest ) {
     BOOST_CHECK( isnan(r[3]) );
     BOOST_CHECK_CLOSE( r[4], z[4], tol );
     BOOST_CHECK_CLOSE( r[5], z[5], tol );
-    x[3] = 0.0; r[3] = 0.0;
+    x.set( 3, 0.0 ); r.set( 3, 0.0 );
     BOOST_CHECK( x == r );
     x = xorg;
     r = x.pwBinaryOp( y, std::divides<Real>() );
@@ -582,10 +584,10 @@ BOOST_AUTO_TEST_CASE( VectorOperationsTest ) {
     BOOST_CHECK( isnan(r[3]) );
     BOOST_CHECK_CLOSE( r[4], z[4], tol );
     BOOST_CHECK_CLOSE( r[5], z[5], tol );
-    x[3] = 0.0; r[3] = 0.0;
+    x.set( 3, 0.0 ); r.set( 3, 0.0 );
     BOOST_CHECK( x == r );
 
-    z[0] = std::sqrt(2.0); z[1] = INFINITY; z[2] = 1.0; z[3] = 1.0; z[4] = 0.25; z[5] = 27.0;
+    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 );
@@ -604,40 +606,40 @@ BOOST_AUTO_TEST_CASE( VectorOperationsTest ) {
 BOOST_AUTO_TEST_CASE( VectorTransformationsTest ) {
     size_t N = 6;
     Prob x(N);
-    x[0] = 2.0; x[1] = 0.0; x[2] = 1.0; x[3] = 0.0; x[4] = 2.0; x[5] = 3.0;
+    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);
-    y[0] = 0.5; y[1] = -1.0; y[2] = 0.0; y[3] = 0.0; y[4] = -2.0; y[5] = 3.0;
+    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);
 
-    z[0] = 2.5; z[1] = -1.0; z[2] = 1.0; z[3] = 0.0; z[4] = 0.0; z[5] = 6.0;
+    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>() );
     BOOST_CHECK( r == z );
 
-    z[0] = 1.5; z[1] = 1.0; z[2] = 1.0; z[3] = 0.0; z[4] = 4.0; z[5] = 0.0;
+    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>() );
     BOOST_CHECK( r == z );
 
-    z[0] = 1.0; z[1] = 0.0; z[2] = 0.0; z[3] = 0.0; z[4] = -4.0; z[5] = 9.0;
+    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>() );
     BOOST_CHECK( r == z );
 
-    z[0] = 4.0; z[1] = 0.0; z[2] = 0.0; z[3] = 0.0; z[4] = -1.0; z[5] = 1.0;
+    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>() );
     BOOST_CHECK( r == z );
 
-    z[0] = 4.0; z[1] = 0.0; z[2] = INFINITY; /*z[3] = INFINITY;*/ z[4] = -1.0; z[5] = 1.0;
+    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 );
@@ -653,7 +655,7 @@ BOOST_AUTO_TEST_CASE( VectorTransformationsTest ) {
     BOOST_CHECK_CLOSE( r[4], z[4], tol );
     BOOST_CHECK_CLOSE( r[5], z[5], tol );
 
-    z[0] = std::sqrt(2.0); z[1] = INFINITY; z[2] = 1.0; z[3] = 1.0; z[4] = 0.25; z[5] = 27.0;
+    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] );
@@ -668,12 +670,12 @@ BOOST_AUTO_TEST_CASE( VectorTransformationsTest ) {
 
 BOOST_AUTO_TEST_CASE( RelatedFunctionsTest ) {
     Prob x(3), y(3), z(3);
-    x[0] = 0.2;
-    x[1] = 0.8;
-    x[2] = 0.0;
-    y[0] = 0.0;
-    y[1] = 0.6;
-    y[2] = 0.4;
+    x.set( 0, 0.2 );
+    x.set( 1, 0.8 );
+    x.set( 2, 0.0 );
+    y.set( 0, 0.0 );
+    y.set( 1, 0.6 );
+    y.set( 2, 0.4 );
 
     z = min( x, y );
     BOOST_CHECK_EQUAL( z[0], 0.0 );
@@ -714,8 +716,8 @@ BOOST_AUTO_TEST_CASE( RelatedFunctionsTest ) {
     BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTHEL ), 0.5 * (0.2 + std::pow(std::sqrt(0.8) - std::sqrt(0.6), 2.0) + 0.4) );
     BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTHEL ), x.innerProduct( y, 0.0, std::plus<Real>(), fo_Hellinger<Real>() ) / 2.0 );
     BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTHEL ), y.innerProduct( x, 0.0, std::plus<Real>(), fo_Hellinger<Real>() ) / 2.0 );
-    x[1] = 0.7; x[2] = 0.1;
-    y[0] = 0.1; y[1] = 0.5;
+    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>() ) );
index 848eac5..7d6180e 100644 (file)
@@ -97,7 +97,7 @@ int main( int argc, char *argv[] ) {
         cout << "Pairwise interactions? " << fg.isPairwise() << endl;
         // Calculate treewidth using various heuristics, if requested
         if( calc_tw ) {
-            std::pair<size_t,size_t> tw;
+            std::pair<size_t,double> tw;
             tw = boundTreewidth(fg, &eliminationCost_MinNeighbors);
             cout << "Treewidth (MinNeighbors):     " << tw.first << " (" << tw.second << " states)" << endl;
             tw = boundTreewidth(fg, &eliminationCost_MinWeight);