Merge branch 'eaton'
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 4 Aug 2009 15:35:49 +0000 (17:35 +0200)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 4 Aug 2009 15:35:49 +0000 (17:35 +0200)
Conflicts:

ChangeLog
Makefile
include/dai/daialg.h

51 files changed:
ChangeLog
Makefile
Makefile.CYGWIN
Makefile.LINUX
Makefile.MACOSX
Makefile.WINDOWS
README
include/dai/alldai.h
include/dai/bipgraph.h
include/dai/bp.h
include/dai/daialg.h
include/dai/doc.h
include/dai/emalg.h [new file with mode: 0644]
include/dai/evidence.h [new file with mode: 0644]
include/dai/exactinf.h
include/dai/exceptions.h
include/dai/factor.h
include/dai/factorgraph.h
include/dai/gibbs.h
include/dai/hak.h
include/dai/jtree.h
include/dai/lc.h
include/dai/mf.h
include/dai/mr.h
include/dai/prob.h
include/dai/regiongraph.h
include/dai/treeep.h
include/dai/util.h
include/dai/weightedgraph.h
src/bp.cpp
src/emalg.cpp [new file with mode: 0644]
src/evidence.cpp [new file with mode: 0644]
src/exceptions.cpp
src/factorgraph.cpp
src/matlab/dai.cpp
src/util.cpp
tests/testdai.cpp
tests/testem/2var.em [new file with mode: 0644]
tests/testem/2var.fg [new file with mode: 0644]
tests/testem/2var_data.tab [new file with mode: 0644]
tests/testem/3var.em [new file with mode: 0644]
tests/testem/3var.fg [new file with mode: 0644]
tests/testem/hoi1_data.tab [new file with mode: 0644]
tests/testem/hoi1_infer_f2.em [new file with mode: 0644]
tests/testem/hoi1_share_f0_f1_f2.em [new file with mode: 0644]
tests/testem/hoi1_share_f0_f2.em [new file with mode: 0644]
tests/testem/runtests [new file with mode: 0755]
tests/testem/runtests.bat [new file with mode: 0755]
tests/testem/testem.cpp [new file with mode: 0644]
tests/testem/testem.out [new file with mode: 0644]
tests/testregression

index 1dc0efe..45e54c5 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,13 @@
-git b6efbfd9e4e76de65c44e4060a9b56f99543b70a (head of branch eaton before merge)
---------------------------------------------------------------------------------
+* Added work-around for bug in Boost Graph Library
+* Improvements of TFactor<T>:
+  - Extended functionality of TFactor<T>::operators +,-,+=,-= to handle different VarSets
+  - Added TFactor<T>::maxMarginal (similar to marginal but with max instead of sum)
+* Added BipartiteGraph::eraseEdge
+* Removed all the virtual default constructors *::create(), as they are not used.
+* Fixed bug in MaxSpanningTree (it wrongly assumed that the tree was not empty).
+* [Charlie Vaske] Added Expectation Maximization code.
+* Added MatLab QuickStart to README.
+* MEX file dai now also returns variable and factor beliefs.
 * Cleanups and extra documentation of the code contributed by Frederik Eaton
 * Removed erroneous 'inline' directives in gibbs.cpp
 * [Frederik Eaton] Added script for automatically generating properties code (used by CBP and BBP)
@@ -57,7 +65,7 @@ git b6efbfd9e4e76de65c44e4060a9b56f99543b70a (head of branch eaton before merge)
   errors (thanks to Dan Preston for reporting this)
 * toc() now returns POSIX system time with maximum accuracy of microseconds
 * Exception objects now remember their error code
-* Added examples/example_springler.cpp
+* Added examples/example_sprinkler.cpp
 
 
 git 065eae35cbfcc36f1a945ae3053c80c23f366306
index d47ce9e..749847b 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -31,7 +31,7 @@ SRC=src
 LIB=lib
 
 # Define build targets
-TARGETS=tests utils lib examples testregression
+TARGETS=tests utils lib examples testregression testem
 ifdef WITH_DOC
   TARGETS:=$(TARGETS) doc 
 endif
@@ -45,7 +45,7 @@ ifdef WITH_MATLAB
 endif
 
 # Define conditional build targets
-OBJECTS:=exactinf$(OE)
+OBJECTS:=exactinf$(OE) evidence$(OE) emalg$(OE)
 ifdef WITH_BP
   CCFLAGS:=$(CCFLAGS) -DDAI_WITH_BP
   OBJECTS:=$(OBJECTS) bp$(OE)
@@ -109,7 +109,7 @@ examples : examples/example$(EE) examples/example_bipgraph$(EE) examples/example
 
 matlabs : matlab/dai$(ME) matlab/dai_readfg$(ME) matlab/dai_writefg$(ME) matlab/dai_potstrength$(ME)
 
-tests : tests/testdai$(EE) tests/testbbp$(EE)
+tests : tests/testdai$(EE) tests/testem/testem$(EE) tests/testbbp$(EE)
 
 utils : utils/createfg$(EE) utils/fg2dot$(EE) utils/fginfo$(EE)
 
@@ -176,6 +176,12 @@ mr$(OE) : $(SRC)/mr.cpp $(INC)/mr.h $(HEADERS)
 gibbs$(OE) : $(SRC)/gibbs.cpp $(INC)/gibbs.h $(HEADERS)
        $(CC) -c $(SRC)/gibbs.cpp
 
+evidence$(OE) : $(SRC)/evidence.cpp $(INC)/evidence.h $(HEADERS)
+       $(CC) -c $(SRC)/evidence.cpp
+
+emalg$(OE) : $(SRC)/emalg.cpp $(INC)/emalg.h $(INC)/evidence.h $(HEADERS)
+       $(CC) -c $(SRC)/emalg.cpp
+
 properties$(OE) : $(SRC)/properties.cpp $(HEADERS)
        $(CC) -c $(SRC)/properties.cpp
 
@@ -207,6 +213,8 @@ examples/example_sprinkler$(EE) : examples/example_sprinkler.cpp $(HEADERS) $(LI
 
 tests/testdai$(EE) : tests/testdai.cpp $(HEADERS) $(LIB)/libdai$(LE)
        $(CC) $(CCO)tests/testdai$(EE) tests/testdai.cpp $(LIBS) $(BOOSTLIBS)
+tests/testem/testem$(EE) : tests/testem/testem.cpp $(HEADERS) $(LIB)/libdai$(LE)
+       $(CC) $(CCO)$@ $< $(LIBS) $(BOOSTLIBS)
 
 tests/testbbp$(EE) : tests/testbbp.cpp $(HEADERS) $(LIB)/libdai$(LE)
        $(CC) $(CCO)tests/testbbp$(EE) tests/testbbp.cpp $(LIBS)
@@ -261,16 +269,22 @@ endif
 # REGRESSION TESTS
 ###################
 
-ifneq ($(OS),WINDOWS)
 testregression : tests/testdai$(EE)
        @echo Starting regression test...this can take a minute or so!
+ifneq ($(OS),WINDOWS)
        cd tests && ./testregression && cd ..
 else
-testregression : tests/testdai$(EE)
-       @echo Starting regression test...this can take a minute or so!
        cd tests && testregression.bat && cd ..
 endif
 
+testem : tests/testem/testem$(EE)
+       @echo Starting EM tests
+ifneq ($(OS),WINDOWS)
+       cd tests/testem && ./runtests && cd ../..
+else
+       cd tests\testem && runtests && cd ..\..
+endif
+
 
 # DOCUMENTATION
 ################
@@ -292,12 +306,12 @@ clean :
        -rm *$(OE)
        -rm matlab/*$(ME)
        -rm examples/example$(EE) examples/example_bipgraph$(EE) examples/example_varset$(EE) examples/example_sprinkler$(EE)
-       -rm tests/testdai$(EE) tests/testbbp$(EE)
+       -rm tests/testdai$(EE) tests/testem/testem$(EE) tests/testbbp$(EE)
        -rm utils/fg2dot$(EE) utils/createfg$(EE) utils/fginfo$(EE)
        -rm -R doc
        -rm -R lib
 else
 .PHONY : clean
 clean :
-       -del *$(OE) *.ilk *.pdb *$(EE) matlab\*$(ME) examples\*$(EE) examples\*.ilk examples\*.pdb tests\testdai$(EE) tests\*.pdb tests\*.ilk utils\*$(EE) utils\*.pdb utils\*.ilk $(LIB)\libdai$(LE)
+       -del *$(OE) *.ilk *.pdb *$(EE) matlab\*$(ME) examples\*$(EE) examples\*.ilk examples\*.pdb tests\testdai$(EE) tests\testem\testem$(EE) tests\*.pdb tests\*.ilk utils\*$(EE) utils\*.pdb utils\*.ilk $(LIB)\libdai$(LE)
 endif
index 4c7676e..3a11de9 100644 (file)
@@ -22,7 +22,7 @@ WITH_DOC=true
 DEBUG=true
 # Build matlab interface?
 WITH_MATLAB=
-# New/old matlab version?
+# MatLab version 7.3 (R2006b) or newer?
 NEW_MATLAB=true
 
 # OPERATING SYSTEM
index ea2814e..2bc8b2f 100644 (file)
@@ -23,7 +23,7 @@ WITH_DOC=true
 DEBUG=true
 # Build matlab interface?
 WITH_MATLAB=
-# New/old matlab version?
+# MatLab version 7.3 (R2006b) or newer?
 NEW_MATLAB=true
 
 # OPERATING SYSTEM
index 1599dcd..b176114 100644 (file)
@@ -22,7 +22,7 @@ WITH_DOC=true
 DEBUG=true
 # Build matlab interface?
 WITH_MATLAB=
-# New/old matlab version?
+# MatLab version 7.3 (R2006b) or newer?
 NEW_MATLAB=true
 
 # OPERATING SYSTEM
index 2c5fdd5..9728988 100644 (file)
@@ -23,7 +23,7 @@ WITH_DOC=
 DEBUG=true
 # Build matlab interface?
 WITH_MATLAB=
-# New/old matlab version?
+# MatLab version 7.3 (R2006b) or newer?
 NEW_MATLAB=true
 
 # OPERATING SYSTEM
diff --git a/README b/README
index 1d85a74..06c144a 100644 (file)
--- a/README
+++ b/README
@@ -13,6 +13,7 @@ with contributions from:
 Martijn Leisink
 Giuseppe Passino
 Frederik Eaton
+Charlie Vaske
 Bastian Wemmenhove
 Christian Wojek
 Claudio Lima
@@ -47,7 +48,7 @@ to the literature, to allow replication; (b) mention this software in the
 Acknowledgements section.  The appropriate citation is: 
 
 J. M. Mooij (2008) "libDAI 0.2.2: A free/open source C++ library for Discrete 
-Approximate Inference methods", http://mloss.org/software/view/77/.
+Approximate Inference methods", http://www.libdai.org
 
 Moreover, as a personal note, I would appreciate it if you would email me with
 citations of papers referencing this work so I can mention them to my funding
@@ -90,17 +91,20 @@ Currently, libDAI supports the following (approximate) inference methods:
     * Gibbs sampler;
     * Conditioned BP [EaG09].
 
+In addition, libDAI supports parameter learning of conditional probability
+tables by Expectation Maximization.
+
 
 Why C++?
 --------
 Because libDAI is implemented in C++, it is very fast compared with
 implementations in MatLab (a factor 1000 faster is not uncommon). libDAI does
-provide a MatLab interface for easy integration with MatLab.
+provide a (limited) MatLab interface for easy integration with MatLab.
 
 
 Releases
 --------
-Releases can be obtained from http://mloss.org/software/view/77/
+Releases can be obtained from www.libdai.org
 License: GNU Public License v2 (or higher).
 
 libDAI-0.2      December 1, 2006
@@ -165,7 +169,7 @@ from http://www.boost.org/ and compile/install it with:
 To build the libDAI source, first copy a template Makefile.* to Makefile.conf
 (for example, copy Makefile.LINUX to Makefile.conf if you use GNU/Linux). 
 Then, edit the Makefile.conf template to adapt it to your local setup.
-Especially directories may change from system to system. Finally, run
+Especially directories may differ from system to system. Finally, run
     
     make
 
@@ -199,3 +203,42 @@ If the build was successful, you can test the example program:
 or the more elaborate test program:
 
     tests\testdai --aliases tests\aliases.conf --filename tests\alarm.fg --methods JTREE_HUGIN BP_SEQMAX
+
+
+Quick start (MatLab)
+--------------------
+You need:
+- MatLab
+- The platform-dependent requirements described above
+
+First, you need to build the libDAI source as described above for your
+platform. By default, the MatLab interface is disabled, so before compiling the
+source, you have to enable it in the Makefile.conf by setting
+"WITH_MATLAB=true". Also, you have to configure the MatLab-specific parts of
+Makefile.conf to match your system (in particular, the Makefile variables ME,
+MATLABDIR and MEX). The MEX file extension depends on your platform; for a
+64-bit linux x86_64 system this would be "ME=.mexa64", for a 32-bit linux x86
+system this would be "ME=.mexglx". If you are unsure about your MEX file
+extension: it needs to be the same as what the MatLab command "mexext" returns.
+The required MEX files are built by issuing
+
+    make
+
+from the command line. The MatLab interface is much less powerful than using
+libDAI from C++. There are two reasons for this: (i) it is boring to write MEX
+files; (ii) the large performance penalty paid when large data structures (like
+factor graphs) have to be converted between their native C++ data structure to
+something that MatLab understands.
+
+A simple example of how to use the MatLab interface is the following (entered
+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');
+    [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]')
+
+where "path_to_libdai" has to be replaced with the directory in which libDAI
+was installed. For other algorithms and their parameters, see
+tests/aliases.conf.
index 61e4547..862bf8d 100644 (file)
@@ -33,6 +33,8 @@
 #include <dai/daialg.h>
 #include <dai/properties.h>
 #include <dai/exactinf.h>
+#include <dai/evidence.h>
+#include <dai/emalg.h>
 #ifdef DAI_WITH_BP
     #include <dai/bp.h>
 #endif
index 2a9d59e..15845de 100644 (file)
@@ -313,6 +313,22 @@ class BipartiteGraph {
         /// Removes node n2 of type 2 and all incident edges.
         void erase2( size_t n2 );
 
+        /// Removes edge between node n1 of type 1 and node n2 of type 2.
+        void eraseEdge( size_t n1, size_t n2 ) {
+            assert( n1 < nr1() );
+            assert( n2 < nr2() );
+            for( Neighbors::iterator i1 = _nb1[n1].begin(); i1 != _nb1[n1].end(); i1++ )
+                if( i1->node == n2 ) {
+                    _nb1[n1].erase( i1 );
+                    break;
+                }
+            for( Neighbors::iterator i2 = _nb2[n2].begin(); i2 != _nb2[n2].end(); i2++ )
+                if( i2->node == n1 ) {
+                    _nb2[n2].erase( i2 );
+                    break;
+                }
+        }
+
         /// Adds an edge between node n1 of type 1 and node n2 of type 2.
         /** If check == true, only adds the edge if it does not exist already.
          */
index 2d0a1f8..9117f9b 100644 (file)
@@ -137,7 +137,6 @@ class BP : public DAIAlgFG {
         /// @name General InfAlg interface
         //@{
         virtual BP* clone() const { return new BP(*this); }
-        virtual BP* create() const { return new BP(); }
         virtual std::string identify() const;
         virtual Factor belief( const Var &n ) const;
         virtual Factor belief( const VarSet &ns ) const;
index 0257311..4dbae3a 100644 (file)
@@ -54,9 +54,6 @@ class InfAlg {
         /// Returns a pointer to a new, cloned copy of *this (i.e., virtual copy constructor)
         virtual InfAlg* clone() const = 0;
 
-        /// Returns a pointer to a newly constructed object *this (i.e., virtual default constructor)
-        virtual InfAlg* create() const { DAI_THROW(NOT_IMPLEMENTED); }
-        
         /// Identifies itself for logging purposes
         virtual std::string identify() const = 0;
 
index 53e7362..e7f9071 100644 (file)
  *  - Gibbs sampler;
  *  - Conditioned BP [\ref EaG09];
  *
+ *  In addition, libDAI supports parameter learning of conditional probability
+ *  tables by Expectation Maximization.
+ *
  *  \section language Why C++?
  *  Because libDAI is implemented in C++, it is very fast compared with
  *  implementations in MatLab (a factor 1000 faster is not uncommon).
  */
 
 
-/** \page biboliography Bibliography
+/** \page bibliography Bibliography
  *  \section Bibliograpy
  *  \anchor KFL01 \ref KFL01
  *  F. R. Kschischang and B. J. Frey and H.-A. Loeliger (2001):
diff --git a/include/dai/emalg.h b/include/dai/emalg.h
new file mode 100644 (file)
index 0000000..021cf7e
--- /dev/null
@@ -0,0 +1,318 @@
+/*  Copyright (C) 2009  Charles Vaske  [cvaske at soe dot ucsc dot edu]
+    University of California Santa Cruz
+
+    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
+*/
+
+
+#ifndef __defined_libdai_emalg_h
+#define __defined_libdai_emalg_h
+
+
+#include <vector>
+#include <map>
+
+#include <dai/factor.h>
+#include <dai/daialg.h>
+#include <dai/evidence.h>
+#include <dai/index.h>
+#include <dai/properties.h>
+
+
+/// \file
+/// \brief Defines classes related to Expectation Maximization: EMAlg, ParameterEstimation, CondProbEstimation and SharedParameters
+/// \todo Describe EM file format
+
+
+namespace dai {
+
+
+/// Interface for a parameter estimation method. 
+/** This parameter estimation interface is based on sufficient statistics. 
+ *  Implementations are responsible for collecting data from a probability 
+ *  vector passed to it from a SharedParameters container object.
+ *
+ *  Implementations of this interface should register a factory function
+ *  via the static ParameterEstimation::registerMethod function.
+ *  The default registry only contains CondProbEstimation, named
+ *  "ConditionalProbEstimation".
+ */
+class ParameterEstimation {
+    public:  
+        /// A pointer to a factory function.
+        typedef ParameterEstimation* (*ParamEstFactory)( const PropertySet& );
+
+        /// Virtual destructor for deleting pointers to derived classes.
+        virtual ~ParameterEstimation() {}
+        /// Virtual copy constructor.
+        virtual ParameterEstimation* clone() const = 0;
+
+        /// General factory method for construction of ParameterEstimation subclasses.
+        static ParameterEstimation* construct( const std::string &method, const PropertySet &p );
+
+        /// Register a subclass so that it can be used with ParameterEstimation::construct.
+        static void registerMethod( const std::string &method, const ParamEstFactory &f ) {
+            if( _registry == NULL )
+                loadDefaultRegistry();
+            (*_registry)[method] = f;
+        }
+
+        /// Estimate the factor using the accumulated sufficient statistics and reset.
+        virtual Prob estimate() = 0;
+
+        /// Accumulate the sufficient statistics for p.
+        virtual void addSufficientStatistics( const Prob &p ) = 0;
+
+        /// Returns the size of the Prob that should be passed to addSufficientStatistics.
+        virtual size_t probSize() const = 0;
+
+    private:
+        /// A static registry containing all methods registered so far.
+        static std::map<std::string, ParamEstFactory> *_registry;
+
+        /// Registers default ParameterEstimation subclasses (currently, only CondProbEstimation).
+        static void loadDefaultRegistry();
+};
+
+
+/// Estimates the parameters of a conditional probability table, using pseudocounts.
+class CondProbEstimation : private ParameterEstimation {
+    private:
+        /// Number of states of the variable of interest
+        size_t _target_dim;
+        /// Current pseudocounts
+        Prob _stats;
+        /// Initial pseudocounts
+        Prob _initial_stats;
+
+    public:
+        /// Constructor
+        /** For a conditional probability \f$ P( X | Y ) \f$, 
+         *  \param target_dimension should equal \f$ | X | \f$
+         *  \param pseudocounts has length \f$ |X| \cdot |Y| \f$
+         */
+        CondProbEstimation( size_t target_dimension, const Prob &pseudocounts );
+
+        /// Virtual constructor, using a PropertySet.
+        /** Some keys in the PropertySet are required.
+         *  For a conditional probability \f$ P( X | Y ) \f$, 
+         *     - target_dimension should be equal to \f$ | X | \f$
+         *     - total_dimension should be equal to \f$ |X| \cdot |Y| \f$
+         *  
+         *  An optional key is:
+         *     - pseudo_count, which specifies the initial counts (defaults to 1)
+         */
+        static ParameterEstimation* factory( const PropertySet &p );
+        
+        /// Virtual copy constructor
+        virtual ParameterEstimation* clone() const { return new CondProbEstimation( _target_dim, _initial_stats ); }
+
+        /// Virtual destructor
+        virtual ~CondProbEstimation() {}
+        
+        /// Returns an estimate of the conditional probability distribution.
+        /** The format of the resulting Prob keeps all the values for 
+         *  \f$ P(X | Y=y) \f$ in sequential order in the array.
+         */
+        virtual Prob estimate();
+        
+        /// Accumulate sufficient statistics from the expectations in p.
+        virtual void addSufficientStatistics( const Prob &p );
+        
+        /// Returns the required size for arguments to addSufficientStatistics.
+        virtual size_t probSize() const { return _stats.size(); }
+};
+
+
+/// A single factor or set of factors whose parameters should be estimated.
+/** To ensure that parameters can be shared between different factors during
+ *  EM learning, each factor's values are reordered to match a desired variable 
+ *  ordering. The ordering of the variables in a factor may therefore differ 
+ *  from the canonical ordering used in libDAI. The SharedParameters
+ *  class couples one or more factors (together with the specified orderings
+ *  of the variables) with a ParameterEstimation object, taking care of the
+ *  necessary permutations of the factor entries / parameters.
+ */
+class SharedParameters {
+    public:
+        /// Convenience label for an index into a factor in a FactorGraph.
+        typedef size_t FactorIndex;
+        /// Convenience label for a grouping of factor orientations.
+        typedef std::map<FactorIndex, std::vector<Var> > FactorOrientations;
+
+    private:
+        /// Maps factor indices to the corresponding VarSets
+        std::map<FactorIndex, VarSet> _varsets;
+        /// Maps factor indices to the corresponding Permute objects that permute the desired ordering into the canonical ordering
+        std::map<FactorIndex, Permute> _perms;
+        /// Maps factor indices to the corresponding desired variable orderings
+        FactorOrientations _varorders;
+        /// Parameter estimation method to be used
+        ParameterEstimation *_estimation;
+        /// Indicates whether the object pointed to by _estimation should be deleted upon destruction
+        bool _deleteEstimation;
+
+        /// Calculates the permutation that permutes the variables in varorder into the canonical ordering
+        static Permute calculatePermutation( const std::vector<Var> &varorder, VarSet &outVS );
+
+        /// Initializes _varsets and _perms from _varorders
+        void setPermsAndVarSetsFromVarOrders();
+
+    public:
+        /// Copy constructor
+        SharedParameters( const SharedParameters &sp );
+
+        /// Constructor 
+        /** \param varorders  all the factor orientations for this parameter
+         *  \param estimation a pointer to the parameter estimation method
+         */ 
+        SharedParameters( const FactorOrientations &varorders, ParameterEstimation *estimation );
+
+        /// Constructor for making an object from a stream and a factor graph
+        SharedParameters( std::istream &is, const FactorGraph &fg_varlookup );
+
+        /// Destructor
+        ~SharedParameters() { 
+            if( _deleteEstimation )
+                delete _estimation;
+        }
+
+        /// Collect the necessary statistics from expected values
+        void collectSufficientStatistics( InfAlg &alg );
+
+        /// Estimate and set the shared parameters
+        void setParameters( FactorGraph &fg );
+};
+
+
+/// A MaximizationStep groups together several parameter estimation tasks into a single unit.
+class MaximizationStep { 
+    private:
+        std::vector<SharedParameters> _params;
+
+    public:
+        /// Default constructor
+        MaximizationStep() : _params() {}
+
+        /// Constructor from a vector of SharedParameters objects
+        MaximizationStep( std::vector<SharedParameters> &maximizations ) : _params(maximizations) {}  
+
+        /// Constructor from an input stream and a corresponding factor graph
+        MaximizationStep( std::istream &is, const FactorGraph &fg_varlookup );
+
+        /// Collect the beliefs from this InfAlg as expectations for the next Maximization step.
+        void addExpectations( InfAlg &alg );
+
+        /// Using all of the currently added expectations, make new factors with maximized parameters and set them in the FactorGraph.
+        void maximize( FactorGraph &fg );
+};
+
+
+/// EMAlg performs Expectation Maximization to learn factor parameters.
+/** This requires specifying:
+ *     - Evidence (instances of observations from the graphical model),
+ *     - InfAlg for performing the E-step, which includes the factor graph,
+ *     - a vector of MaximizationSteps steps to be performed.
+ *
+ *  This implementation can perform incremental EM by using multiple 
+ *  MaximizationSteps.  An expectation step is performed between execution
+ *  of each MaximizationStep.  A call to iterate() will cycle through all
+ *  MaximizationSteps.
+ *
+ *  Having multiple and separate maximization steps allows for maximizing some
+ *  parameters, performing another E step, and then maximizing separate
+ *  parameters, which may result in faster convergence in some cases.
+ */  
+class EMAlg {
+    private:
+        /// All the data samples used during learning
+        const Evidence &_evidence;
+
+        /// How to do the expectation step
+        InfAlg &_estep;
+
+        /// The maximization steps to take
+        std::vector<MaximizationStep> _msteps;
+
+        /// Number of iterations done
+        size_t _iters;
+
+        /// History of likelihoods
+        std::vector<Real> _lastLogZ;
+
+        /// Maximum number of iterations
+        size_t _max_iters;
+
+        /// Convergence tolerance
+        Real _log_z_tol;
+
+    public:
+        /// Key for setting maximum iterations @see setTermConditions
+        static const std::string MAX_ITERS_KEY;
+        /// Default maximum iterations @see setTermConditions
+        static const size_t MAX_ITERS_DEFAULT;
+        /// Key for setting likelihood termination condition @see setTermConditions
+        static const std::string LOG_Z_TOL_KEY;
+        /// Default likelihood tolerance @see setTermConditions
+        static const Real LOG_Z_TOL_DEFAULT;
+
+        /// Construct an EMAlg from all these objects
+        EMAlg( const Evidence &evidence, InfAlg &estep, std::vector<MaximizationStep> &msteps, const PropertySet &termconditions ) 
+          : _evidence(evidence), _estep(estep), _msteps(msteps), _iters(0), _lastLogZ(), _max_iters(MAX_ITERS_DEFAULT), _log_z_tol(LOG_Z_TOL_DEFAULT)
+        { 
+              setTermConditions( termconditions );
+        }
+  
+        /// Construct an EMAlg from an Evidence object, an InfAlg object, and an input stream
+        EMAlg( const Evidence &evidence, InfAlg &estep, std::istream &mstep_file );
+
+        /// Change the conditions for termination
+        /** There are two possible parameters in the PropertySet
+         *    - max_iters maximum number of iterations
+         *    - log_z_tol proportion of increase in logZ
+         *
+         *  \see hasSatisifiedTermConditions()
+         */
+        void setTermConditions( const PropertySet &p );
+
+        /// Determine if the termination conditions have been met.
+        /** There are two sufficient termination conditions:
+         *    -# the maximum number of iterations has been performed
+         *    -# the ratio of logZ increase over previous logZ is less than the 
+         *       tolerance, i.e.,
+         *       \f$ \frac{\log(Z_t) - \log(Z_{t-1})}{| \log(Z_{t-1}) | } < \mathrm{tol} \f$.
+         */
+        bool hasSatisfiedTermConditions() const;
+
+        /// Returns number of iterations done so far
+        size_t getCurrentIters() const { return _iters; }
+
+        /// Perform an iteration over all maximization steps
+        Real iterate();
+
+        /// Perform an iteration over a single MaximizationStep
+        Real iterate( MaximizationStep &mstep );
+
+        /// Iterate until termination conditions are satisfied
+        void run();
+};
+
+
+} // end of namespace dai
+
+
+#endif
diff --git a/include/dai/evidence.h b/include/dai/evidence.h
new file mode 100644 (file)
index 0000000..275c3bf
--- /dev/null
@@ -0,0 +1,108 @@
+/*  Copyright (C) 2009  Charles Vaske  [cvaske at soe dot ucsc dot edu]
+    University of California Santa Cruz
+
+    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
+*/
+
+
+/// \file
+/// \brief Defines classes Evidence and Observation
+/// \todo Describe tabular data file format
+
+
+#ifndef __defined_libdai_evidence_h
+#define __defined_libdai_evidence_h
+
+
+#include <istream>
+#include <dai/daialg.h>
+
+
+namespace dai {
+
+
+/// Stores observed values of a subset of variables
+class Observation {
+    private:
+        /// Used to store the state of some variables
+        std::map<Var, size_t> _obs;
+
+    public:
+        /// Default constructor
+        Observation() : _obs() {}
+
+        /// Get all observations
+        const std::map<Var, size_t>& observations() const { return _obs; }
+        
+        /// Add an observation
+        void addObservation( Var node, size_t setting );
+        
+        /// Clamp variables in the graphical model to their observed values
+        void applyEvidence( InfAlg& alg ) const;
+};
+
+
+/// Stores multiple joint observations of sets of variables.
+/** The Evidence class stores multiple samples, where each sample is the joint
+ *  observation of the states of some variables.
+ */
+class Evidence {
+    private:
+        /// Each sample is the joint observation of the states of some variables
+        std::vector<Observation> _samples;
+
+    public:
+        /// Default constructor
+        Evidence() : _samples() {}
+      
+        /// Read in tabular data from a stream. 
+        /** Each line contains one sample, and the first line is a header line with names.
+         */
+        void addEvidenceTabFile( std::istream& is, std::map<std::string, Var> &varMap );
+
+        /// Read in tabular data from a stream. 
+        /** Each line contains one sample, and the first line is a header line with 
+         *  variable labels which should correspond with a subset of the variables in fg.
+         */
+        void addEvidenceTabFile( std::istream& is, FactorGraph& fg );
+      
+        /// Returns number of stored samples
+        size_t nrSamples() const { return _samples.size(); }
+
+        /// @name Iterator interface
+        //@{
+        /// Iterator over the elements
+        typedef std::vector<Observation>::iterator iterator;
+        /// Constant iterator over the elements
+        typedef std::vector<Observation>::const_iterator const_iterator;
+
+        /// Returns iterator that points to the first element
+        iterator begin() { return _samples.begin(); }
+        /// Returns constant iterator that points to the first element
+        const_iterator begin() const { return _samples.begin(); }
+        /// Returns iterator that points beyond the last element
+        iterator end() { return _samples.end(); }
+        /// Returns constant iterator that points beyond the last element
+        const_iterator end() const { return _samples.end(); }
+        //@}
+};
+
+
+} // end of namespace dai
+
+
+#endif
index 3347337..6094478 100644 (file)
@@ -69,7 +69,6 @@ class ExactInf : public DAIAlgFG {
         /// @name General InfAlg interface
         //@{
         virtual ExactInf* clone() const { return new ExactInf(*this); }
-        virtual ExactInf* create() const { return new ExactInf(); }
         virtual std::string identify() const;
         virtual Factor belief( const Var &n ) const { return beliefV( findVar( n ) ); }
         virtual Factor belief( const VarSet &ns ) const;
index 5f5681c..6023df0 100644 (file)
@@ -72,6 +72,12 @@ class Exception : public std::runtime_error {
                    IMPOSSIBLE_TYPECAST,
                    INTERNAL_ERROR,
                    NOT_NORMALIZABLE,
+                   INVALID_EVIDENCE_FILE,
+                   INVALID_EVIDENCE_LINE,
+                   INVALID_EVIDENCE_OBSERVATION,
+                   INVALID_SHARED_PARAMETERS_ORDER,
+                   INVALID_SHARED_PARAMETERS_INPUT_LINE,
+                   UNKNOWN_PARAMETER_ESTIMATION_METHOD,
                    NUM_ERRORS};  // NUM_ERRORS should be the last entry
 
         /// Constructor
index b441116..42df4ce 100644 (file)
@@ -32,6 +32,7 @@
 
 
 #include <iostream>
+#include <functional>
 #include <cmath>
 #include <dai/prob.h>
 #include <dai/varset.h>
 namespace dai {
 
 
+/// Function object similar to std::divides(), but different in that dividing by zero results in zero
+template<typename T> struct divides0 : public std::binary_function<T, T, T> {
+    T operator()(const T& i, const T& j) const {
+        if( j == (T)0 )
+            return (T)0;
+        else
+            return i / j;
+    }
+};
+
+
 /// Represents a (probability) factor.
 /** Mathematically, a \e factor is a function mapping joint states of some
  *  variables to the nonnegative real numbers.
@@ -222,64 +234,58 @@ template <typename T> class TFactor {
             return *this;
         }
 
+        /// Adds the TFactor f to *this
+        TFactor<T>& operator+= (const TFactor<T>& f) { 
+            if( f._vs == _vs ) // optimize special case
+                _p += f._p;
+            else
+                *this = (*this + f); 
+            return *this;
+        }
+
+        /// Subtracts the TFactor f from *this
+        TFactor<T>& operator-= (const TFactor<T>& f) { 
+            if( f._vs == _vs ) // optimize special case
+                _p -= f._p;
+            else
+                *this = (*this - f); 
+            return *this;
+        }
+
         /// Returns product of *this with the TFactor f
         /** The product of two factors is defined as follows: if 
          *  \f$f : \prod_{l\in L} X_l \to [0,\infty)\f$ and \f$g : \prod_{m\in M} X_m \to [0,\infty)\f$, then
          *  \f[fg : \prod_{l\in L\cup M} X_l \to [0,\infty) : x \mapsto f(x_L) g(x_M).\f]
          */
-        TFactor<T> operator* (const TFactor<T>& f) const;
+        TFactor<T> operator* (const TFactor<T>& f) const {
+            return pointwiseOp(*this,f,std::multiplies<T>());
+        }
 
         /// Returns quotient of *this by the TFactor f
         /** The quotient of two factors is defined as follows: if 
          *  \f$f : \prod_{l\in L} X_l \to [0,\infty)\f$ and \f$g : \prod_{m\in M} X_m \to [0,\infty)\f$, then
          *  \f[\frac{f}{g} : \prod_{l\in L\cup M} X_l \to [0,\infty) : x \mapsto \frac{f(x_L)}{g(x_M)}.\f]
          */
-        TFactor<T> operator/ (const TFactor<T>& f) const;
-
-        /// Adds the TFactor f to *this
-        /** \pre this->vars() == f.vars()
-         */
-        TFactor<T>& operator+= (const TFactor<T>& f) { 
-#ifdef DAI_DEBUG
-            assert( f._vs == _vs );
-#endif
-            _p += f._p;
-            return *this;
-        }
-
-        /// Subtracts the TFactor f from *this
-        /** \pre this->vars() == f.vars()
-         */
-        TFactor<T>& operator-= (const TFactor<T>& f) { 
-#ifdef DAI_DEBUG
-            assert( f._vs == _vs );
-#endif
-            _p -= f._p;
-            return *this;
+        TFactor<T> operator/ (const TFactor<T>& f) const {
+            return pointwiseOp(*this,f,divides0<T>());
         }
 
         /// Returns sum of *this and the TFactor f
-        /** \pre this->vars() == f.vars()
+        /** The sum of two factors is defined as follows: if 
+         *  \f$f : \prod_{l\in L} X_l \to [0,\infty)\f$ and \f$g : \prod_{m\in M} X_m \to [0,\infty)\f$, then
+         *  \f[f+g : \prod_{l\in L\cup M} X_l \to [0,\infty) : x \mapsto f(x_L) + g(x_M).\f]
          */
         TFactor<T> operator+ (const TFactor<T>& f) const {
-#ifdef DAI_DEBUG
-            assert( f._vs == _vs );
-#endif
-            TFactor<T> sum(*this); 
-            sum._p += f._p; 
-            return sum; 
+            return pointwiseOp(*this,f,std::plus<T>());
         }
 
         /// Returns *this minus the TFactor f
-        /** \pre this->vars() == f.vars()
+        /** The difference of two factors is defined as follows: if 
+         *  \f$f : \prod_{l\in L} X_l \to [0,\infty)\f$ and \f$g : \prod_{m\in M} X_m \to [0,\infty)\f$, then
+         *  \f[f-g : \prod_{l\in L\cup M} X_l \to [0,\infty) : x \mapsto f(x_L) - g(x_M).\f]
          */
         TFactor<T> operator- (const TFactor<T>& f) const {
-#ifdef DAI_DEBUG
-            assert( f._vs == _vs );
-#endif
-            TFactor<T> sum(*this); 
-            sum._p -= f._p; 
-            return sum; 
+            return pointwiseOp(*this,f,std::minus<T>());
         }
 
 
@@ -372,6 +378,9 @@ template <typename T> class TFactor {
         /// Returns marginal on ns, obtained by summing out all variables except those in ns, and normalizing the result if normed==true
         TFactor<T> marginal(const VarSet & ns, bool normed=true) const;
 
+        /// Returns max-marginal on ns, obtained by maximizing all variables except those in ns, and normalizing the result if normed==true
+        TFactor<T> maxMarginal(const VarSet & ns, bool normed=true) const;
+
         /// Embeds this factor in a larger VarSet
         /** \pre vars() should be a subset of ns
          *
@@ -383,7 +392,7 @@ template <typename T> class TFactor {
             if( _vs == ns )
                 return *this;
             else
-                return (*this) * TFactor<T>(ns / _vs, 1);
+                return (*this) * TFactor<T>(ns / _vs, (T)1);
         }
 
         /// Returns true if *this has NaN values
@@ -428,40 +437,39 @@ template<typename T> TFactor<T> TFactor<T>::marginal(const VarSet & ns, bool nor
 }
 
 
-template<typename T> TFactor<T> TFactor<T>::operator* (const TFactor<T>& f) const {
-    if( f._vs == _vs ) { // optimizate special case
-        TFactor<T> prod(*this); 
-        prod._p *= f._p; 
-        return prod; 
-    } else {
-        TFactor<T> prod( _vs | f._vs, 0.0 );
+template<typename T> TFactor<T> TFactor<T>::maxMarginal(const VarSet & ns, bool normed) const {
+    VarSet res_ns = ns & _vs;
+
+    TFactor<T> res( res_ns, 0.0 );
 
-        IndexFor i1(_vs, prod._vs);
-        IndexFor i2(f._vs, prod._vs);
+    IndexFor i_res( res_ns, _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];
 
-        for( size_t i = 0; i < prod._p.size(); i++, ++i1, ++i2 )
-            prod._p[i] += _p[i1] * f._p[i2];
+    if( normed )
+        res.normalize( Prob::NORMPROB );
 
-        return prod;
-    }
+    return res;
 }
 
 
-template<typename T> TFactor<T> TFactor<T>::operator/ (const TFactor<T>& f) const {
-    if( f._vs == _vs ) { // optimizate special case
-        TFactor<T> quot(*this); 
-        quot._p /= f._p; 
-        return quot; 
+template<typename T, typename binaryOp> TFactor<T> pointwiseOp( const TFactor<T> &f, const TFactor<T> &g, binaryOp op ) {
+    if( f.vars() == g.vars() ) { // optimizate special case
+        TFactor<T> result(f); 
+        for( size_t i = 0; i < result.states(); i++ )
+            result[i] = op( result[i], g[i] );
+        return result; 
     } else {
-        TFactor<T> quot( _vs | f._vs, 0.0 );
+        TFactor<T> result( f.vars() | g.vars(), 0.0 );
 
-        IndexFor i1(_vs, quot._vs);
-        IndexFor i2(f._vs, quot._vs);
+        IndexFor i1(f.vars(), result.vars());
+        IndexFor i2(g.vars(), result.vars());
 
-        for( size_t i = 0; i < quot._p.size(); i++, ++i1, ++i2 )
-            quot._p[i] += _p[i1] / f._p[i2];
+        for( size_t i = 0; i < result.states(); i++, ++i1, ++i2 )
+            result[i] = op( f[i1], g[i2] );
 
-        return quot;
+        return result;
     }
 }
 
@@ -554,7 +562,7 @@ template<typename T> Real MutualInfo(const TFactor<T> &f) {
     VarSet::const_iterator it = f.vars().begin();
     Var i = *it; it++; Var j = *it;
     TFactor<T> projection = f.marginal(i) * f.marginal(j);
-    return real( dist( f.normalized(), projection, Prob::DISTKL ) );
+    return dist( f.normalized(), projection, Prob::DISTKL );
 }
 
 
index 863b694..5beae95 100644 (file)
@@ -113,9 +113,6 @@ class FactorGraph {
         /// Clone *this (virtual copy constructor)
         virtual FactorGraph* clone() const { return new FactorGraph(); }
 
-        /// Create (virtual default constructor)
-        virtual FactorGraph* create() const { return new FactorGraph(*this); }
-
         /// Returns const reference to i'th variable
         const Var & var(size_t i) const { return _vars[i]; }
         /// Returns const reference to all factors
@@ -255,7 +252,7 @@ class FactorGraph {
         void ReadFromFile(const char *filename);
 
         /// Writes a FactorGraph to a file
-        void WriteToFile(const char *filename) const;
+        void WriteToFile(const char *filename, size_t precision=15) const;
 
         /// Writes a FactorGraph to a GraphViz .dot file
         void printDot( std::ostream& os ) const;
index 2e5ec68..d63d0d4 100644 (file)
@@ -74,7 +74,6 @@ class Gibbs : public DAIAlgFG {
         /// @name General InfAlg interface
         //@{
         virtual Gibbs* clone() const { return new Gibbs(*this); }
-        virtual Gibbs* create() const { return new Gibbs(); }
         virtual std::string identify() const { return std::string(Name) + printProperties(); }
         virtual Factor belief( const Var &n ) const;
         virtual Factor belief( const VarSet &ns ) const;
index d2728e0..caa09f2 100644 (file)
@@ -98,7 +98,6 @@ class HAK : public DAIAlgRG {
         /// @name General InfAlg interface
         //@{
         virtual HAK* clone() const { return new HAK(*this); }
-        virtual HAK* create() const { return new HAK(); }
         virtual std::string identify() const;
         virtual Factor belief( const Var &n ) const;
         virtual Factor belief( const VarSet &ns ) const;
index aaff0e2..5bcc096 100644 (file)
@@ -86,7 +86,6 @@ class JTree : public DAIAlgRG {
         /// @name General InfAlg interface
         //@{
         virtual JTree* clone() const { return new JTree(*this); }
-        virtual JTree* create() const { return new JTree(); }
         virtual std::string identify() const;
         virtual Factor belief( const Var &n ) const;
         virtual Factor belief( const VarSet &ns ) const;
index 3bfb51a..536dc8d 100644 (file)
@@ -107,7 +107,6 @@ class LC : public DAIAlgFG {
         /// @name General InfAlg interface
         //@{
         virtual LC* clone() const { return new LC(*this); }
-        virtual LC* create() const { return new LC(); }
         virtual std::string identify() const;
         virtual Factor belief( const Var &n ) const { return( _beliefs[findVar(n)] ); }
         virtual Factor belief( const VarSet &/*ns*/ ) const { DAI_THROW(NOT_IMPLEMENTED); return Factor(); }
index ebdb577..e983ce7 100644 (file)
@@ -80,7 +80,6 @@ class MF : public DAIAlgFG {
         /// @name General InfAlg interface
         //@{
         virtual MF* clone() const { return new MF(*this); }
-        virtual MF* create() const { return new MF(); }
         virtual std::string identify() const;
         virtual Factor belief( const Var &n ) const;
         virtual Factor belief( const VarSet &ns ) const;
index 39bf995..fa6b7ca 100644 (file)
@@ -104,7 +104,6 @@ class MR : public DAIAlgFG {
         /// @name General InfAlg interface
         //@{
         virtual MR* clone() const { return new MR(*this); }
-        virtual MR* create() const { return new MR(); }
         virtual std::string identify() const;
         virtual Factor belief( const Var &n ) const;
         virtual Factor belief( const VarSet &/*ns*/ ) const { DAI_THROW(NOT_IMPLEMENTED); return Factor(); }
index 73c87e0..4b85374 100644 (file)
@@ -455,7 +455,7 @@ template <typename T> class TProb {
                 arg = i;
               }
             }
-            return make_pair(arg,max);
+            return std::make_pair(arg,max);
         }
 
         /// Normalizes vector using the specified norm
index 8100120..979915d 100644 (file)
@@ -110,9 +110,6 @@ class RegionGraph : public FactorGraph {
         /// Clone *this (virtual copy constructor)
         virtual RegionGraph* clone() const { return new RegionGraph(*this); }
 
-        /// Create (virtual default constructor)
-        virtual RegionGraph* create() const { return new RegionGraph(); }
-
         /// Set the content of the I'th factor and make a backup of its old content if backup == true
         virtual void setFactor( size_t I, const Factor &newFactor, bool backup = false ) {
             FactorGraph::setFactor( I, newFactor, backup ); 
index df3e685..832e382 100644 (file)
@@ -151,7 +151,6 @@ class TreeEP : public JTree {
         /// @name General InfAlg interface
         //@{
         virtual TreeEP* clone() const { return new TreeEP(*this); }
-        virtual TreeEP* create() const { return new TreeEP(); }
         virtual std::string identify() const;
         virtual Real logZ() const;
         virtual void init();
index c939b73..89c1dee 100644 (file)
@@ -240,6 +240,12 @@ class Diffs : public std::vector<double> {
 };
 
 
+/// Split a string into tokens
+void tokenizeString( const std::string& s,
+                   std::vector<std::string>& outTokens,
+                   const std::string& delim="\t\n" );
+
+
 } // end of namespace dai
 
 
index 40ef21d..6cff0bb 100644 (file)
@@ -37,6 +37,7 @@
 #include <set>
 #include <cassert>
 #include <limits>
+#include <climits>   // Work-around for bug in boost graph library
 
 #include <boost/graph/adjacency_list.hpp>
 #include <boost/graph/prim_minimum_spanning_tree.hpp>
@@ -203,17 +204,21 @@ template<typename T> DEdgeVec MinSpanningTreePrims( const WeightedGraph<T> &G )
 /** Uses implementation in Boost Graph Library.
  */
 template<typename T> DEdgeVec MaxSpanningTreePrims( const WeightedGraph<T> & Graph ) {
-    T maxweight = Graph.begin()->second;
-    for( typename WeightedGraph<T>::const_iterator it = Graph.begin(); it != Graph.end(); it++ )
-        if( it->second > maxweight )
-            maxweight = it->second;
-    // make a copy of the graph
-    WeightedGraph<T> gr( Graph );
-    // invoke MinSpanningTreePrims with negative weights
-    // (which have to be shifted to satisfy positivity criterion)
-    for( typename WeightedGraph<T>::iterator it = gr.begin(); it != gr.end(); it++ )
-        it->second = maxweight - it->second;
-    return MinSpanningTreePrims( gr );
+    if( Graph.size() == 0 )
+        return DEdgeVec();
+    else {
+        T maxweight = Graph.begin()->second;
+        for( typename WeightedGraph<T>::const_iterator it = Graph.begin(); it != Graph.end(); it++ )
+            if( it->second > maxweight )
+                maxweight = it->second;
+        // make a copy of the graph
+        WeightedGraph<T> gr( Graph );
+        // invoke MinSpanningTreePrims with negative weights
+        // (which have to be shifted to satisfy positivity criterion)
+        for( typename WeightedGraph<T>::iterator it = gr.begin(); it != gr.end(); it++ )
+            it->second = maxweight - it->second;
+        return MinSpanningTreePrims( gr );
+    }
 }
 
 
index eb851cf..938fa09 100644 (file)
@@ -145,18 +145,6 @@ void BP::init() {
 
 
 void BP::findMaxResidual( size_t &i, size_t &_I ) {
-/*
-    i = 0;
-    _I = 0;
-    double maxres = residual( i, _I );
-    for( size_t j = 0; j < nrVars(); ++j )
-        foreach( const Neighbor &I, nbV(j) )
-            if( residual( j, I.iter ) > maxres ) {
-                i = j;
-                _I = I.iter;
-                maxres = residual( i, _I );
-            }
-*/
     assert( !_lut.empty() );
     LutType::const_iterator largestEl = _lut.end();
     --largestEl;
@@ -169,41 +157,36 @@ void BP::calcNewMessage( size_t i, size_t _I ) {
     // calculate updated message I->i
     size_t I = nbV(i,_I);
 
-    if( !DAI_BP_FAST ) {
-        /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
-        Factor prod( factor( I ) );
-        foreach( const Neighbor &j, nbF(I) )
-            if( j != i ) {     // for all j in I \ i
-                foreach( const Neighbor &J, nbV(j) )
-                    if( J != I ) {     // for all J in nb(j) \ I 
-                        prod *= Factor( var(j), message(j, J.iter) );
-                    }
-            }
-        newMessage(i,_I) = prod.marginal( var(i) ).p();
-    } else {
-        /* OPTIMIZED VERSION */
-        Prob prod( factor(I).p() );
-        if( props.logdomain ) 
-            prod.takeLog();
+    Factor Fprod( factor(I) );
+    Prob &prod = Fprod.p();
+    if( props.logdomain ) 
+        prod.takeLog();
+
+    // Calculate product of incoming messages and factor I
+    foreach( const Neighbor &j, nbF(I) )
+        if( j != i ) { // for all j in I \ i
+            // prod_j will be the product of messages coming into j
+            Prob prod_j( var(j).states(), props.logdomain ? 0.0 : 1.0 ); 
+            foreach( const Neighbor &J, nbV(j) )
+                if( J != I ) { // for all J in nb(j) \ I 
+                    if( props.logdomain )
+                        prod_j += message( j, J.iter );
+                    else
+                        prod_j *= message( j, J.iter );
+                }
 
-        // Calculate product of incoming messages and factor I
-        foreach( const Neighbor &j, nbF(I) ) {
-            if( j != i ) {     // for all j in I \ i
+            // multiply prod with prod_j
+            if( !DAI_BP_FAST ) {
+                /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
+                if( props.logdomain )
+                    Fprod += Factor( var(j), prod_j );
+                else
+                    Fprod *= Factor( var(j), prod_j );
+            } else {
+                /* OPTIMIZED VERSION */
                 size_t _I = j.dual;
                 // 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);
-
-                // prod_j will be the product of messages coming into j
-                Prob prod_j( var(j).states(), props.logdomain ? 0.0 : 1.0 ); 
-                foreach( const Neighbor &J, nbV(j) )
-                    if( J != I ) { // for all J in nb(j) \ I 
-                        if( props.logdomain )
-                            prod_j += message( j, J.iter );
-                        else
-                            prod_j *= message( j, J.iter );
-                    }
-
-                // multiply prod with prod_j
                 for( size_t r = 0; r < prod.size(); ++r )
                     if( props.logdomain )
                         prod[r] += prod_j[ind[r]];
@@ -211,13 +194,23 @@ void BP::calcNewMessage( size_t i, size_t _I ) {
                         prod[r] *= prod_j[ind[r]];
             }
         }
-        if( props.logdomain ) {
-            prod -= prod.max();
-            prod.takeExp();
-        }
 
-        // Marginalize onto i
-        Prob marg( var(i).states(), 0.0 );
+    if( props.logdomain ) {
+        prod -= prod.max();
+        prod.takeExp();
+    }
+
+    // Marginalize onto i
+    Prob marg;
+    if( !DAI_BP_FAST ) {
+        /* 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 */
+        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 ) 
@@ -228,14 +221,14 @@ void BP::calcNewMessage( size_t i, size_t _I ) {
                 if( prod[r] > marg[ind[r]] ) 
                     marg[ind[r]] = prod[r];
         marg.normalize();
-
-        // Store result
-        if( props.logdomain )
-            newMessage(i,_I) = marg.log();
-        else
-            newMessage(i,_I) = marg;
     }
 
+    // Store result
+    if( props.logdomain )
+        newMessage(i,_I) = marg.log();
+    else
+        newMessage(i,_I) = marg;
+
     // Update the residual if necessary
     if( props.updates == Properties::UpdateType::SEQMAX )
         updateResidual( i, _I , dist( newMessage( i, _I ), message( i, _I ), Prob::DISTLINF ) );
@@ -360,34 +353,46 @@ void BP::calcBeliefV( size_t i, Prob &p ) const {
 
 
 void BP::calcBeliefF( size_t I, Prob &p ) const {
-    p = factor(I).p();
-    if( props.logdomain )
-        p.takeLog();
+    Factor Fprod( factor( I ) );
+    Prob &prod = Fprod.p();
 
-    foreach( const Neighbor &j, nbF(I) ) {
-        size_t _I = j.dual;
-        // 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);
+    if( props.logdomain ) 
+        prod.takeLog();
 
+    foreach( const Neighbor &j, nbF(I) ) {
         // prod_j will be the product of messages coming into j
-        Prob prod_j( var(j).states(), props.logdomain ? 0.0 : 1.0 ); 
-        foreach( const Neighbor &J, nbV(j) ) {
+        Prob prod_j( var(j).states(), props.logdomain ? 0.0 : 1.0 );
+        foreach( const Neighbor &J, nbV(j) )
             if( J != I ) { // for all J in nb(j) \ I 
                 if( props.logdomain )
                     prod_j += newMessage( j, J.iter );
                 else
                     prod_j *= newMessage( j, J.iter );
             }
-        }
 
-        // multiply p with prod_j
-        for( size_t r = 0; r < p.size(); ++r ) {
+        // multiply prod with prod_j
+        if( !DAI_BP_FAST ) {
+            /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
             if( props.logdomain )
-                p[r] += prod_j[ind[r]];
+                Fprod += Factor( var(j), prod_j );
             else
-                p[r] *= prod_j[ind[r]];
+                Fprod *= Factor( var(j), prod_j );
+        } else {
+            /* OPTIMIZED VERSION */
+            size_t _I = j.dual;
+            // 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 ) {
+                if( props.logdomain )
+                    prod[r] += prod_j[ind[r]];
+                else
+                    prod[r] *= prod_j[ind[r]];
+            }
         }
     }
+
+    p = prod;
 }
 
 
@@ -399,12 +404,26 @@ Factor BP::beliefV( size_t i ) const {
         p -= p.max();
         p.takeExp();
     }
-
     p.normalize();
+
     return( Factor( var(i), p ) );
 }
 
 
+Factor BP::beliefF( size_t I ) const {
+    Prob p;
+    calcBeliefF( I, p );
+
+    if( props.logdomain ) {
+        p -= p.max();
+        p.takeExp();
+    }
+    p.normalize();
+
+    return( Factor( factor(I).vars(), p ) );
+}
+
+
 Factor BP::belief( const Var &n ) const {
     return( beliefV( findVar( n ) ) );
 }
@@ -434,37 +453,6 @@ Factor BP::belief( const VarSet &ns ) const {
 }
 
 
-Factor BP::beliefF( size_t I ) const {
-    if( !DAI_BP_FAST ) {
-        /*  UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
-
-        Factor prod( factor(I) );
-        foreach( const Neighbor &j, nbF(I) ) {
-            foreach( const Neighbor &J, nbV(j) ) {
-                if( J != I )  // for all J in nb(j) \ I
-                    prod *= Factor( var(j), newMessage(j, J.iter) );
-            }
-        }
-        return prod.normalized();
-    } else {
-        /* OPTIMIZED VERSION */
-
-        Prob prod;
-        calcBeliefF( I, prod );
-
-        if( props.logdomain ) {
-            prod -= prod.max();
-            prod.takeExp();
-        }
-        prod.normalize();
-
-        Factor result( factor(I).vars(), prod );
-
-        return( result );
-    }
-}
-
-
 Real BP::logZ() const {
     Real sum = 0.0;
     for(size_t i = 0; i < nrVars(); ++i )
diff --git a/src/emalg.cpp b/src/emalg.cpp
new file mode 100644 (file)
index 0000000..8e93700
--- /dev/null
@@ -0,0 +1,346 @@
+/*  Copyright (C) 2009  Charles Vaske  [cvaske at soe dot ucsc dot edu]
+    University of California Santa Cruz
+
+    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
+*/
+
+
+#include <dai/util.h>
+#include <dai/emalg.h>
+
+
+namespace dai {
+
+
+std::map<std::string, ParameterEstimation::ParamEstFactory> *ParameterEstimation::_registry = NULL;
+
+
+void ParameterEstimation::loadDefaultRegistry() {
+    _registry = new std::map<std::string, ParamEstFactory>();
+    (*_registry)["ConditionalProbEstimation"] = CondProbEstimation::factory;
+}
+
+
+ParameterEstimation* ParameterEstimation::construct( const std::string &method, const PropertySet &p ) {
+    if( _registry == NULL )
+        loadDefaultRegistry();
+    std::map<std::string, ParamEstFactory>::iterator i = _registry->find(method);
+    if( i == _registry->end() )
+        DAI_THROW(UNKNOWN_PARAMETER_ESTIMATION_METHOD);
+    ParamEstFactory factory = i->second;
+    return factory(p);
+}
+
+
+ParameterEstimation* CondProbEstimation::factory( const PropertySet &p ) {
+    size_t target_dimension = p.getStringAs<size_t>("target_dim");
+    size_t total_dimension = p.getStringAs<size_t>("total_dim");
+    Real pseudo_count = 1;
+    if( p.hasKey("pseudo_count") )
+        pseudo_count = p.getStringAs<Real>("pseudo_count");
+    return new CondProbEstimation( target_dimension, Prob( total_dimension, pseudo_count ) );
+}
+
+
+CondProbEstimation::CondProbEstimation( size_t target_dimension, const Prob &pseudocounts ) 
+  : _target_dim(target_dimension), _stats(pseudocounts), _initial_stats(pseudocounts) 
+{
+    if( _stats.size() % _target_dim )
+        DAI_THROW(MALFORMED_PROPERTY);
+}
+
+
+void CondProbEstimation::addSufficientStatistics( const Prob &p ) {
+    _stats += p;
+}
+
+
+Prob CondProbEstimation::estimate() {
+    // normalize pseudocounts
+    for( size_t parent = 0; parent < _stats.size(); parent += _target_dim ) {
+        // calculate norm
+        Real norm = 0.0;
+        size_t top = parent + _target_dim;
+        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;
+    }
+    // reset _stats to _initial_stats
+    Prob result = _stats;
+    _stats = _initial_stats;
+    return result;
+}
+
+
+Permute SharedParameters::calculatePermutation( const std::vector<Var> &varorder, VarSet &outVS ) {
+    // Collect all labels and dimensions, and order them in vs
+    std::vector<size_t> dims;
+    dims.reserve( varorder.size() );
+    std::vector<long> labels;
+    labels.reserve( varorder.size() );
+    for( size_t i = 0; i < varorder.size(); i++ ) {
+        dims.push_back( varorder[i].states() );
+        labels.push_back( varorder[i].label() );
+        outVS |= varorder[i];
+    }
+  
+    // Construct the sigma array for the permutation object
+    std::vector<size_t> sigma;
+    sigma.reserve( dims.size() );
+    for( VarSet::iterator set_iterator = outVS.begin(); sigma.size() < dims.size(); ++set_iterator )
+        sigma.push_back( find(labels.begin(), labels.end(), set_iterator->label()) - labels.begin() );
+  
+    return Permute( dims, sigma );
+}
+
+
+void SharedParameters::setPermsAndVarSetsFromVarOrders() {
+    if( _varorders.size() == 0 )
+        return;
+    if( _estimation == NULL )
+        DAI_THROW(INVALID_SHARED_PARAMETERS_ORDER);
+
+    // Construct the permutation objects and the varsets
+    for( FactorOrientations::const_iterator foi = _varorders.begin(); foi != _varorders.end(); ++foi ) {
+        VarSet vs;
+        _perms[foi->first] = calculatePermutation( foi->second, vs );
+        _varsets[foi->first] = vs;
+        if( _estimation->probSize() != vs.nrStates() )
+            DAI_THROW(INVALID_SHARED_PARAMETERS_ORDER);
+    }
+}
+
+
+SharedParameters::SharedParameters( std::istream &is, const FactorGraph &fg_varlookup )
+  : _varsets(), _perms(), _varorders(), _estimation(NULL), _deleteEstimation(true) 
+{
+    // Read the desired parameter estimation method from the stream
+    std::string est_method;
+    PropertySet props;
+    is >> est_method;
+    is >> props;
+
+    // Construct a corresponding object
+    _estimation = ParameterEstimation::construct( est_method, props );
+
+    // Read in the factors that are to be estimated
+    size_t num_factors;
+    is >> num_factors;
+    for( size_t sp_i = 0; sp_i < num_factors; ++sp_i ) {
+        std::string line;
+        while( line.size() == 0 && getline(is, line) )
+            ;
+
+        std::vector<std::string> fields;
+        tokenizeString(line, fields, " \t");
+
+        // Lookup the factor in the factorgraph
+        if( fields.size() < 1 )
+            DAI_THROW(INVALID_SHARED_PARAMETERS_INPUT_LINE);
+        std::istringstream iss;
+        iss.str( fields[0] );
+        size_t factor;
+        iss >> factor;
+        const VarSet &vs = fg_varlookup.factor(factor).vars();
+        if( fields.size() != vs.size() + 1 )
+            DAI_THROW(INVALID_SHARED_PARAMETERS_INPUT_LINE);
+
+        // Construct the vector of Vars
+        std::vector<Var> var_order;
+        var_order.reserve( vs.size() );
+        for( size_t fi = 1; fi < fields.size(); ++fi ) {
+            // Lookup a single variable by label
+            long label;
+            std::istringstream labelparse( fields[fi] );
+            labelparse >> label;
+            VarSet::const_iterator vsi = vs.begin();
+            for( ; vsi != vs.end(); ++vsi )
+                if( vsi->label() == label ) 
+                    break;
+            if( vsi == vs.end() )
+                DAI_THROW(INVALID_SHARED_PARAMETERS_INPUT_LINE);
+            var_order.push_back( *vsi );
+        }
+        _varorders[factor] = var_order;
+    }
+
+    // Calculate the necessary permutations
+    setPermsAndVarSetsFromVarOrders();
+}
+
+
+SharedParameters::SharedParameters( const SharedParameters &sp ) 
+  : _varsets(sp._varsets), _perms(sp._perms), _varorders(sp._varorders), _estimation(sp._estimation), _deleteEstimation(sp._deleteEstimation)
+{
+    // If sp owns its _estimation object, we should clone it instead
+    if( _deleteEstimation )
+        _estimation = _estimation->clone();
+}
+
+
+SharedParameters::SharedParameters( const FactorOrientations &varorders, ParameterEstimation *estimation )
+  : _varsets(), _perms(), _varorders(varorders), _estimation(estimation), _deleteEstimation(false) 
+{
+    // Calculate the necessary permutations
+    setPermsAndVarSetsFromVarOrders();
+}
+
+
+void SharedParameters::collectSufficientStatistics( InfAlg &alg ) {
+    for( std::map< FactorIndex, Permute >::iterator i = _perms.begin(); i != _perms.end(); ++i ) {
+        Permute &perm = i->second;
+        VarSet &vs = _varsets[i->first];
+        
+        Factor b = alg.belief(vs);
+        Prob p( b.states(), 0.0 );
+        for( size_t entry = 0; entry < b.states(); ++entry )
+            p[entry] = b[perm.convert_linear_index(entry)];
+        _estimation->addSufficientStatistics( p );
+    }
+}
+
+
+void SharedParameters::setParameters( FactorGraph &fg ) {
+    Prob p = _estimation->estimate();
+    for( std::map<FactorIndex, Permute>::iterator i = _perms.begin(); i != _perms.end(); ++i ) {
+        Permute &perm = i->second;
+        VarSet &vs = _varsets[i->first];
+        
+        Factor f( vs, 0.0 );
+        for( size_t entry = 0; entry < f.states(); ++entry )
+            f[perm.convert_linear_index(entry)] = p[entry];
+
+        fg.setFactor( i->first, f );
+    }
+}
+
+
+MaximizationStep::MaximizationStep( std::istream &is, const FactorGraph &fg_varlookup ) : _params() {
+    size_t num_params = -1;
+    is >> num_params;
+    _params.reserve( num_params );
+    for( size_t i = 0; i < num_params; ++i )
+        _params.push_back( SharedParameters( is, fg_varlookup ) );
+}
+
+
+void MaximizationStep::addExpectations( InfAlg &alg ) {
+    for( size_t i = 0; i < _params.size(); ++i )
+        _params[i].collectSufficientStatistics( alg );
+}
+
+
+void MaximizationStep::maximize( FactorGraph &fg ) {
+    for( size_t i = 0; i < _params.size(); ++i )
+        _params[i].setParameters( fg );
+}
+
+
+const std::string EMAlg::MAX_ITERS_KEY("max_iters");
+const std::string EMAlg::LOG_Z_TOL_KEY("log_z_tol");
+const size_t EMAlg::MAX_ITERS_DEFAULT = 30;
+const Real EMAlg::LOG_Z_TOL_DEFAULT = 0.01;
+
+
+EMAlg::EMAlg( const Evidence &evidence, InfAlg &estep, std::istream &msteps_file )
+  : _evidence(evidence), _estep(estep), _msteps(), _iters(0), _lastLogZ(), _max_iters(MAX_ITERS_DEFAULT), _log_z_tol(LOG_Z_TOL_DEFAULT)
+{
+    msteps_file.exceptions( std::istream::eofbit | std::istream::failbit | std::istream::badbit );
+    size_t num_msteps = -1;
+    msteps_file >> num_msteps;
+    _msteps.reserve(num_msteps);
+    for( size_t i = 0; i < num_msteps; ++i )
+        _msteps.push_back( MaximizationStep( msteps_file, estep.fg() ) );
+}      
+
+
+void EMAlg::setTermConditions( const PropertySet &p ) {
+    if( p.hasKey(MAX_ITERS_KEY) )
+        _max_iters = p.getStringAs<size_t>(MAX_ITERS_KEY);
+    if( p.hasKey(LOG_Z_TOL_KEY) )
+        _log_z_tol = p.getStringAs<Real>(LOG_Z_TOL_KEY);
+}
+
+
+bool EMAlg::hasSatisfiedTermConditions() const {
+    if( _iters >= _max_iters )
+        return true;
+    else if( _lastLogZ.size() < 3 )
+        // need at least 2 to calculate ratio
+        // Also, throw away first iteration, as the parameters may not
+        // have been normalized according to the estimation method
+        return false;
+    else {
+        Real current = _lastLogZ[_lastLogZ.size() - 1];
+        Real previous = _lastLogZ[_lastLogZ.size() - 2];
+        if( previous == 0 )
+            return false;
+        Real diff = current - previous;
+        if( diff < 0 ) {
+            std::cerr << "Error: in EM log-likehood decreased from " << previous << " to " << current << std::endl;
+            return true;
+        }
+        return diff / abs(previous) <= _log_z_tol;
+    }
+}
+
+
+Real EMAlg::iterate( MaximizationStep &mstep ) {
+    Real logZ = 0;
+
+    // Expectation calculation
+    for( Evidence::const_iterator e = _evidence.begin(); e != _evidence.end(); ++e ) {
+        InfAlg* clamped = _estep.clone();
+        e->applyEvidence( *clamped );
+        clamped->init();
+        clamped->run();
+      
+        logZ += clamped->logZ();
+
+        mstep.addExpectations( *clamped );
+
+        delete clamped;
+    }
+    
+    // Maximization of parameters
+    mstep.maximize( _estep.fg() );
+
+    return logZ;
+}
+
+
+Real EMAlg::iterate() {
+    Real likelihood;
+    for( size_t i = 0; i < _msteps.size(); ++i )
+        likelihood = iterate( _msteps[i] );
+    _lastLogZ.push_back( likelihood );
+    ++_iters;
+    return likelihood;
+}
+
+
+void EMAlg::run() {
+    while( !hasSatisfiedTermConditions() )
+        iterate();
+}
+
+
+} // end of namespace dai
diff --git a/src/evidence.cpp b/src/evidence.cpp
new file mode 100644 (file)
index 0000000..e003807
--- /dev/null
@@ -0,0 +1,98 @@
+/*  Copyright (C) 2009  Charles Vaske  [cvaske at soe dot ucsc dot edu]
+    University of California Santa Cruz
+
+    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
+*/
+
+
+#include <sstream>
+#include <string>
+#include <cstdlib>
+
+#include <dai/util.h>
+#include <dai/evidence.h>
+
+
+namespace dai {
+
+
+void Observation::addObservation( Var node, size_t setting ) {
+    _obs[node] = setting;
+}
+
+
+void Observation::applyEvidence( InfAlg &alg ) const {
+    for( std::map<Var, size_t>::const_iterator i = _obs.begin(); i != _obs.end(); ++i )
+        alg.clamp( i->first, i->second );
+}
+  
+
+void Evidence::addEvidenceTabFile( std::istream &is, FactorGraph &fg ) {
+    std::map<std::string, Var> varMap;
+    for( std::vector<Var>::const_iterator v = fg.vars().begin(); v != fg.vars().end(); ++v ) {
+        std::stringstream s;
+        s << v->label();
+        varMap[s.str()] = *v;
+    }
+
+    addEvidenceTabFile( is, varMap );
+}
+
+
+void Evidence::addEvidenceTabFile( std::istream &is, std::map<std::string, Var> &varMap ) {
+    std::string line;
+    getline( is, line );
+    
+    // Parse header
+    std::vector<std::string> header_fields;
+    tokenizeString( line, header_fields );
+    std::vector<std::string>::const_iterator p_field = header_fields.begin();
+    if( p_field == header_fields.end() ) 
+        DAI_THROW(INVALID_EVIDENCE_LINE);
+
+    std::vector<Var> vars;
+    for( ; p_field != header_fields.end(); ++p_field ) {
+        std::map<std::string, Var>::iterator elem = varMap.find( *p_field );
+        if( elem == varMap.end() )
+            DAI_THROW(INVALID_EVIDENCE_FILE);
+        vars.push_back( elem->second );
+    }
+    
+    // Read samples
+    while( getline(is, line) ) {
+        std::vector<std::string> fields;
+        tokenizeString( line, fields );
+        if( fields.size() != vars.size() ) 
+            DAI_THROW(INVALID_EVIDENCE_LINE);
+        
+        Observation sampleData;
+        for( size_t i = 0; i < vars.size(); ++i ) {
+            if( fields[i].size() > 0 ) { // skip if missing observation
+                if( fields[i].find_first_not_of("0123456789") != std::string::npos )
+                    DAI_THROW(INVALID_EVIDENCE_OBSERVATION);
+                size_t state = atoi( fields[i].c_str() );
+                if( state >= vars[i].states() )
+                    DAI_THROW(INVALID_EVIDENCE_OBSERVATION);
+                sampleData.addObservation( vars[i], state );
+            }
+        }
+        _samples.push_back( sampleData );
+    } // finished sample line
+}
+
+
+} // end of namespace dai
index 738a4ad..84b2e22 100644 (file)
@@ -40,7 +40,13 @@ namespace dai {
         "FactorGraph is not connected",
         "Impossible typecast",
         "Internal error",
-        "Quantity not normalizable"
+        "Quantity not normalizable",
+        "Cannot parse Evidence file",
+        "Cannot parse Evidence line",
+        "Invalid observation in Evidence file",
+        "Invalid variable order in SharedParameters",
+        "Variable order line in EM file is invalid",
+        "Unrecognized parameter estimation method"
     }; 
 
 
index 94f7589..d8b9fea 100644 (file)
@@ -21,6 +21,7 @@
 
 
 #include <iostream>
+#include <iomanip>
 #include <iterator>
 #include <map>
 #include <set>
@@ -100,11 +101,8 @@ ostream& operator << (ostream& os, const FactorGraph& fg) {
                 nr_nonzeros++;
         os << nr_nonzeros << endl;
         for( size_t k = 0; k < fg.factor(I).states(); k++ )
-            if( fg.factor(I)[k] != 0.0 ) {
-                char buf[20];
-                sprintf(buf,"%18.14g", fg.factor(I)[k]);
-                os << k << " " << buf << endl;
-            }
+            if( fg.factor(I)[k] != 0.0 )
+                os << k << " " << setw(os.precision()+4) << fg.factor(I)[k] << endl;
     }
 
     return(os);
@@ -268,10 +266,11 @@ void FactorGraph::ReadFromFile( const char *filename ) {
 }
 
 
-void FactorGraph::WriteToFile( const char *filename ) const {
+void FactorGraph::WriteToFile( const char *filename, size_t precision ) const {
     ofstream outfile;
     outfile.open( filename );
     if( outfile.is_open() ) {
+        outfile.precision( precision );
         outfile << *this;
         outfile.close();
     } else
index 149d7b0..1169e35 100644 (file)
@@ -24,7 +24,7 @@
  *                                                                 * 
  * This is a MEX-file for MATLAB.                                  *
  *                                                                 * 
- *   [logZ,q,md] = dai(psi,method,opts);                           *
+ *   [logZ,q,md,qv,qf] = dai(psi,method,opts);                     *
  *                                                                 * 
  *=================================================================*/
 
@@ -45,6 +45,7 @@ using namespace dai;
 #define METHOD_IN       prhs[1]
 #define OPTS_IN         prhs[2]
 #define NR_IN           3
+#define NR_IN_OPT       0
 
 
 /* Output Arguments */
@@ -52,16 +53,19 @@ using namespace dai;
 #define LOGZ_OUT        plhs[0]
 #define Q_OUT           plhs[1]
 #define MD_OUT          plhs[2]
+#define QV_OUT          plhs[3]
+#define QF_OUT          plhs[4]
 #define NR_OUT          3
+#define NR_OUT_OPT      2
 
 
 void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[] )
 { 
     size_t buflen;
     
-    /* Check for proper number of arguments */
-    if( (nrhs != NR_IN) || (nlhs != NR_OUT) ) { 
-        mexErrMsgTxt("Usage: [logZ,q,md] = dai(psi,method,opts)\n\n"
+    // Check for proper number of arguments
+    if( ((nrhs < NR_IN) || (nrhs > NR_IN + NR_IN_OPT)) || ((nlhs < NR_OUT) || (nlhs > NR_OUT + NR_OUT_OPT)) ) {
+        mexErrMsgTxt("Usage: [logZ,q,md,qv,qf] = dai(psi,method,opts)\n\n"
         "\n"
         "INPUT:  psi        = linear cell array containing the factors \n"
         "                     psi{i} should be a structure with a Member field\n"
@@ -71,7 +75,9 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[] )
         "\n"
         "OUTPUT: logZ       = approximation of the logarithm of the partition sum.\n"
         "        q          = linear cell array containing all final beliefs.\n"
-        "        md         = maxdiff (final linf-dist between new and old single node beliefs).\n");
+        "        md         = maxdiff (final linf-dist between new and old single node beliefs).\n"
+        "        qv         = linear cell array containing all variable beliefs.\n"
+        "        qf         = linear cell array containing all factor beliefs.\n");
     } 
     
     char *method;
@@ -119,6 +125,21 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[] )
     MD_OUT = mxCreateDoubleMatrix(1,1,mxREAL);
     *(mxGetPr(MD_OUT)) = maxdiff;
 
+    if( nlhs >= 4 ) {
+        vector<Factor> qv;
+        qv.reserve( fg.nrVars() );
+        for( size_t i = 0; i < fg.nrVars(); i++ )
+            qv.push_back( obj->belief( fg.var(i) ) );
+        QV_OUT = Factors2mx( qv );
+    }
+
+    if( nlhs >= 5 ) {
+        vector<Factor> qf;
+        qf.reserve( fg.nrFactors() );
+        for( size_t I = 0; I < fg.nrFactors(); I++ )
+            qf.push_back( obj->belief( fg.factor(I).vars() ) );
+        QF_OUT = Factors2mx( qf );
+    }
     
     return;
 }
index f15db26..b5e1d34 100644 (file)
@@ -105,5 +105,18 @@ int rnd_int( int min, int max ) {
     return (int)floor(_uni_rnd() * (max + 1 - min) + min);
 }
 
+void tokenizeString(const std::string& s,
+                    std::vector<std::string>& outTokens,
+                    const std::string& delim)
+{
+    size_t start = 0;
+    while (start < s.size()) {
+        size_t end = s.find_first_of(delim, start);
+        if (end > s.size())
+            end = s.size();
+        outTokens.push_back(s.substr(start, end - start));
+        start = end + 1;
+    }
+}
 
 } // end of namespace dai
index 8e47385..918ec52 100644 (file)
@@ -253,124 +253,132 @@ int main( int argc, char *argv[] ) {
         cout << "error and relative logZ error (comparing with the results of" << endl;
         cout << "<method0>, the base method, for which one can use JTREE_HUGIN)." << endl << endl;
         cout << opts_required << opts_optional << endl;
+#ifdef DAI_DEBUG
+        cout << "This is a debugging (unoptimised) build of libDAI." << endl;
+#endif
         return 1;
     }
 
-    // Read aliases
-    map<string,string> Aliases;
-    if( !aliases.empty() ) {
-        ifstream infile;
-        infile.open (aliases.c_str());
-        if (infile.is_open()) {
-            while( true ) {
-                string line;
-                getline(infile,line);
-                if( infile.fail() )
-                    break;
-                if( (!line.empty()) && (line[0] != '#') ) {
-                    string::size_type pos = line.find(':',0);
-                    if( pos == string::npos )
-                        throw "Invalid alias";
-                    else {
-                        string::size_type posl = line.substr(0, pos).find_last_not_of(" \t");
-                        string key = line.substr(0, posl + 1);
-                        string::size_type posr = line.substr(pos + 1, line.length()).find_first_not_of(" \t");
-                        string val = line.substr(pos + 1 + posr, line.length());
-                        Aliases[key] = val;
+    try {
+        // Read aliases
+        map<string,string> Aliases;
+        if( !aliases.empty() ) {
+            ifstream infile;
+            infile.open (aliases.c_str());
+            if (infile.is_open()) {
+                while( true ) {
+                    string line;
+                    getline(infile,line);
+                    if( infile.fail() )
+                        break;
+                    if( (!line.empty()) && (line[0] != '#') ) {
+                        string::size_type pos = line.find(':',0);
+                        if( pos == string::npos )
+                            throw string("Invalid alias");
+                        else {
+                            string::size_type posl = line.substr(0, pos).find_last_not_of(" \t");
+                            string key = line.substr(0, posl + 1);
+                            string::size_type posr = line.substr(pos + 1, line.length()).find_first_not_of(" \t");
+                            string val = line.substr(pos + 1 + posr, line.length());
+                            Aliases[key] = val;
+                        }
                     }
                 }
-            }
-            infile.close();
-        } else
-            throw "Error opening aliases file";
-    }
-
-    FactorGraph fg;
-    fg.ReadFromFile( filename.c_str() );
-
-    vector<Factor> q0;
-    double logZ0 = 0.0;
-
-    cout.setf( ios_base::scientific );
-    cout.precision( 3 );
-
-    cout << "# " << filename << endl;
-    cout.width( 39 );
-    cout << left << "# METHOD" << "\t";
-    if( report_time )
-        cout << right << "SECONDS  " << "\t";
-    if( report_iters )
-        cout << "ITERS" << "\t";
-    cout << "MAX ERROR" << "\t";
-    cout << "AVG ERROR" << "\t";
-    cout << "LOGZ ERROR" << "\t";
-    cout << "MAXDIFF" << "\t";
-    cout << endl;
-
-    for( size_t m = 0; m < methods.size(); m++ ) {
-        pair<string, PropertySet> meth = parseMethod( methods[m], Aliases );
-
-        if( vm.count("tol") )
-            meth.second.Set("tol",tol);
-        if( vm.count("maxiter") )
-            meth.second.Set("maxiter",maxiter);
-        if( vm.count("verbose") )
-            meth.second.Set("verbose",verbose);
-        TestDAI piet(fg, meth.first, meth.second );
-        piet.doDAI();
-        if( m == 0 ) {
-            q0 = piet.q;
-            logZ0 = piet.logZ;
+                infile.close();
+            } else
+                throw string("Error opening aliases file");
         }
-        piet.calcErrs(q0);
 
+        FactorGraph fg;
+        fg.ReadFromFile( filename.c_str() );
+
+        vector<Factor> q0;
+        double logZ0 = 0.0;
+
+        cout.setf( ios_base::scientific );
+        cout.precision( 3 );
+
+        cout << "# " << filename << endl;
         cout.width( 39 );
-        cout << left << methods[m] << "\t";
+        cout << left << "# METHOD" << "\t";
         if( report_time )
-            cout << right << piet.time << "\t";
-        if( report_iters ) {
-            if( piet.has_iters ) {
-                cout << piet.iters << "\t";
-            } else {
-                cout << "N/A  \t";
+            cout << right << "SECONDS  " << "\t";
+        if( report_iters )
+            cout << "ITERS" << "\t";
+        cout << "MAX ERROR" << "\t";
+        cout << "AVG ERROR" << "\t";
+        cout << "LOGZ ERROR" << "\t";
+        cout << "MAXDIFF" << "\t";
+        cout << endl;
+
+        for( size_t m = 0; m < methods.size(); m++ ) {
+            pair<string, PropertySet> meth = parseMethod( methods[m], Aliases );
+
+            if( vm.count("tol") )
+                meth.second.Set("tol",tol);
+            if( vm.count("maxiter") )
+                meth.second.Set("maxiter",maxiter);
+            if( vm.count("verbose") )
+                meth.second.Set("verbose",verbose);
+            TestDAI piet(fg, meth.first, meth.second );
+            piet.doDAI();
+            if( m == 0 ) {
+                q0 = piet.q;
+                logZ0 = piet.logZ;
+            }
+            piet.calcErrs(q0);
+
+            cout.width( 39 );
+            cout << left << methods[m] << "\t";
+            if( report_time )
+                cout << right << piet.time << "\t";
+            if( report_iters ) {
+                if( piet.has_iters ) {
+                    cout << piet.iters << "\t";
+                } else {
+                    cout << "N/A  \t";
+                }
             }
-        }
 
-        if( m > 0 ) {
-            cout.setf( ios_base::scientific );
-            cout.precision( 3 );
-            
-            double me = clipdouble( piet.maxErr(), 1e-9 );
-            cout << me << "\t";
-            
-            double ae = clipdouble( piet.avgErr(), 1e-9 );
-            cout << ae << "\t";
-            
-            if( piet.has_logZ ) {
-                cout.setf( ios::showpos );
-                double le = clipdouble( piet.logZ / logZ0 - 1.0, 1e-9 );
-                cout << le << "\t";
-                cout.unsetf( ios::showpos );
-            } else
-                cout << "N/A       \t";
-
-            if( piet.has_maxdiff ) {
-                double md = clipdouble( piet.maxdiff, 1e-9 );
-                if( isnan( me ) )
-                    md = me;
-                if( isnan( ae ) )
-                    md = ae;
-                cout << md << "\t";
-            } else
-                cout << "N/A    \t";
-        }
-        cout << endl;
+            if( m > 0 ) {
+                cout.setf( ios_base::scientific );
+                cout.precision( 3 );
+                
+                double me = clipdouble( piet.maxErr(), 1e-9 );
+                cout << me << "\t";
+                
+                double ae = clipdouble( piet.avgErr(), 1e-9 );
+                cout << ae << "\t";
+                
+                if( piet.has_logZ ) {
+                    cout.setf( ios::showpos );
+                    double le = clipdouble( piet.logZ / logZ0 - 1.0, 1e-9 );
+                    cout << le << "\t";
+                    cout.unsetf( ios::showpos );
+                } else
+                    cout << "N/A       \t";
+
+                if( piet.has_maxdiff ) {
+                    double md = clipdouble( piet.maxdiff, 1e-9 );
+                    if( isnan( me ) )
+                        md = me;
+                    if( isnan( ae ) )
+                        md = ae;
+                    cout << md << "\t";
+                } else
+                    cout << "N/A    \t";
+            }
+            cout << endl;
 
-        if( marginals ) {
-            for( size_t i = 0; i < piet.q.size(); i++ )
-                cout << "# " << piet.q[i] << endl;
+            if( marginals ) {
+                for( size_t i = 0; i < piet.q.size(); i++ )
+                    cout << "# " << piet.q[i] << endl;
+            }
         }
-    }
 
-    return 0;
+        return 0;
+    } catch( string &s ) {
+        cerr << "Exception: " << s << endl;
+        return 2;
+    }
 }
diff --git a/tests/testem/2var.em b/tests/testem/2var.em
new file mode 100644 (file)
index 0000000..457ebe1
--- /dev/null
@@ -0,0 +1,6 @@
+1
+
+1
+ConditionalProbEstimation [target_dim=2,total_dim=4,pseudo_count=1]
+1
+0 1 0
diff --git a/tests/testem/2var.fg b/tests/testem/2var.fg
new file mode 100644 (file)
index 0000000..279870c
--- /dev/null
@@ -0,0 +1,10 @@
+1
+
+2
+0 1
+2 2
+4
+0      0.5
+1      0.5
+2      0.5
+3      0.5
\ No newline at end of file
diff --git a/tests/testem/2var_data.tab b/tests/testem/2var_data.tab
new file mode 100644 (file)
index 0000000..e59a890
--- /dev/null
@@ -0,0 +1,21 @@
+0      1
+0      1
+0      1
+0      1
+0      1
+0      1
+0      1
+0      1
+0      1
+0      1
+0      0
+1      1
+1      1
+1      1
+1      0
+1      0
+1      0
+1      0
+1      0
+1      0
+1      0
diff --git a/tests/testem/3var.em b/tests/testem/3var.em
new file mode 100644 (file)
index 0000000..dff5410
--- /dev/null
@@ -0,0 +1,6 @@
+1
+
+1
+ConditionalProbEstimation [target_dim=2,total_dim=8,pseudo_count=1]
+1
+0 1 0 2
diff --git a/tests/testem/3var.fg b/tests/testem/3var.fg
new file mode 100644 (file)
index 0000000..69f0077
--- /dev/null
@@ -0,0 +1,14 @@
+1
+
+3
+0 1 2
+2 2 2
+8
+0      0.5
+1      0.5
+2      0.5
+3      0.5
+4      0.5
+5      0.5
+6      0.5
+7      0.5
\ No newline at end of file
diff --git a/tests/testem/hoi1_data.tab b/tests/testem/hoi1_data.tab
new file mode 100644 (file)
index 0000000..ad5f220
--- /dev/null
@@ -0,0 +1,6 @@
+0      1       2       4       6       7
+0      0       1               1       0
+1                      1       1       0
+0      1       0       0       0       1
+0      0       0       1       0       0
+1      0       0               1       0
diff --git a/tests/testem/hoi1_infer_f2.em b/tests/testem/hoi1_infer_f2.em
new file mode 100644 (file)
index 0000000..977eb5b
--- /dev/null
@@ -0,0 +1,6 @@
+1
+
+1
+ConditionalProbEstimation [target_dim=2,total_dim=8,pseudo_count=1]
+1
+2 1 2 4
diff --git a/tests/testem/hoi1_share_f0_f1_f2.em b/tests/testem/hoi1_share_f0_f1_f2.em
new file mode 100644 (file)
index 0000000..1f7756f
--- /dev/null
@@ -0,0 +1,8 @@
+1
+
+1
+ConditionalProbEstimation [target_dim=2,total_dim=8,pseudo_count=1]
+3
+2 1 2 4
+1 0 1 6
+0 6 2 7
diff --git a/tests/testem/hoi1_share_f0_f2.em b/tests/testem/hoi1_share_f0_f2.em
new file mode 100644 (file)
index 0000000..fd656c6
--- /dev/null
@@ -0,0 +1,7 @@
+1
+
+1
+ConditionalProbEstimation [target_dim=2,total_dim=8,pseudo_count=1]
+2
+2 1 2 4
+0 6 2 7
diff --git a/tests/testem/runtests b/tests/testem/runtests
new file mode 100755 (executable)
index 0000000..2442db8
--- /dev/null
@@ -0,0 +1,11 @@
+#!/bin/bash
+TMPFILE1=`mktemp /var/tmp/testem.XXXXXX`
+trap 'rm -f $TMPFILE1' 0 1 15
+
+./testem 2var.fg 2var_data.tab 2var.em > $TMPFILE1
+./testem 3var.fg 2var_data.tab 3var.em >> $TMPFILE1
+./testem ../hoi1.fg hoi1_data.tab hoi1_share_f0_f2.em >> $TMPFILE1
+./testem ../hoi1.fg hoi1_data.tab hoi1_share_f0_f1_f2.em >> $TMPFILE1
+diff -s $TMPFILE1 testem.out || exit 1 
+
+rm -f $TMPFILE1
diff --git a/tests/testem/runtests.bat b/tests/testem/runtests.bat
new file mode 100755 (executable)
index 0000000..44bdf2e
--- /dev/null
@@ -0,0 +1,7 @@
+testem 2var.fg 2var_data.tab 2var.em > testem.out.tmp
+testem 3var.fg 2var_data.tab 3var.em >> testem.out.tmp
+testem ..\hoi1.fg hoi1_data.tab hoi1_share_f0_f2.em >> testem.out.tmp
+testem ..\hoi1.fg hoi1_data.tab hoi1_share_f0_f1_f2.em >> testem.out.tmp
+diff -s testem.out.tmp testem.out
+
+del testem.out.tmp
diff --git a/tests/testem/testem.cpp b/tests/testem/testem.cpp
new file mode 100644 (file)
index 0000000..a77e68c
--- /dev/null
@@ -0,0 +1,79 @@
+/*  Copyright (C) 2009  Charles Vaske  [cvaske at soe dot ucsc dot edu]
+    University of California Santa Cruz
+
+    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
+*/
+
+
+#include <iostream>
+#include <fstream>
+#include <string>
+
+#include <dai/factorgraph.h>
+#include <dai/evidence.h>
+#include <dai/alldai.h>
+
+
+using namespace std;
+using namespace dai;
+
+
+void usage( const string &msg ) {
+    cerr << msg << endl;
+    cerr << "Usage:" << endl;
+    cerr << " testem factorgraph.fg evidence.tab emconfig.em" << endl;
+    exit( 1 );
+}
+
+
+int main( int argc, char** argv ) {
+    if( argc != 4 )
+        usage("Incorrect number of arguments.");
+    
+    FactorGraph fg;
+    fg.ReadFromFile( argv[1] );
+
+    PropertySet infprops;
+    infprops.Set( "verbose", (size_t)1 );
+    infprops.Set( "updates", string("HUGIN") );
+    InfAlg* inf = newInfAlg( "JTREE", fg, infprops );
+    inf->init();
+
+    Evidence e;
+    ifstream estream( argv[2] );
+    e.addEvidenceTabFile( estream, fg );
+
+    cout << "Number of samples: " << e.nrSamples() << endl;
+    for( Evidence::iterator ps = e.begin(); ps != e.end(); ps++ )
+        cout << "Sample #" << (ps - e.begin()) << " has " << ps->observations().size() << " observations." << endl;
+
+    ifstream emstream( argv[3] );
+    EMAlg em(e, *inf, emstream);
+
+    while( !em.hasSatisfiedTermConditions() ) {
+        Real l = em.iterate();
+        cout << "Iteration " << em.getCurrentIters() << " likelihood: " << l <<endl;
+    }
+
+    cout << endl << "Inferred Factor Graph:" << endl << "######################" << endl;
+    cout.precision(12);
+    cout << inf->fg();
+
+    delete inf;
+
+    return 0;
+}
diff --git a/tests/testem/testem.out b/tests/testem/testem.out
new file mode 100644 (file)
index 0000000..5e777cb
--- /dev/null
@@ -0,0 +1,184 @@
+Number of samples: 20
+Sample #0 has 2 observations.
+Sample #1 has 2 observations.
+Sample #2 has 2 observations.
+Sample #3 has 2 observations.
+Sample #4 has 2 observations.
+Sample #5 has 2 observations.
+Sample #6 has 2 observations.
+Sample #7 has 2 observations.
+Sample #8 has 2 observations.
+Sample #9 has 2 observations.
+Sample #10 has 2 observations.
+Sample #11 has 2 observations.
+Sample #12 has 2 observations.
+Sample #13 has 2 observations.
+Sample #14 has 2 observations.
+Sample #15 has 2 observations.
+Sample #16 has 2 observations.
+Sample #17 has 2 observations.
+Sample #18 has 2 observations.
+Sample #19 has 2 observations.
+Iteration 1 likelihood: -13.8629
+Iteration 2 likelihood: -9.56675
+Iteration 3 likelihood: -9.56675
+
+Inferred Factor Graph:
+######################
+1
+
+2
+0 1 
+2 2 
+4
+0   0.166666666667
+1   0.666666666667
+2   0.833333333333
+3   0.333333333333
+Number of samples: 20
+Sample #0 has 2 observations.
+Sample #1 has 2 observations.
+Sample #2 has 2 observations.
+Sample #3 has 2 observations.
+Sample #4 has 2 observations.
+Sample #5 has 2 observations.
+Sample #6 has 2 observations.
+Sample #7 has 2 observations.
+Sample #8 has 2 observations.
+Sample #9 has 2 observations.
+Sample #10 has 2 observations.
+Sample #11 has 2 observations.
+Sample #12 has 2 observations.
+Sample #13 has 2 observations.
+Sample #14 has 2 observations.
+Sample #15 has 2 observations.
+Sample #16 has 2 observations.
+Sample #17 has 2 observations.
+Sample #18 has 2 observations.
+Sample #19 has 2 observations.
+Iteration 1 likelihood: 0
+Iteration 2 likelihood: 3.97035
+Iteration 3 likelihood: 3.97035
+
+Inferred Factor Graph:
+######################
+1
+
+3
+0 1 2 
+2 2 2 
+8
+0   0.214285714286
+1   0.642857142857
+2   0.785714285714
+3   0.357142857143
+4   0.214285714286
+5   0.642857142857
+6   0.785714285714
+7   0.357142857143
+Number of samples: 5
+Sample #0 has 5 observations.
+Sample #1 has 4 observations.
+Sample #2 has 6 observations.
+Sample #3 has 6 observations.
+Sample #4 has 5 observations.
+Iteration 1 likelihood: 11.1646
+Iteration 2 likelihood: 1.53723
+Iteration 3 likelihood: 1.64691
+Iteration 4 likelihood: 1.67497
+Iteration 5 likelihood: 1.68191
+
+Inferred Factor Graph:
+######################
+3
+
+3
+2 6 7 
+2 2 2 
+8
+0   0.398343360806
+1   0.351464146565
+2   0.601656639194
+3   0.648535853435
+4   0.804844687957
+5    0.67374245803
+6   0.195155312043
+7    0.32625754197
+
+3
+0 1 6 
+2 2 2 
+8
+0    1.03521336269
+1    1.54747895212
+2    2.31765218974
+3    1.28041900719
+4      4.922079813
+5    2.52725575019
+6   0.831279296316
+7   0.262805630803
+
+3
+1 2 4 
+2 2 2 
+8
+0   0.398343360806
+1   0.601656639194
+2   0.351464146565
+3   0.648535853435
+4   0.804844687957
+5   0.195155312043
+6    0.67374245803
+7    0.32625754197
+Number of samples: 5
+Sample #0 has 5 observations.
+Sample #1 has 4 observations.
+Sample #2 has 6 observations.
+Sample #3 has 6 observations.
+Sample #4 has 5 observations.
+Iteration 1 likelihood: 11.1646
+Iteration 2 likelihood: -7.29331
+Iteration 3 likelihood: -7.261
+
+Inferred Factor Graph:
+######################
+3
+
+3
+2 6 7 
+2 2 2 
+8
+0   0.495312199726
+1   0.499107942908
+2   0.504687800274
+3   0.500892057092
+4   0.640416549958
+5   0.495122293924
+6   0.359583450042
+7   0.504877706076
+
+3
+0 1 6 
+2 2 2 
+8
+0   0.495312199726
+1   0.504687800274
+2   0.499107942908
+3   0.500892057092
+4   0.640416549958
+5   0.359583450042
+6   0.495122293924
+7   0.504877706076
+
+3
+1 2 4 
+2 2 2 
+8
+0   0.495312199726
+1   0.504687800274
+2   0.499107942908
+3   0.500892057092
+4   0.640416549958
+5   0.359583450042
+6   0.495122293924
+7   0.504877706076
index 55630a1..cdae018 100755 (executable)
@@ -3,6 +3,6 @@ TMPFILE1=`mktemp /var/tmp/testfast.XXXXXX`
 trap 'rm -f $TMP_FILE' 0 1 15
 
 ./testall testfast.fg > $TMPFILE1
-diff -s $TMPFILE1 testfast.out
+diff -s $TMPFILE1 testfast.out || exit 1
 
 rm -f $TMPFILE1