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

34 files changed:
ChangeLog
Makefile
Makefile.CYGWIN
Makefile.LINUX
Makefile.MACOSX
Makefile.WINDOWS
README
doxygen.conf
include/dai/alldai.h
include/dai/bbp.h [new file with mode: 0644]
include/dai/bipgraph.h
include/dai/bp.h
include/dai/bp_dual.h [new file with mode: 0644]
include/dai/cbp.h [new file with mode: 0644]
include/dai/daialg.h
include/dai/doc.h
include/dai/factorgraph.h
include/dai/lc.h
include/dai/mr.h
include/dai/properties.h
include/dai/util.h
scripts/regenerate-properties [new file with mode: 0755]
src/alldai.cpp
src/bbp.cpp [new file with mode: 0644]
src/bp.cpp
src/bp_dual.cpp [new file with mode: 0644]
src/cbp.cpp [new file with mode: 0644]
src/factorgraph.cpp
src/gibbs.cpp
tests/aliases.conf
tests/testall
tests/testbbp.cpp [new file with mode: 0644]
tests/testdai.cpp
tests/testfast.out

index 734e9c3..45e54c5 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,4 +1,3 @@
-Branch master:
 * Added work-around for bug in Boost Graph Library
 * Improvements of TFactor<T>:
   - Extended functionality of TFactor<T>::operators +,-,+=,-= to handle different VarSets
@@ -9,7 +8,26 @@ Branch master:
 * [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)
+* [Frederik Eaton] Updated doxygen.conf to version 1.5.9
+* [Frederik Eaton] Added newInfAlgFromString to alldai.h/cpp
+* [Frederik Eaton] Added FactorGraph::clampVar and FactorGraph::clampFactor
+* [Frederik Eaton] Changed BP to be able to record order in which messages are sent (for use by BBP)
+* [Frederik Eaton] Improved documentation of BipartiteGraph::Neighbor
+* [Frederik Eaton]: Extended InfAlg interface with beliefV() and beliefF() methods;
+* [Frederik Eaton]: Improved PropertySet class:
+  - Added default constructor
+  - Allow construction from string
+  - Added method allKeys() to produce list of keys
+  - In getStringAs(), check for typeid(ValueType) before typeid(std::string)
+    (this fixes a strange bug for empty string property)
+* [Frederik Eaton]: Added some utility macros (DAI_PV, DAI_DMSG, DAI_ACCMUT, DAI_IFVERB)
+  and functions (concat(),rnd()) to include/dai/util.h
+* [Frederik Eaton] Added backwards compatibility layer for edge interface to
+  BipartiteGraph and FactorGraph (which will be obsoleted soon)
+* [Frederik Eaton] Added BP_dual, BBP and CBP algorithms
 * [Frederik Eaton] Added Gibbs::state() accessors/mutators
 * [Frederik Eaton] Fixed Gibbs::getVarDist(size_t) to return uniform 
   distribution if no state is allowed
index 12eff2f..749847b 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -78,9 +78,14 @@ ifdef WITH_GIBBS
   CCFLAGS:=$(CCFLAGS) -DDAI_WITH_GIBBS
   OBJECTS:=$(OBJECTS) gibbs$(OE)
 endif
+ifdef WITH_CBP
+  CCFLAGS:=$(CCFLAGS) -DDAI_WITH_CBP
+  OBJECTS:=$(OBJECTS) bbp$(OE) cbp$(OE) bp_dual$(OE)
+endif
+
 
 # Define standard libDAI header dependencies
-HEADERS=$(INC)/bipgraph.h $(INC)/index.h $(INC)/var.h $(INC)/factor.h $(INC)/varset.h $(INC)/smallset.h $(INC)/prob.h $(INC)/daialg.h $(INC)/properties.h $(INC)/alldai.h $(INC)/enum.h $(INC)/exceptions.h
+HEADERS=$(INC)/bipgraph.h $(INC)/index.h $(INC)/var.h $(INC)/factor.h $(INC)/varset.h $(INC)/smallset.h $(INC)/prob.h $(INC)/daialg.h $(INC)/properties.h $(INC)/alldai.h $(INC)/enum.h $(INC)/exceptions.h $(INC)/util.h
 
 # Setup final command for C++ compiler and MEX
 ifdef DEBUG
@@ -104,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/testem/testem$(EE)
+tests : tests/testdai$(EE) tests/testem/testem$(EE) tests/testbbp$(EE)
 
 utils : utils/createfg$(EE) utils/fg2dot$(EE) utils/fginfo$(EE)
 
@@ -126,6 +131,15 @@ exactinf$(OE) : $(SRC)/exactinf.cpp $(INC)/exactinf.h $(HEADERS)
 bp$(OE) : $(SRC)/bp.cpp $(INC)/bp.h $(HEADERS)
        $(CC) -c $(SRC)/bp.cpp
 
+bp_dual$(OE) : $(SRC)/bp_dual.cpp $(INC)/bp_dual.h $(HEADERS)
+       $(CC) -c $(SRC)/bp_dual.cpp
+
+bbp$(OE) : $(SRC)/bbp.cpp $(INC)/bbp.h $(INC)/bp_dual.h $(HEADERS)
+       $(CC) -c $(SRC)/bbp.cpp
+
+cbp$(OE) : $(SRC)/cbp.cpp $(INC)/cbp.h $(INC)/bbp.h $(INC)/bp_dual.h $(HEADERS)
+       $(CC) -c $(SRC)/cbp.cpp
+
 lc$(OE) : $(SRC)/lc.cpp $(INC)/lc.h $(HEADERS)
        $(CC) -c $(SRC)/lc.cpp
 
@@ -202,6 +216,9 @@ tests/testdai$(EE) : tests/testdai.cpp $(HEADERS) $(LIB)/libdai$(LE)
 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)
+
 
 # MATLAB INTERFACE
 ###################
@@ -289,8 +306,7 @@ 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)
-       -rm tests/testem/testem$(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
index c62a51f..3a11de9 100644 (file)
@@ -15,6 +15,7 @@ WITH_TREEEP=true
 WITH_JTREE=true
 WITH_MR=true
 WITH_GIBBS=true
+WITH_CBP=true
 # Build doxygen documentation?
 WITH_DOC=true
 # Build with debug info?
index 66988e3..2bc8b2f 100644 (file)
@@ -16,6 +16,7 @@ WITH_TREEEP=true
 WITH_JTREE=true
 WITH_MR=true
 WITH_GIBBS=true
+WITH_CBP=true
 # Build doxygen documentation?
 WITH_DOC=true
 # Build with debug info?
index e7fc46f..b176114 100644 (file)
@@ -15,6 +15,7 @@ WITH_TREEEP=true
 WITH_JTREE=true
 WITH_MR=true
 WITH_GIBBS=true
+WITH_CBP=true
 # Build doxygen documentation?
 WITH_DOC=true
 # Build with debug info?
index 788ebd5..9728988 100644 (file)
@@ -16,6 +16,7 @@ WITH_TREEEP=true
 WITH_JTREE=true
 WITH_MR=true
 WITH_GIBBS=true
+WITH_CBP=true
 # Build doxygen documentation?
 WITH_DOC=
 # Build with debug info?
diff --git a/README b/README
index c164cd6..06c144a 100644 (file)
--- a/README
+++ b/README
@@ -88,7 +88,8 @@ Currently, libDAI supports the following (approximate) inference methods:
     * Generalized Belief Propagation [YFW05];
     * Double-loop GBP [HAK03];
     * Various variants of Loop Corrected Belief Propagation [MoK07, MoR05];
-    * Gibbs sampler.
+    * Gibbs sampler;
+    * Conditioned BP [EaG09].
 
 In addition, libDAI supports parameter learning of conditional probability
 tables by Expectation Maximization.
index 83f9a31..c26af9d 100644 (file)
@@ -1,4 +1,4 @@
-# Doxyfile 1.5.5
+# Doxyfile 1.5.9
 
 # This file describes the settings to be used by the documentation system
 # doxygen (www.doxygen.org) for a project
@@ -54,11 +54,11 @@ CREATE_SUBDIRS         = NO
 # information to generate all constant output in the proper language. 
 # The default language is English, other supported languages are: 
 # Afrikaans, Arabic, Brazilian, Catalan, Chinese, Chinese-Traditional, 
-# Croatian, Czech, Danish, Dutch, Farsi, Finnish, French, German, Greek
-# Hungarian, Italian, Japanese, Japanese-en (Japanese with English messages), 
-# Korean, Korean-en, Lithuanian, Norwegian, Macedonian, Persian, Polish
-# Portuguese, Romanian, Russian, Serbian, Slovak, Slovene, Spanish, Swedish
-# and Ukrainian.
+# Croatian, Czech, Danish, Dutch, Esperanto, Farsi, Finnish, French, German
+# Greek, Hungarian, Italian, Japanese, Japanese-en (Japanese with English 
+# messages), Korean, Korean-en, Lithuanian, Norwegian, Macedonian, Persian
+# Polish, Portuguese, Romanian, Russian, Serbian, Serbian-Cyrilic, Slovak
+# Slovene, Spanish, Swedish, Ukrainian, and Vietnamese.
 
 OUTPUT_LANGUAGE        = English
 
@@ -155,13 +155,6 @@ QT_AUTOBRIEF           = NO
 
 MULTILINE_CPP_IS_BRIEF = NO
 
-# If the DETAILS_AT_TOP tag is set to YES then Doxygen 
-# will output the detailed description near the top, like JavaDoc.
-# If set to NO, the detailed description appears after the member 
-# documentation.
-
-DETAILS_AT_TOP         = YES
-
 # If the INHERIT_DOCS tag is set to YES (the default) then an undocumented 
 # member inherits the documentation from any documented member that it 
 # re-implements.
@@ -214,6 +207,17 @@ OPTIMIZE_FOR_FORTRAN   = NO
 
 OPTIMIZE_OUTPUT_VHDL   = NO
 
+# Doxygen selects the parser to use depending on the extension of the files it parses. 
+# With this tag you can assign which parser to use for a given extension. 
+# Doxygen has a built-in mapping, but you can override or extend it using this tag. 
+# The format is ext=language, where ext is a file extension, and language is one of 
+# the parsers supported by doxygen: IDL, Java, Javascript, C#, C, C++, D, PHP, 
+# Objective-C, Python, Fortran, VHDL, C, C++. For instance to make doxygen treat 
+# .inc files as Fortran files (default is PHP), and .f files as C (default is Fortran), 
+# use: inc=Fortran f=C. Note that for custom extensions you also need to set FILE_PATTERNS otherwise the files are not read by doxygen.
+
+EXTENSION_MAPPING      = 
+
 # If you use STL classes (i.e. std::string, std::vector, etc.) but do not want 
 # to include (a tag file for) the STL sources as input, then you should 
 # set this tag to YES in order to let doxygen match functions declarations and 
@@ -223,7 +227,7 @@ OPTIMIZE_OUTPUT_VHDL   = NO
 
 BUILTIN_STL_SUPPORT    = YES
 
-# If you use Microsoft's C++/CLI language, you should set this option to YES to
+# If you use Microsoft's C++/CLI language, you should set this option to YES to 
 # enable parsing support.
 
 CPP_CLI_SUPPORT        = NO
@@ -234,6 +238,15 @@ CPP_CLI_SUPPORT        = NO
 
 SIP_SUPPORT            = NO
 
+# For Microsoft's IDL there are propget and propput attributes to indicate getter 
+# and setter methods for a property. Setting this option to YES (the default) 
+# will make doxygen to replace the get and set methods by a property in the 
+# documentation. This will only work if the methods are indeed getting or 
+# setting a simple type. If this is not the case, or you want to show the 
+# methods anyway, you should set this option to NO.
+
+IDL_PROPERTY_SUPPORT   = YES
+
 # If member grouping is used in the documentation and the DISTRIBUTE_GROUP_DOC 
 # tag is set to YES, then doxygen will reuse the documentation of the first 
 # member in the group (if any) for the other members of the group. By default 
@@ -259,6 +272,22 @@ SUBGROUPING            = YES
 
 TYPEDEF_HIDES_STRUCT   = NO
 
+# The SYMBOL_CACHE_SIZE determines the size of the internal cache use to 
+# determine which symbols to keep in memory and which to flush to disk. 
+# When the cache is full, less often used symbols will be written to disk. 
+# For small to medium size projects (<1000 input files) the default value is 
+# probably good enough. For larger projects a too small cache size can cause 
+# doxygen to be busy swapping symbols to and from disk most of the time 
+# causing a significant performance penality. 
+# If the system has enough physical memory increasing the cache will improve the 
+# performance by keeping more symbols in memory. Note that the value works on 
+# a logarithmic scale so increasing the size by one will rougly double the 
+# memory usage. The cache size is given by this formula: 
+# 2^(16+SYMBOL_CACHE_SIZE). The valid range is 0..9, the default is 0, 
+# corresponding to a cache size of 2^16 = 65536 symbols
+
+SYMBOL_CACHE_SIZE      = 0
+
 #---------------------------------------------------------------------------
 # Build related configuration options
 #---------------------------------------------------------------------------
@@ -386,7 +415,7 @@ SORT_GROUP_NAMES       = NO
 # sorted by fully-qualified names, including namespaces. If set to 
 # NO (the default), the class list will be sorted only by class name, 
 # not including the namespace part. 
-# Note: This option is not very useful if HIDE_SCOPE_NAMES is set to YES.
+# Note: This option is not very useful if HIDE_SCOPE_NAMES is set to YES. 
 # Note: This option applies only to the class list, not to the 
 # alphabetical list.
 
@@ -443,6 +472,19 @@ SHOW_USED_FILES        = YES
 
 SHOW_DIRECTORIES       = YES
 
+# Set the SHOW_FILES tag to NO to disable the generation of the Files page. 
+# This will remove the Files entry from the Quick Index and from the 
+# Folder Tree View (if specified). The default is YES.
+
+SHOW_FILES             = YES
+
+# Set the SHOW_NAMESPACES tag to NO to disable the generation of the 
+# Namespaces page. 
+# This will remove the Namespaces entry from the Quick Index 
+# and from the Folder Tree View (if specified). The default is YES.
+
+SHOW_NAMESPACES        = YES
+
 # The FILE_VERSION_FILTER tag can be used to specify a program or script that 
 # doxygen should invoke to get the current version for each file (typically from 
 # the version control system). Doxygen will invoke the program by executing (via 
@@ -453,6 +495,15 @@ SHOW_DIRECTORIES       = YES
 
 FILE_VERSION_FILTER    = 
 
+# The LAYOUT_FILE tag can be used to specify a layout file which will be parsed by 
+# doxygen. The layout file controls the global structure of the generated output files 
+# in an output format independent way. The create the layout file that represents 
+# doxygen's defaults, run doxygen with the -l option. You can optionally specify a 
+# file name after the option, if omitted DoxygenLayout.xml will be used as the name 
+# of the layout file.
+
+LAYOUT_FILE            = 
+
 #---------------------------------------------------------------------------
 # configuration options related to warning and progress messages
 #---------------------------------------------------------------------------
@@ -513,7 +564,8 @@ WARN_LOGFILE           =
 # directories like "/usr/src/myproject". Separate the files or directories 
 # with spaces.
 
-INPUT                  = src include/dai
+INPUT                  = src \
+                         include/dai
 
 # This tag can be used to specify the character encoding of the source files 
 # that doxygen parses. Internally doxygen uses the UTF-8 encoding, which is 
@@ -597,14 +649,17 @@ IMAGE_PATH             =
 # by executing (via popen()) the command <filter> <input-file>, where <filter> 
 # is the value of the INPUT_FILTER tag, and <input-file> is the name of an 
 # input file. Doxygen will then use the output that the filter program writes 
-# to standard output.  If FILTER_PATTERNS is specified, this tag will be 
+# to standard output. 
+# If FILTER_PATTERNS is specified, this tag will be 
 # ignored.
 
 INPUT_FILTER           = 
 
 # The FILTER_PATTERNS tag can be used to specify filters on a per file pattern 
-# basis.  Doxygen will compare the file name with each pattern and apply the 
-# filter if there is a match.  The filters are a list of the form: 
+# basis. 
+# Doxygen will compare the file name with each pattern and apply the 
+# filter if there is a match. 
+# The filters are a list of the form: 
 # pattern=filter (like *.cpp=my_cpp_filter). See INPUT_FILTER for further 
 # info on how filters are used. If FILTER_PATTERNS is empty, INPUT_FILTER 
 # is applied to all files.
@@ -639,22 +694,23 @@ INLINE_SOURCES         = NO
 
 STRIP_CODE_COMMENTS    = YES
 
-# If the REFERENCED_BY_RELATION tag is set to YES (the default) 
+# If the REFERENCED_BY_RELATION tag is set to YES 
 # then for each documented function all documented 
 # functions referencing it will be listed.
 
 REFERENCED_BY_RELATION = NO
 
-# If the REFERENCES_RELATION tag is set to YES (the default) 
+# If the REFERENCES_RELATION tag is set to YES 
 # then for each documented function all documented entities 
 # called/used by that function will be listed.
 
 REFERENCES_RELATION    = NO
 
-# If the REFERENCES_LINK_SOURCE tag is set to YES (the default)
-# and SOURCE_BROWSER tag is set to YES, then the hyperlinks from
-# functions in REFERENCES_RELATION and REFERENCED_BY_RELATION lists will
-# link to the source code.  Otherwise they will link to the documentstion.
+# If the REFERENCES_LINK_SOURCE tag is set to YES (the default) 
+# and SOURCE_BROWSER tag is set to YES, then the hyperlinks from 
+# functions in REFERENCES_RELATION and REFERENCED_BY_RELATION lists will 
+# link to the source code. 
+# Otherwise they will link to the documentation.
 
 REFERENCES_LINK_SOURCE = YES
 
@@ -743,12 +799,13 @@ HTML_STYLESHEET        =
 
 HTML_ALIGN_MEMBERS     = YES
 
-# If the GENERATE_HTMLHELP tag is set to YES, additional index files 
-# will be generated that can be used as input for tools like the 
-# Microsoft HTML help workshop to generate a compiled HTML help file (.chm) 
-# of the generated HTML documentation.
+# If the HTML_DYNAMIC_SECTIONS tag is set to YES then the generated HTML 
+# documentation will contain sections that can be hidden and shown after the 
+# page has loaded. For this to work a browser that supports 
+# JavaScript and DHTML is required (for instance Mozilla 1.0+, Firefox 
+# Netscape 6.0+, Internet explorer 5.0+, Konqueror, or Safari).
 
-GENERATE_HTMLHELP      = NO
+HTML_DYNAMIC_SECTIONS  = NO
 
 # If the GENERATE_DOCSET tag is set to YES, additional index files 
 # will be generated that can be used as input for Apple's Xcode 3 
@@ -757,7 +814,8 @@ GENERATE_HTMLHELP      = NO
 # HTML output directory. Running make will produce the docset in that 
 # directory and running "make install" will install the docset in 
 # ~/Library/Developer/Shared/Documentation/DocSets so that Xcode will find 
-# it at startup.
+# it at startup. 
+# See http://developer.apple.com/tools/creatingdocsetswithdoxygen.html for more information.
 
 GENERATE_DOCSET        = NO
 
@@ -775,13 +833,12 @@ DOCSET_FEEDNAME        = "Doxygen generated docs"
 
 DOCSET_BUNDLE_ID       = org.doxygen.Project
 
-# If the HTML_DYNAMIC_SECTIONS tag is set to YES then the generated HTML 
-# documentation will contain sections that can be hidden and shown after the 
-# page has loaded. For this to work a browser that supports 
-# JavaScript and DHTML is required (for instance Mozilla 1.0+, Firefox 
-# Netscape 6.0+, Internet explorer 5.0+, Konqueror, or Safari).
+# If the GENERATE_HTMLHELP tag is set to YES, additional index files 
+# will be generated that can be used as input for tools like the 
+# Microsoft HTML help workshop to generate a compiled HTML help file (.chm) 
+# of the generated HTML documentation.
 
-HTML_DYNAMIC_SECTIONS  = NO
+GENERATE_HTMLHELP      = NO
 
 # If the GENERATE_HTMLHELP tag is set to YES, the CHM_FILE tag can 
 # be used to specify the file name of the resulting .chm file. You 
@@ -803,6 +860,12 @@ HHC_LOCATION           =
 
 GENERATE_CHI           = NO
 
+# If the GENERATE_HTMLHELP tag is set to YES, the CHM_INDEX_ENCODING 
+# is used to encode HtmlHelp index (hhk), content (hhc) and project file 
+# content.
+
+CHM_INDEX_ENCODING     = 
+
 # If the GENERATE_HTMLHELP tag is set to YES, the BINARY_TOC flag 
 # controls whether a binary table of contents is generated (YES) or a 
 # normal table of contents (NO) in the .chm file.
@@ -814,6 +877,55 @@ BINARY_TOC             = NO
 
 TOC_EXPAND             = NO
 
+# If the GENERATE_QHP tag is set to YES and both QHP_NAMESPACE and QHP_VIRTUAL_FOLDER 
+# are set, an additional index file will be generated that can be used as input for 
+# Qt's qhelpgenerator to generate a Qt Compressed Help (.qch) of the generated 
+# HTML documentation.
+
+GENERATE_QHP           = NO
+
+# If the QHG_LOCATION tag is specified, the QCH_FILE tag can 
+# be used to specify the file name of the resulting .qch file. 
+# The path specified is relative to the HTML output folder.
+
+QCH_FILE               = 
+
+# The QHP_NAMESPACE tag specifies the namespace to use when generating 
+# Qt Help Project output. For more information please see 
+# http://doc.trolltech.com/qthelpproject.html#namespace
+
+QHP_NAMESPACE          = 
+
+# The QHP_VIRTUAL_FOLDER tag specifies the namespace to use when generating 
+# Qt Help Project output. For more information please see 
+# http://doc.trolltech.com/qthelpproject.html#virtual-folders
+
+QHP_VIRTUAL_FOLDER     = doc
+
+# If QHP_CUST_FILTER_NAME is set, it specifies the name of a custom filter to add. 
+# For more information please see 
+# http://doc.trolltech.com/qthelpproject.html#custom-filters
+
+QHP_CUST_FILTER_NAME   = 
+
+# The QHP_CUST_FILT_ATTRS tag specifies the list of the attributes of the custom filter to add.For more information please see 
+# <a href="http://doc.trolltech.com/qthelpproject.html#custom-filters">Qt Help Project / Custom Filters</a>.
+
+QHP_CUST_FILTER_ATTRS  = 
+
+# The QHP_SECT_FILTER_ATTRS tag specifies the list of the attributes this project's 
+# filter section matches. 
+# <a href="http://doc.trolltech.com/qthelpproject.html#filter-attributes">Qt Help Project / Filter Attributes</a>.
+
+QHP_SECT_FILTER_ATTRS  = 
+
+# If the GENERATE_QHP tag is set to YES, the QHG_LOCATION tag can 
+# be used to specify the location of Qt's qhelpgenerator. 
+# If non-empty doxygen will try to run qhelpgenerator on the generated 
+# .qhp file.
+
+QHG_LOCATION           = 
+
 # The DISABLE_INDEX tag can be used to turn on/off the condensed index at 
 # top of each HTML page. The value NO (the default) enables the index and 
 # the value YES disables it.
@@ -825,12 +937,20 @@ DISABLE_INDEX          = NO
 
 ENUM_VALUES_PER_LINE   = 4
 
-# If the GENERATE_TREEVIEW tag is set to YES, a side panel will be
-# generated containing a tree-like index structure (just like the one that 
+# The GENERATE_TREEVIEW tag is used to specify whether a tree-like index 
+# structure should be generated to display hierarchical information. 
+# If the tag value is set to FRAME, a side panel will be generated 
+# containing a tree-like index structure (just like the one that 
 # is generated for HTML Help). For this to work a browser that supports 
 # JavaScript, DHTML, CSS and frames is required (for instance Mozilla 1.0+, 
 # Netscape 6.0+, Internet explorer 5.0+, or Konqueror). Windows users are 
-# probably better off using the HTML help feature.
+# probably better off using the HTML help feature. Other possible values 
+# for this tag are: HIERARCHIES, which will generate the Groups, Directories, 
+# and Class Hierarchy pages using a tree view instead of an ordered list; 
+# ALL, which combines the behavior of FRAME and HIERARCHIES; and NONE, which 
+# disables this behavior completely. For backwards compatibility with previous 
+# releases of Doxygen, the values YES and NO are equivalent to FRAME and NONE 
+# respectively.
 
 GENERATE_TREEVIEW      = NO
 
@@ -840,6 +960,14 @@ GENERATE_TREEVIEW      = NO
 
 TREEVIEW_WIDTH         = 250
 
+# Use this tag to change the font size of Latex formulas included 
+# as images in the HTML documentation. The default is 10. Note that 
+# when you change the font size after a successful doxygen run you need 
+# to manually remove any form_*.png images from the HTML output directory 
+# to force them to be regenerated.
+
+FORMULA_FONTSIZE       = 10
+
 #---------------------------------------------------------------------------
 # configuration options related to the LaTeX output
 #---------------------------------------------------------------------------
@@ -916,6 +1044,10 @@ LATEX_BATCHMODE        = NO
 
 LATEX_HIDE_INDICES     = NO
 
+# If LATEX_SOURCE_CODE is set to YES then doxygen will include source code with syntax highlighting in the LaTeX output. Note that which sources are shown also depends on other settings such as SOURCE_BROWSER.
+
+LATEX_SOURCE_CODE      = NO
+
 #---------------------------------------------------------------------------
 # configuration options related to the RTF output
 #---------------------------------------------------------------------------
@@ -1052,8 +1184,10 @@ GENERATE_PERLMOD       = NO
 PERLMOD_LATEX          = NO
 
 # If the PERLMOD_PRETTY tag is set to YES the Perl module output will be 
-# nicely formatted so it can be parsed by a human reader.  This is useful 
-# if you want to understand what is going on.  On the other hand, if this 
+# nicely formatted so it can be parsed by a human reader. 
+# This is useful 
+# if you want to understand what is going on. 
+# On the other hand, if this 
 # tag is set to NO the size of the Perl module output will be much smaller 
 # and Perl will parse it just the same.
 
@@ -1081,13 +1215,13 @@ ENABLE_PREPROCESSING   = YES
 # compilation will be performed. Macro expansion can be done in a controlled 
 # way by setting EXPAND_ONLY_PREDEF to YES.
 
-MACRO_EXPANSION        = NO
+MACRO_EXPANSION        = YES
 
 # If the EXPAND_ONLY_PREDEF and MACRO_EXPANSION tags are both set to YES 
 # then the macro expansion is limited to the macros specified with the 
 # PREDEFINED and EXPAND_AS_DEFINED tags.
 
-EXPAND_ONLY_PREDEF     = NO
+EXPAND_ONLY_PREDEF     = YES
 
 # If the SEARCH_INCLUDES tag is set to YES (the default) the includes files 
 # in the INCLUDE_PATH (see below) will be search if a #include is found.
@@ -1115,14 +1249,23 @@ INCLUDE_FILE_PATTERNS  =
 # undefined via #undef or recursively expanded use the := operator 
 # instead of the = operator.
 
-PREDEFINED             = DAI_WITH_BP DAI_WITH_MF DAI_WITH_HAK DAI_WITH_LC DAI_WITH_TREEEP DAI_WITH_JTREE DAI_WITH_MR
+PREDEFINED             = DAI_WITH_BP \
+                         DAI_WITH_MF \
+                         DAI_WITH_HAK \
+                         DAI_WITH_LC \
+                         DAI_WITH_TREEEP \
+                         DAI_WITH_JTREE \
+                         DAI_WITH_MR \
+                         DAI_WITH_CBP \
+                         DAI_ACCMUT \
+                         DAI_DEBUG
 
 # If the MACRO_EXPANSION and EXPAND_ONLY_PREDEF tags are set to YES then 
 # this tag can be used to specify a list of macro names that should be expanded. 
 # The macro definition that is found in the sources will be used. 
 # Use the PREDEFINED tag if you want to use a different macro definition.
 
-EXPAND_AS_DEFINED      = 
+EXPAND_AS_DEFINED      = DAI_ACCMUT
 
 # If the SKIP_FUNCTION_MACROS tag is set to YES (the default) then 
 # doxygen's preprocessor will remove all function-like macros that are alone 
@@ -1130,7 +1273,7 @@ EXPAND_AS_DEFINED      =
 # function macros are typically used for boiler-plate code, and will confuse 
 # the parser if not removed.
 
-SKIP_FUNCTION_MACROS   = YES
+SKIP_FUNCTION_MACROS   = NO
 
 #---------------------------------------------------------------------------
 # Configuration::additions related to external references   
@@ -1140,14 +1283,16 @@ SKIP_FUNCTION_MACROS   = YES
 # Optionally an initial location of the external documentation 
 # can be added for each tagfile. The format of a tag file without 
 # this location is as follows: 
-#   TAGFILES = file1 file2 ... 
+#  
+# TAGFILES = file1 file2 ... 
 # Adding location for the tag files is done as follows: 
-#   TAGFILES = file1=loc1 "file2 = loc2" ... 
+#  
+# TAGFILES = file1=loc1 "file2 = loc2" ... 
 # where "loc1" and "loc2" can be relative or absolute paths or 
 # URLs. If a location is present for each tag, the installdox tool 
-# does not have to be run to correct the links.
-# Note that each tag file must have a unique name
-# (where the name does NOT include the path)
+# does not have to be run to correct the links. 
+# Note that each tag file must have a unique name 
+# (where the name does NOT include the path) 
 # If a tag file is not located in the directory in which doxygen 
 # is run, you must also specify the path to the tagfile here.
 
@@ -1210,6 +1355,29 @@ HIDE_UNDOC_RELATIONS   = YES
 
 HAVE_DOT               = NO
 
+# By default doxygen will write a font called FreeSans.ttf to the output 
+# directory and reference it in all dot files that doxygen generates. This 
+# font does not include all possible unicode characters however, so when you need 
+# these (or just want a differently looking font) you can specify the font name 
+# using DOT_FONTNAME. You need need to make sure dot is able to find the font, 
+# which can be done by putting it in a standard location or by setting the 
+# DOTFONTPATH environment variable or by setting DOT_FONTPATH to the directory 
+# containing the font.
+
+DOT_FONTNAME           = FreeSans
+
+# The DOT_FONTSIZE tag can be used to set the size of the font of dot graphs. 
+# The default size is 10pt.
+
+DOT_FONTSIZE           = 10
+
+# By default doxygen will tell dot to use the output directory to look for the 
+# FreeSans.ttf font (which doxygen will put there itself). If you specify a 
+# different font using DOT_FONTNAME you can set the path where dot 
+# can find it using this tag.
+
+DOT_FONTPATH           = 
+
 # If the CLASS_GRAPH and HAVE_DOT tags are set to YES then doxygen 
 # will generate a graph for each documented class showing the direct and 
 # indirect inheritance relations. Setting this tag to YES will force the 
@@ -1277,13 +1445,13 @@ GRAPHICAL_HIERARCHY    = YES
 
 # If the DIRECTORY_GRAPH, SHOW_DIRECTORIES and HAVE_DOT tags are set to YES 
 # then doxygen will show the dependencies a directory has on other directories 
-# in a graphical way. The dependency relations are determined by the #include
+# in a graphical way. The dependency relations are determined by the #include 
 # relations between the files in the directories.
 
 DIRECTORY_GRAPH        = YES
 
 # The DOT_IMAGE_FORMAT tag can be used to set the image format of the images 
-# generated by dot. Possible values are png, jpg, or gif
+# generated by dot. Possible values are png, jpg, or gif 
 # If left blank png will be used.
 
 DOT_IMAGE_FORMAT       = png
@@ -1299,7 +1467,7 @@ DOT_PATH               =
 
 DOTFILE_DIRS           = 
 
-# The MAX_DOT_GRAPH_MAX_NODES tag can be used to set the maximum number of 
+# The DOT_GRAPH_MAX_NODES tag can be used to set the maximum number of 
 # nodes that will be shown in the graph. If the number of nodes in a graph 
 # becomes larger than this value, doxygen will truncate the graph, which is 
 # visualized by representing a node as a red box. Note that doxygen if the 
@@ -1320,10 +1488,10 @@ DOT_GRAPH_MAX_NODES    = 50
 MAX_DOT_GRAPH_DEPTH    = 0
 
 # Set the DOT_TRANSPARENT tag to YES to generate images with a transparent 
-# background. This is enabled by default, which results in a transparen
-# background. Warning: Depending on the platform used, enabling this option 
-# may lead to badly anti-aliased labels on the edges of a graph (i.e. they 
-# become hard to read).
+# background. This is disabled by default, because dot on Windows does no
+# seem to support this out of the box. Warning: Depending on the platform used, 
+# enabling this option may lead to badly anti-aliased labels on the edges of 
+# a graph (i.e. they become hard to read).
 
 DOT_TRANSPARENT        = NO
 
@@ -1347,7 +1515,7 @@ GENERATE_LEGEND        = YES
 DOT_CLEANUP            = YES
 
 #---------------------------------------------------------------------------
-# Configuration::additions related to the search engine   
+# Options related to the search engine
 #---------------------------------------------------------------------------
 
 # The SEARCHENGINE tag specifies whether or not a search engine should be 
index c0010ca..862bf8d 100644 (file)
@@ -59,6 +59,9 @@
 #ifdef DAI_WITH_GIBBS
     #include <dai/gibbs.h>
 #endif
+#ifdef DAI_WITH_CBP
+    #include <dai/cbp.h>
+#endif
 
 
 /// Namespace for libDAI
@@ -74,6 +77,14 @@ namespace dai {
 InfAlg *newInfAlg( const std::string &name, const FactorGraph &fg, const PropertySet &opts );
 
 
+/// Constructs a new approximate inference algorithm.
+/** \param nameOpts The name and options of the approximate inference algorithm (should be in the format "name[opts]").
+ *  \param fg The FactorGraph that the algorithm should be applied to.
+ *  \return Returns a pointer to the new InfAlg object; it is the responsibility of the caller to delete it later.
+ */
+InfAlg *newInfAlgFromString( const std::string &nameOpts, const FactorGraph &fg );
+
+
 /// Contains the names of all approximate inference algorithms compiled into libDAI.
 static const char* DAINames[] = {
     ExactInf::Name,
@@ -100,6 +111,9 @@ static const char* DAINames[] = {
 #endif
 #ifdef DAI_WITH_GIBBS
     Gibbs::Name,
+#endif
+#ifdef DAI_WITH_CBP
+    CBP::Name,
 #endif
     ""
 };
diff --git a/include/dai/bbp.h b/include/dai/bbp.h
new file mode 100644 (file)
index 0000000..2b41568
--- /dev/null
@@ -0,0 +1,354 @@
+/*  Copyright (C) 2009  Frederik Eaton [frederik at ofb dot net]
+
+    This file is part of libDAI.
+
+    libDAI is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    libDAI is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with libDAI; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+
+/// \file
+/// \brief Defines class BBP [\ref EaG09]
+/// \todo Improve documentation
+/// \todo Debug clean_updates
+
+
+#ifndef ___defined_libdai_bbp_h
+#define ___defined_libdai_bbp_h
+
+
+#include <vector>
+#include <utility>
+
+#include <dai/prob.h>
+#include <dai/daialg.h>
+#include <dai/factorgraph.h>
+#include <dai/enum.h>
+#include <dai/bp_dual.h>
+
+
+namespace dai {
+
+
+/// Computes the adjoint of the unnormed probability vector from the normalizer and the adjoint of the normalized probability vector @see eqn. (13) in [\ref EaG09]
+Prob unnormAdjoint( const Prob &w, Real Z_w, const Prob &adj_w );
+
+/// Runs Gibbs sampling for \a iters iterations on ia.fg(), and returns state
+std::vector<size_t> getGibbsState( const InfAlg &ia, size_t iters );
+
+
+/// Implements BBP (Back-Belief-Propagation) [\ref EaG09]
+class BBP {
+    protected:
+        /// @name Inputs
+        //@{
+        BP_dual _bp_dual;
+        const FactorGraph *_fg;
+        const InfAlg *_ia;
+        //@}
+
+        /// Number of iterations done
+        size_t _iters;
+
+        /// @name Outputs
+        //@{
+        /// Variable factor adjoints
+        std::vector<Prob> _adj_psi_V;
+        /// Factor adjoints
+        std::vector<Prob> _adj_psi_F;
+        /// Variable->factor message adjoints (indexed [i][_I])
+        std::vector<std::vector<Prob> > _adj_n;
+        /// Factor->variable message adjoints (indexed [i][_I])
+        std::vector<std::vector<Prob> > _adj_m;
+        /// Normalized variable belief adjoints
+        std::vector<Prob> _adj_b_V;
+        /// Normalized factor belief adjoints
+        std::vector<Prob> _adj_b_F;
+        //@}
+
+        /// @name Helper quantities computed from the BP messages
+        //@{
+        /// _T[i][_I] (see eqn. (41) in [\ref EaG09])
+        std::vector<std::vector<Prob > > _T;
+        /// _U[I][_i] (see eqn. (42) in [\ref EaG09])
+        std::vector<std::vector<Prob > > _U;
+        /// _S[i][_I][_j] (see eqn. (43) in [\ref EaG09])
+        std::vector<std::vector<std::vector<Prob > > > _S;
+        /// _R[I][_i][_J] (see eqn. (44) in [\ref EaG09])
+        std::vector<std::vector<std::vector<Prob > > > _R;
+        //@}
+
+        /// Unnormalized variable belief adjoints
+        std::vector<Prob> _adj_b_V_unnorm;
+        /// Unnormalized factor belief adjoints
+        std::vector<Prob> _adj_b_F_unnorm;
+
+        /// Initial variable factor adjoints
+        std::vector<Prob> _init_adj_psi_V;
+        /// Initial factor adjoints
+        std::vector<Prob> _init_adj_psi_F;
+
+        /// Unnormalized variable->factor message adjoint (indexed [i][_I])
+        std::vector<std::vector<Prob> > _adj_n_unnorm;
+        /// Unnormalized factor->variable message adjoint (indexed [i][_I])
+        std::vector<std::vector<Prob> > _adj_m_unnorm;
+        /// Updated normalized variable->factor message adjoint (indexed [i][_I])
+        std::vector<std::vector<Prob> > _new_adj_n;
+        /// Updated normalized factor->variable message adjoint (indexed [i][_I])
+        std::vector<std::vector<Prob> > _new_adj_m;
+
+        /// @name Optimized indexing (for performance)
+        //@{
+        /// Calculates _indices, which is a cache of IndexFor @see bp.cpp
+        void RegenerateInds();
+        
+        /// Index type
+        typedef std::vector<size_t>  _ind_t;
+        /// Cached indices (indexed [i][_I])
+        std::vector<std::vector<_ind_t> >  _indices; 
+        /// Returns an index from the cache
+        const _ind_t& _index(size_t i, size_t _I) const { return _indices[i][_I]; }
+        //@}
+
+        /// @name Initialization
+        //@{
+        /// Calculate T values; see eqn. (41) in [\ref EaG09]
+        void RegenerateT();
+        /// Calculate U values; see eqn. (42) in [\ref EaG09]
+        void RegenerateU();
+        /// Calculate S values; see eqn. (43) in [\ref EaG09]
+        void RegenerateS();
+        /// Calculate R values; see eqn. (44) in [\ref EaG09]
+        void RegenerateR();
+        /// Calculate _adj_b_V_unnorm and _adj_b_F_unnorm from _adj_b_V and _adj_b_F
+        void RegenerateInputs();
+        /// Initialise members for factor adjoints (call after RegenerateInputs)
+        void RegeneratePsiAdjoints();
+        /// Initialise members for message adjoints (call after RegenerateInputs) for parallel algorithm
+        void RegenerateParMessageAdjoints();
+        /// Initialise members for message adjoints (call after RegenerateInputs) for sequential algorithm
+        /** Same as RegenerateMessageAdjoints, but calls sendSeqMsgN rather
+         *  than updating _adj_n (and friends) which are unused in the sequential algorithm.
+         */
+        void RegenerateSeqMessageAdjoints();
+        //@}
+
+        /// Returns T value; see eqn. (41) in [\ref EaG09]
+        DAI_ACCMUT(Prob & T(size_t i, size_t _I), { return _T[i][_I]; });
+        /// Retunrs U value; see eqn. (42) in [\ref EaG09]
+        DAI_ACCMUT(Prob & U(size_t I, size_t _i), { return _U[I][_i]; });
+        /// Returns S value; see eqn. (43) in [\ref EaG09]
+        DAI_ACCMUT(Prob & S(size_t i, size_t _I, size_t _j), { return _S[i][_I][_j]; });
+        /// Returns R value; see eqn. (44) in [\ref EaG09]
+        DAI_ACCMUT(Prob & R(size_t I, size_t _i, size_t _J), { return _R[I][_i][_J]; });
+
+        /// @name Parallel algorithm
+        //@{
+        /// Calculates new variable->factor message adjoint
+        /** Increases variable factor adjoint according to eqn. (27) in [\ref EaG09] and
+         *  calculates the new variable->factor message adjoint according to eqn. (29) in [\ref EaG09].
+         */
+        void calcNewN( size_t i, size_t _I );
+        /// Calculates new factor->variable message adjoint
+        /** Increases factor adjoint according to eqn. (28) in [\ref EaG09] and
+         *  calculates the new factor->variable message adjoint according to the r.h.s. of eqn. (30) in [\ref EaG09].
+         */
+        void calcNewM( size_t i, size_t _I );
+        /// Calculates unnormalized variable->factor message adjoint from the normalized one
+        void calcUnnormMsgN( size_t i, size_t _I );
+        /// Calculates unnormalized factor->variable message adjoint from the normalized one
+        void calcUnnormMsgM( size_t i, size_t _I );
+        /// Updates (un)normalized variable->factor message adjoints
+        void upMsgN( size_t i, size_t _I );
+        /// Updates (un)normalized factor->variable message adjoints
+        void upMsgM( size_t i, size_t _I );
+        /// Do one parallel update of all message adjoints
+        void doParUpdate();
+        //@}
+
+        /// @name Sequential algorithm
+        //@{
+        /// Helper function for sendSeqMsgM: increases factor->variable message adjoint by p and calculates the corresponding unnormalized adjoint
+        void incrSeqMsgM( size_t i, size_t _I, const Prob& p );
+        void updateSeqMsgM( size_t i, size_t _I );
+        /// Sets normalized factor->variable message adjoint and calculates the corresponding unnormalized adjoint
+        void setSeqMsgM( size_t i, size_t _I, const Prob &p ); 
+        /// Implements routine Send-n in Figure 5 in [\ref EaG09]
+        void sendSeqMsgN( size_t i, size_t _I, const Prob &f );
+        /// Implements routine Send-m in Figure 5 in [\ref EaG09]
+        void sendSeqMsgM( size_t i, size_t _I );
+        //@}
+
+        /// Calculates averaged L-1 norm of unnormalized message adjoints
+        Real getUnMsgMag();
+        /// Calculates averaged L-1 norms of current and new normalized message adjoints
+        void getMsgMags( Real &s, Real &new_s );
+
+        /// Sets all vectors _adj_b_F to zero
+        void zero_adj_b_F() {
+            _adj_b_F.clear();
+            _adj_b_F.reserve( _fg->nrFactors() );
+            for( size_t I = 0; I < _fg->nrFactors(); I++ )
+                _adj_b_F.push_back( Prob( _fg->factor(I).states(), Real( 0.0 ) ) );
+        }
+
+        /// Returns indices and magnitude of the largest normalized factor->variable message adjoint
+        void getArgmaxMsgM( size_t &i, size_t &_I, Real &mag );
+        /// Returns magnitude of the largest (in L1-norm) normalized factor->variable message adjoint
+        Real getMaxMsgM();
+        /// Calculates sum of L1 norms of all normalized factor->variable message adjoints
+        Real getTotalMsgM();
+        /// Calculates sum of L1 norms of all updated normalized factor->variable message adjoints
+        Real getTotalNewMsgM();
+        /// Calculates sum of L1 norms of all normalized variable->factor message adjoints
+        Real getTotalMsgN();
+
+    public:
+        /// Called by \a init, recalculates intermediate values
+        void Regenerate();
+
+        /// Constructor
+        BBP( const InfAlg *ia, const PropertySet &opts ) : _bp_dual(ia), _fg(&(ia->fg())), _ia(ia) {
+            props.set(opts);
+        }
+
+        /// Returns a vector of Probs (filled with zeroes) with state spaces corresponding to the factors in the factor graph fg
+        std::vector<Prob> getZeroAdjF( const FactorGraph &fg );
+        /// Returns a vector of Probs (filled with zeroes) with state spaces corresponding to the variables in the factor graph fg
+        std::vector<Prob> getZeroAdjV( const FactorGraph &fg );
+
+        /// Initializes belief adjoints and initial factor adjoints and regenerates
+        void init( const std::vector<Prob> &adj_b_V, const std::vector<Prob> &adj_b_F, const std::vector<Prob> &adj_psi_V, const std::vector<Prob> &adj_psi_F ) {
+            _adj_b_V = adj_b_V;
+            _adj_b_F = adj_b_F;
+            _init_adj_psi_V = adj_psi_V;
+            _init_adj_psi_F = adj_psi_F;
+            Regenerate(); 
+        }
+
+        /// Initializes belief adjoints and with zero initial factor adjoints and regenerates
+        void init( const std::vector<Prob> &adj_b_V, const std::vector<Prob> &adj_b_F ) {
+            init( adj_b_V, adj_b_F, getZeroAdjV(*_fg), getZeroAdjF(*_fg) );
+        }
+
+        /// Initializes variable belief adjoints (and sets factor belief adjoints to zero) and with zero initial factor adjoints and regenerates
+        void init( const std::vector<Prob> &adj_b_V ) {
+            init(adj_b_V, getZeroAdjF(*_fg));
+        }
+
+        /// Run until change is less than given tolerance
+        void run();
+
+        /// Return number of iterations done so far
+        size_t doneIters() { return _iters; }
+
+        /// Returns variable factor adjoint
+        DAI_ACCMUT(Prob& adj_psi_V(size_t i), { return _adj_psi_V[i]; });
+        /// Returns factor adjoint
+        DAI_ACCMUT(Prob& adj_psi_F(size_t I), { return _adj_psi_F[I]; });
+        /// Returns variable belief adjoint
+        DAI_ACCMUT(Prob& adj_b_V(size_t i), { return _adj_b_V[i]; });
+        /// Returns factor belief adjoint
+        DAI_ACCMUT(Prob& adj_b_F(size_t I), { return _adj_b_F[I]; });
+
+     protected:
+        /// Returns variable->factor message adjoint
+        DAI_ACCMUT(Prob& adj_n(size_t i, size_t _I), { return _adj_n[i][_I]; });
+        /// Returns factor->variable message adjoint
+        DAI_ACCMUT(Prob& adj_m(size_t i, size_t _I), { return _adj_m[i][_I]; });
+
+     public: 
+        /// Parameters of this algorithm
+        /* PROPERTIES(props,BBP) {
+           /// Enumeration of possible update schedules
+           DAI_ENUM(UpdateType,SEQ_FIX,SEQ_MAX,SEQ_BP_REV,SEQ_BP_FWD,PAR);
+
+           /// Verbosity
+           size_t verbose;
+
+           /// Maximum number of iterations
+           size_t maxiter;
+
+           /// Tolerance (not used for updates = SEQ_BP_REV, SEQ_BP_FWD)
+           double tol;
+
+           /// Damping constant (0 for none); damping = 1 - lambda where lambda is the damping constant used in [\ref EaG09]
+           double damping;
+
+           /// Update schedule
+           UpdateType updates;
+
+           bool clean_updates;
+        } 
+        */
+/* {{{ GENERATED CODE: DO NOT EDIT. Created by 
+    ./scripts/regenerate-properties include/dai/bbp.h src/bbp.cpp 
+*/
+        struct Properties {
+            /// Enumeration of possible update schedules
+            DAI_ENUM(UpdateType,SEQ_FIX,SEQ_MAX,SEQ_BP_REV,SEQ_BP_FWD,PAR);
+            /// Verbosity
+            size_t verbose;
+            /// Maximum number of iterations
+            size_t maxiter;
+            /// Tolerance (not used for updates = SEQ_BP_REV, SEQ_BP_FWD)
+            double tol;
+            /// Damping constant (0 for none); damping = 1 - lambda where lambda is the damping constant used in [\ref EaG09]
+            double damping;
+            /// Update schedule
+            UpdateType updates;
+            bool clean_updates;
+
+            /// Set members from PropertySet
+            void set(const PropertySet &opts);
+            /// Get members into PropertySet
+            PropertySet get() const;
+            /// Convert to a string which can be parsed as a PropertySet
+            std::string toString() const;
+        } props;
+/* }}} END OF GENERATED CODE */
+};
+
+
+/// Enumeration of several cost functions that can be used with BBP.
+DAI_ENUM(bbp_cfn_t,CFN_GIBBS_B,CFN_GIBBS_B2,CFN_GIBBS_EXP,CFN_GIBBS_B_FACTOR,CFN_GIBBS_B2_FACTOR,CFN_GIBBS_EXP_FACTOR,CFN_VAR_ENT,CFN_FACTOR_ENT,CFN_BETHE_ENT);
+
+/// Initialise BBP using InfAlg, cost function, and stateP
+/** Calls bbp.init with adjoints calculated from ia.beliefV and
+ *  ia.beliefF. stateP is a Gibbs state and can be NULL, it will be
+ *  initialised using a Gibbs run of 2*fg.Iterations() iterations.
+ */
+void initBBPCostFnAdj( BBP &bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const std::vector<size_t> *stateP );
+
+/// Answers question: does the given cost function depend on having a Gibbs state?
+bool needGibbsState( bbp_cfn_t cfn );
+
+/// Calculate actual value of cost function (cfn_type, stateP)
+/** This function returns the actual value of the cost function whose
+ *  gradient with respect to singleton beliefs is given by
+ *  gibbsToB1Adj on the same arguments
+ */
+Real getCostFn( const InfAlg &fg, bbp_cfn_t cfn_type, const std::vector<size_t> *stateP );
+
+/// Function to test the validity of adjoints computed by BBP given a state for each variable using numerical derivatives.
+/** Factors containing a variable are multiplied by psi_1 adjustments to verify accuracy of _adj_psi_V.
+ *  \param h controls size of perturbation.
+ */
+double numericBBPTest( const InfAlg &bp, const std::vector<size_t> *state, const PropertySet &bbp_props, bbp_cfn_t cfn, double h );
+
+
+} // end of namespace dai
+
+
+#endif
index 17b39a7..15845de 100644 (file)
@@ -56,14 +56,59 @@ namespace dai {
  */
 class BipartiteGraph {
     public:
-        /// Describes a neighboring node of some other node in a BipartiteGraph.
-        /** A Neighbor structure has three members: \a iter, \a node and \a dual. The \a 
-         *  node member is the most important member: it contains the index of the neighboring node. The \a iter 
-         *  member is useful for iterating over neighbors, and contains the index of this Neighbor entry in the
-         *  corresponding BipartiteGraph::Neighbors variable. The \a dual member is useful to find the dual Neighbor 
-         *  element: a pair of neighboring nodes can be either specified as a node of type 1 and a neighbor of type 
-         *  2, or as a node of type 2 and a neighbor of type 1; the \a dual member contains the index of the dual 
-         *  Neighbor element (see the example for another explanation of the dual member).
+        /// Describes the neighbor relationship of two nodes in a BipartiteGraph.
+        /** Sometimes we want to do an action, such as sending a
+         *  message, for all edges in a graph. However, most graphs
+         *  will be sparse, so we need some way of storing a set of
+         *  the neighbors of a node, which is both fast and
+         *  memory-efficient. We also need to be able to go between
+         *  viewing node \a A as a neighbor of node \a B, and node \a
+         *  B as a neighbor of node \a A. The Neighbor struct solves
+         *  both of these problems. Each node has a list of neighbors,
+         *  stored as a vector<Neighbor>, and extra information is
+         *  included in the Neighbor struct which allows us to access
+         *  a node as a neighbor of its neighbor (the \a dual member).
+         *
+         *  By convention, variable identifiers naming indices into a
+         *  vector of neighbors are prefixed with an underscore ("_"). 
+         *  The neighbor list which they point into is then understood
+         *  from context. For example:
+         *
+         *  \code
+         *  void BP::calcNewMessage( size_t i, size_t _I )
+         *  \endcode
+         *
+         *  Here, \a i is the "absolute" index of node i, but \a _I is
+         *  understood as a "relative" index, giving node I's entry in
+         *  nb1(i). The corresponding Neighbor structure can be
+         *  accessed as nb1(i,_I) or nb1(i)[_I]. The absolute index of
+         *  \a _I, which would be called \a I, can be recovered from
+         *  the \a node member: nb1(i,_I).node. The \a iter member
+         *  gives the relative index \a _I, and the \a dual member
+         *  gives the "dual" relative index, i.e. the index of \a i in
+         *  \a I's neighbor list.
+         *
+         *  \code
+         *  Neighbor n = nb1(i,_I);
+         *  n.node == I &&
+         *  n.iter == _I &&
+         *  nb2(n.node,n.dual).node == i
+         *  \endcode
+         *
+         *  In a FactorGraph, nb1 is called nbV, and nb2 is called
+         *  nbF.
+         *
+         *  There is no easy way to transform a pair of absolute node
+         *  indices \a i and \a I into a Neighbor structure relative
+         *  to one of the nodes. Such a feature has never yet been
+         *  found to be necessary. Iteration over edges can always be
+         *  accomplished using the Neighbor lists, and by writing
+         *  functions that accept relative indices:
+         *  \code
+         *  for( size_t i = 0; i < nrVars(); ++i )
+         *      foreach( const Neighbor &I, nbV(i) )
+         *          calcNewMessage( i, I.iter );
+         *  \endcode
          */
         struct Neighbor {
             /// Corresponds to the index of this Neighbor entry in the vector of neighbors
@@ -101,9 +146,19 @@ class BipartiteGraph {
             std::vector<size_t> ind2;       // indices of nodes of type 2
         };
 
+        /// @name Backwards compatibility layer (to be removed soon)
+        //@{
+        /// Enable backwards compatibility layer?
+        bool _edge_indexed;
+        /// Call indexEdges() first to initialize these members
+        std::vector<Edge> _edges;
+        /// Call indexEdges() first to initialize these members
+        hash_map<Edge,size_t> _vv2e;
+        //}@
+
     public:
         /// Default constructor (creates an empty bipartite graph)
-        BipartiteGraph() : _nb1(), _nb2() {}
+        BipartiteGraph() : _nb1(), _nb2(), _edge_indexed(false) {}
 
         /// Constructs BipartiteGraph from a range of edges.
         /** \tparam EdgeInputIterator Iterator that iterates over instances of BipartiteGraph::Edge.
@@ -113,7 +168,7 @@ class BipartiteGraph {
          *  \param end Points just beyond the last edge.
          */
         template<typename EdgeInputIterator>
-        BipartiteGraph( size_t nr1, size_t nr2, EdgeInputIterator begin, EdgeInputIterator end ) : _nb1( nr1 ), _nb2( nr2 ) {
+        BipartiteGraph( size_t nr1, size_t nr2, EdgeInputIterator begin, EdgeInputIterator end ) : _nb1( nr1 ), _nb2( nr2 ), _edge_indexed(false) {
             construct( nr1, nr2, begin, end );
         }
 
@@ -318,6 +373,53 @@ class BipartiteGraph {
         /// Writes this BipartiteGraph to an output stream in GraphViz .dot syntax
         void printDot( std::ostream& os ) const;
 
+        /// @name Backwards compatibility layer (to be removed soon)
+        //@{
+        void indexEdges() {
+            std::cerr << "Warning: this BipartiteGraph edge interface is obsolete!" << std::endl;
+            _edges.clear();
+            _vv2e.clear();
+            size_t i=0;
+            foreach(const Neighbors &nb1s, _nb1) {
+                foreach(const Neighbor &n2, nb1s) {
+                    Edge e(i, n2.node);
+                    _edges.push_back(e);
+                }
+                i++;
+            }
+            sort(_edges.begin(), _edges.end()); // unnecessary?
+
+            i=0;
+            foreach(const Edge& e, _edges) {
+                _vv2e[e] = i++;
+            }
+
+            _edge_indexed = true;
+        }
+        
+        const Edge& edge(size_t e) const {
+            assert(_edge_indexed);
+            return _edges[e];
+        }
+        
+        const std::vector<Edge>& edges() const {
+            return _edges;
+        }
+        
+        size_t VV2E(size_t n1, size_t n2) const {
+            assert(_edge_indexed);
+            Edge e(n1,n2);
+            hash_map<Edge,size_t>::const_iterator i = _vv2e.find(e);
+            assert(i != _vv2e.end());
+            return i->second;
+        }
+
+        size_t nr_edges() const {
+            assert(_edge_indexed);
+            return _edges.size();
+        }
+        //}@
+
     private:
         /// Checks internal consistency
         void check() const;
index 1b1e577..9117f9b 100644 (file)
@@ -57,15 +57,17 @@ class BP : public DAIAlgFG {
         double _maxdiff;
         /// Number of iterations needed
         size_t _iters;
+        /// The history of message updates (only recorded if recordSentMessages is true)
+        std::vector<std::pair<std::size_t, std::size_t> > _sentMessages;
     
     public:
         /// Parameters of this inference algorithm
         struct Properties {
             /// Enumeration of possible update schedules
-            DAI_ENUM(UpdateType,SEQFIX,SEQRND,SEQMAX,PARALL)
+            DAI_ENUM(UpdateType,SEQFIX,SEQRND,SEQMAX,PARALL);
 
             /// Enumeration of inference variants
-            DAI_ENUM(InfType,SUMPROD,MAXPROD)
+            DAI_ENUM(InfType,SUMPROD,MAXPROD);
         
             /// Verbosity
             size_t verbose;
@@ -92,12 +94,18 @@ class BP : public DAIAlgFG {
         /// Name of this inference algorithm
         static const char *Name;
 
+        /// Specifies whether the history of message updates should be recorded
+        bool recordSentMessages;
+
     public:
         /// Default constructor
-        BP() : DAIAlgFG(), _edges(), _edge2lut(), _lut(), _maxdiff(0.0), _iters(0U), props() {}
+        BP() : DAIAlgFG(), _edges(), _edge2lut(), _lut(), _maxdiff(0.0), _iters(0U), _sentMessages(), props(), recordSentMessages(false) {}
 
         /// Copy constructor
-        BP( const BP &x ) : DAIAlgFG(x), _edges(x._edges), _edge2lut(x._edge2lut), _lut(x._lut), _maxdiff(x._maxdiff), _iters(x._iters), props(x.props) {
+        BP( const BP &x ) : DAIAlgFG(x), _edges(x._edges), _edge2lut(x._edge2lut),
+            _lut(x._lut), _maxdiff(x._maxdiff), _iters(x._iters), _sentMessages(x._sentMessages), 
+            props(x.props), recordSentMessages(x.recordSentMessages) 
+        {
             for( LutType::iterator l = _lut.begin(); l != _lut.end(); ++l )
                 _edge2lut[l->second.first][l->second.second] = l;
         }
@@ -112,13 +120,15 @@ class BP : public DAIAlgFG {
                     _edge2lut[l->second.first][l->second.second] = l;
                 _maxdiff = x._maxdiff;
                 _iters = x._iters;
+                _sentMessages = x._sentMessages;
                 props = x.props;
+                recordSentMessages = x.recordSentMessages;
             }
             return *this;
         }
 
         /// Construct from FactorGraph fg and PropertySet opts
-        BP( const FactorGraph & fg, const PropertySet &opts ) : DAIAlgFG(fg), _edges(), _maxdiff(0.0), _iters(0U), props() {
+        BP( const FactorGraph & fg, const PropertySet &opts ) : DAIAlgFG(fg), _edges(), _maxdiff(0.0), _iters(0U), _sentMessages(), props(), recordSentMessages(false) {
             setProperties( opts );
             construct();
         }
@@ -151,6 +161,16 @@ class BP : public DAIAlgFG {
          */
         std::vector<std::size_t> findMaximum() const;
 
+        /// Returns history of sent messages
+        const std::vector<std::pair<std::size_t, std::size_t> >& getSentMessages() const {
+            return _sentMessages;
+        }
+
+        /// Clears history of sent messages
+        void clearSentMessages() {
+            _sentMessages.clear();
+        }
+
     private:
         const Prob & message(size_t i, size_t _I) const { return _edges[i][_I].message; }
         Prob & message(size_t i, size_t _I) { return _edges[i][_I].message; }
diff --git a/include/dai/bp_dual.h b/include/dai/bp_dual.h
new file mode 100644 (file)
index 0000000..0164635
--- /dev/null
@@ -0,0 +1,137 @@
+/*  Copyright (C) 2009  Frederik Eaton [frederik at ofb dot net]
+
+    This file is part of libDAI.
+
+    libDAI is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    libDAI is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with libDAI; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+
+/// \file
+/// \brief Defines class BP_dual
+/// \todo This replicates a large part of the functionality of BP; would it not be shorter to adapt BP instead?
+
+
+#ifndef __defined_libdai_bp_dual_h
+#define __defined_libdai_bp_dual_h
+
+
+#include <dai/daialg.h>
+#include <dai/factorgraph.h>
+#include <dai/enum.h>
+
+
+namespace dai {
+
+
+/// Calculates both types of BP messages and their normalizers from an InfAlg.
+/** BP_dual calculates "dual" versions of BP messages (both messages from factors
+ *  to variables and messages from variables to factors), and normalizers, given an InfAlg. 
+ *  These are computed from the variable and factor beliefs of the InfAlg.
+ *  This class is used primarily by BBP.
+ */
+class BP_dual {
+
+    protected:
+        /// Convenience label for storing edge properties
+        template<class T>
+        struct _edges_t : public std::vector<std::vector<T> > {};
+
+        /// Groups together the data structures for storing the two types of messages and their normalizers
+        struct messages {            
+            /// Unnormalized variable->factor messages
+            _edges_t<Prob> n;
+            /// Normalizers of variable->factor messages
+            _edges_t<Real> Zn;
+            /// Unnormalized Factor->variable messages
+            _edges_t<Prob> m;
+            /// Normalizers of factor->variable messages
+            _edges_t<Real> Zm;
+        };
+        /// Stores all messages
+        messages _msgs;
+
+        /// Groups together the data structures for storing the two types of beliefs and their normalizers
+        struct beliefs {
+            /// Unnormalized variable beliefs
+            std::vector<Prob> b1;
+            /// Normalizers of variable beliefs
+            std::vector<Real> Zb1;
+            /// Unnormalized factor beliefs
+            std::vector<Prob> b2;
+            /// Normalizers of factor beliefs
+            std::vector<Real> Zb2;
+        };
+        /// Stores all beliefs
+        beliefs _beliefs;
+
+        /// Pointer to the InfAlg object
+        const InfAlg *_ia;
+        
+        /// Does all necessary preprocessing
+        void init();
+        /// Allocates space for _msgs
+        void regenerateMessages();
+        /// Allocates space for _beliefs
+        void regenerateBeliefs();
+
+        /// Calculate all messages from InfAlg beliefs
+        void calcMessages();
+        /// Update factor->variable message (i->I)
+        void calcNewM(size_t i, size_t _I);
+        /// Update variable->factor message (I->i)
+        void calcNewN(size_t i, size_t _I);
+
+        /// Calculate all variable and factor beliefs from messages
+        void calcBeliefs();
+        /// Calculate variable belief
+        void calcBeliefV(size_t i);
+        /// Calculate factor belief
+        void calcBeliefF(size_t I);
+
+    public:
+        /// Construct BP_dual object from (converged) InfAlg object's beliefs and factors. 
+        /*  A pointer to the the InfAlg object is stored, 
+         *  so the object must not be destroyed before the BP_dual is destroyed.
+         */
+        BP_dual( const InfAlg *ia ) : _ia(ia) { init(); }
+
+        /// Returns the underlying FactorGraph
+        const FactorGraph& fg() const { return _ia->fg(); }
+
+        /// Returns factor -> var message (I->i)
+        DAI_ACCMUT(Prob & msgM( size_t i, size_t _I ), { return _msgs.m[i][_I]; });
+        /// Returns var -> factor message (i->I)
+        DAI_ACCMUT(Prob & msgN( size_t i, size_t _I ), { return _msgs.n[i][_I]; });
+        /// Returns normalizer for msgM
+        DAI_ACCMUT(Real & zM( size_t i, size_t _I ), { return _msgs.Zm[i][_I]; });
+        /// Returns normalizer for msgN
+        DAI_ACCMUT(Real & zN( size_t i, size_t _I ), { return _msgs.Zn[i][_I]; });
+
+        /// Returns variable belief
+        Factor beliefV( size_t i ) const { return Factor( _ia->fg().var(i), _beliefs.b1[i] ); }
+        /// Returns factor belief
+        Factor beliefF( size_t I ) const { return Factor( _ia->fg().factor(I).vars(), _beliefs.b2[I] ); }
+
+        /// Returns normalizer for variable belief
+        Real beliefVZ( size_t i ) const { return _beliefs.Zb1[i]; }
+        /// Returns normalizer for factor belief
+        Real beliefFZ( size_t I ) const { return _beliefs.Zb2[I]; }
+};
+
+
+} // end of namespace dai
+
+
+#endif
diff --git a/include/dai/cbp.h b/include/dai/cbp.h
new file mode 100644 (file)
index 0000000..2a562e1
--- /dev/null
@@ -0,0 +1,269 @@
+/*  Copyright (C) 2009  Frederik Eaton [frederik at ofb dot net]
+
+    This file is part of libDAI.
+
+    libDAI is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    libDAI is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with libDAI; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+
+/// \file
+/// \brief Defines class CBP [\ref EaG09]
+/// \todo Improve documentation
+
+
+#ifndef __defined_libdai_cbp_h
+#define __defined_libdai_cbp_h
+
+
+#include <fstream>
+#include <boost/shared_ptr.hpp>
+
+#include <dai/daialg.h>
+#include <dai/cbp.h>
+#include <dai/bbp.h>
+
+
+namespace dai {
+
+
+/// Find a variable to clamp using BBP (goes with maximum adjoint) 
+/// \see BBP
+std::pair<size_t, size_t> bbpFindClampVar( const InfAlg &in_bp, bool clampingVar, const PropertySet &bbp_props, bbp_cfn_t cfn, Real *maxVarOut );
+
+
+/// Class for CBP (Clamped Belief Propagation)
+/** This algorithm uses configurable heuristics to choose a variable
+ *  x_i and a state x_i*. Inference is done with x_i "clamped" to x_i*
+ *  (i.e., conditional on x_i == x_i*), and also with the negation of this
+ *  condition. Clamping is done recursively up to a fixed number of
+ *  levels (other stopping criteria are also implemented, see
+ *  \a recursion property). The resulting approximate marginals are
+ *  combined using logZ estimates.
+ */
+class CBP : public DAIAlgFG {
+    private:
+        /// Variable beliefs
+        std::vector<Factor> _beliefsV;
+        /// Factor beliefs
+        std::vector<Factor> _beliefsF;
+        /// Log-partition sum
+        double _logZ;
+
+        /// Counts number of clampings at each leaf node
+        double _sum_level;
+
+        /// Number of leaves of recursion tree
+        size_t _num_leaves;
+
+        /// Output stream where information about the clampings is written
+        boost::shared_ptr<std::ofstream> _clamp_ofstream;
+
+        /// Returns BBP cost function used
+        bbp_cfn_t BBP_cost_function() { return props.bbp_cfn; }
+
+        /// Prints beliefs, variables and partition sum, in case of a debugging build
+        void printDebugInfo();
+
+        /// Called by 'run', and by itself. Implements the main algorithm. 
+        /** Chooses a variable to clamp, recurses, combines the logZ and
+         *  beliefs estimates of the children, and returns the improved
+         *  estimates in \a lz_out and \a beliefs_out to its parent
+         */
+        void runRecurse( InfAlg *bp, double orig_logZ, std::vector<size_t> clamped_vars_list, size_t &num_leaves,
+                         size_t &choose_count, double &sum_level, Real &lz_out, std::vector<Factor> &beliefs_out );
+
+        /// Choose the next variable to clamp
+        /** Choose the next variable to clamp, given a converged InfAlg (\a bp),
+         *  and a vector of variables that are already clamped (\a
+         *  clamped_vars_list). Returns the chosen variable in \a i, and
+         *  the set of states in \a xis. If \a maxVarOut is non-NULL and
+         *  props.choose==CHOOSE_BBP then it is used to store the
+         *  adjoint of the chosen variable
+         */
+        virtual bool chooseNextClampVar( InfAlg* bp, std::vector<size_t> &clamped_vars_list, size_t &i, std::vector<size_t> &xis, Real *maxVarOut );
+
+        /// Return the InfAlg to use at each step of the recursion. 
+        /// \todo At present, only returns a BP instance
+        InfAlg* getInfAlg();
+
+        /// Numer of iterations needed
+        size_t _iters;
+        /// Maximum difference encountered so far
+        double _maxdiff;
+
+        /// Sets variable beliefs, factor beliefs and logZ
+        /** \param bs should be a concatenation of the variable beliefs followed by the factor beliefs
+         */
+        void setBeliefs( const std::vector<Factor> &bs, double logZ );
+
+        /// Constructor helper function
+        void construct();
+
+    public:
+        /// Construct CBP object from FactorGraph fg and PropertySet opts
+        CBP( const FactorGraph &fg, const PropertySet &opts ) : DAIAlgFG(fg) {
+            props.set( opts );
+            construct();
+        }
+
+        /// Name of this inference algorithm
+        static const char *Name;
+
+        /// @name General InfAlg interface
+        //@{
+        virtual CBP* clone() const { return new CBP(*this); }
+        virtual std::string identify() const { return std::string(Name) + props.toString(); }
+        virtual Factor belief (const Var &n) const { return _beliefsV[findVar(n)]; }
+        virtual Factor belief (const VarSet &) const { DAI_THROW(NOT_IMPLEMENTED); }
+        virtual std::vector<Factor> beliefs() const { return concat(_beliefsV, _beliefsF); }
+        virtual Real logZ() const { return _logZ; }
+        virtual void init() {};
+        virtual void init( const VarSet & ) {};
+        virtual double run();
+        virtual double maxDiff() const { return _maxdiff; }
+        virtual size_t Iterations() const { return _iters; }
+        //@}
+
+        Factor beliefV( size_t i ) const { return _beliefsV[i]; }
+        Factor beliefF( size_t I ) const { return _beliefsF[I]; }
+
+        //----------------------------------------------------------------
+
+        /// Parameters of this inference algorithm
+        /* PROPERTIES(props,CBP) {
+            /// Enumeration of possible update schedules
+            typedef BP::Properties::UpdateType UpdateType;
+            /// Enumeration of possible methods for deciding when to stop recursing
+            DAI_ENUM(RecurseType,REC_FIXED,REC_LOGZ,REC_BDIFF);
+            /// Enumeration of possible heuristics for choosing clamping variable
+            DAI_ENUM(ChooseMethodType,CHOOSE_RANDOM,CHOOSE_MAXENT,CHOOSE_BBP,CHOOSE_BP_L1,CHOOSE_BP_CFN);
+            /// Enumeration of possible clampings: variables or factors
+            DAI_ENUM(ClampType,CLAMP_VAR,CLAMP_FACTOR);
+            
+            /// Verbosity
+            size_t verbose = 0;
+
+            /// Tolerance to use in BP
+            double tol;
+            /// Update style for BP
+            UpdateType updates;
+            /// Maximum number of iterations for BP
+            size_t maxiter;
+
+            /// Tolerance to use for controlling recursion depth (\a recurse is REC_LOGZ or REC_BDIFF)
+            double rec_tol;
+            /// Maximum number of levels of recursion (\a recurse is REC_FIXED)
+            size_t max_levels = 10;
+            /// If choose==CHOOSE_BBP and maximum adjoint is less than this value, don't recurse
+            double min_max_adj;
+            /// Heuristic for choosing clamping variable
+            ChooseMethodType choose;
+            /// Method for deciding when to stop recursing
+            RecurseType recursion;
+            /// Whether to clamp variables or factors
+            ClampType clamp;
+            /// Properties to pass to BBP
+            PropertySet bbp_props;
+            /// Cost function to use for BBP
+            bbp_cfn_t bbp_cfn;
+            /// Random seed
+            size_t rand_seed = 0;
+
+            /// If non-empty, write clamping choices to this file
+            std::string clamp_outfile = "";
+        }
+        */
+/* {{{ GENERATED CODE: DO NOT EDIT. Created by 
+    ./scripts/regenerate-properties include/dai/cbp.h src/cbp.cpp 
+*/
+        struct Properties {
+            /// Enumeration of possible update schedules
+            typedef BP::Properties::UpdateType UpdateType;
+            /// Enumeration of possible methods for deciding when to stop recursing
+            DAI_ENUM(RecurseType,REC_FIXED,REC_LOGZ,REC_BDIFF);
+            /// Enumeration of possible heuristics for choosing clamping variable
+            DAI_ENUM(ChooseMethodType,CHOOSE_RANDOM,CHOOSE_MAXENT,CHOOSE_BBP,CHOOSE_BP_L1,CHOOSE_BP_CFN);
+            /// Enumeration of possible clampings: variables or factors
+            DAI_ENUM(ClampType,CLAMP_VAR,CLAMP_FACTOR);
+            /// Verbosity
+            size_t verbose;
+            /// Tolerance to use in BP
+            double tol;
+            /// Update style for BP
+            UpdateType updates;
+            /// Maximum number of iterations for BP
+            size_t maxiter;
+            /// Tolerance to use for controlling recursion depth (\a recurse is REC_LOGZ or REC_BDIFF)
+            double rec_tol;
+            /// Maximum number of levels of recursion (\a recurse is REC_FIXED)
+            size_t max_levels;
+            /// If choose=CHOOSE_BBP and maximum adjoint is less than this value, don't recurse
+            double min_max_adj;
+            /// Heuristic for choosing clamping variable
+            ChooseMethodType choose;
+            /// Method for deciding when to stop recursing
+            RecurseType recursion;
+            /// Whether to clamp variables or factors
+            ClampType clamp;
+            /// Properties to pass to BBP
+            PropertySet bbp_props;
+            /// Cost function to use for BBP
+            bbp_cfn_t bbp_cfn;
+            /// Random seed
+            size_t rand_seed;
+            /// If non-empty, write clamping choices to this file
+            std::string clamp_outfile;
+
+            /// Set members from PropertySet
+            void set(const PropertySet &opts);
+            /// Get members into PropertySet
+            PropertySet get() const;
+            /// Convert to a string which can be parsed as a PropertySet
+            std::string toString() const;
+        } props;
+/* }}} END OF GENERATED CODE */
+
+        /// Returns heuristic used for clamping variable
+        Properties::ChooseMethodType ChooseMethod() { return props.choose; }
+        /// Returns method used for deciding when to stop recursing
+        Properties::RecurseType Recursion() { return props.recursion; }
+        /// Returns clamping type used
+        Properties::ClampType Clamping() { return props.clamp; }
+        /// Returns maximum number of levels of recursion
+        size_t maxClampLevel() { return props.max_levels; }
+        /// Returns props.min_max_adj @see CBP::Properties::min_max_adj
+        double minMaxAdj() { return props.min_max_adj; }
+        /// Returns tolerance used for controlling recursion depth
+        double recTol() { return props.rec_tol; }
+};
+
+
+/// Given a sorted vector of states \a xis and total state count \a n_states, return a vector of states not in \a xis
+std::vector<size_t> complement( std::vector<size_t>& xis, size_t n_states );
+
+/// Computes \f$\frac{\exp(a)}{\exp(a)+\exp(b)}\f$
+Real unSoftMax( Real a, Real b );
+
+/// Computes log of sum of exponents, i.e., \f$\log\left(\exp(a) + \exp(b)\right)\f$
+Real logSumExp( Real a, Real b );
+
+/// Compute sum of pairwise L-infinity distances of the first \a nv factors in each vector
+Real dist( const std::vector<Factor>& b1, const std::vector<Factor>& b2, size_t nv );
+
+
+} // end of namespace dai
+
+
+#endif
index 4b1001e..4dbae3a 100644 (file)
@@ -63,6 +63,18 @@ class InfAlg {
         /// Returns the "belief" (i.e., approximate marginal probability distribution) of a set of variables
         virtual Factor belief( const VarSet &n ) const = 0;
 
+        /// Returns marginal for a variable.
+        /** Sometimes preferred to belief() for performance reasons.
+          * Faster implementations exist in e.g. BP.
+          */
+        virtual Factor beliefV( size_t i ) const { return belief( fg().var(i) ); }
+
+        /// Returns marginal for a factor.
+        /** Sometimes preferred to belief() for performance reasons.
+          * Faster implementations exist in e.g. BP.
+          */
+        virtual Factor beliefF( size_t I ) const { return belief( fg().factor(I).vars() ); }
+
         /// Returns all "beliefs" (i.e., approximate marginal probability distribution) calculated by the algorithm
         virtual std::vector<Factor> beliefs() const = 0;
 
index a87c5ff..e7f9071 100644 (file)
  *  - Double-loop GBP [\ref HAK03];
  *  - Various variants of Loop Corrected Belief Propagation
  *    [\ref MoK07, \ref MoR05];
- *  - Gibbs sampler.
+ *  - Gibbs sampler;
+ *  - Conditioned BP [\ref EaG09];
  *
  *  In addition, libDAI supports parameter learning of conditional probability
  *  tables by Expectation Maximization.
  *  "Sufficient Conditions for Convergence of the Sum-Product Algorithm",
  *  <em>IEEE Transactions on Information Theory</em> 53(12):4422-4437.
  *  http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=4385778
+ *
+ *  \anchor EaG09 \ref EaG09
+ *  F. Eaton and Z. Ghahramani (2009):
+ *  "Choosing a Variable to Clamp",
+ *  <em>Proceedings of the Twelfth International Conference on Artificial Intelligence and Statistics (AISTATS 2009)</em> 5:145-152
+ *  http://jmlr.csail.mit.edu/proceedings/papers/v5/eaton09a/eaton09a.pdf
  */
 
 
index 0c351c7..5beae95 100644 (file)
@@ -212,10 +212,21 @@ class FactorGraph {
             }
         }
 
-        /// Clamp variable n to value i (i.e. multiply with a Kronecker delta \f$\delta_{x_n, i}\f$);
-        /// If backup == true, make a backup of all factors that are changed
+        /// Clamp variable n to value i (i.e. multiply with a Kronecker delta \f$\delta_{x_n, i}\f$)
+        /** If backup == true, make a backup of all factors that are changed
+         */
         virtual void clamp( const Var & n, size_t i, bool backup = false );
 
+        /// Clamp a variable in a factor graph to have one out of a list of values
+        /** If backup == true, make a backup of all factors that are changed
+         */
+        void clampVar( size_t i, const std::vector<size_t> &xis, bool backup = false );
+
+        /// Clamp a factor in a factor graph to have one out of a list of values
+        /** If backup == true, make a backup of all factors that are changed
+         */
+        void clampFactor( size_t I, const std::vector<size_t> &xIs, bool backup = false );
+
         /// Set all factors interacting with the i'th variable 1
         virtual void makeCavity( unsigned i, bool backup = false );
 
@@ -273,6 +284,15 @@ class FactorGraph {
         friend std::ostream& operator << (std::ostream& os, const FactorGraph& fg);
         friend std::istream& operator >> (std::istream& is, FactorGraph& fg);
 
+        /// @name Backwards compatibility layer (to be removed soon)
+        //@{
+        size_t VV2E(size_t n1, size_t n2) const { return G.VV2E(n1,n2); }
+        const Edge& edge(size_t e) const { return G.edge(e); }
+        void indexEdges() { G.indexEdges(); }
+        size_t nr_edges() const { return G.nr_edges(); }
+        const std::vector<Edge>& edges() const { return G.edges(); }
+        //}@
+
     private:
         /// Part of constructors (creates edges, neighbors and adjacency matrix)
         void constructGraph( size_t nrEdges );
index d0e508b..536dc8d 100644 (file)
@@ -60,10 +60,10 @@ class LC : public DAIAlgFG {
         /// Parameters of this inference algorithm
         struct Properties {
             /// Enumeration of possible ways to initialize the cavities
-            DAI_ENUM(CavityType,FULL,PAIR,PAIR2,UNIFORM)
+            DAI_ENUM(CavityType,FULL,PAIR,PAIR2,UNIFORM);
 
             /// Enumeration of different update schedules
-            DAI_ENUM(UpdateType,SEQFIX,SEQRND,NONE)
+            DAI_ENUM(UpdateType,SEQFIX,SEQRND,NONE);
 
             /// Verbosity
             size_t verbose;
index 89a9f81..fa6b7ca 100644 (file)
@@ -72,10 +72,10 @@ class MR : public DAIAlgFG {
         /// Parameters of this inference algorithm
         struct Properties {
             /// Enumeration of different types of update equations
-            DAI_ENUM(UpdateType,FULL,LINEAR)
+            DAI_ENUM(UpdateType,FULL,LINEAR);
 
             /// Enumeration of different ways of initializing the cavity correlations
-            DAI_ENUM(InitType,RESPPROP,CLAMPING,EXACT)
+            DAI_ENUM(InitType,RESPPROP,CLAMPING,EXACT);
 
             /// Verbosity
             size_t verbose;
index 196fcd4..6ad5f2a 100644 (file)
@@ -59,6 +59,16 @@ std::ostream& operator<< (std::ostream & os, const Property & p);
 /// Represents a set of properties, mapping keys (of type PropertyKey) to values (of type PropertyValue)
 class PropertySet : private std::map<PropertyKey, PropertyValue> {
     public:
+        /// Default constructor
+        PropertySet() {}
+
+        /// Construct PropertySet from a string
+        PropertySet( const std::string &s ) {
+            std::stringstream ss;
+            ss << s;
+            ss >> *this;
+        }
+
         /// Gets a property
         const PropertyValue & Get(const PropertyKey &key) const { 
             PropertySet::const_iterator x = find(key); 
@@ -114,14 +124,14 @@ class PropertySet : private std::map<PropertyKey, PropertyValue> {
         template<typename ValueType>
         ValueType getStringAs(const PropertyKey &key) const { 
             PropertyValue val = Get(key);
-            if( val.type() == typeid(std::string) ) {
+            if( val.type() == typeid(ValueType) ) {
+                return boost::any_cast<ValueType>(val);
+            } else if( val.type() == typeid(std::string) ) {
                 std::stringstream ss;
                 ss << GetAs<std::string>(key);
                 ValueType result;
                 ss >> result;
                 return result;
-            } else if( val.type() == typeid(ValueType) ) {
-                return boost::any_cast<ValueType>(val);
             } else {
                 DAI_THROW(IMPOSSIBLE_TYPECAST);
                 return ValueType();
@@ -146,6 +156,15 @@ class PropertySet : private std::map<PropertyKey, PropertyValue> {
         /// Check if a property with the given key exists
         bool hasKey(const PropertyKey &key) const { PropertySet::const_iterator x = find(key); return (x != this->end()); }
 
+        /// Returns a set containing all keys
+        std::set<PropertyKey> allKeys() const {
+            std::set<PropertyKey> res;
+            const_iterator i;
+            for( i = begin(); i != end(); i++ )
+                res.insert( i->first );
+            return res;
+        }
+
         // Friends
         friend std::ostream& operator<< (std::ostream & os, const PropertySet & ps);
         friend std::istream& operator>> (std::istream& is, PropertySet & ps);
index ffcc96f..89c1dee 100644 (file)
@@ -35,6 +35,7 @@
 #include <iostream>
 #include <cstdio>
 #include <boost/foreach.hpp>
+#include <boost/functional/hash.hpp>
 #include <algorithm>
 
 
 /// An alias to the BOOST_FOREACH macro from the boost::foreach library
 #define foreach BOOST_FOREACH
 
+#ifdef DAI_DEBUG
+/// \brief "Print variable". Prints the text of an expression, followed by its value (only if DAI_DEBUG is defined)
+/**
+ *  Useful debugging macro to see what your code is doing. 
+ *  Example: \code DAI_PV(3+4) \endcode
+ *  Output: \code 3+4= 7 \endcode
+ */
+#define DAI_PV(x) do {std::cerr << #x "= " << (x) << std::endl;} while(0)
+/// "Debugging message": Prints a message (only if DAI_DEBUG is defined)
+#define DAI_DMSG(str) do {std::cerr << str << std::endl;} while(0)
+#else
+#define DAI_PV(x) do {} while(0)
+#define DAI_DMSG(str) do {} while(0)
+#endif
+
+/// Produces accessor and mutator methods according to the common pattern.
+/** Example:
+ *  \code DAI_ACCMUT(size_t& maxIter(), { return props.maxiter; }); \endcode
+ *  \todo At the moment, only the mutator appears in doxygen documentation.
+ */
+#define DAI_ACCMUT(x,y)                     \
+      x y;                                  \
+      const x const y;
+
+/// Macro to give error message \a stmt if props.verbose>=\a n
+#define DAI_IFVERB(n, stmt) if(props.verbose>=n) { cerr << stmt; }
+
 
 /// Real number (alias for double, which could be changed to long double if necessary)
 typedef double Real;
@@ -77,14 +105,14 @@ namespace dai {
     /// hash_map is an alias for std::map.
     /** Since there is no TR1 unordered_map implementation available yet, we fall back on std::map.
      */
-    template <typename T, typename U>
+    template <typename T, typename U, typename H = boost::hash<T> >
         class hash_map : public std::map<T,U> {};
 #else
     /// hash_map is an alias for std::tr1::unordered_map.
     /** We use the (experimental) TR1 unordered_map implementation included in modern GCC distributions.
      */
-    template <typename T, typename U>
-        class hash_map : public std::tr1::unordered_map<T,U> {};
+    template <typename T, typename U, typename H = boost::hash<T> >
+        class hash_map : public std::tr1::unordered_map<T,U,H> {};
 #endif
 
 
@@ -104,6 +132,11 @@ double rnd_stdnormal();
 /// Returns a random integer in interval [min, max]
 int rnd_int( int min, int max );
 
+/// Returns a random integer in the half-open interval \f$[0,n)\f$
+inline int rnd( int n) {
+    return rnd_int( 0, n-1 );
+}
+
 
 /// Writes a std::vector to a std::ostream
 template<class T> 
@@ -142,6 +175,17 @@ std::ostream& operator << (std::ostream& os, const std::pair<T1,T2> & x) {
     return os;
 }
 
+/// Concatenate two vectors
+template<class T>
+std::vector<T> concat( const std::vector<T>& u, const std::vector<T>& v ) {
+    std::vector<T> w;
+    w.reserve( u.size() + v.size() );
+    for( size_t i = 0; i < u.size(); i++ )
+        w.push_back( u[i] );
+    for( size_t i = 0; i < v.size(); i++ )
+        w.push_back( v[i] );
+    return w;
+}
 
 /// Used to keep track of the progress made by iterative algorithms
 class Diffs : public std::vector<double> {
diff --git a/scripts/regenerate-properties b/scripts/regenerate-properties
new file mode 100755 (executable)
index 0000000..1c3b45a
--- /dev/null
@@ -0,0 +1,342 @@
+#!/usr/bin/perl
+use warnings;
+use strict;
+
+use Cwd 'abs_path';
+
+@ARGV == 2 or die "Need 2 arguments";
+
+my ($header_file,$source_file) = @ARGV;
+@ARGV=();
+
+# Regular expressions for nested brackets (uses perl 5.10 features)
+my ($nested_angle_brackets) = qr/(<(?:[^<>]++|(?1))*>)/;
+my ($nested_curly_brackets) = qr/({(?:[^{}]++|(?1))*})/;
+my ($nested_parentheses) = qr/(\((?:[^()]++|(?1))*\))/;
+
+# Patterns to match generated code blocks
+my ($gen_code_start_pat, $gen_code_end_pat, $props_start_pat) =
+  (qr(/\*.*GENERATED CODE: DO NOT EDIT.*),
+   qr(/\*.*END OF GENERATED CODE.*\*/),
+   qr(/\*.*PROPERTIES));
+# Actual delimiters to use for generated code blocks
+my ($gen_code_start, $gen_code_end) =
+  ("/* {{{ GENERATED CODE: DO NOT EDIT. Created by \n".
+   "    $0 $header_file $source_file \n*/\n",
+   "/* }}} END OF GENERATED CODE */\n");
+
+# Strings to hold text of files
+my $header_buffer="";
+my $source_buffer="";
+
+# ----------------------------------------------------------------
+# Useful routines
+
+# remove leading and trailing whitespace
+sub stripws ($) {
+  my ($s) = @_;
+  $s =~ s/^\s*//s;
+  $s =~ s/\s*$//s;
+  return $s;
+}
+
+# split comments and non-comments, returning 2-element array
+# of (uncommented part, comments)
+my $comment_pat = qr#(/\*[^*]*\*+(?:[^/*][^*]*\*+)*/|//[^\n]*)|("(?:\\.|[^"\\])*"|'(?:\\.|[^'\\])*'|.[^/"'\\]*)#;
+sub sepcomment ($) {
+  my ($c) = @_;
+  my ($u) = $c;
+  $u =~ s#$comment_pat#defined $2 ? $2 : ""#gse;
+  $c =~ s#$comment_pat#defined $1 ? $1 : ""#gse;
+  return ($u, $c);
+}
+
+# Run diff, return output
+sub getdiff ($$) {
+  my ($a,$b) = @_;
+my $diff="";
+open DIFF, "diff -u \Q$a\E \Q$b\E |"
+  or die "Couldn't run diff";
+while(<DIFF>) {
+  $diff .= $_;
+}
+close DIFF;
+  return $diff;
+}
+
+sub myrename ($$) {
+  my ($a,$b) = @_;
+  $b = abs_path $b;
+  rename($a, $b) or die "Couldn't rename $a to $b";
+}
+
+sub writefile ($$) {
+  my ($buf, $f) = @_;
+  open OUT, ">", $f or die "Couldn't open $f for writing";
+  print OUT $buf;
+  close OUT;
+}
+
+our ($buffer);
+
+# Step through file, appending lines to buffer
+# - Lines which match generated code blocks are omitted
+# - Other lines are passed to the subroutine. They are added to buffer
+# unless the subroutine returns 1 (e.g. if it has already added them)
+sub process_file (&$) {
+  my ($sub, $file) = @_;
+  local ($buffer) = "";
+  # step through lines of file, appending each one to buffer
+  open IN, "<", $file or die "Couldn't open $file for reading";
+  while (<IN>) {
+    # delete anything between GENERATED CODE etc.
+    if (/$gen_code_end_pat/) {
+      die "Unmatched generated code block end at $file line $.";
+    } elsif (/$gen_code_start_pat/) {
+      my $s = $.;
+      my $found=0;
+      while (<IN>) {
+        chomp;
+        if (/$gen_code_end_pat/) {
+          $found=1;
+          last;
+        }
+      }
+      die "Unterminated generated code block at $file line $s"
+        unless $found;
+      next;
+    } else {
+      my ($res) = &$sub;
+      next if $res;
+    }
+    $buffer .= $_;
+  }
+  close IN;
+  return $buffer;
+}
+
+# Parse a PROPERTIES() block
+sub process_properties ($$$) {
+  my ($file, $start_line, $props_text) = @_;
+  my ($errline)="$file:$start_line";
+  $props_text =~ s/^.*PROPERTIES\s*($nested_parentheses)//s
+    or die "Malformed PROPERTIES, $errline: expected PROPERTIES(args)";
+  my ($args) = $1;
+  $args =~ s/^\(//g;
+  $args =~ s/\)$//g;
+  my (@args) = split /,/, $args;
+  @args == 2 or die "PROPERTIES needs two arguments at $errline";
+  my ($struct_name, $class) = @args;
+
+  $props_text =~ m/^\s*($nested_curly_brackets)\s*$/s
+    or die "Malformed PROPERTIES, $errline: expected {} after PROPERTIES";
+  my ($body) = $1;
+  $body =~ s/^{(.*)}$/$1/s; # get rid of curly brackets
+  my (@stmts) = split /;\s*\n/s, $body; # split on ";"
+  my (@typedecls) = ();
+  my (@vars) = ();
+  foreach my $s (@stmts) {
+    my ($s, $cmt) = sepcomment $s;
+    $cmt = stripws $cmt;
+    if($s =~ /^\s*$/s) { # extra ";"
+      next;
+    } elsif($s =~ /DAI_ENUM|typedef/) {
+      push @typedecls, [stripws $s, $cmt];
+    } else {
+      $s =~ /^\s*[a-zA-Z_][\w:]+\s*$nested_angle_brackets?/
+        or die "Couldn't match statement $s in PROPERTIES at $errline";
+      my $type = stripws $&;
+      my $name = stripws $';
+      my $default = undef;
+      if($name =~ /^(.*)=(.*)$/) {
+        $name = stripws $1;
+        $default = stripws $2;
+      }
+      push @vars, [$type, $name, $default, $cmt];
+    }
+  }
+  
+  my ($stext) = "";
+  my ($text) = <<EOF;
+        struct Properties {
+EOF
+  my $indent = (' 'x12);
+  for my $d (@typedecls) {
+    my ($decl,$cmt) = @$d;
+    if($cmt ne '') {
+      $text .= join("\n", map { "$indent$_" } split /\n/s, $cmt)."\n";
+    }
+    $text .= "$indent$decl;\n"
+  }
+  for my $v (@vars) {
+    my ($type,$name,$default,$cmt) = @$v;
+    if($cmt ne '') {
+      $text .= join("\n", map { "$indent$_" } split /\n/s, $cmt)."\n"
+    }
+    $text .= "$indent$type $name;\n";
+  }
+
+  $text .= <<EOF;
+
+            /// Set members from PropertySet
+            void set(const PropertySet &opts);
+            /// Get members into PropertySet
+            PropertySet get() const;
+            /// Convert to a string which can be parsed as a PropertySet
+            std::string toString() const;
+        } $struct_name;
+EOF
+
+  $stext .= <<EOF;
+namespace dai {
+
+void ${class}::Properties::set(const PropertySet &opts)
+{
+    const std::set<PropertyKey> &keys = opts.allKeys();
+    std::set<PropertyKey>::const_iterator i;
+    bool die=false;
+    for(i=keys.begin(); i!=keys.end(); i++) {
+EOF
+  for my $v (@vars) {
+    my ($type,$name,$default,$cmt) = @$v;
+    $stext .= <<EOF;
+        if(*i == "$name") continue;
+EOF
+  }
+  $stext .= <<EOF;
+        cerr << "$class: Unknown property " << *i << endl;
+        die=true;
+    }
+    if(die) {
+        DAI_THROW(UNKNOWN_PROPERTY_TYPE);
+    }
+EOF
+  for my $v (@vars) {
+    my ($type,$name,$default,$cmt) = @$v;
+    if(!defined $default) {
+      $stext .= <<EOF;
+    if(!opts.hasKey("$name")) {
+        cerr << "$class: Missing property \\\"$name\\\" for method \\\"\Q$class\E\\\"" << endl;
+        die=true;
+    }
+EOF
+    }
+
+  }
+  $stext .= <<EOF;
+    if(die) {
+        DAI_THROW(NOT_ALL_PROPERTIES_SPECIFIED);
+    }
+EOF
+  for my $v (@vars) {
+    my ($type,$name,$default,$cmt) = @$v;
+    if(defined $default) {
+      $stext .= <<EOF;
+    if(opts.hasKey("$name")) {
+        $name = opts.getStringAs<$type>("$name");
+    } else {
+        $name = $default;
+    }
+EOF
+    } else {
+      $stext .= <<EOF;
+    $name = opts.getStringAs<$type>("$name");
+EOF
+    }
+  }
+
+  $stext .= <<EOF;
+}
+PropertySet ${class}::Properties::get() const {
+    PropertySet opts;
+EOF
+
+  for my $v (@vars) {
+    my ($type,$name,$default,$cmt) = @$v;
+#     $text .= qq{opts.set("$name", $name);};
+    $stext .= <<EOF;
+    opts.Set("$name", $name);
+EOF
+
+  }
+  $stext .= <<EOF;
+    return opts;
+}
+string ${class}::Properties::toString() const {
+    stringstream s(stringstream::out);
+    s << "[";
+EOF
+
+  my ($i)=0;
+  for my $v (@vars) {
+    my ($type,$name,$default,$cmt) = @$v;
+    $i++;
+    my $c = ($i<@vars)?" << \",\"":"";
+      $stext .= <<EOF;
+    s << "$name=" << $name$c;
+EOF
+  }
+
+  $stext .= <<EOF;
+    s << "]";
+    return s.str();
+}
+} // end of namespace dai
+EOF
+  return ($gen_code_start.$text.$gen_code_end, $gen_code_start.$stext.$gen_code_end);
+}
+
+# ----------------------------------------------------------------
+# Main loop
+
+$source_buffer = process_file { return 0; } $source_file;
+$source_buffer =~ s/\n+$//s;
+
+$header_buffer = process_file { 
+  if (/$props_start_pat/) {
+    # when we see something resembling properties, record it, and when we
+    # get to the end, process and emit the generated code
+    my $start_line = $.;
+    my $props_text = $_;
+    while (<IN>) {
+      last if(m(\*/));
+      $props_text .= $_;
+    }
+    $buffer .= $props_text;
+    $buffer .= $_;
+    my ($htext, $stext) = process_properties($header_file, $start_line, $props_text);
+    $buffer .= $htext;
+    $source_buffer .= "\n\n\n".$stext;
+    return 1;
+  }
+  return 0;
+} $header_file;
+
+# Write new contents of files to temporary locations
+my $header_tmp = "$header_file.new";
+my $source_tmp = "$source_file.new";
+writefile ($header_buffer, $header_tmp);
+writefile ($source_buffer, $source_tmp);
+
+# Get a diff of changes, show it to the user, and ask for confirmation
+my ($diff);
+$diff = getdiff ($header_file, $header_tmp);
+$diff .= getdiff ($source_file, $source_tmp);
+
+if($diff eq '') {
+  warn "No differences\n";
+} else {
+  my $pager = $ENV{PAGER} || "less";
+  open PAGER, "|$pager";
+  print PAGER $diff;
+  close PAGER;
+
+  print STDERR "Replace old with new version? (y|N) ";
+  $_=<>;
+  if (/^y/i) {
+    myrename($header_tmp, $header_file);
+    myrename($source_tmp, $source_file);
+  } else {
+    warn "Aborted\n";
+  }
+}
index 643c917..3c5fbba 100644 (file)
@@ -32,7 +32,7 @@ namespace dai {
 using namespace std;
 
 
-InfAlg *newInfAlg( const string &name, const FactorGraph &fg, const PropertySet &opts ) {
+InfAlg *newInfAlg( const std::string &name, const FactorGraph &fg, const PropertySet &opts ) {
     if( name == ExactInf::Name )
         return new ExactInf (fg, opts);
 #ifdef DAI_WITH_BP
@@ -66,9 +66,31 @@ InfAlg *newInfAlg( const string &name, const FactorGraph &fg, const PropertySet
 #ifdef DAI_WITH_GIBBS
     if( name == Gibbs::Name )
         return new Gibbs (fg, opts);
+#endif
+#ifdef DAI_WITH_CBP
+    if( name == CBP::Name )
+        return new CBP (fg, opts);
 #endif
     DAI_THROW(UNKNOWN_DAI_ALGORITHM);
 }
 
 
+/// \todo Make alias file non-testdai-specific, and use it in newInfAlgFromString
+InfAlg *newInfAlgFromString( const std::string &nameOpts, const FactorGraph &fg ) {
+    string::size_type pos = nameOpts.find_first_of('[');
+    string name;
+    PropertySet opts;
+    if( pos == string::npos ) {
+        name = nameOpts;
+    } else {
+        name = nameOpts.substr(0,pos);
+
+        stringstream ss;
+        ss << nameOpts.substr(pos,nameOpts.length());
+        ss >> opts;
+    }
+    return newInfAlg(name,fg,opts);
+}
+
+
 } // end of namespace dai
diff --git a/src/bbp.cpp b/src/bbp.cpp
new file mode 100644 (file)
index 0000000..1e55c2e
--- /dev/null
@@ -0,0 +1,1237 @@
+/*  Copyright (C) 2009  Frederik Eaton [frederik at ofb dot net]
+
+    This file is part of libDAI.
+
+    libDAI is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    libDAI is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with libDAI; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+
+#include <dai/bp.h>
+#include <dai/bbp.h>
+#include <dai/gibbs.h>
+#include <dai/util.h>
+#include <dai/bipgraph.h>
+
+
+namespace dai {
+
+
+using namespace std;
+
+
+typedef BipartiteGraph::Neighbor Neighbor;
+
+
+Prob unnormAdjoint( const Prob &w, Real Z_w, const Prob &adj_w ) {
+    assert( w.size() == adj_w.size() );
+    Prob adj_w_unnorm( w.size(), 0.0 );
+    Real s = 0.0;
+    for( size_t i = 0; i < w.size(); i++ )
+        s += w[i] * adj_w[i];
+    for( size_t i = 0; i < w.size(); i++ )
+        adj_w_unnorm[i] = (adj_w[i] - s) / Z_w;
+    return adj_w_unnorm;
+//  THIS WOULD BE ABOUT 50% SLOWER:  return (adj_w - (w * adj_w).sum()) / Z_w;
+}
+
+
+std::vector<size_t> getGibbsState( const InfAlg &ia, size_t iters ) {
+    PropertySet gibbsProps;
+    gibbsProps.Set("iters", iters);
+    gibbsProps.Set("verbose", size_t(0));
+    Gibbs gibbs( ia.fg(), gibbsProps );
+    gibbs.run();
+    return gibbs.state();
+}
+
+
+size_t getFactorEntryForState( const FactorGraph &fg, size_t I, const vector<size_t> &state ) {
+    size_t f_entry = 0;
+    for( int _j = fg.nbF(I).size() - 1; _j >= 0; _j-- ) {
+        // note that iterating over nbF(I) yields the same ordering
+        // of variables as iterating over factor(I).vars()
+        size_t j = fg.nbF(I)[_j];
+        f_entry *= fg.var(j).states();
+        f_entry += state[j];
+    }
+    return f_entry;
+}
+
+
+#define LOOP_ij(body) {                             \
+    size_t i_states = _fg->var(i).states();         \
+    size_t j_states = _fg->var(j).states();         \
+    if(_fg->var(i) > _fg->var(j)) {                 \
+        size_t xij=0;                               \
+        for(size_t xi=0; xi<i_states; xi++) {       \
+            for(size_t xj=0; xj<j_states; xj++) {   \
+                body;                               \
+                xij++;                              \
+            }                                       \
+        }                                           \
+    } else {                                        \
+        size_t xij=0;                               \
+        for(size_t xj=0; xj<j_states; xj++) {       \
+            for(size_t xi=0; xi<i_states; xi++) {   \
+                body;                               \
+                xij++;                              \
+            }                                       \
+        }                                           \
+    }                                               \
+}
+
+
+void BBP::RegenerateInds() {
+    // initialise _indices
+    //     typedef std::vector<size_t>        _ind_t;
+    //     std::vector<std::vector<_ind_t> >  _indices; 
+    _indices.resize( _fg->nrVars() );
+    for( size_t i = 0; i < _fg->nrVars(); i++ ) {
+        _indices[i].clear();
+        _indices[i].reserve( _fg->nbV(i).size() );
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
+            _ind_t index;
+            index.reserve( _fg->factor(I).states() );
+            for( IndexFor k(_fg->var(i), _fg->factor(I).vars()); k >= 0; ++k )
+                index.push_back( k );
+            _indices[i].push_back( index );
+        }
+    }
+}
+
+
+void BBP::RegenerateT() {
+    // _T[i][_I]
+    _T.resize( _fg->nrVars() );
+    for( size_t i = 0; i < _fg->nrVars(); i++ ) {
+        _T[i].resize( _fg->nbV(i).size() );
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
+            Prob prod( _fg->var(i).states(), 1.0 );
+            foreach( const Neighbor &J, _fg->nbV(i) )
+                if( J.node != I.node )
+                    prod *= _bp_dual.msgM( i, J.iter );
+            _T[i][I.iter] = prod;
+        }
+    }
+}
+
+
+void BBP::RegenerateU() {
+    // _U[I][_i]
+    _U.resize( _fg->nrFactors() );
+    for( size_t I = 0; I < _fg->nrFactors(); I++ ) {
+        _U[I].resize( _fg->nbF(I).size() );
+        foreach( const Neighbor &i, _fg->nbF(I) ) {
+            Prob prod( _fg->factor(I).states(), 1.0 );
+            foreach( const Neighbor &j, _fg->nbF(I) )
+                if( i.node != j.node ) {
+                    Prob n_jI( _bp_dual.msgN( j, j.dual ) );
+                    const _ind_t &ind = _index( j, j.dual );
+                    // multiply prod by n_jI
+                    for( size_t x_I = 0; x_I < prod.size(); x_I++ )
+                        prod[x_I] *= n_jI[ind[x_I]];
+                }
+            _U[I][i.iter] = prod;
+        }
+    }
+}
+
+
+void BBP::RegenerateS() {
+    // _S[i][_I][_j]
+    _S.resize( _fg->nrVars() );
+    for( size_t i = 0; i < _fg->nrVars(); i++ ) {
+        _S[i].resize( _fg->nbV(i).size() );
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
+            _S[i][I.iter].resize( _fg->nbF(I).size() );
+            foreach( const Neighbor &j, _fg->nbF(I) )
+                if( i != j ) {
+                    Factor prod( _fg->factor(I) );
+                    foreach( const Neighbor &k, _fg->nbF(I) ) {
+                        if( k != i && k.node != j.node ) {
+                            const _ind_t &ind = _index( k, k.dual );
+                            Prob p( _bp_dual.msgN( k, k.dual ) );
+                            for( size_t x_I = 0; x_I < prod.states(); x_I++ )
+                                prod.p()[x_I] *= p[ind[x_I]];
+                        }
+                    }
+                    // "Marginalize" onto i|j (unnormalized)
+                    Prob marg;
+                    marg = prod.marginal( VarSet(_fg->var(i), _fg->var(j)), false ).p();
+                    _S[i][I.iter][j.iter] = marg;
+                }
+        }
+    }
+}
+
+
+void BBP::RegenerateR() {
+    // _R[I][_i][_J]
+    _R.resize( _fg->nrFactors() );
+    for( size_t I = 0; I < _fg->nrFactors(); I++ ) {
+        _R[I].resize( _fg->nbF(I).size() );
+        foreach( const Neighbor &i, _fg->nbF(I) ) {
+            _R[I][i.iter].resize( _fg->nbV(i).size() );
+            foreach( const Neighbor &J, _fg->nbV(i) ) {
+                if( I != J ) {
+                    Prob prod( _fg->var(i).states(), 1.0 );
+                    foreach( const Neighbor &K, _fg->nbV(i) )
+                        if( K.node != I && K.node != J.node ) 
+                            prod *= _bp_dual.msgM( i, K.iter );
+                    _R[I][i.iter][J.iter] = prod;
+                }
+            }
+        }
+    }
+}
+
+
+void BBP::RegenerateInputs() {
+    _adj_b_V_unnorm.clear();
+    _adj_b_V_unnorm.reserve( _fg->nrVars() );
+    for( size_t i = 0; i < _fg->nrVars(); i++ )
+        _adj_b_V_unnorm.push_back( unnormAdjoint( _bp_dual.beliefV(i).p(), _bp_dual.beliefVZ(i), _adj_b_V[i] ) );
+    _adj_b_F_unnorm.clear();
+    _adj_b_F_unnorm.reserve( _fg->nrFactors() );
+    for( size_t I = 0; I < _fg->nrFactors(); I++ )
+        _adj_b_F_unnorm.push_back( unnormAdjoint( _bp_dual.beliefF(I).p(), _bp_dual.beliefFZ(I), _adj_b_F[I] ) );
+}
+
+
+void BBP::RegeneratePsiAdjoints() {
+    _adj_psi_V.clear();
+    _adj_psi_V.reserve( _fg->nrVars() );
+    for( size_t i = 0; i < _fg->nrVars(); i++ ) {
+        Prob p( _adj_b_V_unnorm[i] );
+        assert( p.size() == _fg->var(i).states() );
+        foreach( const Neighbor &I, _fg->nbV(i) )
+            p *= _bp_dual.msgM( i, I.iter );
+        p += _init_adj_psi_V[i];
+        _adj_psi_V.push_back( p );
+    }
+    _adj_psi_F.clear();
+    _adj_psi_F.reserve( _fg->nrFactors() );
+    for( size_t I = 0; I < _fg->nrFactors(); I++ ) {
+        Prob p( _adj_b_F_unnorm[I] );
+        assert( p.size() == _fg->factor(I).states() );
+        foreach( const Neighbor &i, _fg->nbF(I) ) {
+            Prob n_iI( _bp_dual.msgN( i, i.dual ) );
+            const _ind_t& ind = _index( i, i.dual );
+            // multiply prod with n_jI
+            for( size_t x_I = 0; x_I < p.size(); x_I++ )
+                p[x_I] *= n_iI[ind[x_I]];
+        }
+        p += _init_adj_psi_F[I];
+        _adj_psi_F.push_back( p );
+    }
+}
+
+
+void BBP::RegenerateParMessageAdjoints() {
+    size_t nv = _fg->nrVars();
+    _adj_n.resize( nv );
+    _adj_m.resize( nv );
+    _adj_n_unnorm.resize( nv );
+    _adj_m_unnorm.resize( nv );
+    _new_adj_n.resize( nv );
+    _new_adj_m.resize( nv );
+    for( size_t i = 0; i < _fg->nrVars(); i++ ) {
+        size_t n_i = _fg->nbV(i).size();
+        _adj_n[i].resize( n_i );
+        _adj_m[i].resize( n_i );
+        _adj_n_unnorm[i].resize( n_i );
+        _adj_m_unnorm[i].resize( n_i );
+        _new_adj_n[i].resize( n_i );
+        _new_adj_m[i].resize( n_i );
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
+            { // calculate adj_n
+                Prob prod( _fg->factor(I).p() );
+                prod *= _adj_b_F_unnorm[I];
+                foreach( const Neighbor &j, _fg->nbF(I) )
+                    if( i != j ) {
+                        Prob n_jI( _bp_dual.msgN( j, j.dual ) );
+                        const _ind_t &ind = _index( j, j.dual );
+                        // multiply prod with n_jI
+                        for( size_t x_I = 0; x_I < prod.size(); x_I++ )
+                            prod[x_I] *= n_jI[ind[x_I]];
+                    }
+                Prob marg( _fg->var(i).states(), 0.0 );
+                const _ind_t &ind = _index( i, I.iter );
+                for( size_t r = 0; r < prod.size(); r++ )
+                    marg[ind[r]] += prod[r];
+                _new_adj_n[i][I.iter] = marg;
+                upMsgN( i, I.iter );
+            }
+
+            { // calculate adj_m
+                Prob prod( _adj_b_V_unnorm[i] );
+                assert( prod.size() == _fg->var(i).states() );
+                foreach( const Neighbor &J, _fg->nbV(i) )
+                    if( J.node != I.node )
+                        prod *= _bp_dual.msgM(i,J.iter);
+                _new_adj_m[i][I.iter] = prod;
+                upMsgM( i, I.iter );
+            }
+        }
+    }
+}
+
+
+void BBP::RegenerateSeqMessageAdjoints() {
+    size_t nv = _fg->nrVars();
+    _adj_m.resize( nv );
+    _adj_m_unnorm.resize( nv );
+    _new_adj_m.resize( nv );
+    for( size_t i = 0; i < _fg->nrVars(); i++ ) {
+        size_t n_i = _fg->nbV(i).size();
+        _adj_m[i].resize( n_i );
+        _adj_m_unnorm[i].resize( n_i );
+        _new_adj_m[i].resize( n_i );
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
+            // calculate adj_m
+            Prob prod( _adj_b_V_unnorm[i] );
+            assert( prod.size() == _fg->var(i).states() );
+            foreach( const Neighbor &J, _fg->nbV(i) )
+                if( J.node != I.node )
+                    prod *= _bp_dual.msgM( i, J.iter );
+            _adj_m[i][I.iter] = prod;
+            calcUnnormMsgM( i, I.iter );
+            _new_adj_m[i][I.iter] = Prob( _fg->var(i).states(), 0.0 );
+        }
+    }
+    for( size_t i = 0; i < _fg->nrVars(); i++ ) {
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
+            // calculate adj_n
+            Prob prod( _fg->factor(I).p() );
+            prod *= _adj_b_F_unnorm[I];
+            foreach( const Neighbor &j, _fg->nbF(I) )
+                if( i != j ) {
+                    Prob n_jI( _bp_dual.msgN( j, j.dual) );
+                    const _ind_t& ind = _index( j, j.dual );
+                    // multiply prod with n_jI
+                    for( size_t x_I = 0; x_I < prod.size(); x_I++ )
+                        prod[x_I] *= n_jI[ind[x_I]];
+                }
+            Prob marg( _fg->var(i).states(), 0.0 );
+            const _ind_t &ind = _index( i, I.iter );
+            for( size_t r = 0; r < prod.size(); r++ )
+                marg[ind[r]] += prod[r];
+            sendSeqMsgN( i, I.iter,marg );
+        }
+    }
+}
+
+
+void BBP::calcNewN( size_t i, size_t _I ) {
+    _adj_psi_V[i] += T(i,_I) * _adj_n_unnorm[i][_I];
+    Prob &new_adj_n_iI = _new_adj_n[i][_I];
+    new_adj_n_iI = Prob( _fg->var(i).states(), 0.0 );
+    size_t I = _fg->nbV(i)[_I];
+    foreach( const Neighbor &j, _fg->nbF(I) )
+        if( j != i ) {
+            const Prob &p = _S[i][_I][j.iter];
+            const Prob &_adj_m_unnorm_jI = _adj_m_unnorm[j][j.dual];
+            LOOP_ij(
+                new_adj_n_iI[xi] += p[xij] * _adj_m_unnorm_jI[xj];
+            );
+            /* THE FOLLOWING WOULD BE ABOUT TWICE AS SLOW:
+            Var vi = _fg->var(i);
+            Var vj = _fg->var(j);
+            new_adj_n_iI = (Factor(VarSet(vi, vj), p) * Factor(vj,_adj_m_unnorm_jI)).marginal(vi,false).p();
+            */
+        }
+}
+
+
+void BBP::calcNewM( size_t i, size_t _I ) {
+    const Neighbor &I = _fg->nbV(i)[_I];
+    Prob p( U(I, I.dual) );
+    const Prob &adj = _adj_m_unnorm[i][_I];
+    const _ind_t &ind = _index(i,_I);
+    for( size_t x_I = 0; x_I < p.size(); x_I++ )
+        p[x_I] *= adj[ind[x_I]];
+    _adj_psi_F[I] += p;
+    /* THE FOLLOWING WOULD BE SLIGHTLY SLOWER:
+    _adj_psi_F[I] += (Factor( _fg->factor(I).vars(), U(I, I.dual) ) * Factor( _fg->var(i), _adj_m_unnorm[i][_I] )).p();
+    */
+
+    _new_adj_m[i][_I] = Prob( _fg->var(i).states(), 0.0 );
+    foreach( const Neighbor &J, _fg->nbV(i) )
+        if( J != I )
+            _new_adj_m[i][_I] += _R[I][I.dual][J.iter] * _adj_n_unnorm[i][J.iter];
+}
+
+
+void BBP::calcUnnormMsgN( size_t i, size_t _I ) {
+    _adj_n_unnorm[i][_I] = unnormAdjoint( _bp_dual.msgN(i,_I), _bp_dual.zN(i,_I), _adj_n[i][_I] );
+}
+
+
+void BBP::calcUnnormMsgM( size_t i, size_t _I ) {
+    _adj_m_unnorm[i][_I] = unnormAdjoint( _bp_dual.msgM(i,_I), _bp_dual.zM(i,_I), _adj_m[i][_I] );
+}
+
+
+void BBP::upMsgN( size_t i, size_t _I ) {
+    _adj_n[i][_I] = _new_adj_n[i][_I];
+    calcUnnormMsgN( i, _I );
+}
+
+
+void BBP::upMsgM( size_t i, size_t _I ) {
+    _adj_m[i][_I] = _new_adj_m[i][_I];
+    calcUnnormMsgM( i, _I );
+}
+
+
+void BBP::doParUpdate() {
+    for( size_t i = 0; i < _fg->nrVars(); i++ )
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
+            calcNewM( i, I.iter );
+            calcNewN( i, I.iter );
+        }
+    for( size_t i = 0; i < _fg->nrVars(); i++ )
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
+            upMsgM( i, I.iter );
+            upMsgN( i, I.iter );
+        }
+}
+
+
+void BBP::incrSeqMsgM( size_t i, size_t _I, const Prob &p ) {
+    if( props.clean_updates )
+        _new_adj_m[i][_I] += p;
+    else {
+        _adj_m[i][_I] += p;
+        calcUnnormMsgM(i, _I);
+    }
+}
+
+
+#if 0
+Real pv_thresh=1000;
+#endif
+
+
+void BBP::updateSeqMsgM( size_t i, size_t _I ) {
+    if( props.clean_updates ) {
+#if 0
+        if(_new_adj_m[i][_I].sumAbs() > pv_thresh ||
+           _adj_m[i][_I].sumAbs() > pv_thresh) {
+
+            DAI_DMSG("in updateSeqMsgM");
+            DAI_PV(i);
+            DAI_PV(_I);
+            DAI_PV(_adj_m[i][_I]);
+            DAI_PV(_new_adj_m[i][_I]);
+        }
+#endif
+        _adj_m[i][_I] += _new_adj_m[i][_I];
+        calcUnnormMsgM( i, _I );
+        _new_adj_m[i][_I].fill( 0.0 );
+    }
+}
+
+
+void BBP::setSeqMsgM( size_t i, size_t _I, const Prob &p ) {
+    _adj_m[i][_I] = p;
+    calcUnnormMsgM( i, _I );
+}
+
+
+void BBP::sendSeqMsgN( size_t i, size_t _I, const Prob &f ) {
+    Prob f_unnorm = unnormAdjoint( _bp_dual.msgN(i,_I), _bp_dual.zN(i,_I), f );
+    const Neighbor &I = _fg->nbV(i)[_I];
+    assert( I.iter == _I );
+    _adj_psi_V[i] += f_unnorm * T( i, _I );
+#if 0
+    if(f_unnorm.sumAbs() > pv_thresh) {
+        DAI_DMSG("in sendSeqMsgN");
+        DAI_PV(i);
+        DAI_PV(I);
+        DAI_PV(_I);
+        DAI_PV(_bp_dual.msgN(i,_I));
+        DAI_PV(_bp_dual.zN(i,_I));
+        DAI_PV(_bp_dual.msgM(i,_I));
+        DAI_PV(_bp_dual.zM(i,_I));
+        DAI_PV(_fg->factor(I).p());
+    }
+#endif
+    foreach( const Neighbor &J, _fg->nbV(i) ) {
+        if( J.node != I.node ) {
+#if 0
+            if(f_unnorm.sumAbs() > pv_thresh) {
+                DAI_DMSG("in sendSeqMsgN loop");
+                DAI_PV(J);
+                DAI_PV(f_unnorm);
+                DAI_PV(_R[J][J.dual][_I]);
+                DAI_PV(f_unnorm * _R[J][J.dual][_I]);
+            }
+#endif
+            incrSeqMsgM( i, J.iter, f_unnorm * R(J, J.dual, _I) );
+        }
+    }
+}
+
+
+void BBP::sendSeqMsgM( size_t j, size_t _I ) {
+    const Neighbor &I = _fg->nbV(j)[_I];
+
+//     DAI_PV(j);
+//     DAI_PV(I);
+//     DAI_PV(_adj_m_unnorm_jI);
+//     DAI_PV(_adj_m[j][_I]);
+//     DAI_PV(_bp_dual.zM(j,_I));
+
+    size_t _j = I.dual;
+    const Prob &_adj_m_unnorm_jI = _adj_m_unnorm[j][_I];
+    Prob um( U(I, _j) );
+    const _ind_t &ind = _index(j, _I);
+    for( size_t x_I = 0; x_I < um.size(); x_I++ )
+        um[x_I] *= _adj_m_unnorm_jI[ind[x_I]];
+    um *= 1 - props.damping;
+    _adj_psi_F[I] += um;
+    
+    /* THE FOLLOWING WOULD BE SLIGHTLY SLOWER:
+    _adj_psi_F[I] += (Factor( _fg->factor(I).vars(), U(I, _j) ) * Factor( _fg->var(j), _adj_m_unnorm[j][_I] )).p() * (1.0 - props.damping);
+    */
+
+//     DAI_DMSG("in sendSeqMsgM");
+//     DAI_PV(j);
+//     DAI_PV(I);
+//     DAI_PV(_I);
+//     DAI_PV(_fg->nbF(I).size());
+    foreach( const Neighbor &i, _fg->nbF(I) ) {
+        if( i.node != j ) {
+            const Prob &S = _S[i][i.dual][_j];
+            Prob msg( _fg->var(i).states(), 0.0 );
+            LOOP_ij(
+                msg[xi] += S[xij] * _adj_m_unnorm_jI[xj];
+            );
+            msg *= 1.0 - props.damping;
+            /* THE FOLLOWING WOULD BE ABOUT TWICE AS SLOW:
+            Var vi = _fg->var(i);
+            Var vj = _fg->var(j);
+            msg = (Factor(VarSet(vi,vj), S) * Factor(vj,_adj_m_unnorm_jI)).marginal(vi,false).p() * (1.0 - props.damping);
+            */
+#if 0
+            if(msg.sumAbs() > pv_thresh) {
+                DAI_DMSG("in sendSeqMsgM loop");
+
+                DAI_PV(j);
+                DAI_PV(I);
+                DAI_PV(_I);
+                DAI_PV(_fg->nbF(I).size());
+                DAI_PV(_fg->factor(I).p());
+                DAI_PV(_S[i][i.dual][_j]);
+
+                DAI_PV(i);
+                DAI_PV(i.dual);
+                DAI_PV(msg);
+                DAI_PV(_fg->nbV(i).size());
+            }
+#endif
+            assert( _fg->nbV(i)[i.dual].node == I );
+            sendSeqMsgN( i, i.dual, msg );
+        }
+    }
+    setSeqMsgM( j, _I, _adj_m[j][_I] * props.damping );
+}
+
+
+Real BBP::getUnMsgMag() {
+    Real s = 0.0;
+    size_t e = 0;
+    for( size_t i = 0; i < _fg->nrVars(); i++ )
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
+            s += _adj_m_unnorm[i][I.iter].sumAbs();
+            s += _adj_n_unnorm[i][I.iter].sumAbs();
+            e++;
+        }
+    return s / e;
+}
+
+
+void BBP::getMsgMags( Real &s, Real &new_s ) {
+    s = 0.0; 
+    new_s = 0.0;
+    size_t e = 0;
+    for( size_t i = 0; i < _fg->nrVars(); i++ )
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
+            s += _adj_m[i][I.iter].sumAbs();
+            s += _adj_n[i][I.iter].sumAbs();
+            new_s += _new_adj_m[i][I.iter].sumAbs();
+            new_s += _new_adj_n[i][I.iter].sumAbs();
+            e++;
+        }
+    s /= e; 
+    new_s /= e;
+}
+
+// tuple<size_t,size_t,Real> BBP::getArgMaxPsi1Adj() {
+//     size_t argmax_var=0;
+//     size_t argmax_var_state=0;
+//     Real max_var=0;
+//     for( size_t i = 0; i < _fg->nrVars(); i++ ) {
+//         pair<size_t,Real> argmax_state = adj_psi_V(i).argmax();
+//         if(i==0 || argmax_state.second>max_var) {
+//             argmax_var = i;
+//             max_var = argmax_state.second;
+//             argmax_var_state = argmax_state.first;
+//         }
+//     }
+//     assert(/*0 <= argmax_var_state &&*/
+//            argmax_var_state < _fg->var(argmax_var).states());
+//     return tuple<size_t,size_t,Real>(argmax_var,argmax_var_state,max_var);
+// }
+
+
+void BBP::getArgmaxMsgM( size_t &out_i, size_t &out__I, Real &mag ) {
+    bool found = false;
+    for( size_t i = 0; i < _fg->nrVars(); i++ )
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
+            Real thisMag = _adj_m[i][I.iter].sumAbs();
+            if( !found || mag < thisMag ) {
+                found = true;
+                mag = thisMag;
+                out_i = i;
+                out__I = I.iter;
+            }
+        }
+    assert( found );
+}
+
+
+Real BBP::getMaxMsgM() {
+    size_t dummy;
+    Real mag;
+    getArgmaxMsgM( dummy, dummy, mag );
+    return mag;
+}
+
+
+Real BBP::getTotalMsgM() {
+    Real mag = 0.0;
+    for( size_t i = 0; i < _fg->nrVars(); i++ )
+        foreach( const Neighbor &I, _fg->nbV(i) )
+            mag += _adj_m[i][I.iter].sumAbs();
+    return mag;
+}
+
+
+Real BBP::getTotalNewMsgM() {
+    Real mag = 0.0;
+    for( size_t i = 0; i < _fg->nrVars(); i++ )
+        foreach( const Neighbor &I, _fg->nbV(i) )
+            mag += _new_adj_m[i][I.iter].sumAbs();
+    return mag;
+}
+
+
+Real BBP::getTotalMsgN() {
+    Real mag = 0.0;
+    for( size_t i = 0; i < _fg->nrVars(); i++ )
+        foreach( const Neighbor &I, _fg->nbV(i) )
+            mag += _adj_n[i][I.iter].sumAbs();
+    return mag;
+}
+
+
+void BBP::Regenerate() {
+    RegenerateInds();
+    RegenerateT();
+    RegenerateU();
+    RegenerateS();
+    RegenerateR();
+    RegenerateInputs();
+    RegeneratePsiAdjoints();
+    if( props.updates == Properties::UpdateType::PAR )
+        RegenerateParMessageAdjoints();
+    else
+        RegenerateSeqMessageAdjoints();
+    _iters = 0;
+}
+
+
+std::vector<Prob> BBP::getZeroAdjF( const FactorGraph &fg ) {
+    vector<Prob> adj_2;
+    adj_2.reserve( fg.nrFactors() );
+    for( size_t I = 0; I < fg.nrFactors(); I++ )
+        adj_2.push_back( Prob( fg.factor(I).states(), 0.0 ) );
+    return adj_2;
+}
+
+
+std::vector<Prob> BBP::getZeroAdjV( const FactorGraph &fg ) {
+    vector<Prob> adj_1;
+    adj_1.reserve( fg.nrVars() );
+    for( size_t i = 0; i < fg.nrVars(); i++ )
+        adj_1.push_back( Prob( fg.var(i).states(), 0.0 ) );
+    return adj_1;
+}
+
+
+void BBP::run() {
+    typedef BBP::Properties::UpdateType UT;
+    Real &tol = props.tol;
+    UT &updates = props.updates;
+
+    Real tic = toc();
+    switch( (size_t)updates ) {
+        case UT::SEQ_MAX: {
+            size_t i, _I;
+            Real mag;
+            do {
+                _iters++;
+                getArgmaxMsgM( i, _I, mag );
+                sendSeqMsgM( i, _I );
+            } while( mag > tol && _iters < props.maxiter );
+
+            if( _iters >= props.maxiter )
+                if( props.verbose >= 1 )
+                    cerr << "Warning: BBP didn't converge in " << _iters << " iterations (greatest message magnitude = " << mag << ")" << endl;
+            break;
+        } case UT::SEQ_FIX: {
+            Real mag;
+            do {
+                _iters++;
+                mag = getTotalMsgM();
+                if( mag < tol ) 
+                    break;
+
+                for( size_t i = 0; i < _fg->nrVars(); i++ )
+                    foreach( const Neighbor &I, _fg->nbV(i) )
+                        sendSeqMsgM( i, I.iter );
+                for( size_t i = 0; i < _fg->nrVars(); i++ )
+                    foreach( const Neighbor &I, _fg->nbV(i) )
+                        updateSeqMsgM( i, I.iter );
+            } while( mag > tol && _iters < props.maxiter );
+
+            if( _iters >= props.maxiter )
+                if( props.verbose >= 1 )
+                    cerr << "Warning: BBP didn't converge in " << _iters << " iterations (greatest message magnitude = " << mag << ")" << endl;
+            break;
+        } case UT::SEQ_BP_REV:
+          case UT::SEQ_BP_FWD: {
+            const BP *bp = static_cast<const BP*>(_ia);
+            vector<pair<size_t, size_t> > sentMessages = bp->getSentMessages();
+            size_t totalMessages = sentMessages.size();
+            if( totalMessages == 0 ) {
+                cerr << "Asked for updates = " << updates << " but no BP messages; did you forget to set recordSentMessages?" << endl;
+                DAI_THROW(INTERNAL_ERROR);
+            }
+            if( updates==UT::SEQ_BP_FWD )
+                reverse( sentMessages.begin(), sentMessages.end() );
+//          DAI_PV(sentMessages.size());
+//          DAI_PV(_iters);
+//          DAI_PV(props.maxiter);
+            while( sentMessages.size() > 0 && _iters < props.maxiter ) {
+//              DAI_PV(sentMessages.size());
+//              DAI_PV(_iters);
+                _iters++;
+                pair<size_t, size_t> e = sentMessages.back();
+                sentMessages.pop_back();
+                size_t i = e.first, _I = e.second;
+                sendSeqMsgM( i, _I );
+            }
+            if( _iters >= props.maxiter )
+                if( props.verbose >= 1 )
+                    cerr << "Warning: BBP updates limited to " << props.maxiter << " iterations, but using UpdateType " << updates << " with " << totalMessages << " messages" << endl;
+            break;
+        } case UT::PAR: {
+            do {
+                _iters++;
+                doParUpdate();
+            } while( (_iters < 2 || getUnMsgMag() > tol) && _iters < props.maxiter );
+            if( _iters == props.maxiter ) {
+                Real s, new_s;
+                getMsgMags( s, new_s );
+                if( props.verbose >= 1 )
+                    cerr << "Warning: BBP didn't converge in " << _iters << " iterations (unnorm message magnitude = " << getUnMsgMag() << ", norm message mags = " << s << " -> " << new_s << ")" << endl;
+            }
+            break;
+        }
+    }
+    if( props.verbose >= 3 )
+        cerr << "BBP::run() took " << toc()-tic << " seconds " << doneIters() << " iterations" << endl;
+}
+
+
+double numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const PropertySet &bbp_props, bbp_cfn_t cfn, double h ) {
+    // calculate the value of the unperturbed cost function
+    Real cf0 = getCostFn( bp, cfn, state );
+
+    // run BBP to estimate adjoints
+    BBP bbp( &bp, bbp_props );
+    initBBPCostFnAdj( bbp, bp, cfn, state );
+    bbp.run();
+
+    Real d = 0;
+    const FactorGraph& fg = bp.fg();
+
+    if( 1 ) {
+        // verify bbp.adj_psi_V
+
+        // for each variable i
+        for( size_t i = 0; i < fg.nrVars(); i++ ) {
+            vector<double> adj_est;
+            // for each value xi
+            for( size_t xi = 0; xi < fg.var(i).states(); xi++ ) {
+                // Clone 'bp' (which may be any InfAlg)
+                InfAlg *bp_prb = bp.clone();
+
+                // perturb it
+                size_t n = bp_prb->fg().var(i).states();
+                Prob psi_1_prb( n, 1.0 );
+                psi_1_prb[xi] += h;
+//                 psi_1_prb.normalize();
+                size_t I = bp_prb->fg().nbV(i)[0]; // use first factor in list of neighbors of i
+                bp_prb->fg().factor(I) *= Factor( bp_prb->fg().var(i), psi_1_prb );
+                
+                // call 'init' on the perturbed variables
+                bp_prb->init( bp_prb->fg().var(i) );
+                
+                // run copy to convergence
+                bp_prb->run();
+
+                // calculate new value of cost function
+                Real cf_prb = getCostFn( *bp_prb, cfn, state );
+
+                // use to estimate adjoint for i
+                adj_est.push_back( (cf_prb - cf0) / h );
+                
+                // free cloned InfAlg
+                delete bp_prb;
+            }
+            Prob p_adj_est( adj_est.begin(), adj_est.end() );
+            // compare this numerical estimate to the BBP estimate; sum the distances
+            cout << "i: " << i
+                 << ", p_adj_est: " << p_adj_est
+                 << ", bbp.adj_psi_V(i): " << bbp.adj_psi_V(i) << endl;
+            d += dist( p_adj_est, bbp.adj_psi_V(i), Prob::DISTL1 );
+        }
+    }
+    /*    if(1) {
+        // verify bbp.adj_n and bbp.adj_m
+
+        // We actually want to check the responsiveness of objective
+        // function to changes in the final messages. But at the end of a
+        // BBP run, the message adjoints are for the initial messages (and
+        // they should be close to zero, see paper). So this resets the
+        // BBP adjoints to the refer to the desired final messages
+        bbp.RegenerateMessageAdjoints();
+
+        // for each variable i
+        for(size_t i=0; i<bp_dual.nrVars(); i++) {
+            // for each factor I ~ i
+            foreach(size_t I, bp_dual.nbV(i)) {
+                vector<double> adj_n_est;
+                // for each value xi
+                for(size_t xi=0; xi<bp_dual.var(i).states(); xi++) {
+                    BP_dual bp_dual_prb(bp_dual);
+                    // make h-sized change to newMsgN
+                    bp_dual_prb.newMsgN(i,I)[xi] += h;
+                    // recalculate beliefs
+                    bp_dual_prb.CalcBeliefs();
+                    // get cost function value
+                    Real cf_prb = getCostFn(bp_dual_prb, cfn, &state);
+                    // add it to list of adjoints
+                    adj_n_est.push_back((cf_prb-cf0)/h);
+                }
+        
+                vector<double> adj_m_est;
+                // for each value xi
+                for(size_t xi=0; xi<bp_dual.var(i).states(); xi++) {
+                    BP_dual bp_dual_prb(bp_dual);
+                    // make h-sized change to newMsgM
+                    bp_dual_prb.newMsgM(I,i)[xi] += h;
+                    // recalculate beliefs
+                    bp_dual_prb.CalcBeliefs();
+                    // get cost function value
+                    Real cf_prb = getCostFn(bp_dual_prb, cfn, &state);
+                    // add it to list of adjoints
+                    adj_m_est.push_back((cf_prb-cf0)/h);
+                }
+
+                Prob p_adj_n_est(adj_n_est.begin(), adj_n_est.end());
+                // compare this numerical estimate to the BBP estimate; sum the distances
+                cerr << "i: " << i << ", I: " << I
+                     << ", adj_n_est: " << p_adj_n_est
+                     << ", bbp.adj_n(i,I): " << bbp.adj_n(i,I) << endl;
+                d += dist(p_adj_n_est, bbp.adj_n(i,I), Prob::DISTL1);
+
+                Prob p_adj_m_est(adj_m_est.begin(), adj_m_est.end());
+                // compare this numerical estimate to the BBP estimate; sum the distances
+                cerr << "i: " << i << ", I: " << I
+                     << ", adj_m_est: " << p_adj_m_est
+                     << ", bbp.adj_m(I,i): " << bbp.adj_m(I,i) << endl;
+                d += dist(p_adj_m_est, bbp.adj_m(I,i), Prob::DISTL1);
+            }
+        }
+    }
+    */
+    /*    if(0) {
+        // verify bbp.adj_b_V
+        for(size_t i=0; i<bp_dual.nrVars(); i++) {
+            vector<double> adj_b_V_est;
+            // for each value xi
+            for(size_t xi=0; xi<bp_dual.var(i).states(); xi++) {
+                BP_dual bp_dual_prb(bp_dual);
+
+                // make h-sized change to b_1(i)[x_i]
+                bp_dual_prb._beliefs.b1[i][xi] += h;
+
+                // get cost function value
+                Real cf_prb = getCostFn(bp_dual_prb, cfn, &state);
+
+                // add it to list of adjoints
+                adj_b_V_est.push_back((cf_prb-cf0)/h);
+            }
+            Prob p_adj_b_V_est(adj_b_V_est.begin(), adj_b_V_est.end());
+            // compare this numerical estimate to the BBP estimate; sum the distances
+            cerr << "i: " << i
+                 << ", adj_b_V_est: " << p_adj_b_V_est
+                 << ", bbp.adj_b_V(i): " << bbp.adj_b_V(i) << endl;
+            d += dist(p_adj_b_V_est, bbp.adj_b_V(i), Prob::DISTL1);
+        }
+    }
+    */
+
+    // return total of distances
+    return d;
+}
+
+
+bool needGibbsState( bbp_cfn_t cfn ) {
+    switch( (size_t)cfn ) {
+        case bbp_cfn_t::CFN_GIBBS_B:
+        case bbp_cfn_t::CFN_GIBBS_B2:
+        case bbp_cfn_t::CFN_GIBBS_EXP:
+        case bbp_cfn_t::CFN_GIBBS_B_FACTOR:
+        case bbp_cfn_t::CFN_GIBBS_B2_FACTOR:
+        case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR:
+            return true;
+        default:
+            return false;
+    }
+}
+
+
+void initBBPCostFnAdj( BBP &bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const vector<size_t> *stateP ) {
+    const FactorGraph &fg = ia.fg();
+    
+    switch( (size_t)cfn_type ) {
+        case bbp_cfn_t::CFN_BETHE_ENT: {
+            vector<Prob> b1_adj;
+            vector<Prob> b2_adj;
+            vector<Prob> psi1_adj;
+            vector<Prob> psi2_adj;
+            b1_adj.reserve( fg.nrVars() );
+            psi1_adj.reserve( fg.nrVars() );
+            b2_adj.reserve( fg.nrFactors() );
+            psi2_adj.reserve( fg.nrFactors() );
+            for( size_t i = 0; i < fg.nrVars(); i++ ) {
+                size_t dim = fg.var(i).states();
+                int c = fg.nbV(i).size();
+                Prob p(dim,0.0);
+                for( size_t xi = 0; xi < dim; xi++ )
+                    p[xi] = (1 - c) * (1 + log( ia.beliefV(i)[xi] ));
+                b1_adj.push_back( p );
+
+                for( size_t xi = 0; xi < dim; xi++ )
+                    p[xi] = -ia.beliefV(i)[xi];
+                psi1_adj.push_back( p );
+            }
+            for( size_t I = 0; I < fg.nrFactors(); I++ ) {
+                size_t dim = fg.factor(I).states();
+                Prob p( dim, 0.0 );
+                for( size_t xI = 0; xI < dim; xI++ )
+                    p[xI] = 1 + log( ia.beliefF(I)[xI] / fg.factor(I).p()[xI] );
+                b2_adj.push_back( p );
+
+                for( size_t xI = 0; xI < dim; xI++ )
+                    p[xI] = -ia.beliefF(I)[xI] / fg.factor(I).p()[xI];
+                psi2_adj.push_back( p );
+            }
+            bbp.init( b1_adj, b2_adj, psi1_adj, psi2_adj );
+            break;
+        } case bbp_cfn_t::CFN_FACTOR_ENT: {
+            vector<Prob> b2_adj;
+            b2_adj.reserve( fg.nrFactors() );
+            for( size_t I = 0; I < fg.nrFactors(); I++ ) {
+                size_t dim = fg.factor(I).states();
+                Prob p( dim, 0.0 );
+                for( size_t xI = 0; xI < dim; xI++ ) {
+                    double bIxI = ia.beliefF(I)[xI];
+                    if( bIxI < 1.0e-15 )
+                        p[xI] = -1.0e10;
+                    else
+                        p[xI] = 1 + log( bIxI );
+                }
+                b2_adj.push_back(p);
+            }
+            bbp.init( bbp.getZeroAdjV(fg), b2_adj );
+            break;
+        } case bbp_cfn_t::CFN_VAR_ENT: {
+            vector<Prob> b1_adj;
+            b1_adj.reserve( fg.nrVars() );
+            for( size_t i = 0; i < fg.nrVars(); i++ ) {
+                size_t dim = fg.var(i).states();
+                Prob p( dim, 0.0 );
+                for( size_t xi = 0; xi < fg.var(i).states(); xi++ ) {
+                    double bixi = ia.beliefV(i)[xi];
+                    if( bixi < 1.0e-15 )
+                        p[xi] = -1.0e10;
+                    else
+                        p[xi] = 1 + log( bixi );
+                }
+                b1_adj.push_back( p );
+            }
+            bbp.init( b1_adj );
+            break;
+        } case bbp_cfn_t::CFN_GIBBS_B:
+          case bbp_cfn_t::CFN_GIBBS_B2:
+          case bbp_cfn_t::CFN_GIBBS_EXP: {
+            // cost functions that use Gibbs sample, summing over variable marginals
+            vector<size_t> state;
+            if( stateP == NULL )
+                state = getGibbsState( ia, 2*ia.Iterations() );
+            else
+                state = *stateP;
+            assert( state.size() == fg.nrVars() );
+
+            vector<Prob> b1_adj;
+            b1_adj.reserve(fg.nrVars());
+            for( size_t i = 0; i < state.size(); i++ ) {
+                size_t n = fg.var(i).states();
+                Prob delta( n, 0.0 );
+                assert(/*0<=state[i] &&*/ state[i] < n);
+                double b = ia.beliefV(i)[state[i]];
+                switch( (size_t)cfn_type ) {
+                    case bbp_cfn_t::CFN_GIBBS_B:
+                        delta[state[i]] = 1.0; 
+                        break;
+                    case bbp_cfn_t::CFN_GIBBS_B2:
+                        delta[state[i]] = b; 
+                        break;
+                    case bbp_cfn_t::CFN_GIBBS_EXP:
+                        delta[state[i]] = exp(b); 
+                        break;
+                    default: 
+                        DAI_THROW(INTERNAL_ERROR);
+                }
+                b1_adj.push_back( delta );
+            }
+            bbp.init( b1_adj );
+            break;
+        } case bbp_cfn_t::CFN_GIBBS_B_FACTOR:
+          case bbp_cfn_t::CFN_GIBBS_B2_FACTOR:
+          case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR: {
+            // cost functions that use Gibbs sample, summing over factor marginals
+            vector<size_t> state;
+            if( stateP == NULL )
+                state = getGibbsState( ia, 2*ia.Iterations() );
+            else
+                state = *stateP;
+            assert( state.size() == fg.nrVars() );
+
+            vector<Prob> b2_adj;
+            b2_adj.reserve( fg.nrVars() );
+            for( size_t I = 0; I <  fg.nrFactors(); I++ ) {
+                size_t n = fg.factor(I).states();
+                Prob delta( n, 0.0 );
+
+                size_t x_I = getFactorEntryForState( fg, I, state );
+                assert(/*0<=x_I &&*/ x_I < n);
+
+                double b = ia.beliefF(I)[x_I];
+                switch( (size_t)cfn_type ) {
+                    case bbp_cfn_t::CFN_GIBBS_B_FACTOR:
+                        delta[x_I] = 1.0; 
+                        break;
+                    case bbp_cfn_t::CFN_GIBBS_B2_FACTOR:
+                        delta[x_I] = b; 
+                        break;
+                    case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR:
+                        delta[x_I] = exp( b );
+                        break;
+                    default: 
+                        DAI_THROW(INTERNAL_ERROR);
+                }
+                b2_adj.push_back( delta );
+            }
+            bbp.init( bbp.getZeroAdjV(fg), b2_adj );
+            break;
+        } default: 
+            DAI_THROW(UNKNOWN_ENUM_VALUE);
+    }
+}
+
+
+Real getCostFn( const InfAlg &ia, bbp_cfn_t cfn_type, const vector<size_t> *stateP ) {
+    double cf = 0.0;
+    const FactorGraph &fg = ia.fg();
+
+    switch( (size_t)cfn_type ) {
+        case bbp_cfn_t::CFN_BETHE_ENT: // ignores state
+            cf = -ia.logZ();
+            break;
+        case bbp_cfn_t::CFN_VAR_ENT: // ignores state
+            for( size_t i = 0; i < fg.nrVars(); i++ )
+                cf += -ia.beliefV(i).entropy();
+            break;
+        case bbp_cfn_t::CFN_FACTOR_ENT: // ignores state
+            for( size_t I = 0; I < fg.nrFactors(); I++ )
+                cf += -ia.beliefF(I).entropy();
+            break;
+        case bbp_cfn_t::CFN_GIBBS_B:
+        case bbp_cfn_t::CFN_GIBBS_B2:
+        case bbp_cfn_t::CFN_GIBBS_EXP: {
+            assert( stateP != NULL );
+            vector<size_t> state = *stateP;
+            assert( state.size() == fg.nrVars() );
+            for( size_t i = 0; i < fg.nrVars(); i++ ) {
+                double b = ia.beliefV(i)[state[i]];
+                switch( (size_t)cfn_type ) {
+                    case bbp_cfn_t::CFN_GIBBS_B:
+                        cf += b;
+                        break;
+                    case bbp_cfn_t::CFN_GIBBS_B2:
+                        cf += b * b / 2.0;
+                        break;
+                    case bbp_cfn_t::CFN_GIBBS_EXP:
+                        cf += exp( b );
+                        break;
+                    default:
+                        DAI_THROW(INTERNAL_ERROR);
+                }
+            }
+            break;
+        } case bbp_cfn_t::CFN_GIBBS_B_FACTOR:
+          case bbp_cfn_t::CFN_GIBBS_B2_FACTOR:
+          case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR: {
+            assert( stateP != NULL );
+            vector<size_t> state = *stateP;
+            assert( state.size() == fg.nrVars() );
+            for( size_t I = 0; I < fg.nrFactors(); I++ ) {
+                size_t x_I = getFactorEntryForState( fg, I, state );
+                double b = ia.beliefF(I)[x_I];
+                switch( (size_t)cfn_type ) {
+                    case bbp_cfn_t::CFN_GIBBS_B_FACTOR:
+                        cf += b;
+                        break;
+                    case bbp_cfn_t::CFN_GIBBS_B2_FACTOR:
+                        cf += b * b / 2.0;
+                        break;
+                    case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR:
+                        cf += exp( b );
+                        break;
+                    default:
+                        DAI_THROW(INTERNAL_ERROR);
+                }
+            }
+            break;
+        } default: 
+            cerr << "Unknown cost function " << cfn_type << endl;
+            DAI_THROW(UNKNOWN_ENUM_VALUE);
+    }
+    return cf;
+}
+
+
+} // end of namespace dai
+
+
+/* {{{ GENERATED CODE: DO NOT EDIT. Created by 
+    ./scripts/regenerate-properties include/dai/bbp.h src/bbp.cpp 
+*/
+namespace dai {
+
+void BBP::Properties::set(const PropertySet &opts)
+{
+    const std::set<PropertyKey> &keys = opts.allKeys();
+    std::set<PropertyKey>::const_iterator i;
+    bool die=false;
+    for(i=keys.begin(); i!=keys.end(); i++) {
+        if(*i == "verbose") continue;
+        if(*i == "maxiter") continue;
+        if(*i == "tol") continue;
+        if(*i == "damping") continue;
+        if(*i == "updates") continue;
+        if(*i == "clean_updates") continue;
+        cerr << "BBP: Unknown property " << *i << endl;
+        die=true;
+    }
+    if(die) {
+        DAI_THROW(UNKNOWN_PROPERTY_TYPE);
+    }
+    if(!opts.hasKey("verbose")) {
+        cerr << "BBP: Missing property \"verbose\" for method \"BBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("maxiter")) {
+        cerr << "BBP: Missing property \"maxiter\" for method \"BBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("tol")) {
+        cerr << "BBP: Missing property \"tol\" for method \"BBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("damping")) {
+        cerr << "BBP: Missing property \"damping\" for method \"BBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("updates")) {
+        cerr << "BBP: Missing property \"updates\" for method \"BBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("clean_updates")) {
+        cerr << "BBP: Missing property \"clean_updates\" for method \"BBP\"" << endl;
+        die=true;
+    }
+    if(die) {
+        DAI_THROW(NOT_ALL_PROPERTIES_SPECIFIED);
+    }
+    verbose = opts.getStringAs<size_t>("verbose");
+    maxiter = opts.getStringAs<size_t>("maxiter");
+    tol = opts.getStringAs<double>("tol");
+    damping = opts.getStringAs<double>("damping");
+    updates = opts.getStringAs<UpdateType>("updates");
+    clean_updates = opts.getStringAs<bool>("clean_updates");
+}
+PropertySet BBP::Properties::get() const {
+    PropertySet opts;
+    opts.Set("verbose", verbose);
+    opts.Set("maxiter", maxiter);
+    opts.Set("tol", tol);
+    opts.Set("damping", damping);
+    opts.Set("updates", updates);
+    opts.Set("clean_updates", clean_updates);
+    return opts;
+}
+string BBP::Properties::toString() const {
+    stringstream s(stringstream::out);
+    s << "[";
+    s << "verbose=" << verbose << ",";
+    s << "maxiter=" << maxiter << ",";
+    s << "tol=" << tol << ",";
+    s << "damping=" << damping << ",";
+    s << "updates=" << updates << ",";
+    s << "clean_updates=" << clean_updates;
+    s << "]";
+    return s.str();
+}
+} // end of namespace dai
+/* }}} END OF GENERATED CODE */
index 17fc9d8..938fa09 100644 (file)
@@ -483,6 +483,8 @@ void BP::init( const VarSet &ns ) {
 
 
 void BP::updateMessage( size_t i, size_t _I ) {
+    if( recordSentMessages )
+        _sentMessages.push_back(make_pair(i,_I));
     if( props.damping == 0.0 ) {
         message(i,_I) = newMessage(i,_I);
         if( props.updates == Properties::UpdateType::SEQMAX )
diff --git a/src/bp_dual.cpp b/src/bp_dual.cpp
new file mode 100644 (file)
index 0000000..44014fe
--- /dev/null
@@ -0,0 +1,162 @@
+/*  Copyright (C) 2009  Frederik Eaton [frederik at ofb dot net]
+
+    This file is part of libDAI.
+
+    libDAI is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    libDAI is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with libDAI; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+
+#include <iostream>
+#include <sstream>
+#include <algorithm>
+
+#include <dai/bp_dual.h>
+#include <dai/util.h>
+#include <dai/bipgraph.h>
+
+
+namespace dai {
+
+
+using namespace std;
+
+
+typedef BipartiteGraph::Neighbor Neighbor;
+
+
+void BP_dual::init() {
+    regenerateMessages();
+    regenerateBeliefs();
+    calcMessages();
+    calcBeliefs();
+}
+
+
+void BP_dual::regenerateMessages() {
+    size_t nv = fg().nrVars();
+    _msgs.Zn.resize(nv);
+    _msgs.Zm.resize(nv);
+    _msgs.m.resize(nv);
+    _msgs.n.resize(nv);
+    for( size_t i = 0; i < nv; i++ ) {
+        size_t nvf = fg().nbV(i).size();
+        _msgs.Zn[i].resize(nvf, 1.0);
+        _msgs.Zm[i].resize(nvf, 1.0);
+        size_t states = fg().var(i).states();
+        _msgs.n[i].resize(nvf, Prob(states));
+        _msgs.m[i].resize(nvf, Prob(states));
+    }
+}
+
+
+void BP_dual::regenerateBeliefs() {
+    _beliefs.b1.clear();
+    _beliefs.b1.reserve(fg().nrVars());
+    _beliefs.Zb1.resize(fg().nrVars(), 1.0);
+    _beliefs.b2.clear();
+    _beliefs.b2.reserve(fg().nrFactors());
+    _beliefs.Zb2.resize(fg().nrFactors(), 1.0);
+
+    for( size_t i = 0; i < fg().nrVars(); i++ )
+        _beliefs.b1.push_back( Prob( fg().var(i).states() ) );
+    for( size_t I = 0; I < fg().nrFactors(); I++ )
+        _beliefs.b2.push_back( Prob( fg().factor(I).states() ) );
+}
+
+
+void BP_dual::calcMessages() {
+    // calculate 'n' messages from "factor marginal / factor"
+    for( size_t I = 0; I < fg().nrFactors(); I++ ) {
+        Factor f = _ia->beliefF(I) / fg().factor(I);
+        foreach( const Neighbor &i, fg().nbF(I) )
+            msgN(i, i.dual) = f.marginal( fg().var(i) ).p();
+    }
+    // calculate 'm' messages and normalizers from 'n' messages
+    for( size_t i = 0; i < fg().nrVars(); i++ )
+        foreach( const Neighbor &I, fg().nbV(i) )
+            calcNewM( i, I.iter );
+    // recalculate 'n' messages and normalizers from 'm' messages
+    for( size_t i = 0; i < fg().nrVars(); i++ )
+        foreach( const Neighbor &I, fg().nbV(i) )
+            calcNewN(i, I.iter);
+}
+
+
+void BP_dual::calcNewM( size_t i, size_t _I ) {
+    // calculate updated message I->i
+    const Neighbor &I = fg().nbV(i)[_I];
+    Prob prod( fg().factor(I).p() );
+    foreach( const Neighbor &j, fg().nbF(I) )
+        if( j != i ) { // for all j in I \ i
+            Prob &n = msgN(j,j.dual);
+            IndexFor ind( fg().var(j), fg().factor(I).vars() );
+            for( size_t x = 0; ind >= 0; x++, ++ind )
+                prod[x] *= n[ind];
+        }
+    // Marginalize onto i
+    Prob marg( fg().var(i).states(), 0.0 );
+    // ind is the precalculated Index(i,I) i.e. to x_I == k corresponds x_i == ind[k]
+    IndexFor ind( fg().var(i), fg().factor(I).vars() );
+    for( size_t x = 0; ind >= 0; x++, ++ind )
+        marg[ind] += prod[x];
+    
+    _msgs.Zm[i][_I] = marg.normalize();
+    _msgs.m[i][_I] = marg;
+}
+
+
+void BP_dual::calcNewN( size_t i, size_t _I ) {
+    // calculate updated message i->I
+    const Neighbor &I = fg().nbV(i)[_I];
+    Prob prod( fg().var(i).states(), 1.0 );
+    foreach( const Neighbor &J, fg().nbV(i) )
+        if( J.node != I.node ) // for all J in i \ I
+            prod *= msgM(i,J.iter);
+    _msgs.Zn[i][_I] = prod.normalize();
+    _msgs.n[i][_I] = prod;
+}
+
+
+void BP_dual::calcBeliefs() {
+    for( size_t i = 0; i < fg().nrVars(); i++ )
+        calcBeliefV(i);  // calculate b_i
+    for( size_t I = 0; I < fg().nrFactors(); I++ )
+        calcBeliefF(I);  // calculate b_I
+}
+
+
+void BP_dual::calcBeliefV( size_t i ) {
+    Prob prod( fg().var(i).states(), 1.0 );
+    foreach( const Neighbor &I, fg().nbV(i) )
+        prod *= msgM(i,I.iter);
+    _beliefs.Zb1[i] = prod.normalize();
+    _beliefs.b1[i] = prod;
+}
+
+
+void BP_dual::calcBeliefF( size_t I ) {
+    Prob prod( fg().factor(I).p() );
+    foreach( const Neighbor &j, fg().nbF(I) ) {
+        IndexFor ind( fg().var(j), fg().factor(I).vars() );
+        Prob n( msgN(j,j.dual) );
+        for( size_t x = 0; ind >= 0; x++, ++ind )
+            prod[x] *= n[ind];
+    }
+    _beliefs.Zb2[I] = prod.normalize();
+    _beliefs.b2[I] = prod;
+}
+
+
+} // end of namespace dai
diff --git a/src/cbp.cpp b/src/cbp.cpp
new file mode 100644 (file)
index 0000000..a5b82f0
--- /dev/null
@@ -0,0 +1,673 @@
+/*  Copyright (C) 2009  Frederik Eaton [frederik at ofb dot net]
+
+    This file is part of libDAI.
+
+    libDAI is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    libDAI is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with libDAI; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+
+#include <iostream>
+#include <sstream>
+#include <map>
+#include <set>
+#include <algorithm>
+
+#include <dai/util.h>
+#include <dai/properties.h>
+
+#include <dai/bp.h>
+#include <dai/cbp.h>
+#include <dai/bbp.h>
+
+
+namespace dai {
+
+
+using namespace std;
+using boost::shared_ptr;
+
+
+const char *CBP::Name = "CBP";
+
+
+void CBP::setBeliefs( const std::vector<Factor> &bs, double logZ ) {
+    size_t i = 0;
+    _beliefsV.clear(); 
+    _beliefsV.reserve( nrVars() );
+    _beliefsF.clear(); 
+    _beliefsF.reserve( nrFactors() );
+    for( i = 0; i < nrVars(); i++ ) 
+        _beliefsV.push_back( bs[i] );
+    for( ; i < nrVars() + nrFactors(); i++ ) 
+        _beliefsF.push_back( bs[i] );
+    _logZ = logZ;
+}
+
+
+void CBP::construct() {
+    _beliefsV.clear();
+    _beliefsV.reserve(nrVars());
+    for( size_t i = 0; i < nrVars(); i++ )
+        _beliefsV.push_back( Factor(var(i)).normalized() );
+
+    _beliefsF.clear();
+    _beliefsF.reserve(nrFactors());
+    for( size_t I = 0; I < nrFactors(); I++ ) {
+        Factor f = factor(I);
+        f.fill(1); f.normalize();
+        _beliefsF.push_back( f );
+    }
+
+    // to compute average level
+    _sum_level = 0;
+    _num_leaves = 0;
+
+    _maxdiff = 0;
+    _iters = 0;
+
+    if( props.clamp_outfile.length() > 0 ) {
+        _clamp_ofstream = shared_ptr<ofstream>(new ofstream( props.clamp_outfile.c_str(), ios_base::out|ios_base::trunc ));
+        *_clamp_ofstream << "# COUNT LEVEL VAR STATE" << endl;
+    }
+}
+
+
+/// Calculates a vector of mixtures p * b + (1-p) * c
+static vector<Factor> mixBeliefs( Real p, const vector<Factor> &b, const vector<Factor> &c ) {
+    vector<Factor> out;
+    assert( b.size() == c.size() );
+    out.reserve( b.size() );
+    Real pc = 1 - p;
+    for( size_t i = 0; i < b.size(); i++ )
+        // probably already normalized, but do it again just in case
+        out.push_back( b[i].normalized() * p + c[i].normalized() * pc );
+    return out;
+}
+
+
+double CBP::run() {
+    size_t seed = props.rand_seed;
+    if( seed > 0) 
+        rnd_seed( seed );
+
+    InfAlg *bp = getInfAlg();
+    bp->init();
+    bp->run();
+    _iters += bp->Iterations();
+
+    vector<Factor> beliefs_out;
+    Real lz_out;
+    size_t choose_count=0;
+    runRecurse( bp, bp->logZ(), vector<size_t>(0), _num_leaves, choose_count, _sum_level, lz_out, beliefs_out );
+    if( props.verbose >= 1 )
+        cerr << "CBP average levels = " << (_sum_level / _num_leaves) << ", leaves = " << _num_leaves << endl;
+    setBeliefs( beliefs_out, lz_out );
+    return 0.0;
+}
+
+
+/// \todo Eventually this allow other inference algorithms that BP to be selected, via a property
+InfAlg* CBP::getInfAlg() {
+    PropertySet bpProps;
+    bpProps.Set("updates", props.updates);
+    bpProps.Set("tol", props.tol);
+    bpProps.Set("maxiter", props.maxiter);
+    bpProps.Set("verbose", props.verbose);
+    bpProps.Set("logdomain", false);
+    bpProps.Set("damping", 0.0);
+    BP *bp = new BP( *this, bpProps );
+    bp->recordSentMessages = true;
+    bp->init();
+    return bp;
+}
+
+
+void CBP::runRecurse( InfAlg *bp, double orig_logZ, vector<size_t> clamped_vars_list, size_t &num_leaves,
+                      size_t &choose_count, double &sum_level, Real &lz_out, vector<Factor>& beliefs_out) {
+    // choose a variable/states to clamp:
+    size_t i;
+    vector<size_t> xis;
+    Real maxVar = 0.0;
+    bool found;
+    bool clampingVar = (Clamping() == Properties::ClampType::CLAMP_VAR);
+
+    if( Recursion() == Properties::RecurseType::REC_LOGZ && recTol() > 0 && exp( bp->logZ() - orig_logZ ) < recTol() )
+        found = false;
+    else
+        found = chooseNextClampVar( bp, clamped_vars_list, i, xis, &maxVar );
+
+    if( !found ) {
+        num_leaves++;
+        sum_level += clamped_vars_list.size();
+        beliefs_out = bp->beliefs();
+        lz_out = bp->logZ();
+        return;
+    }
+
+    choose_count++;
+    if( props.clamp_outfile.length() > 0 )
+        *_clamp_ofstream << choose_count << "\t" << clamped_vars_list.size() << "\t" << i << "\t" << xis[0] << endl;
+    
+    if( clampingVar )
+        foreach( size_t xi, xis ) 
+            assert(/*0<=xi &&*/ xi < var(i).states() );
+    else
+        foreach( size_t xI, xis ) 
+            assert(/*0<=xI &&*/ xI < factor(i).states() );
+    // - otherwise, clamp and recurse, saving margin estimates for each
+    // clamp setting. afterwards, combine estimates.
+
+    // compute complement of 'xis'
+    vector<size_t> cmp_xis = complement( xis, clampingVar ? var(i).states() : factor(i).states() );
+
+    /// \todo could do this more efficiently with a nesting version of saveProbs/undoProbs
+    Real lz; 
+    vector<Factor> b;
+    InfAlg *bp_c = bp->clone();
+    if( clampingVar ) {
+        bp_c->fg().clampVar( i, xis );
+        bp_c->init( var(i) );
+    } else {
+        bp_c->fg().clampFactor( i, xis );
+        bp_c->init( factor(i).vars() );
+    }
+    bp_c->run();
+    _iters += bp_c->Iterations();
+
+    lz = bp_c->logZ();
+    b = bp_c->beliefs();
+
+    Real cmp_lz; 
+    vector<Factor> cmp_b;
+    InfAlg *cmp_bp_c = bp->clone();
+    if( clampingVar ) {
+        cmp_bp_c->fg().clampVar( i, cmp_xis );
+        cmp_bp_c->init(var(i));
+    } else {
+        cmp_bp_c->fg().clampFactor( i, cmp_xis );
+        cmp_bp_c->init( factor(i).vars() );
+    }
+    cmp_bp_c->run();
+    _iters += cmp_bp_c->Iterations();
+
+    cmp_lz = cmp_bp_c->logZ();
+    cmp_b = cmp_bp_c->beliefs();
+
+    double p = unSoftMax( lz, cmp_lz );
+    Real bp__d = 0.0;
+    
+    if( Recursion() == Properties::RecurseType::REC_BDIFF && recTol() > 0 ) {
+        vector<Factor> combined_b( mixBeliefs( p, b, cmp_b ) );
+        Real new_lz = logSumExp( lz,cmp_lz );
+        bp__d = dist( bp->beliefs(), combined_b, nrVars() );
+        if( exp( new_lz - orig_logZ) * bp__d < recTol() ) {
+            num_leaves++;
+            sum_level += clamped_vars_list.size();
+            beliefs_out = combined_b;
+            lz_out = new_lz;
+            return;
+        }
+    }
+
+    // either we are not doing REC_BDIFF or the distance was large
+    // enough to recurse:
+    runRecurse( bp_c, orig_logZ, clamped_vars_list, num_leaves, choose_count, sum_level, lz, b );
+    runRecurse( cmp_bp_c, orig_logZ, clamped_vars_list, num_leaves, choose_count, sum_level, cmp_lz, cmp_b );
+
+    p = unSoftMax( lz, cmp_lz );
+
+    beliefs_out = mixBeliefs( p, b, cmp_b );
+    lz_out = logSumExp( lz, cmp_lz );
+
+    if( props.verbose >= 2 ) {
+        Real d = dist( bp->beliefs(), beliefs_out, nrVars() );
+        cerr << "Distance (clamping " << i << "): " << d;
+        if( Recursion() == Properties::RecurseType::REC_BDIFF )
+            cerr << "; bp_dual predicted " << bp__d;
+        cerr << "; max_adjoint = " << maxVar << "; logZ = " << lz_out << " (in " << bp->logZ() << ") (orig " << orig_logZ << "); p = " << p << "; level = " << clamped_vars_list.size() << endl;
+    }
+
+    delete bp_c;
+    delete cmp_bp_c;
+}
+
+
+// 'xis' must be sorted
+bool CBP::chooseNextClampVar( InfAlg *bp, vector<size_t> &clamped_vars_list, size_t &i, vector<size_t> &xis, Real *maxVarOut ) {
+    Real tiny = 1.0e-14;
+    if( props.verbose >= 3 )
+        cerr << "clamped_vars_list" << clamped_vars_list << endl;
+    if( clamped_vars_list.size() >= maxClampLevel() )
+        return false;
+    if( ChooseMethod() == Properties::ChooseMethodType::CHOOSE_RANDOM ) {
+        if( Clamping() == Properties::ClampType::CLAMP_VAR ) {
+            int t = 0, t1 = 100;
+            do {
+                i = rnd( nrVars() );
+                t++;
+            } while( abs( bp->beliefV(i).p().max() - 1) < tiny && t < t1 );
+            if( t == t1 ) {
+                return false;
+                // die("Too many levels requested in CBP");
+            }
+            // only pick probable values for variable
+            size_t xi;
+            do {
+                xi = rnd( var(i).states() );
+                t++;
+            } while( bp->beliefV(i).p()[xi] < tiny && t < t1 );
+            assert( t < t1 );
+            xis.resize( 1, xi );
+            //         assert(!_clamped_vars.count(i)); // not true for >2-ary variables
+            DAI_IFVERB(2, "CHOOSE_RANDOM at level " << clamped_vars_list.size() << " chose variable " << i << " state " << xis[0] << endl);
+        } else {
+            int t = 0, t1 = 100;
+            do {
+                i = rnd( nrFactors() );
+                t++;
+            } while( abs( bp->beliefF(i).p().max() - 1) < tiny && t < t1 );
+            if( t == t1 )
+                return false;
+                // die("Too many levels requested in CBP");
+            // only pick probable values for variable
+            size_t xi;
+            do {
+                xi = rnd( factor(i).states() );
+                t++;
+            } while( bp->beliefF(i).p()[xi] < tiny && t < t1 );
+            assert( t < t1 );
+            xis.resize( 1, xi );
+            //         assert(!_clamped_vars.count(i)); // not true for >2-ary variables
+            DAI_IFVERB(2, endl<<"CHOOSE_RANDOM chose factor "<<i<<" state "<<xis[0]<<endl);
+        }
+    } else if( ChooseMethod() == Properties::ChooseMethodType::CHOOSE_MAXENT ) {
+        if( Clamping() == Properties::ClampType::CLAMP_VAR ) {
+            Real max_ent = -1.0;
+            int win_k = -1, win_xk = -1;
+            for( size_t k = 0; k < nrVars(); k++ ) {
+                Real ent=bp->beliefV(k).entropy();
+                if( max_ent < ent ) {
+                    max_ent = ent;
+                    win_k = k;
+                    win_xk = bp->beliefV(k).p().argmax().first;
+                }
+            }
+            assert( win_k >= 0 ); 
+            assert( win_xk >= 0 );
+            i = win_k; 
+            xis.resize( 1, win_xk );
+            DAI_IFVERB(2, endl<<"CHOOSE_MAXENT chose variable "<<i<<" state "<<xis[0]<<endl);
+            if( bp->beliefV(i).p()[xis[0]] < tiny ) {
+                DAI_IFVERB(2, "Warning: CHOOSE_MAXENT found unlikely state, not recursing");
+                return false;
+            }
+        } else {
+            Real max_ent = -1.0;
+            int win_k = -1, win_xk = -1;
+            for( size_t k = 0; k < nrFactors(); k++ ) {
+                Real ent = bp->beliefF(k).entropy();
+                if( max_ent < ent ) {
+                    max_ent = ent;
+                    win_k = k;
+                    win_xk = bp->beliefF(k).p().argmax().first;
+                }
+            }
+            assert( win_k >= 0 ); 
+            assert( win_xk >= 0 );
+            i = win_k; 
+            xis.resize( 1, win_xk );
+            DAI_IFVERB(2, endl<<"CHOOSE_MAXENT chose factor "<<i<<" state "<<xis[0]<<endl);
+            if( bp->beliefF(i).p()[xis[0]] < tiny ) {
+                DAI_IFVERB(2, "Warning: CHOOSE_MAXENT found unlikely state, not recursing");
+                return false;
+            }
+        }
+    } else if( ChooseMethod()==Properties::ChooseMethodType::CHOOSE_BP_L1 ||
+               ChooseMethod()==Properties::ChooseMethodType::CHOOSE_BP_CFN ) {
+        bool doL1 = (ChooseMethod() == Properties::ChooseMethodType::CHOOSE_BP_L1);
+        vector<size_t> state;
+        if( !doL1 && needGibbsState( BBP_cost_function() ) )
+            state = getGibbsState( *bp, 2*bp->Iterations() );
+        // try clamping each variable manually
+        assert( Clamping() == Properties::ClampType::CLAMP_VAR );
+        Real max_cost = 0.0;
+        int win_k = -1, win_xk = -1;
+        for( size_t k = 0; k < nrVars(); k++ ) {
+            for( size_t xk = 0; xk < var(k).states(); xk++ ) {
+                if( bp->beliefV(k)[xk] < tiny ) 
+                    continue;
+                InfAlg *bp1 = bp->clone();
+                bp1->clamp( var(k), xk );
+                bp1->init( var(k) );
+                bp1->run();
+                Real cost = 0;
+                if( doL1 )
+                    for( size_t j = 0; j < nrVars(); j++ )
+                        cost += dist( bp->beliefV(j), bp1->beliefV(j), Prob::DISTL1 );
+                else
+                    cost = getCostFn( *bp1, BBP_cost_function(), &state );
+                if( cost > max_cost || win_k == -1 ) {
+                    max_cost = cost;
+                    win_k = k;
+                    win_xk = xk;
+                }
+                delete bp1;
+            }
+        }
+        assert( win_k >= 0 ); 
+        assert( win_xk >= 0 );
+        i = win_k; 
+        xis.resize( 1, win_xk );
+    } else if( ChooseMethod() == Properties::ChooseMethodType::CHOOSE_BBP ) {
+        Real mvo; 
+        if( !maxVarOut ) 
+            maxVarOut = &mvo;
+        bool clampingVar = (Clamping() == Properties::ClampType::CLAMP_VAR);
+        pair<size_t, size_t> cv = bbpFindClampVar( *bp, clampingVar, props.bbp_props, BBP_cost_function(), &mvo );
+
+        // if slope isn't big enough then don't clamp
+        if( mvo < minMaxAdj() )
+            return false;
+
+        size_t xi = cv.second;
+        i = cv.first;
+#define VAR_INFO (clampingVar?"variable ":"factor ")                       \
+            << i << " state " << xi                                     \
+            << " (p=" << (clampingVar?bp->beliefV(i)[xi]:bp->beliefF(i)[xi]) \
+            << ", entropy = " << (clampingVar?bp->beliefV(i):bp->beliefF(i)).entropy() \
+            << ", maxVar = "<< mvo << ")" 
+        Prob b = ( clampingVar ? bp->beliefV(i).p() : bp->beliefF(i).p());
+        if( b[xi] < tiny ) {
+            cerr << "Warning, at level " << clamped_vars_list.size() << ", bbpFindClampVar found unlikely " << VAR_INFO << endl;
+            return false;
+        }
+        if( abs(b[xi] - 1) < tiny ) {
+            cerr << "Warning at level " << clamped_vars_list.size() << ", bbpFindClampVar found overly likely " << VAR_INFO << endl;
+            return false;
+        }
+
+        xis.resize( 1, xi );
+        if( clampingVar )
+            assert(/*0<=xi &&*/ xi < var(i).states() );
+        else
+            assert(/*0<=xi &&*/ xi < factor(i).states() );
+        DAI_IFVERB(2, "CHOOSE_BBP (num clamped = " << clamped_vars_list.size() << ") chose " << i << " state " << xi << endl);
+    } else
+        DAI_THROW(INTERNAL_ERROR);
+    clamped_vars_list.push_back( i );
+    return true;
+}
+
+
+void CBP::printDebugInfo() {
+    DAI_PV(_beliefsV);
+    DAI_PV(_beliefsF);
+    DAI_PV(_logZ);
+}
+
+
+//----------------------------------------------------------------
+
+
+/** Function which takes a factor graph as input, runs Gibbs and BP_dual,
+ *  creates and runs a BBP object, finds best variable, returns
+ *  (variable,state) pair for clamping
+ */
+pair<size_t, size_t> bbpFindClampVar( const InfAlg &in_bp, bool clampingVar, const PropertySet &bbp_props, bbp_cfn_t cfn, Real *maxVarOut ) {
+    BBP bbp( &in_bp, bbp_props );
+    initBBPCostFnAdj( bbp, in_bp, cfn, NULL );
+    bbp.run();
+  
+    // find and return the (variable,state) with the largest adj_psi_V
+    size_t argmax_var = 0;
+    size_t argmax_var_state = 0;
+    Real max_var = 0;
+    if( clampingVar ) {
+        for( size_t i = 0; i < in_bp.fg().nrVars(); i++ ) {
+            Prob adj_psi_V = bbp.adj_psi_V(i);
+            if(0) {
+                // helps to account for amount of movement possible in variable
+                // i's beliefs? seems not..
+                adj_psi_V *= in_bp.beliefV(i).entropy();
+            }
+            if(0){
+//                 adj_psi_V *= Prob(in_bp.fg().var(i).states(),1.0)-in_bp.beliefV(i).p();
+                adj_psi_V *= in_bp.beliefV(i).p();
+            }
+            // try to compensate for effect on same variable (doesn't work)
+            //     adj_psi_V[gibbs.state()[i]] -= bp_dual.beliefV(i)[gibbs.state()[i]]/10;
+            pair<size_t,Real> argmax_state = adj_psi_V.argmax();
+
+            if( i == 0 || argmax_state.second > max_var ) {
+                argmax_var = i;
+                max_var = argmax_state.second;
+                argmax_var_state = argmax_state.first;
+            }
+        }
+        assert(/*0 <= argmax_var_state &&*/
+               argmax_var_state < in_bp.fg().var(argmax_var).states() );
+    } else {
+        for( size_t I = 0; I < in_bp.fg().nrFactors(); I++ ) {
+            Prob adj_psi_F = bbp.adj_psi_F(I);
+            if(0) {
+                // helps to account for amount of movement possible in variable
+                // i's beliefs? seems not..
+                adj_psi_F *= in_bp.beliefF(I).entropy();
+            }
+            // try to compensate for effect on same variable (doesn't work)
+            //     adj_psi_V[gibbs.state()[i]] -= bp_dual.beliefV(i)[gibbs.state()[i]]/10;
+            pair<size_t,Real> argmax_state = adj_psi_F.argmax();
+
+            if( I == 0 || argmax_state.second > max_var ) {
+                argmax_var = I;
+                max_var = argmax_state.second;
+                argmax_var_state = argmax_state.first;
+            }
+        }
+        assert(/*0 <= argmax_var_state &&*/
+               argmax_var_state < in_bp.fg().factor(argmax_var).states() );
+    }
+    if( maxVarOut ) 
+        *maxVarOut = max_var;
+    return make_pair( argmax_var, argmax_var_state );
+}
+
+
+vector<size_t> complement( vector<size_t> &xis, size_t n_states ) {
+    vector<size_t> cmp_xis( 0 );
+    size_t j = 0;
+    for( size_t xi = 0; xi < n_states; xi++ ) {
+        while( j < xis.size() && xis[j] < xi )
+            j++;
+        if( j >= xis.size() || xis[j] > xi )
+            cmp_xis.push_back(xi);
+    }
+    assert( xis.size()+cmp_xis.size() == n_states );
+    return cmp_xis;
+}
+
+
+Real unSoftMax( Real a, Real b ) {
+    if( a > b )
+        return 1.0 / (1.0 + exp(b-a));
+    else
+        return exp(a-b) / (exp(a-b) + 1.0);
+}
+
+
+Real logSumExp( Real a, Real b ) {
+    if( a > b )
+        return a + log1p( exp( b-a ) );
+    else
+        return b + log1p( exp( a-b ) );
+}
+
+
+Real dist( const vector<Factor> &b1, const vector<Factor> &b2, size_t nv ) {
+    Real d = 0.0;
+    for( size_t k = 0; k < nv; k++ )
+        d += dist( b1[k], b2[k], Prob::DISTLINF );
+    return d;
+}
+
+
+} // end of namespace dai
+
+
+/* {{{ GENERATED CODE: DO NOT EDIT. Created by 
+    ./scripts/regenerate-properties include/dai/cbp.h src/cbp.cpp 
+*/
+namespace dai {
+
+void CBP::Properties::set(const PropertySet &opts)
+{
+    const std::set<PropertyKey> &keys = opts.allKeys();
+    std::set<PropertyKey>::const_iterator i;
+    bool die=false;
+    for(i=keys.begin(); i!=keys.end(); i++) {
+        if(*i == "verbose") continue;
+        if(*i == "tol") continue;
+        if(*i == "updates") continue;
+        if(*i == "maxiter") continue;
+        if(*i == "rec_tol") continue;
+        if(*i == "max_levels") continue;
+        if(*i == "min_max_adj") continue;
+        if(*i == "choose") continue;
+        if(*i == "recursion") continue;
+        if(*i == "clamp") continue;
+        if(*i == "bbp_props") continue;
+        if(*i == "bbp_cfn") continue;
+        if(*i == "rand_seed") continue;
+        if(*i == "clamp_outfile") continue;
+        cerr << "CBP: Unknown property " << *i << endl;
+        die=true;
+    }
+    if(die) {
+        DAI_THROW(UNKNOWN_PROPERTY_TYPE);
+    }
+    if(!opts.hasKey("tol")) {
+        cerr << "CBP: Missing property \"tol\" for method \"CBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("updates")) {
+        cerr << "CBP: Missing property \"updates\" for method \"CBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("maxiter")) {
+        cerr << "CBP: Missing property \"maxiter\" for method \"CBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("rec_tol")) {
+        cerr << "CBP: Missing property \"rec_tol\" for method \"CBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("min_max_adj")) {
+        cerr << "CBP: Missing property \"min_max_adj\" for method \"CBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("choose")) {
+        cerr << "CBP: Missing property \"choose\" for method \"CBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("recursion")) {
+        cerr << "CBP: Missing property \"recursion\" for method \"CBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("clamp")) {
+        cerr << "CBP: Missing property \"clamp\" for method \"CBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("bbp_props")) {
+        cerr << "CBP: Missing property \"bbp_props\" for method \"CBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("bbp_cfn")) {
+        cerr << "CBP: Missing property \"bbp_cfn\" for method \"CBP\"" << endl;
+        die=true;
+    }
+    if(die) {
+        DAI_THROW(NOT_ALL_PROPERTIES_SPECIFIED);
+    }
+    if(opts.hasKey("verbose")) {
+        verbose = opts.getStringAs<size_t>("verbose");
+    } else {
+        verbose = 0;
+    }
+    tol = opts.getStringAs<double>("tol");
+    updates = opts.getStringAs<UpdateType>("updates");
+    maxiter = opts.getStringAs<size_t>("maxiter");
+    rec_tol = opts.getStringAs<double>("rec_tol");
+    if(opts.hasKey("max_levels")) {
+        max_levels = opts.getStringAs<size_t>("max_levels");
+    } else {
+        max_levels = 10;
+    }
+    min_max_adj = opts.getStringAs<double>("min_max_adj");
+    choose = opts.getStringAs<ChooseMethodType>("choose");
+    recursion = opts.getStringAs<RecurseType>("recursion");
+    clamp = opts.getStringAs<ClampType>("clamp");
+    bbp_props = opts.getStringAs<PropertySet>("bbp_props");
+    bbp_cfn = opts.getStringAs<bbp_cfn_t>("bbp_cfn");
+    if(opts.hasKey("rand_seed")) {
+        rand_seed = opts.getStringAs<size_t>("rand_seed");
+    } else {
+        rand_seed = 0;
+    }
+    if(opts.hasKey("clamp_outfile")) {
+        clamp_outfile = opts.getStringAs<std::string>("clamp_outfile");
+    } else {
+        clamp_outfile = "";
+    }
+}
+PropertySet CBP::Properties::get() const {
+    PropertySet opts;
+    opts.Set("verbose", verbose);
+    opts.Set("tol", tol);
+    opts.Set("updates", updates);
+    opts.Set("maxiter", maxiter);
+    opts.Set("rec_tol", rec_tol);
+    opts.Set("max_levels", max_levels);
+    opts.Set("min_max_adj", min_max_adj);
+    opts.Set("choose", choose);
+    opts.Set("recursion", recursion);
+    opts.Set("clamp", clamp);
+    opts.Set("bbp_props", bbp_props);
+    opts.Set("bbp_cfn", bbp_cfn);
+    opts.Set("rand_seed", rand_seed);
+    opts.Set("clamp_outfile", clamp_outfile);
+    return opts;
+}
+string CBP::Properties::toString() const {
+    stringstream s(stringstream::out);
+    s << "[";
+    s << "verbose=" << verbose << ",";
+    s << "tol=" << tol << ",";
+    s << "updates=" << updates << ",";
+    s << "maxiter=" << maxiter << ",";
+    s << "rec_tol=" << rec_tol << ",";
+    s << "max_levels=" << max_levels << ",";
+    s << "min_max_adj=" << min_max_adj << ",";
+    s << "choose=" << choose << ",";
+    s << "recursion=" << recursion << ",";
+    s << "clamp=" << clamp << ",";
+    s << "bbp_props=" << bbp_props << ",";
+    s << "bbp_cfn=" << bbp_cfn << ",";
+    s << "rand_seed=" << rand_seed << ",";
+    s << "clamp_outfile=" << clamp_outfile;
+    s << "]";
+    return s.str();
+}
+} // end of namespace dai
+/* }}} END OF GENERATED CODE */
index 87bbe47..d8b9fea 100644 (file)
@@ -330,6 +330,37 @@ void FactorGraph::clamp( const Var & n, size_t i, bool backup ) {
 }
 
 
+void FactorGraph::clampVar( size_t i, const vector<size_t> &is, bool backup ) {
+    Var n = var(i);
+    Factor mask_n( n, 0.0 );
+
+    foreach( size_t i, is ) {
+        assert( i <= n.states() );
+        mask_n[i] = 1.0;
+    }
+
+    map<size_t, Factor> newFacs;
+    for( size_t I = 0; I < nrFactors(); I++ ) 
+        if( factor(I).vars().contains( n ) ) {
+            newFacs[I] = factor(I) * mask_n;
+        }
+    setFactors( newFacs, backup );
+}
+
+
+void FactorGraph::clampFactor( size_t I, const vector<size_t> &is, bool backup ) {
+    size_t st = factor(I).states();
+    Factor newF( factor(I).vars(), 0.0 );
+
+    foreach( size_t i, is ) { 
+        assert( i <= st ); 
+        newF[i] = factor(I)[i];
+    }
+
+    setFactor( I, newF, backup );
+}
+
+
 void FactorGraph::backupFactor( size_t I ) {
     map<size_t,Factor>::iterator it = _backup.find( I );
     if( it != _backup.end() )
@@ -372,6 +403,7 @@ void FactorGraph::restoreFactors() {
     _backup.clear();
 }
 
+
 void FactorGraph::backupFactors( const std::set<size_t> & facs ) {
     for( std::set<size_t>::const_iterator fac = facs.begin(); fac != facs.end(); fac++ )
         backupFactor( *fac );
index df55257..3bb07fd 100644 (file)
@@ -198,12 +198,12 @@ double Gibbs::run() {
 }
 
 
-inline Factor Gibbs::beliefV( size_t i ) const {
+Factor Gibbs::beliefV( size_t i ) const {
     return Factor( var(i), _var_counts[i].begin() ).normalized();
 }
 
 
-inline Factor Gibbs::beliefF( size_t I ) const {
+Factor Gibbs::beliefF( size_t I ) const {
     return Factor( factor(I).vars(), _factor_counts[I].begin() ).normalized();
 }
 
index f031109..dc9a71b 100644 (file)
@@ -1,5 +1,6 @@
 # --- BP ----------------------
 
+BP:                             BP[updates=SEQFIX,tol=1e-9,maxiter=10000,logdomain=0]
 BP_SEQFIX:                      BP[updates=SEQFIX,tol=1e-9,maxiter=10000,logdomain=0]
 BP_SEQRND:                      BP[updates=SEQRND,tol=1e-9,maxiter=10000,logdomain=0]
 BP_SEQMAX:                      BP[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0]
@@ -100,3 +101,8 @@ GIBBS_1e6:                      GIBBS[iters=1000000]
 GIBBS_1e7:                      GIBBS[iters=10000000]
 GIBBS_1e8:                      GIBBS[iters=100000000]
 GIBBS_1e9:                      GIBBS[iters=1000000000]
+
+# --- CBP ---------------------
+
+CBP:                            CBP[max_levels=12,updates=SEQMAX,tol=1e-9,rec_tol=1e-9,maxiter=500,choose=CHOOSE_RANDOM,recursion=REC_FIXED,clamp=CLAMP_VAR,min_max_adj=1.0e-9,bbp_cfn=CFN_FACTOR_ENT,verbose=0,rand_seed=0,bbp_props=[verbose=0,tol=1.0e-9,maxiter=10000,damping=0,updates=SEQ_BP_REV,clean_updates=0],clamp_outfile=]
+BBP:                            CBP[choose=CHOOSE_BBP]
index e3ada79..7c76ea9 100755 (executable)
@@ -1,2 +1,2 @@
 #!/bin/bash
-./testdai --report-iters false --report-time false --marginals true --aliases aliases.conf --filename $1 --methods EXACT[verbose=0] JTREE_HUGIN JTREE_SHSH BP_SEQFIX BP_SEQRND BP_SEQMAX BP_PARALL BP_SEQFIX_LOG BP_SEQRND_LOG BP_SEQMAX_LOG BP_PARALL_LOG MF_SEQRND TREEEP TREEEPWC GBP_MIN GBP_DELTA GBP_LOOP3 GBP_LOOP4 GBP_LOOP5 GBP_LOOP6 GBP_LOOP7 HAK_MIN HAK_DELTA HAK_LOOP3 HAK_LOOP4 HAK_LOOP5 MR_RESPPROP_FULL MR_CLAMPING_FULL MR_EXACT_FULL MR_RESPPROP_LINEAR MR_CLAMPING_LINEAR MR_EXACT_LINEAR LCBP_FULLCAV_SEQFIX LCBP_FULLCAVin_SEQFIX LCBP_FULLCAV_SEQRND LCBP_FULLCAVin_SEQRND LCBP_FULLCAV_NONE LCBP_FULLCAVin_NONE LCBP_PAIRCAV_SEQFIX LCBP_PAIRCAVin_SEQFIX LCBP_PAIRCAV_SEQRND LCBP_PAIRCAVin_SEQRND LCBP_PAIRCAV_NONE LCBP_PAIRCAVin_NONE LCBP_PAIR2CAV_SEQFIX LCBP_PAIR2CAVin_SEQFIX LCBP_PAIR2CAV_SEQRND LCBP_PAIR2CAVin_SEQRND LCBP_PAIR2CAV_NONE LCBP_PAIR2CAVin_NONE LCBP_UNICAV_SEQFIX LCBP_UNICAV_SEQRND
+./testdai --report-iters false --report-time false --marginals true --aliases aliases.conf --filename $1 --methods EXACT[verbose=0] JTREE_HUGIN JTREE_SHSH BP_SEQFIX BP_SEQRND BP_SEQMAX BP_PARALL BP_SEQFIX_LOG BP_SEQRND_LOG BP_SEQMAX_LOG BP_PARALL_LOG MF_SEQRND TREEEP TREEEPWC GBP_MIN GBP_DELTA GBP_LOOP3 GBP_LOOP4 GBP_LOOP5 GBP_LOOP6 GBP_LOOP7 HAK_MIN HAK_DELTA HAK_LOOP3 HAK_LOOP4 HAK_LOOP5 MR_RESPPROP_FULL MR_CLAMPING_FULL MR_EXACT_FULL MR_RESPPROP_LINEAR MR_CLAMPING_LINEAR MR_EXACT_LINEAR LCBP_FULLCAV_SEQFIX LCBP_FULLCAVin_SEQFIX LCBP_FULLCAV_SEQRND LCBP_FULLCAVin_SEQRND LCBP_FULLCAV_NONE LCBP_FULLCAVin_NONE LCBP_PAIRCAV_SEQFIX LCBP_PAIRCAVin_SEQFIX LCBP_PAIRCAV_SEQRND LCBP_PAIRCAVin_SEQRND LCBP_PAIRCAV_NONE LCBP_PAIRCAVin_NONE LCBP_PAIR2CAV_SEQFIX LCBP_PAIR2CAVin_SEQFIX LCBP_PAIR2CAV_SEQRND LCBP_PAIR2CAVin_SEQRND LCBP_PAIR2CAV_NONE LCBP_PAIR2CAVin_NONE LCBP_UNICAV_SEQFIX LCBP_UNICAV_SEQRND BBP CBP
diff --git a/tests/testbbp.cpp b/tests/testbbp.cpp
new file mode 100644 (file)
index 0000000..8c061cb
--- /dev/null
@@ -0,0 +1,123 @@
+/*  Copyright (C) 2009  Joris Mooij  [joris dot mooij at tuebingen dot mpg dot de]
+    Max Planck Institute for Biological Cybernetics, Germany
+
+    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 <dai/alldai.h>
+#include <dai/bbp.h>
+
+
+using namespace dai;
+using namespace std;
+
+
+int main( int argc, char *argv[] ) {
+    if ( argc != 2 ) {
+        cout << "Usage: " << argv[0] << " <filename.fg>" << endl << endl;
+        cout << "Reads factor graph <filename.fg> and verifies" << endl;
+        cout << "whether BBP works correctly on it." << endl << endl;
+        return 1;
+    } else {
+        // Read FactorGraph from the file specified by the first command line argument
+        FactorGraph fg;
+        fg.ReadFromFile(argv[1]);
+        
+        // Set some constants
+        size_t verbose = 0;
+        double tol = 1.0e-9;
+        size_t maxiter = 10000;
+        double damping = 0.0;
+        BBP::Properties::UpdateType updates = BBP::Properties::UpdateType::PAR;
+        bool   clean_updates = false;
+
+        // Store the constants in a PropertySet object
+        PropertySet opts;
+        opts.Set("verbose",verbose);  // Verbosity (amount of output generated)
+        opts.Set("tol",tol);          // Tolerance for convergence
+        opts.Set("maxiter",maxiter);  // Maximum number of iterations
+        opts.Set("damping",damping);  // Amount of damping applied
+
+        // Construct a BP (belief propagation) object from the FactorGraph fg
+        BP bp(fg, opts("updates",string("SEQFIX"))("logdomain",false));
+        bp.recordSentMessages = true;
+        bp.init();
+        bp.run();
+
+        vector<size_t> state( fg.nrVars(), 0 );
+
+        for( size_t t = 0; t < 90; t++ ) {
+            clean_updates = t % 2;
+            BBP::Properties::UpdateType updates;
+            switch( (t / 2) % 5 ) {
+                case BBP::Properties::UpdateType::SEQ_FIX:
+                    updates = BBP::Properties::UpdateType::SEQ_FIX;
+                    break;
+                case BBP::Properties::UpdateType::SEQ_MAX:
+                    updates = BBP::Properties::UpdateType::SEQ_MAX;
+                    break;
+                case BBP::Properties::UpdateType::SEQ_BP_REV:
+                    updates = BBP::Properties::UpdateType::SEQ_BP_REV;
+                    break;
+                case BBP::Properties::UpdateType::SEQ_BP_FWD:
+                    updates = BBP::Properties::UpdateType::SEQ_BP_FWD;
+                    break;
+                case BBP::Properties::UpdateType::PAR:
+                    updates = BBP::Properties::UpdateType::PAR;
+                    break;
+            }
+            bbp_cfn_t cfn;
+            switch( (t / 10) % 9 ) {
+                case 0:
+                    cfn = bbp_cfn_t::CFN_GIBBS_B;
+                    break;
+                case 1:
+                    cfn = bbp_cfn_t::CFN_GIBBS_B2;
+                    break;
+                case 2:
+                    cfn = bbp_cfn_t::CFN_GIBBS_EXP;
+                    break;
+                case 3:
+                    cfn = bbp_cfn_t::CFN_GIBBS_B_FACTOR;
+                    break;
+                case 4:
+                    cfn = bbp_cfn_t::CFN_GIBBS_B2_FACTOR;
+                    break;
+                case 5:
+                    cfn = bbp_cfn_t::CFN_GIBBS_EXP_FACTOR;
+                    break;
+                case 6:
+                    cfn = bbp_cfn_t::CFN_VAR_ENT;
+                    break;
+                case 7:
+                    cfn = bbp_cfn_t::CFN_FACTOR_ENT;
+                    break;
+                case 8:
+                    cfn = bbp_cfn_t::CFN_BETHE_ENT;
+                    break;
+            }
+
+            double h = 1e-4;
+            double result = numericBBPTest( bp, &state, opts("updates",updates)("clean_updates",clean_updates), cfn, h );
+            cout << "clean_updates=" << clean_updates << ", updates=" << updates << ", cfn=" << cfn << ", result: " << result << endl;
+        }
+    }
+
+    return 0;
+}
index d3d340f..918ec52 100644 (file)
@@ -197,7 +197,7 @@ pair<string, PropertySet> parseMethod( const string &_s, const map<string,string
         if( ps.first == DAINames[n] )
             break;
     if( strlen( DAINames[n] ) == 0 && (ps.first != "LDPC") )
-        throw string("Unknown DAI algorithm \"") + ps.first + string("\" in \"") + _s + string("\"");
+        throw std::runtime_error(string("Unknown DAI algorithm \"") + ps.first + string("\" in \"") + _s + string("\""));
 
     return ps;
 }
index eb632d8..61e0e54 100644 (file)
@@ -884,3 +884,37 @@ LCBP_UNICAV_SEQRND                         9.483e-02       3.078e-02       N/A             1.000e-09
 # ({x13}, (5.266e-01, 4.734e-01))
 # ({x14}, (6.033e-01, 3.967e-01))
 # ({x15}, (1.558e-01, 8.442e-01))
+BBP                                            7.049e-04       2.319e-04       +7.075e-05      1.000e-09       
+# ({x0}, (3.890e-01, 6.110e-01))
+# ({x1}, (5.555e-01, 4.445e-01))
+# ({x2}, (4.587e-01, 5.413e-01))
+# ({x3}, (5.479e-01, 4.521e-01))
+# ({x4}, (6.663e-01, 3.337e-01))
+# ({x5}, (2.105e-01, 7.895e-01))
+# ({x6}, (8.180e-01, 1.820e-01))
+# ({x7}, (2.324e-01, 7.676e-01))
+# ({x8}, (2.169e-01, 7.831e-01))
+# ({x9}, (2.050e-01, 7.950e-01))
+# ({x10}, (7.666e-01, 2.334e-01))
+# ({x11}, (1.216e-01, 8.784e-01))
+# ({x12}, (4.210e-01, 5.790e-01))
+# ({x13}, (5.351e-01, 4.649e-01))
+# ({x14}, (6.284e-01, 3.716e-01))
+# ({x15}, (1.355e-01, 8.645e-01))
+CBP                                            2.160e-05       1.082e-05       +1.720e-06      1.000e-09       
+# ({x0}, (3.888e-01, 6.112e-01))
+# ({x1}, (5.556e-01, 4.444e-01))
+# ({x2}, (4.587e-01, 5.413e-01))
+# ({x3}, (5.480e-01, 4.520e-01))
+# ({x4}, (6.660e-01, 3.340e-01))
+# ({x5}, (2.107e-01, 7.893e-01))
+# ({x6}, (8.178e-01, 1.822e-01))
+# ({x7}, (2.327e-01, 7.673e-01))
+# ({x8}, (2.171e-01, 7.829e-01))
+# ({x9}, (2.052e-01, 7.948e-01))
+# ({x10}, (7.664e-01, 2.336e-01))
+# ({x11}, (1.217e-01, 8.783e-01))
+# ({x12}, (4.214e-01, 5.786e-01))
+# ({x13}, (5.348e-01, 4.652e-01))
+# ({x14}, (6.291e-01, 3.709e-01))
+# ({x15}, (1.357e-01, 8.643e-01))