Merged Generalized Loop Correction code kindly provided by Siamak Ravanbakhsh
authorJoris Mooij <j.mooij@cs.ru.nl>
Fri, 5 Oct 2012 11:03:34 +0000 (13:03 +0200)
committerJoris Mooij <j.mooij@cs.ru.nl>
Fri, 5 Oct 2012 11:03:34 +0000 (13:03 +0200)
20 files changed:
AUTHORS
ChangeLog
Makefile
Makefile.ALL
README
include/dai/alldai.h
include/dai/cobwebgraph.h [new file with mode: 0644]
include/dai/daialg.h
include/dai/doc.h
include/dai/factorgraph.h
include/dai/glc.h [new file with mode: 0644]
include/dai/regiongraph.h
src/alldai.cpp
src/cobwebgraph.cpp [new file with mode: 0644]
src/factorgraph.cpp
src/glc.cpp [new file with mode: 0644]
tests/aliases.conf
tests/testall
tests/testall.bat
tests/testfast.out

diff --git a/AUTHORS b/AUTHORS
index 5e9ba50..bfd1203 100644 (file)
--- a/AUTHORS
+++ b/AUTHORS
@@ -3,6 +3,7 @@ libDAI was written by Joris M. Mooij with the help of the following people:
 Martijn Leisink (laid down the foundations for the library),
 Frederik Eaton (contributed the Gibbs sampler, conditioned BP, fractional BP and various other improvements),
 Charles Vaske (contributed the parameter learning algorithms),
+Siamak Ravanbakhsh (contributed the generalized loop corrections code),
 Giuseppe Passino (contributed the findMaximum function, and various other improvements),
 Bastian Wemmenhove (contributed the MR algorithm),
 Patrick Pletscher (contributed the SWIG interface).
index f92ed71..664db95 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+git HEAD
+--------
+* Merged Generalized Loop Correction code kindly provided by Siamak Ravanbakhsh
+
+
 libDAI-0.3.1 (2012-09-17)
 -------------------------
 
@@ -13,7 +18,7 @@ Build system improvements:
 
 Bug fixes:
 * Fixed bug (found by Sameh Khamis): define NAN (Windows)
-* Fixed bug (found by Same Khamis): add some necessary casts of size_t to BigInt (Windows)
+* Fixed bug (found by Sameh Khamis): add some necessary casts of size_t to BigInt (Windows)
 * Fixed bug (found by Yan): replaced GNU extension __PRETTY_FUNCTION__ by __FUNCTION__ (Visual Studio) or __func__ (other compilers)
 * Fixed bug (found by Andy Mueller): added GMP library invocations to swig Makefile
 * Fixed bug (found by cax): when building MatLab MEX files, GMP/MPIR libraries were not linked
index 76c69c3..d380833 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -12,8 +12,8 @@ include Makefile.ALL
 include Makefile.conf
 
 # Set version and date
-DAI_VERSION="0.3.1"
-DAI_DATE="September 17, 2012"
+DAI_VERSION="git HEAD"
+DAI_DATE="September 17, 2012 - or later"
 
 # Directories of libDAI sources
 # Location of libDAI headers
@@ -41,7 +41,7 @@ ifdef WITH_DOC
 endif
 
 # Define conditional build targets
-NAMES:=graph dag bipgraph varset daialg alldai clustergraph factor factorgraph properties regiongraph util weightedgraph exceptions exactinf evidence emalg io
+NAMES:=graph dag bipgraph varset daialg alldai clustergraph factor factorgraph properties regiongraph cobwebgraph util weightedgraph exceptions exactinf evidence emalg io
 ifdef WITH_BP
   WITHFLAGS:=$(WITHFLAGS) -DDAI_WITH_BP
   NAMES:=$(NAMES) bp
@@ -90,6 +90,11 @@ ifdef WITH_DECMAP
   WITHFLAGS:=$(WITHFLAGS) -DDAI_WITH_DECMAP
   NAMES:=$(NAMES) decmap
 endif
+ifdef WITH_GLC
+  WITHFLAGS:=$(WITHFLAGS) -DDAI_WITH_GLC
+  NAMES:=$(NAMES) glc 
+endif
+
 
 # Define standard libDAI header dependencies, source file names and object file names
 HEADERS=$(foreach name,graph dag bipgraph index var factor varset smallset prob daialg properties alldai enum exceptions util,$(INC)/$(name).h)
@@ -191,6 +196,9 @@ emalg$(OE) : $(SRC)/emalg.cpp $(INC)/emalg.h $(INC)/evidence.h $(HEADERS)
 decmap$(OE) : $(SRC)/decmap.cpp $(INC)/decmap.h $(HEADERS)
        $(CC) -c $<
 
+glc$(OE) : $(SRC)/glc.cpp $(INC)/glc.h $(HEADERS) $(INC)/cobwebgraph.h
+       $(CC) -c $<
+
 
 # EXAMPLES
 ###########
index a26a8c9..cb1c001 100644 (file)
@@ -30,6 +30,7 @@ WITH_MR=true
 WITH_GIBBS=true
 WITH_CBP=true
 WITH_DECMAP=true
+WITH_GLC=true
 
 # Build with debug info? (slower but safer)
 DEBUG=true
diff --git a/README b/README
index 743f1a7..6d39c57 100644 (file)
--- a/README
+++ b/README
@@ -2,8 +2,8 @@ libDAI  -  A free/open source C++ library for Discrete Approximate Inference
 
 -------------------------------------------------------------------------------
 
-Version:  0.3.1
-Date:     September 17, 2012
+Version:  git HEAD
+Date:     September 17, 2012 - or later
 See also: http://www.libdai.org
 
 -------------------------------------------------------------------------------
index f28fc6d..882f8d7 100644 (file)
@@ -61,6 +61,9 @@
 #ifdef DAI_WITH_DECMAP
     #include <dai/decmap.h>
 #endif
+#ifdef DAI_WITH_GLC
+    #include <dai/glc.h>
+#endif
 
 
 /// Namespace for libDAI
diff --git a/include/dai/cobwebgraph.h b/include/dai/cobwebgraph.h
new file mode 100644 (file)
index 0000000..fc68e1b
--- /dev/null
@@ -0,0 +1,279 @@
+/*  This file is part of libDAI - http://www.libdai.org/
+ *
+ *  Copyright (c) 2006-2012, The libDAI authors. All rights reserved.
+ *
+ *  Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
+ */
+
+
+/// \file
+/// \brief Defines class CobwebGraph, which implements a type of region graph used by GLC.
+
+
+#ifndef __defined_libdai_cobwebgraph_h
+#define __defined_libdai_cobwebgraph_h
+
+
+#include <iostream>
+#include <dai/factorgraph.h>
+#include <dai/weightedgraph.h>
+#include <dai/smallset.h>
+#include <dai/enum.h>
+#include <algorithm>
+#include <map>
+#include <set>
+
+
+namespace dai {
+
+
+/// A CobwebGraph is a special type of region graph used by the GLC algorithm
+/** \author Siamak Ravanbakhsh
+ */
+class CobwebGraph : public FactorGraph {
+    protected:
+        /// Vector of variable indices internal to each region (r)
+        std::vector<SmallSet<size_t> >        _INRs;
+        /// Vector of variable indices on the boundary of each region (\ominus r)
+        std::vector<SmallSet<size_t> >        _EXRs;
+        /// Index of factors in each region
+        std::vector<SmallSet<size_t> >        _Rfs;
+        /// Index of factors internal to each region, i.e., all its variables are internal to the region
+        std::vector<SmallSet<size_t> >        _Rifs;
+        /// Index of factors that bridge each region, i.e., not all its variables are internal to the region
+        std::vector<SmallSet<size_t> >        _Rxfs;
+        /// The vector of domain of messages leaving each region (\ominus r_{p,q})
+        std::vector<std::vector<VarSet> >     _outM;
+
+        /// The information in connection between two regions
+        struct Connection {
+            /// Index of the first region (p)
+            size_t my;
+            /// Index of the second region (q)
+            size_t his;
+            /// Index of this connection in the connections of the first region
+            size_t iter;
+            /// Index of the mirror of this connection in the connections of the second region. (reference of this message in _outM)
+            size_t dual;
+            /// The message sent from region q (his) to p (my)
+            Factor msg;
+            /// "Index" of variables in the message
+            std::vector<size_t> varinds;
+            /// Used as a temporary factor only
+            /** \todo Remove CobwebGraph::Connection::newmsg
+             */
+            Factor newmsg;
+            /// Index of factors in common
+            std::vector<size_t> fc;
+            /// Regions rho that are descendents of msg from q to p in p's msg-region graph (not used in partitioning case)
+            std::vector<VarSet> subregions;
+        };
+
+        /** Indicates what happens if a subset of variables in the boundary of a region (\ominus r_p) is shared by
+         *  some neighbors such that one (\ominus r_{p,q1}) is a subset of another (\ominus r_{p,q2}).
+         *  - ALL     all such neighbors are included in the the updates.
+         *  - TOP     only (\ominus r_{p,q2}) is included unless (\ominus r_{p,q2} = \ominus r_{p,q1}) in which case both are included
+         *  - CLOSEST (default): similar to TOP but in case of a tie the the region r_q with largest r_q \cap r_p is considered
+         *  \note Not important in perfomance!
+         */
+        DAI_ENUM(NeighborType,ALL,TOP,CLOSEST);
+
+        /// Vector of all connections to each region
+        std::vector<std::vector<Connection> > _M;
+
+        /// For each i the index of (cobweb) regions that contain variable i
+        std::vector<std::vector<size_t> > _var2CW;
+
+        /** For each region r_p a mapping from all variables rho in its msg-region graph to a pair:
+         *  first: counting number
+         *  second: the index of top-regions that contain rho
+         */
+        std::vector<std::map<VarSet, std::pair<int,std::vector<size_t> > > > _cn;
+
+        /// Whether a given set of region vars makes a partitioning or not
+        bool isPartition;
+
+
+    public:
+    /// \name Constructors and destructors
+    //@{
+        /// Default constructor
+        CobwebGraph() : FactorGraph(), _INRs(), _EXRs(),_Rfs(),_Rifs(),_Rxfs(), _M(), _var2CW(), _cn(), isPartition(true) {}
+
+        /// Constructs a cobweb graph from a factor graph
+        CobwebGraph( const FactorGraph& fg ): FactorGraph(), _INRs(), _EXRs(),_Rfs(), _M(), _var2CW(), _cn(), isPartition(true) {
+            // Copy factor graph structure
+            FactorGraph::operator=( fg );
+        }
+
+        /// Clone \c *this (virtual copy constructor)
+        virtual CobwebGraph* clone() const { return new CobwebGraph(*this); }
+    //@}
+
+    /// \name Accessors and mutators
+    //@{
+        /// Returns the number of regions
+        size_t nrCWs() const { return _INRs.size(); }
+
+        /// Returns constant reference to the _cn for region \a R
+        const std::map<VarSet, std::pair<int,std::vector<size_t> > >& cn( size_t R ) const {
+            DAI_DEBASSERT( R < nrCWs() );
+            return _cn[R];
+        }
+
+        /// Returns a reference to the _cn for region \a R
+        std::map<VarSet, std::pair<int,std::vector<size_t> > >& cn( size_t R ) {
+            DAI_DEBASSERT( R < nrCWs() );
+            return _cn[R];
+        }
+
+        /// Returns constant reference the vector of domain of all outgoing messages from region \a R
+        const std::vector<VarSet>& outM( size_t R ) const {
+            DAI_DEBASSERT( R < _outM.size() );
+            return _outM[R];
+        }
+
+        /// Returns reference the vector of domain of all outgoing messages from region \a R
+        std::vector<VarSet>& outM( size_t R ) {
+            DAI_DEBASSERT( R < _outM.size() );
+            return _outM[R];
+        }
+
+        /// Returns constant reference to the variables in the outgoing message from region \a R to its \a j'th neighbor 
+        /** \a j corresponds to dual in the connection construct
+         */
+        const VarSet& outM( size_t R, size_t j ) const {
+            DAI_DEBASSERT( R < _outM.size() );
+            DAI_DEBASSERT( j < _outM[R].size() );
+            return _outM[R][j];
+        }
+
+        /// Returns a reference to the variables in the outgoing message from region \a R to its \a j'th neighbor 
+        /** \a j corresponds to dual in the connection construct
+         */
+        VarSet& outM( size_t R, size_t j ) {
+            DAI_DEBASSERT( R < _outM.size() );
+            DAI_DEBASSERT( j < _outM[R].size() );
+            return _outM[R][j];
+        }
+
+        /// Returns constant reference to the index of factors of region \a R
+        const SmallSet<size_t>& Rfs( size_t R ) const {
+            DAI_DEBASSERT( R < _Rfs.size() );
+            return _Rfs[R];
+        }
+
+        /// Returns reference to the index of factors of region \a R
+        SmallSet<size_t>& Rfs( size_t R ) {
+            DAI_DEBASSERT( R < _Rfs.size() );
+            return _Rfs[R];
+        }
+
+        /// Returns constant reference to the index of variables on the boundary of region \a R (\ominus r)
+        const SmallSet<size_t>& EXRs( size_t R ) const {
+            DAI_DEBASSERT( R < _EXRs.size() );
+            return _EXRs[R];
+        }
+
+        /// Returns reference to the index of variables on the boundary of region \a R (\ominus r)
+        SmallSet<size_t>& EXRs( size_t R ) {
+            DAI_DEBASSERT( R < _EXRs.size() );
+            return _EXRs[R];
+        }
+
+        /// Returns constant reference to the index of variables inside region \a R (r)
+        const SmallSet<size_t>& INRs( size_t R ) const {
+            DAI_DEBASSERT( R < _INRs.size() );
+            return _INRs[R];
+        }
+
+        /// Returns reference to the index of variables inside region \a R (r)
+        SmallSet<size_t>& INRs( size_t R ) {
+            DAI_DEBASSERT( R < _INRs.size() );
+            return _INRs[R];
+        }
+
+        /// Returns constant reference to the connection structure from region \a R to its \a i'th neighbour
+        const Connection& M( size_t R, size_t i  ) const {
+            DAI_DEBASSERT(R < _M.size());
+            DAI_DEBASSERT(i < _M[R].size());
+            return _M[R][i]; }
+
+        /// Returns reference to the connection structure from region \a R to its \a i'th neighbour
+        Connection& M( size_t R, size_t i  ) {
+            DAI_DEBASSERT(R < _M.size());
+            DAI_DEBASSERT(i < _M[R].size());
+            return _M[R][i]; }
+
+        /// Returns constant reference to the vector of connection structure from region \a R to all its neighbours
+        const std::vector<Connection>& M( size_t R ) const {
+            DAI_DEBASSERT(R < _M.size());
+            return _M[R]; 
+        }
+
+        /// Returns vector of all connections to region \a R
+        std::vector<Connection>& M( size_t R ) { return _M[R]; }
+
+        /// Returns the vector of region indices that contain \a i as internal variable
+        const std::vector<size_t>& var2CW( size_t i ) const { return _var2CW[i]; }
+    //@}
+
+    /// \name Operations
+    //@{
+        /// Sets up all the regions and messages
+        void setRgn( std::vector<SmallSet<size_t> > regions, NeighborType neighbors, bool debugging = false );
+    //@}
+
+    /// \name Input/output
+    //@{
+        /// Reads a cobweb graph from a file
+        /** \note Not implemented yet
+         */
+        virtual void ReadFromFile( const char* /*filename*/ ) {
+            DAI_THROW(NOT_IMPLEMENTED);
+        }
+
+        /// Writes a cobweb graph to a file
+        /** \note Not implemented yet
+         */
+        virtual void WriteToFile( const char* /*filename*/, size_t /*precision*/=15 ) const {
+            DAI_THROW(NOT_IMPLEMENTED);
+        }
+
+        /// Writes a cobweb graph to an output stream
+        friend std::ostream& operator<< ( std::ostream& os, const CobwebGraph& rg );
+
+        /// Writes a cobweb graph to a GraphViz .dot file
+        /** \note Not implemented yet
+         */
+        virtual void printDot( std::ostream& /*os*/ ) const {
+            DAI_THROW(NOT_IMPLEMENTED);
+        }
+    //@}
+
+
+    protected:
+        /// The function to check for partitioning
+        bool checkPartition( const std::vector<SmallSet<size_t> >& regions ) const;
+
+        /// Helper function that sets the regions containing each variable using the values in _INRs
+        void setVar2CW();
+
+        /// Setup the _INRs, _EXRs and all the factor indices (e.g. Rifs)
+        void setExtnFact();
+
+        /// Helper function that setups the msgs (_M, _outM) using the values in _INRs and _EXRs
+        void setMSGs( NeighborType neighbors );
+
+        /// Sets _cn
+        void setCountingNumbers( bool debugging = false );
+
+        /// For the given set of variables for each region, removes the regions that are non-maximal
+        void eraseNonMaximal( std::vector<SmallSet<size_t> >& regions );
+};
+
+
+} // end of namespace dai
+
+
+#endif
index 975084c..e58a830 100644 (file)
@@ -19,6 +19,7 @@
 #include <vector>
 #include <dai/factorgraph.h>
 #include <dai/regiongraph.h>
+#include <dai/cobwebgraph.h>
 #include <dai/properties.h>
 
 
@@ -152,6 +153,11 @@ class InfAlg {
         /** If \a backup == \c true, make a backup of all factors that are changed.
          */
         virtual void makeCavity( size_t i, bool backup = false ) = 0;
+
+        /// Sets all factors indicated by \a facInds to one.
+        /** If \a backup == \c true, make a backup of all factors that are changed.
+         */
+        virtual void makeRegionCavity( std::vector<size_t> facInds, bool backup = false ) = 0;
     //@}
 
     /// \name Backup/restore mechanism for factors
@@ -229,6 +235,11 @@ class DAIAlg : public InfAlg, public GRM {
         /** If \a backup == \c true, make a backup of all factors that are changed.
          */
         void makeCavity( size_t i, bool backup = false ) { GRM::makeCavity( i, backup ); }
+
+        /// Sets all factors indicated by \a facInds to one.
+        /** If \a backup == \c true, make a backup of all factors that are changed.
+         */
+        void makeRegionCavity( std::vector<size_t> facInds, bool backup ){ GRM::makeRegionCavity( facInds, backup ); }
     //@}
 
     /// \name Backup/restore mechanism for factors
@@ -242,6 +253,8 @@ class DAIAlg : public InfAlg, public GRM {
         void restoreFactor( size_t I ) { GRM::restoreFactor( I ); }
         /// Restore the factors involving the variables in \a vs from their backup copies
         void restoreFactors( const VarSet &vs ) { GRM::restoreFactors( vs ); }
+        /// Restore all factors from their backup copies
+        void restoreFactors() { GRM::restoreFactors(); }
     //@}
 };
 
@@ -252,6 +265,9 @@ typedef DAIAlg<FactorGraph> DAIAlgFG;
 /// Base class for inference algorithms that operate on a RegionGraph
 typedef DAIAlg<RegionGraph> DAIAlgRG;
 
+/// Base class for GLC that operates on CobwebGraph
+typedef DAIAlg<CobwebGraph> DAIAlgCG;
+
 
 /// Calculates the marginal probability distribution for \a vs using inference algorithm \a obj.
 /** calcMarginal() works by clamping all variables in \a vs and calculating the partition sum for each clamped state.
index 851e07f..eb4adc8 100644 (file)
@@ -31,8 +31,8 @@
 
 /** \mainpage Reference manual for libDAI - A free/open source C++ library for Discrete Approximate Inference methods
  *  \author Joris Mooij (with contributions of Frederik Eaton)
- *  \version 0.3.1
- *  \date September 17, 2012
+ *  \version git HEAD
+ *  \date September 17, 2012 - or later
  *
  *  <hr size="1">
  *  \section about About libDAI
  *  <em>Journal of Statistical Mechanics: Theory and Experiment</em> 2005(10)-P10011,
  *  http://stacks.iop.org/1742-5468/2005/P10011
  *
+ *  \anchor RYG12 \ref RYG12
+ *  S. Ravanbakhsh, C.-N. Yu, R. Greiner (2012):
+ *  "A Generalized Loop Correction Method for Approximate Inference in Graphical Models",
+ *  <em>29th International Conference on Machine Learning (ICML 2012)</em>,
+ *  http://www.icml.cc/2012/papers/#paper-304
+ *
  *  \anchor StW99 \ref StW99
  *  A. Steger and N. C. Wormald (1999):
  *  "Generating Random Regular Graphs Quickly",
index a88956b..5aff6a5 100644 (file)
@@ -175,6 +175,14 @@ class FactorGraph {
             return I;
         }
 
+        /// Returns the VarSet corresponding to a vector of variable indices
+        VarSet inds2vars( const std::vector<size_t>& inds ) const {
+               VarSet vs;
+               for( std::vector<size_t>::const_iterator it = inds.begin(); it != inds.end(); it++ )
+                       vs.insert( var(*it) );
+               return vs;
+        }
+
         /// Return all variables that occur in a factor involving the \a i 'th variable, itself included
         VarSet Delta( size_t i ) const;
 
@@ -191,6 +199,14 @@ class FactorGraph {
             return Delta( vs ) / vs;
         }
 
+        /// Return all the indices of variables that occur in a factor involving the \a i 'th variable, itself included
+        SmallSet<size_t> Deltai( size_t i ) const;
+
+        /// Return all the indices of variables that occur in a factor involving the \a i 'th variable, itself excluded
+        SmallSet<size_t> deltai( size_t i ) const {
+            return( Deltai( i ) / i );
+        }
+
         /// Returns \c true if the factor graph is connected
         bool isConnected() const { return _G.isConnected(); }
 
@@ -310,6 +326,11 @@ class FactorGraph {
         /** If \a backup == \c true, make a backup of all factors that are changed
          */
         virtual void makeCavity( size_t i, bool backup = false );
+
+        /// Set all factors indexed by \a facInds to 1
+        /** If \a backup == \c true, make a backup of all factors that are changed
+         */
+        virtual void makeRegionCavity( std::vector<size_t> facInds, bool backup );
     //@}
 
     /// \name Input/Output
diff --git a/include/dai/glc.h b/include/dai/glc.h
new file mode 100644 (file)
index 0000000..f22d27c
--- /dev/null
@@ -0,0 +1,366 @@
+/*  This file is part of libDAI - http://www.libdai.org/
+ *
+ *  Copyright (c) 2006-2012, The libDAI authors. All rights reserved.
+ *
+ *  Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
+ */
+
+
+/// \file
+/// \brief Defines classes GLC and Cobweb, which implement the "Generalized Loop Correction method"
+/// \todo Fix the init of GLC
+
+
+#ifndef __defined_libdai_glc_h
+#define __defined_libdai_glc_h
+
+
+#include <algorithm>
+#include <set>
+#include <iostream>
+#include <string>
+#include <dai/util.h>
+#include <dai/daialg.h>
+#include <dai/enum.h>
+#include <dai/cobwebgraph.h>
+#include <dai/properties.h>
+#include <dai/exceptions.h>
+
+
+namespace dai {
+
+
+/// Represents a region for the GLC algorithm.
+/** This region contains all the proper factors from the factor graph. It also contains a factor for incoming messages and the
+ *  cavity distribution over this region. The outgoing messages also receive a factor mainly so that we can easily calculate the marginal
+ *  of the region over the outgoing message variables.
+ *
+ *  \author Siamak Ravanbakhsh
+ */
+class Cobweb {
+    public:
+        /// Pointer to the approximate inference algorithm and factor graph over region R (provided in parameters of GLC)
+        DAIAlgFG *cav;
+
+        /// Global to local factor index conversion.
+        /** Factor(I) in original graph corresponds to cav->factor(_g2l[I]) in this region
+         *  the j'th incoming message to this region M(R,j) corresponds to cav->factor(_g2l[-j-1])
+         *  the j'th cavity factor (for pairwise case),_cavitydists[R][j] corresponds to cav->factor(_g2l[-m-j-1])
+         *      when m = M(R).size()
+         *      note that for full cavity the _cavitydists[R] has a single factor and for uniform cavity it is empty
+         *  the j'th outgoing message from this region outM(R,j) corresponds to cav->factor(_g2l[-c-m-j-1])
+         *      when m  = M(R).size() and c = _cavitydists[R].size()
+         */
+        std::map<int, size_t> _g2l;
+
+        /// Default constructor
+        Cobweb(){}
+
+        /// Construct from global-to-local factor index conversion structure
+        Cobweb( const std::map<int, size_t>& g2l ): cav(0), _g2l(g2l) {}
+
+        /// Default destructor
+        virtual ~Cobweb(){ delete cav; }
+
+        /// Sets the factor graph and inference algorithm for this region
+        void setInfAlg( DAIAlgFG* alg ) {
+            cav = alg;
+        }
+
+        /// Returns the factor of this region globaly indexed by \a I (mainly used to access messages to and from this region)
+        Factor factor( size_t I ) {
+            return cav->factor( _g2l[I] );
+        }
+
+        /// Initializes the inference algorithm over this region
+        void initialize() {
+            cav->init();
+        }
+
+        /// Calculates the marginal over \a ns after temporarily removing the factors indexed by \a rmgind
+        Factor marginal( const VarSet& ns, const std::vector<size_t>& rmgind ) {
+            // Set the factors in rmgind to one
+            VarSet vs;
+            std::vector<size_t> rmlinds; // local index
+            rmlinds.reserve( rmgind.size() );
+            for( std::vector<size_t>::const_iterator it = rmgind.begin(); it != rmgind.end(); it++ ) {
+                vs |= cav->factor( _g2l[*it] ).vars();
+                rmlinds.push_back( _g2l[*it] );
+            }
+            cav->makeRegionCavity( rmlinds, true );
+            // Initialize if necessary
+            cav->init( vs );
+            cav->run();
+            Factor result;
+            result = cav->belief( ns );
+            cav->restoreFactors();
+            return result;
+        }
+
+        /// Calculates the marginal over the variables in outgoing message indexed by \a outmsgind after temporarily removing the factors indexed by \a rmgind
+        /** This function is used to calculate outgoing messages from this region
+         */
+        Factor marginal( const int& outmsgind, const std::vector<size_t>& rmgind ) {
+            // set the factors in rmgind to one
+            VarSet vs;
+            std::vector<size_t> rmlinds; // local index
+            rmlinds.reserve( rmgind.size() );
+            for( std::vector<size_t>::const_iterator it = rmgind.begin(); it != rmgind.end(); it++ ) {
+                vs |= cav->factor( _g2l[*it] ).vars();
+                rmlinds.push_back( _g2l[*it] );
+            }
+            cav->makeRegionCavity( rmlinds, true );
+            // Initialize if necessary
+            cav->init( vs );
+            cav->run();
+            Factor result;
+            result = cav->beliefF( _g2l[outmsgind] );
+            cav->restoreFactors();
+            return result;
+        }
+
+        /// Sets (or multiplies, if \a multiply == \c true) the factor indexed by \a gind by factor \a msg
+        void updateFactor( int gind, const Factor& msg, bool multiply=false ) {
+            if( !multiply )
+                cav->setFactor( _g2l[gind], msg, false );
+            else
+                cav->setFactor( _g2l[gind], (msg*(cav->factor( _g2l[gind] ))).normalized(), false );
+        }
+
+        /// Runs inference and returns belief of variable \a v
+        Factor belief( const Var v ) {
+            cav->run();
+            return cav->belief( v );
+        }
+
+        /// Runs inference and returns belief of variables \a vs
+        Factor belief( const VarSet& vs ) {
+            cav->run();
+            return cav->belief( vs );
+        }
+};
+
+
+/// Approximate inference algorithm "Generalized Loop Correction" [\ref RYG12]
+/** Inherits from a CobwebGraph which in turn inherits from a FactorGraph.
+ *  Although a CobwebGraph represents the graph with cobweb clusters (for design reasons?), 
+ *  the actual vector of cobwebs is contained in GLC class. All the other functionalities 
+ *  needed to represent the cobwebgraph is in its own class. The cobweb provides the function 
+ *  of marginalization using the provided inference method. This can be beneficial when using 
+ *  uniform cavities (or a subset of pairwise cavities, but this is not implemented yet). 
+ *  However, in general the performance does not improve using Junction Tree instead of exact marginalization...
+ *
+ *  \author Siamak Ravanbakhsh
+ */
+class GLC : public DAIAlgCG {
+    private:
+        /// Stores for each variable the approximate cavity distribution
+        std::vector<Cobweb> _CWs;
+
+        /// Gives the starting global index for the factors of outgoing messages of each region
+        std::vector<int> _outOffset; 
+
+        /// Cavity distributions
+        /** For each region it contains:
+         *  - Nothing if using UniCAV
+         *  - One Factor if using FullCAV
+         *  - All pairwise caivity factors if using PairCAV or Pair2CAV
+         */
+        std::vector<std::vector<Factor> > _cavitydists;
+
+        /// Single variable beliefs
+        std::vector<Factor> _beliefs;
+
+        /// Beliefs over factors
+        std::map<VarSet, Factor> _factorBeliefs;
+
+        /// Maximum difference encountered so far
+        Real _maxdiff;
+
+        /// Number of iterations needed
+        size_t _iters;
+
+    public:
+        /// Parameters for GLC
+        struct Properties {
+            /// Enumeration of possible ways to initialize the cavities
+            /** The following initialization methods are defined:
+             *  - FULL calculates the marginal using calcMarginal()
+             *  - PAIR calculates only second order interactions using calcPairBeliefs() with \a accurate == \c false
+             *  - PAIR2 calculates only second order interactions using calcPairBeliefs() with \a accurate == \c true
+             *  - UNIFORM uses a uniform distribution
+             */
+            DAI_ENUM(CavityType,FULL,PAIR,PAIR2,UNIFORM);
+
+            /// Enumeration of different possible types of region (\ominus r in the paper)
+            /** The following initialization methods are defined:
+             *  - SINGLE using single variable
+             *  - FACTOR partitioning regions of that form corresponding to GLC (in the paper)
+             *  - LOOP partitioning regions of that form corresponding to GLC (in the paper)
+             *  - DELTA partitioning regions of that form corresponding to GLC (in the paper)
+             *  - OVFACTOR, OVLOOP, OVDELTA the corresponding form of region in the overlapping form corresponding to GLC+ (in the paper)
+             *  For example, using OVFACTOR each factor becomes a region (\ominus r), while using FACTOR the regions contain only
+             *  non-overlapping factors and if some variables can't be covered by any remaining factors single variable regions
+             *  are constructed for them.
+             *  With this convention OVDELTA here corresponds to DELTA using CVM (HAK).
+             */
+            DAI_ENUM(RegionType,SINGLE,FACTOR,OVFACTOR,LOOP,OVLOOP,DELTA,OVDELTA)
+
+            /// Enumeration of different update schedules
+            /** The following update schedules are defined:
+             *  - SEQFIX sequential fixed schedule
+             *  - SEQRND sequential random schedule
+             */
+            DAI_ENUM(UpdateType,SEQFIX,SEQRND);
+
+            /// Verbosity (amount of output sent to stderr)
+            size_t verbose;
+
+            // Maximum length of the loops to use in construction of loopy regions
+            size_t loopdepth;
+
+            /// Maximum number of iterations
+            size_t maxiter;
+
+            /// The algorithm will quit and report the current result at reaching this runtime
+            Real maxtime;
+
+            /// Tolerance for convergence test
+            Real tol;
+
+            /// Complete or partial reinitialization of the messages in calculating the cavity distribution?
+            bool reinit;
+
+            /// Whether or not to include auxiliary factors in a region. 
+            /** The auxilary factors are the ones that only depend on variables on the boundary of a region (\oplus r - \ominus r).
+             *  This sometimes produces orders of magnitude improvement (e.g., in grids), but it is not discussed in the paper.
+             */
+            bool auxfactors;
+
+            /// How to initialize the cavities
+            CavityType cavity;
+
+            // Type of region
+            RegionType rgntype;
+
+            /// What update schedule to use
+            UpdateType updates;
+
+            /** Indicates what happens if a subset of variables in the boundary of a region (\ominus r_p) is shared by
+             *  some neighbors such that one (\ominus r_{p,q1}) is a subset of another (\ominus r_{p,q2}).
+             *  - ALL     all such neighbors are included in the the updates.
+             *  - TOP     only (\ominus r_{p,q2}) is included unless (\ominus r_{p,q2} = \ominus r_{p,q1}) in which case both are included
+             *  - CLOSEST (default): similar to TOP but in case of a tie the the region r_q with largest r_q \cap r_p is considered
+             *  \note Not important in perfomance!
+             */
+            CobwebGraph::NeighborType neighbors;
+
+            /// Name of the algorithm used to initialize the cavity distributions
+            std::string cavainame;
+
+            /// Parameters for the algorithm used to initialize the cavity distributions
+            PropertySet cavaiopts;
+
+            // Name of the algorithm used for marginalization over a region (EXACT is fine)
+            std::string inainame;
+
+            // Parameters for algorithm used for marginalization over a region
+            PropertySet inaiopts;
+        } props;
+
+
+    public:
+        /// Default constructor
+        GLC() : DAIAlgCG(), _CWs(),_outOffset(), _cavitydists(), _beliefs(), _factorBeliefs(), _maxdiff(), _iters(), props() {}
+
+        /// Construct from FactorGraph \a fg and PropertySet \a opts
+        /** \param fg Factor graph.
+         *  \param opts Parameters @see Properties
+         */
+        GLC( const FactorGraph &fg, const PropertySet &opts );
+
+    /// \name General InfAlg interface
+    //@{
+        virtual GLC* clone() const { return new GLC(*this); }
+        virtual GLC* construct( const FactorGraph &fg, const PropertySet &opts ) const { return new GLC( fg, opts ); }
+        virtual std::string name() const { return "GLC"; }
+        virtual Factor belief( const Var &v ) const { return beliefV( findVar( v ) ); }
+        virtual Factor belief( const VarSet &vs ) const;
+        virtual Factor beliefV( size_t i ) const { return _beliefs[i]; }
+        virtual std::vector<Factor> beliefs() const;
+        virtual Real logZ() const { DAI_THROW(NOT_IMPLEMENTED); return 0.0; }
+        virtual void init(){ initCWs(); }
+        virtual void init( const VarSet &/*ns*/ ) { init(); }
+        virtual Real run();
+        virtual Real maxDiff() const { return _maxdiff; }
+        virtual size_t Iterations() const { return _iters; }
+        virtual void setMaxIter( size_t maxiter ) { props.maxiter = maxiter; }
+        virtual void setProperties( const PropertySet &opts );
+        virtual PropertySet getProperties() const;
+        virtual std::string printProperties() const;
+    //@}
+
+    /// \name Additional interface specific for GLC
+    //@{
+        /// Returns constant reference to outer region \a alpha
+        const Cobweb& CW( size_t alpha ) const {
+            DAI_DEBASSERT( alpha < nrCWs() );
+            return _CWs[alpha];
+        }
+
+        /// Returns reference to outer region \a alpha
+        Cobweb& CW( size_t alpha ) {
+            DAI_DEBASSERT( alpha < nrCWs() );
+            return _CWs[alpha];
+        }
+
+        /// Approximates the cavity distribution of variable \a i, using the inference algorithm \a name with parameters \a opts
+        Real CalcCavityDist( size_t i, const std::string &name, const PropertySet &opts );
+
+        /// Approximates all cavity distributions using inference algorithm \a name with parameters \a opts
+        Real InitCavityDists( const std::string &name, const PropertySet &opts );
+
+        /// Updates the belief of the Markov blanket of region \a R using its intersection with \a _R2'th neighbor (partitioning case).
+        void NewPancake( size_t R, size_t _R2 );
+
+        /// Updates the belief of the Markov blanket of region \a R using all its overlapping neighbors (non-partitioning case).
+        /** \note The implementation uses a more compact form than the message passing over neighborhood region-graph explained in the paper.
+         */
+        void OVNewPancake( size_t R );
+
+        /// Calculates the belief of variable \a i: if \a isFinal == \c true, it returns the geometric average given by all regions that include \a i
+        void CalcBelief( size_t i, bool isFinal = false );
+
+        /// Calculates the belief of factor \a I
+        void CalcFactorBelief( size_t I );
+
+        /// Designates the variables in each region of the graph based on the rgntype
+        std::vector<SmallSet<size_t> > calcRegions();
+
+        /// Initializes the cobweb regions and all the messages between regions.
+        void initCWs();
+
+        /// Initializes all the cobweb regions.
+        /** Defines their factors using the cavity distribution (as factor(s))
+         *  and all the anticipated messages to be exchanged with other regions (treating a message as a factor),
+         *  and also it builds a mapping from cobweb factors to factors of the original factor graph and messages.
+         */
+        void setCWs( const std::string &name, const PropertySet &opts );
+
+        /// Helper function for partitioning case to build the regions by finding loops of maximum given size.
+        /** If a variable is contained in a loop it is not considered for inclusion in any other loop.
+         */
+        void findLoopClusters( SmallSet<size_t>& remaining, std::set<SmallSet<size_t> > &allcl, SmallSet<size_t> newcl, const size_t& root, size_t length, SmallSet<size_t> vars );
+
+        /// Helper function for general case to build the regions by finding loops of maximum given size.
+        /** All loops of max-size are found and variables that are not covered are contained in smaller loops or in factor-sized or single-variable regions.
+         */
+        void findOVLoopClusters( SmallSet<size_t>& remaining, std::set<SmallSet<size_t> > &allcl, SmallSet<size_t> newcl, const size_t& root, size_t length, SmallSet<size_t> vars );
+    //@}
+};
+
+
+} // end of namespace dai
+
+
+#endif
index 525e9a7..aef3ed3 100644 (file)
@@ -234,14 +234,14 @@ class RegionGraph : public FactorGraph {
             DAI_THROW(NOT_IMPLEMENTED);
         }
 
-        /// Writes a factor graph to a file
+        /// Writes a region graph to a file
         /** \note Not implemented yet
          */
         virtual void WriteToFile( const char* /*filename*/, size_t /*precision*/=15 ) const {
             DAI_THROW(NOT_IMPLEMENTED);
         }
 
-        /// Writes a RegionGraph to an output stream
+        /// Writes a region graph to an output stream
         friend std::ostream& operator<< ( std::ostream& os, const RegionGraph& rg );
 
         /// Writes a region graph to a GraphViz .dot file
index ae0b720..8109cdd 100644 (file)
@@ -59,6 +59,10 @@ class _builtinInfAlgs : public std::map<std::string, InfAlg *> {
 #ifdef DAI_WITH_DECMAP
             operator[]( DecMAP().name() ) = new DecMAP;
 #endif
+#ifdef DAI_WITH_GLC
+            operator[]( GLC().name() ) = new GLC;
+#endif
+
         }
 
         ~_builtinInfAlgs() {
diff --git a/src/cobwebgraph.cpp b/src/cobwebgraph.cpp
new file mode 100644 (file)
index 0000000..73c368d
--- /dev/null
@@ -0,0 +1,272 @@
+/*  This file is part of libDAI - http://www.libdai.org/
+ *
+ *  Copyright (c) 2006-2012, The libDAI authors. All rights reserved.
+ *
+ *  Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
+ */
+
+
+#include <dai/cobwebgraph.h>
+#include <algorithm>
+#include <map>
+#include <set>
+
+
+namespace dai {
+
+
+using namespace std;
+
+
+void CobwebGraph::setRgn( vector<SmallSet<size_t> > partition, NeighborType neighbors, bool debugging ) {
+    if( checkPartition( partition ) )
+        isPartition = true;
+    else
+        isPartition = false;
+
+    eraseNonMaximal( partition );
+    if( debugging )
+        cerr << "isPartition:" << isPartition << endl;
+    _INRs = partition;
+    if( debugging )
+        cerr << "setting external factors and regional factors" << endl;
+    setExtnFact();
+    if( debugging )
+        cerr << "setting msgs" << endl;
+    setMSGs( neighbors );
+//    if(!isPartition)
+    if( debugging )
+        cerr << "setting the counting numbers" << endl;
+    if( !isPartition )
+        setCountingNumbers( debugging );
+    if( debugging )
+        cerr << "setting var2dr" << endl;
+    setVar2CW();
+}
+
+
+void CobwebGraph::eraseNonMaximal( vector<SmallSet<size_t> >& partition ) {
+    vector<SmallSet<size_t> >::iterator del = partition.begin();
+    while( del != partition.end() ){
+        bool isNonMax = false;
+        for( vector<SmallSet<size_t> >::iterator it = partition.begin(); it != partition.end(); it++ )
+            if( (*it) >> (*del) && it != del ) {
+                isNonMax = true;
+                break;
+            }
+        if( isNonMax )
+            del = partition.erase( del );
+        else
+            del++;
+    }
+}
+
+
+void CobwebGraph::setCountingNumbers( bool debugging ) {
+    // should handle the case when one region in the msg-region graph is subset of another or even to top regions are the same
+    _cn.reserve( nrCWs() );
+    for( size_t R = 0; R < _INRs.size(); R++ ) {
+        vector<VarSet> topR;
+        // finding the intersection of messages
+        SmallSet<VarSet> betas;
+        for( size_t m1 = 0; m1 < M(R).size(); m1++ ) {
+            topR.push_back( M(R,m1).msg.vars() );
+            for( size_t m2 = 0; m2 < M(R).size(); m2++ ) {
+                if( m1 == m2 )
+                    continue;
+                VarSet b = M(R,m1).msg.vars() & M(R,m2).msg.vars();
+                if( b.size() == 0 )
+                    continue;
+                betas.insert( b );
+            }
+        }
+        if( debugging )
+            cerr << "topR: " << topR << endl;
+
+        // finding the intersections of intersections
+        bool somechange = true;
+        while( somechange ) {
+            SmallSet<VarSet> newbetas;
+            somechange = false;
+            bforeach( const VarSet &b1, betas ) {
+                bforeach( const VarSet &b2, betas ) {
+                    if( b1 == b2 )
+                        continue;
+                    VarSet b3 = b1 & b2;
+                    if( betas.contains(b3) || b3.size() == 0 )
+                        continue;
+                    newbetas |= b3;
+                    somechange = true;
+                }
+            }
+            betas |= newbetas;
+        }
+        if( debugging )
+            cerr << "betas: " << betas << endl;
+
+        // set the counting number
+        _cn.push_back( map<VarSet, pair<int, vector<size_t> > >() );
+        // adding sub-regions of every message
+        for( size_t i = 0; i < topR.size(); i++ ) {
+            M(R,i).subregions.clear();
+            bforeach( const VarSet& b, betas )
+                if( b << topR[i] )
+                    M(R,i).subregions.push_back( b );
+        }
+        SmallSet<VarSet> subVisited;
+        SmallSet<VarSet> topRSet( topR.begin(), topR.end(), topR.size() );
+        // looks to see if all parents of a sub-region got their counting number
+        while( !betas.empty() ) {
+            bforeach( const VarSet &beta, betas ) {
+                bool allparentsset = true;
+                bforeach( const VarSet &beta2, betas )
+                    if( beta2 >> beta && beta2 != beta ) {
+                        allparentsset = false;
+                        break;
+                    }
+                if( allparentsset ) {
+                    // the first in the pair is cn and the second the index of top regions containing it
+                    _cn[R][beta] = make_pair( 1, vector<size_t>() );
+                    for( size_t TR = 0; TR < topR.size(); TR++ ) {
+                        if( topR[TR] >> beta ) {
+                            _cn[R][beta].first--;
+                            _cn[R][beta].second.push_back( TR );
+                        }
+                    }
+                    bforeach( const VarSet& possibleparent, subVisited )
+                        if( possibleparent >> beta )
+                            _cn[R][beta].first -= _cn[R][possibleparent].first;
+
+                    if( debugging )
+                        cerr << "cn[" << R << "][" << beta << "] <- " << _cn[R][beta] << endl;
+                    subVisited.insert( beta );
+                    betas /= beta;
+                    break; // since betas has changed we need to enter the loop again
+                }
+            }
+        }
+        if( debugging )
+            cerr << "cn[" << R << "] " << _cn[R] << endl;
+    }
+}
+
+
+bool CobwebGraph::checkPartition( const vector<SmallSet<size_t> >& regions ) const {
+    for( size_t i = 0; i < regions.size(); i++ )
+        for( size_t j = 0; j < regions.size(); j++ ) {
+            if( j == i )
+                continue;
+            if( regions[i].intersects( regions[j] ) )
+                return false;
+        }
+    return true;
+}
+
+
+void CobwebGraph::setVar2CW() {
+    _var2CW.clear();
+    _var2CW.resize( nrVars() );
+    for( size_t R = 0; R < _INRs.size(); R++ ) {
+        bforeach( size_t i, _INRs[R] )
+            _var2CW[i].push_back( R );
+    }
+}
+
+
+// assumes availablity of externals and internals
+void CobwebGraph::setMSGs( NeighborType neighbors ) {
+    // initialize
+    _outM.clear();
+    _outM.reserve( nrCWs() );
+    for( size_t R = 0; R < nrCWs(); R++ )
+        _outM.push_back( vector<VarSet>() );
+    _M.clear();
+    _M.reserve( _INRs.size() );
+    for( size_t R = 0; R < _INRs.size(); R++ ) {
+        _M.push_back( vector<Connection>() );
+        for( size_t R2 = 0; R2 < _INRs.size(); R2++ ) {
+            if( R2 == R )
+                continue;
+            SmallSet<size_t> tmpSet = _EXRs[R] & (_INRs[R2]); // if R2 has boundary of R1 inside it
+            if( tmpSet.size() == 0 )
+                continue;
+
+            Connection tmp;
+            tmp.fc = (_Rfs[R] & _Rfs[R2]).elements();
+            tmp.my = R;
+            tmp.iter = _M[R].size();
+            tmp.his = R2;
+            tmp.dual = _outM[R2].size();
+            VarSet tmpVarSet;
+            bforeach( const size_t &i, tmpSet )
+                tmpVarSet.insert( var(i) );
+            tmp.msg = Factor( tmpVarSet, (Real)1.0 ); //*
+            _outM[R2].push_back( tmpVarSet ); //*
+            // tmp.varinds = tmpSet.elements();
+            _M[R].push_back( tmp );
+        }
+        if( neighbors == NeighborType::CLOSEST || neighbors == NeighborType::TOP ) {
+            bool foundone = true;
+            while( foundone ) {
+                for( size_t Rfirst = 0; Rfirst < _M[R].size(); Rfirst++ ) {
+                    foundone = false;
+                    for( size_t Rsec = 0; Rsec < _M[R].size(); Rsec++ ) {
+                        if( Rfirst == Rsec )
+                            continue;
+                        if( _M[R][Rfirst].msg.vars() >> _M[R][Rsec].msg.vars() ) { // one message is subseteq of the other
+                            // only if one exclusively contains the other remove it
+                            if( ((neighbors == NeighborType::TOP || neighbors == NeighborType::CLOSEST) && _M[R][Rfirst].msg.vars().size() > _M[R][Rsec].msg.vars().size()) ) {
+                                _M[R].erase( _M[R].begin() + Rsec );
+                                foundone = true;
+                                break;
+                            } else if( neighbors == NeighborType::CLOSEST && // if both have the same size then in case of CLOSEST delete the second one if it has smaller inner region overlap
+                                    (INRs(_M[R][Rfirst].his) & INRs(R)).size() > (INRs(_M[R][Rsec].his) & INRs(R)).size() ) {
+                                _M[R].erase( _M[R].begin() + Rsec );
+                                foundone = true;
+                                break;
+                            }
+                        }
+                    }
+                    if( foundone )
+                        break;
+                }
+            }
+        }
+    }
+}
+
+
+void CobwebGraph::setExtnFact() {
+    // finds and sets all the externals using _INRs
+    SmallSet<size_t> allpots;
+    _EXRs.clear();
+    _Rfs.clear();
+    _EXRs.reserve( _INRs.size() );
+    _Rfs.reserve( _INRs.size() );
+    _Rifs.reserve( _INRs.size() );
+    _Rxfs.reserve( _INRs.size() );
+    for( size_t R = 0; R < _INRs.size(); R++ ) {
+        _EXRs.push_back( SmallSet<size_t>() );
+        _Rfs.push_back( SmallSet<size_t>() );
+        _Rifs.push_back( SmallSet<size_t>() );
+        _Rxfs.push_back( SmallSet<size_t>() );
+        bforeach( const size_t &i, _INRs[R] ) {
+            _EXRs[R] |= deltai(i) / _INRs[R];
+            bforeach( const Neighbor &I, nbV(i) ) {
+                _Rfs[R] |= I.node;
+                if( factor(I).vars() << inds2vars(INRs(R).elements()) )
+                    _Rifs[R] |= I.node;
+                else
+                    _Rxfs[R] |= I.node;
+            }
+        }
+    }
+}
+
+
+ostream& operator<< ( ostream& os, const CobwebGraph& rg ) {
+    return os;
+}
+
+
+} // end of namespace dai
index 045a46a..3cff7eb 100644 (file)
@@ -219,6 +219,17 @@ VarSet FactorGraph::Delta( const VarSet &ns ) const {
 }
 
 
+SmallSet<size_t> FactorGraph::Deltai( size_t i ) const {
+    // calculate Markov Blanket
+    SmallSet<size_t> Del;
+    bforeach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
+        bforeach( const Neighbor &j, nbF(I) ) // for all neighboring variables j of I
+            Del |= j;
+
+    return Del;
+}
+
+
 void FactorGraph::makeCavity( size_t i, bool backup ) {
     // fills all Factors that include var(i) with ones
     map<size_t,Factor> newFacs;
@@ -228,6 +239,14 @@ void FactorGraph::makeCavity( size_t i, bool backup ) {
 }
 
 
+void FactorGraph::makeRegionCavity( std::vector<size_t> facInds, bool backup ) {
+    map<size_t,Factor> newFacs;
+    for( size_t I = 0; I < facInds.size(); I++ )
+       newFacs[facInds[I]] = Factor(factor(facInds[I]).vars(), (Real)1);
+    setFactors( newFacs, backup );
+}
+
+
 void FactorGraph::ReadFromFile( const char *filename ) {
     ifstream infile;
     infile.open( filename );
diff --git a/src/glc.cpp b/src/glc.cpp
new file mode 100644 (file)
index 0000000..977a563
--- /dev/null
@@ -0,0 +1,591 @@
+/*  This file is part of libDAI - http://www.libdai.org/
+ *
+ *  Copyright (c) 2006-2012, The libDAI authors. All rights reserved.
+ *
+ *  Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
+ */
+
+
+#include <iostream>
+#include <algorithm>
+#include <map>
+#include <set>
+#include <dai/glc.h>
+#include <dai/util.h>
+#include <dai/alldai.h>
+
+
+namespace dai {
+
+
+using namespace std;
+
+
+void GLC::setProperties( const PropertySet &opts ) {
+    DAI_ASSERT( opts.hasKey("tol") );
+    DAI_ASSERT( opts.hasKey("maxiter") );
+    DAI_ASSERT( opts.hasKey("cavity") );
+    DAI_ASSERT( opts.hasKey("updates") );
+    DAI_ASSERT( opts.hasKey("rgntype") );
+    DAI_ASSERT( opts.hasKey("cavainame") );
+    DAI_ASSERT( opts.hasKey("cavaiopts") );
+    DAI_ASSERT( opts.hasKey("inainame") );
+    DAI_ASSERT( opts.hasKey("inaiopts") );
+
+    if( opts.hasKey("maxtime") )
+        props.maxtime = opts.getStringAs<Real>("maxtime");
+    else
+        props.maxtime = INFINITY;
+
+    props.tol = opts.getStringAs<Real>("tol");
+    props.maxiter = opts.getStringAs<size_t>("maxiter");
+    props.cavity = opts.getStringAs<Properties::CavityType>("cavity");
+    props.updates = opts.getStringAs<Properties::UpdateType>("updates");
+    props.rgntype = opts.getStringAs<Properties::RegionType>("rgntype");
+    if(opts.hasKey("neighbors"))
+        props.neighbors = opts.getStringAs<CobwebGraph::NeighborType>("neighbors");
+    else
+        props.neighbors = CobwebGraph::NeighborType::CLOSEST;
+    if( opts.hasKey("verbose") )
+        props.verbose = opts.getStringAs<size_t>("verbose");
+    else
+        props.verbose = 0;
+    if( opts.hasKey("loopdepth") )
+        props.loopdepth = opts.getStringAs<size_t>("loopdepth");
+    else
+        props.loopdepth = 4;
+    if( opts.hasKey("cavainame") )
+        props.cavainame = opts.getStringAs<string>("cavainame");
+    if( opts.hasKey("cavaiopts") )
+        props.cavaiopts = opts.getStringAs<PropertySet>("cavaiopts");
+    if( opts.hasKey("inainame") )
+        props.inainame = opts.getStringAs<string>("inainame");
+    if( opts.hasKey("inaiopts") )
+        props.inaiopts = opts.getStringAs<PropertySet>("inaiopts");
+    if( opts.hasKey("reinit") )
+        props.reinit = opts.getStringAs<bool>("reinit");
+}
+
+
+PropertySet GLC::getProperties() const {
+    PropertySet opts;
+    opts.set( "tol", props.tol );
+    opts.set( "maxiter", props.maxiter );
+    opts.set( "verbose", props.verbose );
+    opts.set( "loopdepth", props.loopdepth );
+    opts.set( "cavity", props.cavity );
+    opts.set( "neighbors", props.neighbors );
+    opts.set( "updates", props.updates );
+    opts.set( "rgntype", props.rgntype );
+    opts.set( "cavainame", props.cavainame );
+    opts.set( "cavaiopts", props.cavaiopts );
+    opts.set( "inainame", props.inainame );
+    opts.set( "inaiopts", props.inaiopts );
+    opts.set( "reinit", props.reinit );
+    opts.set( "maxtime", props.maxtime );
+
+    return opts;
+}
+
+
+string GLC::printProperties() const {
+    stringstream s( stringstream::out );
+    s << "[";
+    s << "tol=" << props.tol << ",";
+    s << "maxiter=" << props.maxiter << ",";
+    s << "verbose=" << props.verbose << ",";
+    s << "loopdepth=" << props.loopdepth << ",";
+    s << "cavity=" << props.cavity << ",";
+    s << "updates=" << props.updates << ",";
+    s << "neighbors=" << props.neighbors << ",";
+    s << "rgntype=" << props.rgntype << ",";
+    s << "cavainame=" << props.cavainame << ",";
+    s << "cavaiopts=" << props.cavaiopts << ",";
+    s << "inainame=" << props.inainame << ",";
+    s << "inaiopts=" << props.inaiopts << ",";
+    s << "reinit=" << props.reinit << ",";
+    s << "maxtime=" << props.maxtime << "]";
+    return s.str();
+}
+
+
+GLC::GLC( const FactorGraph& fg, const PropertySet &opts ) : DAIAlgCG(fg), _CWs(),_outOffset(), _cavitydists(), _beliefs(), _factorBeliefs(), _maxdiff(0.0), _iters(0), props() {
+    setProperties( opts );
+    setRgn( calcRegions(), props.neighbors, false );
+    if( props.verbose >= 3 ){
+        cerr << "#regions: " << nrCWs() << endl;
+        for( size_t i = 0; i < nrCWs(); i++ )
+            cerr << INRs(i).elements() << endl;
+    }
+    if( props.verbose >= 3 )
+        cerr << "initializing cavity dists" << endl;
+    _cavitydists.resize( nrCWs() );
+    _maxdiff = InitCavityDists( props.cavainame, props.cavaiopts );
+    // build the vector of Cobweb Regions
+    setCWs( props.inainame, props.inaiopts );
+    if( props.verbose >= 3 )
+        cerr << "Regions are built" << endl;
+    // create beliefs
+    _beliefs.reserve( nrVars() );
+    for( size_t i = 0; i < nrVars(); i++ )
+        _beliefs.push_back( Factor(var(i)) );
+}
+
+
+vector<SmallSet<size_t> > GLC::calcRegions() {
+    // contains the variable indices for each region
+    vector<SmallSet<size_t> > regions;
+    // single variables
+    if( props.rgntype == Properties::RegionType::SINGLE ) {
+        regions.resize( nrVars() );
+        for( size_t i = 0; i < nrVars(); i++ )
+            regions[i] = SmallSet<size_t>(i);
+    // partitioning factors
+    } else if( props.rgntype == Properties::RegionType::FACTOR ) {
+        SmallSet<size_t> remainingVars;
+        for( size_t i = 0; i < nrVars(); i++ )
+            remainingVars.insert( i );
+        vector<SmallSet<size_t> > facVars( nrFactors() );
+        for( size_t I = 0; I < nrFactors(); I++ )
+            facVars[I] = findVars( factor(I).vars() );
+        bool canadd = true;
+        while( canadd ) {
+            canadd = false;
+            for( size_t I = 0; I < nrFactors(); I++ ) {
+                if( facVars[I] << remainingVars ) {
+                    regions.push_back( facVars[I] );
+                    remainingVars /= facVars[I];
+                    canadd = true;
+                    break;
+                }
+            }
+        }
+        // add the variables that are not covered as single variable regions
+        bforeach(size_t i, remainingVars)
+            regions.push_back( SmallSet<size_t>(i) );
+    // overlapping factors (non-partitioning)
+    } else if( props.rgntype == Properties::RegionType::OVFACTOR ) {
+        SmallSet<size_t> remainingVars;
+        for( size_t I = 0; I < nrFactors(); I++ )
+            regions.push_back( findVars( factor(I).vars() ) );
+    // partitioning: adds loops over variables that are not covered recursively
+    } else if( props.rgntype == Properties::RegionType::LOOP ) {
+        set<SmallSet<size_t> > scl; // to contain clusters
+        SmallSet<size_t> remaining; // variables not in a cluster yet
+        // build the set of all var inds
+        for( size_t v = 0; v < nrVars(); v++ )
+            remaining.insert( v );
+        // while there for a remaining var we can find a loop
+        bool changing = true;
+        while( changing ) {
+            changing = false;
+            for( size_t i0 = 0; i0 < nrVars(); i0++ ) {
+                if( !remaining.contains(i0) )
+                    continue;
+                size_t currentsize = scl.size();
+                if( props.loopdepth > 1 )
+                    findLoopClusters( remaining, scl, SmallSet<size_t>(i0),i0, props.loopdepth - 1, deltai(i0) & remaining );
+                if(scl.size() > currentsize)
+                    changing = true;
+            }
+        }
+        copy( scl.begin(), scl.end(), back_inserter(regions) );
+        // add the vars that are not covered as single variable regions
+        bforeach( size_t i, remaining )
+            regions.push_back( SmallSet<size_t>(i) );
+    // adds overlapping loops of maximum size loopdepth recursively
+    } else if( props.rgntype == Properties::RegionType::OVLOOP ) {
+        for( size_t I = 0; I < nrFactors(); I++ )
+            regions.push_back( findVars( factor(I).vars() ) );
+
+        SmallSet<size_t> remaining; // variables not in a cluster yet
+        // build the set of all var inds
+        for( size_t v = 0; v < nrVars(); v++ )
+            remaining.insert( v );
+
+        set<SmallSet<size_t> > scl; // to contain clusters
+        // while there for a remaining var we can find a loop
+
+        for( size_t i0 = 0; i0 < nrVars(); i0++ )
+            findOVLoopClusters( remaining, scl, SmallSet<size_t>(i0),i0, props.loopdepth - 1, deltai(i0) );
+        copy( scl.begin(), scl.end(), back_inserter(regions) );
+    // non-partitioning: overlapping (variable+markov blanket)s
+    } else if( props.rgntype == Properties::RegionType::OVDELTA ) {
+        for( size_t I = 0; I < nrFactors(); I++ )
+            regions.push_back( findVars(factor(I).vars()) );
+        SmallSet<size_t> remaining; // variables not in a cluster yet
+        // build the set of all var inds
+        for( size_t v = 0; v < nrVars(); v++ )
+            remaining.insert( v );
+        set<SmallSet<size_t> > scl; // to contain clusters
+        // while there for a remaining var we can find a loop
+        for( size_t i0 = 0; i0 < nrVars(); i0++ )
+            findOVLoopClusters( remaining, scl, SmallSet<size_t>(i0),i0, props.loopdepth - 1, deltai(i0) );
+        copy( scl.begin(), scl.end(), back_inserter(regions) );
+    } else if( props.rgntype == Properties::RegionType::DELTA ) {
+        SmallSet<size_t> remaining; // variables not in a cluster yet
+        // build the set of all var inds
+        for( size_t v = 0; v < nrVars(); v++ )
+            remaining.insert( v );
+        while( remaining.size() > 0 ) {
+            size_t cv = remaining.elements()[ rand() % remaining.size() ];
+            regions.push_back( Deltai(cv) & remaining );
+            remaining /= Deltai(cv);
+        }
+    }
+    return regions;
+}
+
+
+void GLC::findOVLoopClusters( SmallSet<size_t>& remaining, set<SmallSet<size_t> > &allcl, SmallSet<size_t> newcl, const size_t& root, size_t length, SmallSet<size_t> vars ) {
+    for( SmallSet<size_t>::const_iterator in = vars.begin(); in != vars.end(); in++ ) {
+        SmallSet<size_t> ind = deltai( *in );
+        if( (newcl.size()) >= 2 && ind.contains( root ) ) {
+            allcl.insert( newcl | *in );
+            remaining /= (newcl | *in);
+        } else if( length > 1 )
+            findOVLoopClusters( remaining, allcl, newcl | *in, root, length - 1, ind / newcl );
+    }
+}
+
+
+void GLC::findLoopClusters( SmallSet<size_t>& remaining, set<SmallSet<size_t> > &allcl, SmallSet<size_t> newcl, const size_t& root, size_t length, SmallSet<size_t> vars ) {
+    for( SmallSet<size_t>::const_iterator in = vars.begin(); in != vars.end(); in++ ) {
+        SmallSet<size_t> ind = deltai( *in ) & remaining;
+        if( (newcl.size()) >= 2 && ind.contains( root ) ) {
+            allcl.insert( newcl | *in );
+            remaining /= (newcl | *in);
+            break;
+        } else if( length > 1 )
+            findLoopClusters( remaining, allcl, newcl | *in, root, length - 1, ind / newcl );
+    }
+}
+
+
+void GLC::CalcBelief( size_t i, bool isFinal ) {
+    if( !isFinal )
+        _beliefs[i] = CW( var2CW(i)[0] ).belief( var(i) );
+    else{ // take geometric average
+        _beliefs[i] = Factor(var(i), (Real)1);
+        for( size_t rind = 0; rind < var2CW(i).size(); rind++ )
+            _beliefs[i] *= CW( var2CW(i)[rind] ).belief( var(i) );
+        _beliefs[i] ^= (Real)1 / var2CW(i).size();
+    }
+}
+
+
+void GLC::CalcFactorBelief( size_t I ) {
+    VarSet ns = factor(I).vars();
+    for( size_t R = 0; R < nrCWs(); R++ )
+        if( ns << inds2vars( INRs(R).elements() ) )
+            _factorBeliefs[ns] = CW(R).belief(ns);
+}
+
+
+Factor GLC::belief( const VarSet &ns ) const {
+    if( ns.size() == 0 )
+        return Factor();
+    else if( ns.size() == 1 )
+        return beliefV( findVar( *(ns.begin()) ) );
+    else {
+        map<VarSet, Factor>::const_iterator it = _factorBeliefs.find(ns);
+        if( it != _factorBeliefs.end() )
+            return it->second;
+    }
+    DAI_THROW(BELIEF_NOT_AVAILABLE);
+    return Factor();
+}
+
+
+vector<Factor> GLC::beliefs() const {
+    vector<Factor> result = _beliefs;
+    for( map<VarSet, Factor>::const_iterator it = _factorBeliefs.begin(); it != _factorBeliefs.end(); it++ )
+        result.push_back( it->second );
+    return result;
+}
+
+
+Real GLC::CalcCavityDist( size_t R, const std::string &name, const PropertySet& opts ) {
+    vector<Factor> pairbeliefs;
+    Real maxdiff = 0;
+
+    if( props.cavity != Properties::CavityType::UNIFORM ) {
+        InfAlg *cav = newInfAlg( name, *this, opts );
+        cav->makeRegionCavity( Rfs(R).elements() );
+
+        if( props.cavity == Properties::CavityType::FULL )
+            pairbeliefs.push_back( calcMarginal( *cav, inds2vars(EXRs(R).elements()), props.reinit ) );
+        if( props.cavity == Properties::CavityType::PAIR )
+            pairbeliefs = calcPairBeliefs( *cav, inds2vars(EXRs(R).elements()), props.reinit, false );
+        else if( props.cavity == Properties::CavityType::PAIR2 )
+            pairbeliefs = calcPairBeliefs( *cav, inds2vars(EXRs(R).elements()), props.reinit, true );
+        maxdiff = cav->maxDiff();
+        delete cav;
+    }
+    if( props.verbose >= 3 )
+        cerr << "R:" << R << "cavity of size " << pairbeliefs.size() << endl;
+    _cavitydists[R] = pairbeliefs;
+    return maxdiff;
+}
+
+
+Real GLC::InitCavityDists( const std::string &name, const PropertySet &opts ) {
+    double tic = toc();
+
+    if( props.verbose >= 2 ) {
+        cerr << this->name() << "::InitCavityDists:  ";
+        if( props.cavity == Properties::CavityType::UNIFORM )
+            cerr << "Using uniform initial cavity distributions" << endl;
+        else if( props.cavity == Properties::CavityType::FULL )
+            cerr << "Using full " << name << opts << "...";
+        else if( props.cavity == Properties::CavityType::PAIR )
+            cerr << "Using pairwise " << name << opts << "...";
+        else if( props.cavity == Properties::CavityType::PAIR2 )
+            cerr << "Using pairwise(new) " << name << opts << "...";
+    }
+
+    Real maxdiff = 0.0;
+    for( size_t R = 0; R < nrCWs(); R++ ) {
+        Real md = CalcCavityDist(R, name, opts);
+        if( md > maxdiff )
+            maxdiff = md;
+    }
+    if( props.verbose >= 2 )
+        cerr << this->name() << "::InitCavityDists used " << toc() - tic << " seconds." << endl;
+    return maxdiff;
+}
+
+
+void GLC::NewPancake( size_t R, size_t _R2 ) {
+    Factor newmsg;
+    // calculate the new msg
+    newmsg = CW(M(R, _R2).his).marginal( M(R, _R2).msg.vars(), M(R, _R2).fc) / ((CW(R).marginal( M(R, _R2).msg.vars(), M(R, _R2).fc))*M(R, _R2).msg.inverse());
+    newmsg.normalize();
+    // update the cobweb-graph with this new msg
+    M(R, _R2).msg = newmsg;
+    // update the corresponding factor in the cobweb region so to be used in future marginalizations (-_R2 - 1) is the index (see cobweb class documents)
+    CW(R).updateFactor( -_R2 - 1 ,newmsg);
+}
+
+
+void GLC::OVNewPancake( size_t R ) {
+    Factor newmsg, allmsgs, marg;
+    for( size_t _R2 = 0; _R2 < M(R).size(); _R2++ ) { // for all _R2 that send messages to R
+        // calculate the message
+        newmsg = (CW(M(R, _R2).his).marginal( M(R, _R2).msg.vars(), M(R, _R2).fc) / ((CW(R).marginal( M(R, _R2).msg.vars(), M(R, _R2).fc))))*M(R, _R2).msg;
+        newmsg.normalize();
+        // allmsgs takes care of the double-counts the way using the region-graph does in the paper. It is more compact but not as efficient
+        allmsgs = newmsg;
+        for( size_t sub = 0; sub < M(R,_R2).subregions.size(); sub++ ) { //for all descendants of rho of (\ominus r_{R, R2})
+            Real c = (Real)cn(R)[M(R, _R2).subregions[sub]].first; // counting number for rho
+            size_t nrParents = cn(R)[M(R, _R2).subregions[sub]].second.size(); // number of rho's parents
+            marg = newmsg.marginal(M(R, _R2).subregions[sub], false); // marginal of the message from R2 to R over rho
+            allmsgs *= (marg^(c / nrParents)); // discount the double counting of rho from the message
+            for( size_t pind = 0; pind < nrParents; pind++ ) // this operation corresponds to the upward pass
+                if( cn(R)[M(R, _R2).subregions[sub]].second[pind] != _R2 )
+                    M(R, cn(R)[M(R, _R2).subregions[sub]].second[pind]).newmsg *= marg^(-c / (nrParents * (nrParents - (Real)1)));
+        }
+        // set the message after taking double-countings into account
+        M(R, _R2).msg = allmsgs.normalized();
+    }
+
+    for( size_t _R2 = 0; _R2 < M(R).size(); _R2++ ) {
+        M(R, _R2).newmsg *= M(R, _R2).msg;
+        M(R, _R2).newmsg.normalize();
+        // update the corresponding factor in the cobweb region: to be used in future marginalizations
+        CW(R).updateFactor( -_R2 - 1 , M(R, _R2).msg, false);
+        M(R, _R2).msg = M(R, _R2).newmsg;
+        M(R, _R2).newmsg.fill( (Real)1 );
+    }
+}
+
+
+Real GLC::run() {
+    if( props.verbose >= 2 )
+        cerr << "Starting " << identify() << "...";
+    if( props.verbose >= 2 )
+        cerr << endl;
+
+    double tic = toc();
+
+    if( props.verbose >= 3 )
+        cerr << "Initializing older beliefs" << endl;
+    vector<vector<Factor> > oldBeliefsM;
+    vector<Factor> oldBeliefsP;
+    if( isPartition ) {
+        oldBeliefsM.reserve( nrCWs() );
+        for( size_t R = 0; R < nrCWs(); R++ ) {
+            oldBeliefsM.push_back( vector<Factor>() );
+            oldBeliefsM[R].reserve( M(R).size() );
+            for( size_t m = 0; m < M(R).size(); m++ )
+                oldBeliefsM[R].push_back( M(R,m).msg );
+        }
+    } else {
+        oldBeliefsP.reserve( nrVars() );
+        for( size_t i = 0; i < nrVars(); i++ ) {
+            CalcBelief( i, true );
+            oldBeliefsP.push_back( _beliefs[i] );
+        }
+    }
+    if( props.verbose >= 3 )
+        cerr << "Setting the update sequence" <<endl;
+
+    size_t nredges = 0;
+    vector<Edge> update_seq;
+    update_seq.reserve( nredges );
+    for( size_t R = 0; R < nrCWs(); ++R )
+        for( size_t R2 = 0; R2 < M(R).size(); R2++ ) {
+            update_seq.push_back( Edge( R, R2 ) );
+            nredges++;
+        }
+
+    vector<size_t> update_seq_node;
+    update_seq_node.reserve( nrCWs() );
+    for( size_t R = 0; R < nrCWs(); ++R )
+        update_seq_node.push_back( R );
+    if( props.verbose >= 3 )
+        cerr << "Going into the main loop" << endl;
+
+    Real maxDiff = INFINITY;
+    for( _iters = 0; _iters < props.maxiter && maxDiff > props.tol && (toc() - tic) < props.maxtime; _iters++ ) {
+        if( props.updates == Properties::UpdateType::SEQRND ) {
+            random_shuffle( update_seq.begin(), update_seq.end() );
+            random_shuffle( update_seq_node.begin(), update_seq_node.end() );
+        }
+        if( isPartition ) {
+            for( size_t t = 0; t < nredges; t++ ) {
+                size_t R = update_seq[t].first;
+                size_t _R2 = update_seq[t].second;
+                NewPancake( R, _R2);
+            }
+        } else {
+            for( size_t t = 0; t < nrCWs(); t++ )
+                OVNewPancake( update_seq_node[t] );
+        }
+
+        if( !isPartition )
+            for( size_t i = 0; i < nrVars(); i++ )
+                CalcBelief( i , true );
+
+        maxDiff = -INFINITY;
+        if( isPartition ) {
+            for( size_t R = 0; R < nrCWs(); R++ ) {
+                for( size_t m = 0; m < M(R).size(); m++ ) {
+                    maxDiff = std::max( maxDiff, dist( M(R, m).msg, oldBeliefsM[R][m], DISTLINF ) );
+                    oldBeliefsM[R][m] = M(R,m).msg;
+                }
+            }
+        } else {
+            for( size_t i = 0; i < nrVars(); i++ ) {
+                maxDiff = std::max( maxDiff, dist( _beliefs[i], oldBeliefsP[i], DISTLINF ) );
+                oldBeliefsP[i] = _beliefs[i];
+            }
+        }
+
+        if( props.verbose >= 2 )
+            cerr << this->name() << "::run:  maxdiff " << maxDiff << " after " << _iters+1 << " passes" << endl;
+    }
+
+    if( isPartition )
+        for( size_t i = 0; i < nrVars(); i++ )
+            CalcBelief( i, true );
+
+    if( maxDiff > _maxdiff )
+        _maxdiff = maxDiff;
+
+    if( props.verbose >= 2 ) {
+        if( maxDiff > props.tol ) {
+            if( props.verbose == 2 )
+                cerr << endl;
+                cerr << this->name() << "::run:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << maxDiff << endl;
+        } else {
+            if( props.verbose >= 2 )
+                cerr << this->name() << "::run:  ";
+                cerr << "converged in " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
+        }
+    }
+
+    for( size_t I = 0; I < nrFactors(); I++ )
+        CalcFactorBelief( I );
+
+    return maxDiff;
+}
+
+
+void GLC::setCWs( const std::string &name, const PropertySet &opts ) {
+    // needs to generate a vector of relevant factors, and a dictionary for mapping between local and global factors
+    _CWs.clear();
+    _CWs.reserve( _INRs.size() );
+    _outOffset.clear();
+    _outOffset.reserve( _INRs.size() );
+    // main factors
+    if( props.verbose >= 3 )
+        cerr << "Setting CWs..." << endl;
+
+    for( size_t R = 0; R < nrCWs(); R++ ) { // for each region
+        map<int, size_t> g2l;
+        vector<Factor> Rfactors; // contains all factors
+        size_t potsize = Rfs(R).size() + M(R).size() + _cavitydists[R].size();
+        Rfactors.reserve( potsize );
+        size_t _I = 0;
+        // adding main factors
+        for( SmallSet<size_t>::const_iterator it = Rfs(R).begin(); it != Rfs(R).end(); it++ ) {
+            Rfactors.push_back( factor(*it) );
+            g2l[(int)(*it)] = _I;
+            _I++;
+        }
+
+        // adding factors for incoming messages
+        size_t m = 0;
+        size_t offset = Rfs(R).size();
+        for( vector<Connection>::const_iterator it = M(R).begin(); it != M(R).end(); it++ ) {
+            Rfactors.push_back( (*it).msg );
+            g2l[(int)-m-1] = m + offset;
+            m++;
+        }
+        // cavities
+        offset = Rfs(R).size() + M(R).size();
+        for( size_t c = 0; c < _cavitydists[R].size(); c++ ) {
+            Rfactors.push_back( _cavitydists[R][c] );
+            g2l[(int)-M(R).size() - c - 1] = offset + c;
+        }
+        // outgoing msgs
+        offset = Rfs(R).size() + M(R).size() + _cavitydists[R].size();
+//      _outOffset.push_back(-M(R).size() - _cavitydists[R].size() - 1);
+        size_t counter = 0;
+        if( props.verbose >= 3 )
+            cerr << "going to the loop!" << endl;
+        for( size_t c = 0; c < outM(R).size(); c++ ) {
+            bool isAvailable = false;
+            for( size_t I = 0; I < Rfactors.size(); I++ )
+                if( outM(R,c) << Rfactors[I].vars() ) {
+                    isAvailable = true;
+                    break;
+                }
+            if( !isAvailable ) {
+                Rfactors.push_back( Factor( outM(R)[c], (Real)1 ) );
+                g2l[(int)-M(R).size() - _cavitydists[R].size() - 1 - counter] = offset + counter;
+                counter++;
+            }
+        }
+
+        _CWs.push_back( Cobweb(g2l) );
+        InfAlg *alg = newInfAlg( name, FactorGraph(Rfactors), opts );
+        DAIAlgFG *dalg = static_cast<DAIAlgFG*> (alg);
+        _CWs[R].setInfAlg( dalg );
+    }
+}
+
+
+void GLC::initCWs(){
+    for( size_t R = 0; R < nrCWs(); R++ ) {
+        _CWs[R].initialize();
+        for( size_t m = 0; m < M(R).size(); m++ ) {
+            M(R)[m].msg.fill( (Real)1 );
+            M(R)[m].newmsg.fill( (Real)1 );
+        }
+    }
+}
+
+
+} // end of namespace dai
index 4df765f..2e71762 100644 (file)
@@ -181,3 +181,19 @@ BBP:                            CBP[choose=CHOOSE_BBP]
 # --- DECMAP ------------------
 
 DECMAP:                                DECMAP[ianame=BP,iaopts=[inference=MAXPROD,updates=SEQRND,logdomain=1,tol=1e-9,maxiter=10000,damping=0.1,verbose=0],reinit=1,verbose=0]
+
+#-------GLC(+)--------------
+GLC_UNICAV_SINGLE:              GLC[rgntype=SINGLE,cavity=UNIFORM,updates=SEQRND,maxiter=10000,tol=1e-9,cavainame=BP,cavaiopts=[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0],inainame=EXACT,inaiopts=[],tol=1e-9,verbose=0]
+GLC_PAIRCAV_SINGLE:             GLC[rgntype=SINGLE,cavity=PAIR,updates=SEQRND,maxiter=10000,tol=1e-9,cavainame=BP,cavaiopts=[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0],inainame=EXACT,inaiopts=[],tol=1e-9]
+GLC_PAIR2CAV_SINGLE:            GLC[rgntype=SINGLE,cavity=PAIR2,updates=SEQRND,maxiter=10000,tol=1e-9,cavainame=BP,cavaiopts=[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0],inainame=EXACT,inaiopts=[],tol=1e-9]
+GLC_FULLCAV_SINGLE:             GLC[rgntype=SINGLE,cavity=FULL,updates=SEQRND,maxiter=10000,tol=1e-9,cavainame=BP,cavaiopts=[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0],inainame=EXACT,inaiopts=[],tol=1e-9]
+GLC+_UNICAV_FACTOR:             GLC[rgntype=OVFACTOR,cavity=UNIFORM,updates=SEQRND,maxiter=10000,tol=1e-9,cavainame=BP,cavaiopts=[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0],inainame=JTREE,inaiopts=[inference=SUMPROD,updates=HUGIN],tol=1e-9]
+GLC+_FULLCAV_FACTOR:            GLC[rgntype=OVFACTOR,cavity=FULL,updates=SEQRND,maxiter=10000,tol=1e-9,cavainame=BP,cavaiopts=[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0],inainame=EXACT,inaiopts=[],tol=1e-9]
+GLC+_UNICAV_LOOP3:              GLC[rgntype=OVLOOP,loopdepth=3,cavity=UNIFORM,updates=SEQRND,maxiter=10000,tol=1e-9,cavainame=BP,cavaiopts=[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0],inainame=JTREE,inaiopts=[inference=SUMPROD,updates=HUGIN],tol=1e-9]
+GLC+_FULLCAV_LOOP3:             GLC[rgntype=OVLOOP,loopdepth=3,cavity=FULL,updates=SEQRND,maxiter=10000,tol=1e-9,cavainame=BP,cavaiopts=[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0],inainame=EXACT,inaiopts=[],tol=1e-9]
+GLC_UNICAV_LOOP4:               GLC[rgntype=LOOP,loopdepth=4,cavity=UNIFORM,updates=SEQRND,maxiter=10000,tol=1e-9,cavainame=BP,cavaiopts=[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0],inainame=JTREE,inaiopts=[inference=SUMPROD,updates=HUGIN],tol=1e-9]
+GLC_FULLCAV_LOOP4:              GLC[rgntype=LOOP,loopdepth=4,cavity=UNIFORM,updates=SEQRND,maxiter=10000,tol=1e-9,cavainame=BP,cavaiopts=[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0],inainame=EXACT,inaiopts=[],tol=1e-9]
+GLC+_UNICAV_LOOP4:              GLC[rgntype=OVLOOP,loopdepth=4,cavity=UNIFORM,updates=SEQRND,maxiter=10000,tol=1e-9,cavainame=BP,cavaiopts=[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0],inainame=JTREE,inaiopts=[inference=SUMPROD,updates=HUGIN],tol=1e-9]
+GLC+_FULLCAV_LOOP4:             GLC[rgntype=OVLOOP,loopdepth=4,cavity=FULL,updates=SEQRND,maxiter=10000,tol=1e-9,cavainame=BP,cavaiopts=[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0],inainame=EXACT,inaiopts=[],tol=1e-9]
+GLC+_UNICAV_LOOP5:              GLC[rgntype=OVLOOP,loopdepth=5,cavity=UNIFORM,updates=SEQRND,maxiter=10000,tol=1e-9,cavainame=BP,cavaiopts=[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0],inainame=JTREE,inaiopts=[inference=SUMPROD,updates=HUGIN],tol=1e-9]
+GLC+_FULLCAV_LOOP5:             GLC[rgntype=OVLOOP,loopdepth=5,cavity=FULL,updates=SEQRND,maxiter=10000,tol=1e-9,cavainame=BP,cavaiopts=[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0],inainame=EXACT,inaiopts=[],tol=1e-9]
index 5607188..8cbbf07 100755 (executable)
@@ -1,6 +1,6 @@
 #!/bin/bash
 # Marginal inference
-./testdai --report-iters false --report-time false --marginals VAR --aliases aliases.conf --filename $1 --methods EXACT JTREE_MINFILL_HUGIN JTREE_MINFILL_SHSH JTREE_WEIGHTEDMINFILL_HUGIN JTREE_WEIGHTEDMINFILL_SHSH JTREE_MINWEIGHT_HUGIN JTREE_MINWEIGHT_SHSH JTREE_MINNEIGHBORS_HUGIN JTREE_MINNEIGHBORS_SHSH BP BP_SEQFIX BP_SEQRND BP_SEQMAX BP_PARALL BP_SEQFIX_LOG BP_SEQRND_LOG BP_SEQMAX_LOG BP_PARALL_LOG FBP FBP_SEQFIX FBP_SEQRND FBP_SEQMAX FBP_PARALL FBP_SEQFIX_LOG FBP_SEQRND_LOG FBP_SEQMAX_LOG FBP_PARALL_LOG TRWBP TRWBP_SEQFIX TRWBP_SEQRND TRWBP_SEQMAX TRWBP_PARALL TRWBP_SEQFIX_LOG TRWBP_SEQRND_LOG TRWBP_SEQMAX_LOG TRWBP_PARALL_LOG MF MF_NAIVE_UNI MF_NAIVE_RND MF_HARDSPIN_UNI MF_HARDSPIN_RND TREEEP TREEEPWC GBP_MIN GBP_BETHE GBP_LOOP3 HAK_MIN HAK_BETHE 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 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 LCTREEEP BBP
+./testdai --report-iters false --report-time false --marginals VAR --aliases aliases.conf --filename $1 --methods EXACT JTREE_MINFILL_HUGIN JTREE_MINFILL_SHSH JTREE_WEIGHTEDMINFILL_HUGIN JTREE_WEIGHTEDMINFILL_SHSH JTREE_MINWEIGHT_HUGIN JTREE_MINWEIGHT_SHSH JTREE_MINNEIGHBORS_HUGIN JTREE_MINNEIGHBORS_SHSH BP BP_SEQFIX BP_SEQRND BP_SEQMAX BP_PARALL BP_SEQFIX_LOG BP_SEQRND_LOG BP_SEQMAX_LOG BP_PARALL_LOG FBP FBP_SEQFIX FBP_SEQRND FBP_SEQMAX FBP_PARALL FBP_SEQFIX_LOG FBP_SEQRND_LOG FBP_SEQMAX_LOG FBP_PARALL_LOG TRWBP TRWBP_SEQFIX TRWBP_SEQRND TRWBP_SEQMAX TRWBP_PARALL TRWBP_SEQFIX_LOG TRWBP_SEQRND_LOG TRWBP_SEQMAX_LOG TRWBP_PARALL_LOG MF MF_NAIVE_UNI MF_NAIVE_RND MF_HARDSPIN_UNI MF_HARDSPIN_RND TREEEP TREEEPWC GBP_MIN GBP_BETHE GBP_LOOP3 HAK_MIN HAK_BETHE 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 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 LCTREEEP BBP GLC_UNICAV_SINGLE GLC_PAIRCAV_SINGLE GLC_PAIR2CAV_SINGLE GLC_FULLCAV_SINGLE GLC+_UNICAV_FACTOR GLC+_FULLCAV_FACTOR GLC+_UNICAV_LOOP3 GLC+_FULLCAV_LOOP3 GLC_UNICAV_LOOP4 GLC_FULLCAV_LOOP4 GLC+_UNICAV_LOOP4 GLC+_FULLCAV_LOOP4 GLC+_UNICAV_LOOP5 GLC+_FULLCAV_LOOP5
 # GBP_DELTA, GBP_LOOP4, GBP_LOOP5, GBP_LOOP6, GBP_LOOP7 misbehave
 # MAP inference
 ./testdai --report-iters false --report-time false --marginals VAR --aliases aliases.conf --filename $1 --methods JTREE_MINFILL_HUGIN_MAP JTREE_MINFILL_SHSH_MAP JTREE_WEIGHTEDMINFILL_HUGIN_MAP JTREE_WEIGHTEDMINFILL_SHSH_MAP JTREE_MINWEIGHT_HUGIN_MAP JTREE_MINWEIGHT_SHSH_MAP JTREE_MINNEIGHBORS_HUGIN_MAP JTREE_MINNEIGHBORS_SHSH_MAP MP_SEQFIX MP_SEQRND MP_PARALL MP_SEQFIX_LOG MP_SEQRND_LOG MP_PARALL_LOG FMP_SEQFIX FMP_SEQRND FMP_PARALL FMP_SEQFIX_LOG FMP_SEQRND_LOG FMP_PARALL_LOG TRWMP_SEQFIX TRWMP_SEQRND TRWMP_PARALL TRWMP_SEQFIX_LOG TRWMP_SEQRND_LOG TRWMP_PARALL_LOG DECMAP
index f2931aa..e983c99 100755 (executable)
@@ -1,6 +1,6 @@
 @ECHO OFF
 REM Marginal inference
-@testdai --report-iters false --report-time false --marginals VAR --aliases aliases.conf --filename %1 --methods EXACT JTREE_MINFILL_HUGIN JTREE_MINFILL_SHSH JTREE_WEIGHTEDMINFILL_HUGIN JTREE_WEIGHTEDMINFILL_SHSH JTREE_MINWEIGHT_HUGIN JTREE_MINWEIGHT_SHSH JTREE_MINNEIGHBORS_HUGIN JTREE_MINNEIGHBORS_SHSH BP BP_SEQFIX BP_SEQRND BP_SEQMAX BP_PARALL BP_SEQFIX_LOG BP_SEQRND_LOG BP_SEQMAX_LOG BP_PARALL_LOG FBP FBP_SEQFIX FBP_SEQRND FBP_SEQMAX FBP_PARALL FBP_SEQFIX_LOG FBP_SEQRND_LOG FBP_SEQMAX_LOG FBP_PARALL_LOG TRWBP TRWBP_SEQFIX TRWBP_SEQRND TRWBP_SEQMAX TRWBP_PARALL TRWBP_SEQFIX_LOG TRWBP_SEQRND_LOG TRWBP_SEQMAX_LOG TRWBP_PARALL_LOG MF MF_NAIVE_UNI MF_NAIVE_RND MF_HARDSPIN_UNI MF_HARDSPIN_RND TREEEP TREEEPWC GBP_MIN GBP_BETHE GBP_LOOP3 HAK_MIN HAK_BETHE 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 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 LCTREEEP BBP
+@testdai --report-iters false --report-time false --marginals VAR --aliases aliases.conf --filename %1 --methods EXACT JTREE_MINFILL_HUGIN JTREE_MINFILL_SHSH JTREE_WEIGHTEDMINFILL_HUGIN JTREE_WEIGHTEDMINFILL_SHSH JTREE_MINWEIGHT_HUGIN JTREE_MINWEIGHT_SHSH JTREE_MINNEIGHBORS_HUGIN JTREE_MINNEIGHBORS_SHSH BP BP_SEQFIX BP_SEQRND BP_SEQMAX BP_PARALL BP_SEQFIX_LOG BP_SEQRND_LOG BP_SEQMAX_LOG BP_PARALL_LOG FBP FBP_SEQFIX FBP_SEQRND FBP_SEQMAX FBP_PARALL FBP_SEQFIX_LOG FBP_SEQRND_LOG FBP_SEQMAX_LOG FBP_PARALL_LOG TRWBP TRWBP_SEQFIX TRWBP_SEQRND TRWBP_SEQMAX TRWBP_PARALL TRWBP_SEQFIX_LOG TRWBP_SEQRND_LOG TRWBP_SEQMAX_LOG TRWBP_PARALL_LOG MF MF_NAIVE_UNI MF_NAIVE_RND MF_HARDSPIN_UNI MF_HARDSPIN_RND TREEEP TREEEPWC GBP_MIN GBP_BETHE GBP_LOOP3 HAK_MIN HAK_BETHE 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 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 LCTREEEP BBP GLC_UNICAV_SINGLE GLC_PAIRCAV_SINGLE GLC_PAIR2CAV_SINGLE GLC_FULLCAV_SINGLE GLC+_UNICAV_FACTOR GLC+_FULLCAV_FACTOR GLC+_UNICAV_LOOP3 GLC+_FULLCAV_LOOP3 GLC_UNICAV_LOOP4 GLC_FULLCAV_LOOP4 GLC+_UNICAV_LOOP4 GLC+_FULLCAV_LOOP4 GLC+_UNICAV_LOOP5 GLC+_FULLCAV_LOOP5
 REM GBP_DELTA, GBP_LOOP4, GBP_LOOP5, GBP_LOOP6, GBP_LOOP7 misbehave
 
 REM MAP inference
index a157b49..51b1500 100644 (file)
@@ -1377,6 +1377,244 @@ BBP                                     1.000e-09       1.000e-09       1.000e-09       1.000e-09
 # ({x13}, (9.038e-01, 9.617e-02))
 # ({x14}, (2.408e-01, 7.592e-01))
 # ({x15}, (6.910e-01, 3.090e-01))
+GLC_UNICAV_SINGLE                              8.924e-03       3.480e-03       N/A             N/A             N/A             1.000e-09       
+# ({x0}, (3.486e-01, 6.514e-01))
+# ({x1}, (6.432e-01, 3.568e-01))
+# ({x2}, (5.007e-01, 4.993e-01))
+# ({x3}, (3.027e-01, 6.973e-01))
+# ({x4}, (3.661e-01, 6.339e-01))
+# ({x5}, (6.415e-01, 3.585e-01))
+# ({x6}, (5.819e-01, 4.181e-01))
+# ({x7}, (5.445e-01, 4.555e-01))
+# ({x8}, (2.718e-01, 7.282e-01))
+# ({x9}, (7.144e-01, 2.856e-01))
+# ({x10}, (5.711e-01, 4.289e-01))
+# ({x11}, (5.339e-01, 4.661e-01))
+# ({x12}, (3.515e-01, 6.485e-01))
+# ({x13}, (9.038e-01, 9.620e-02))
+# ({x14}, (2.497e-01, 7.503e-01))
+# ({x15}, (6.859e-01, 3.141e-01))
+GLC_PAIRCAV_SINGLE                             1.355e-03       5.520e-04       N/A             N/A             N/A             1.000e-09       
+# ({x0}, (3.498e-01, 6.502e-01))
+# ({x1}, (6.455e-01, 3.545e-01))
+# ({x2}, (5.000e-01, 5.000e-01))
+# ({x3}, (3.051e-01, 6.949e-01))
+# ({x4}, (3.694e-01, 6.306e-01))
+# ({x5}, (6.414e-01, 3.586e-01))
+# ({x6}, (5.804e-01, 4.196e-01))
+# ({x7}, (5.433e-01, 4.567e-01))
+# ({x8}, (2.795e-01, 7.205e-01))
+# ({x9}, (7.093e-01, 2.907e-01))
+# ({x10}, (5.788e-01, 4.212e-01))
+# ({x11}, (5.380e-01, 4.620e-01))
+# ({x12}, (3.540e-01, 6.460e-01))
+# ({x13}, (9.041e-01, 9.593e-02))
+# ({x14}, (2.406e-01, 7.594e-01))
+# ({x15}, (6.909e-01, 3.091e-01))
+GLC_PAIR2CAV_SINGLE                            1.254e-03       5.210e-04       N/A             N/A             N/A             1.000e-09       
+# ({x0}, (3.499e-01, 6.501e-01))
+# ({x1}, (6.454e-01, 3.546e-01))
+# ({x2}, (5.000e-01, 5.000e-01))
+# ({x3}, (3.051e-01, 6.949e-01))
+# ({x4}, (3.695e-01, 6.305e-01))
+# ({x5}, (6.413e-01, 3.587e-01))
+# ({x6}, (5.804e-01, 4.196e-01))
+# ({x7}, (5.433e-01, 4.567e-01))
+# ({x8}, (2.796e-01, 7.204e-01))
+# ({x9}, (7.093e-01, 2.907e-01))
+# ({x10}, (5.788e-01, 4.212e-01))
+# ({x11}, (5.380e-01, 4.620e-01))
+# ({x12}, (3.540e-01, 6.460e-01))
+# ({x13}, (9.041e-01, 9.592e-02))
+# ({x14}, (2.405e-01, 7.595e-01))
+# ({x15}, (6.909e-01, 3.091e-01))
+GLC_FULLCAV_SINGLE                             1.671e-04       3.175e-05       N/A             N/A             N/A             1.000e-09       
+# ({x0}, (3.500e-01, 6.500e-01))
+# ({x1}, (6.447e-01, 3.553e-01))
+# ({x2}, (4.997e-01, 5.003e-01))
+# ({x3}, (3.049e-01, 6.951e-01))
+# ({x4}, (3.698e-01, 6.302e-01))
+# ({x5}, (6.401e-01, 3.599e-01))
+# ({x6}, (5.793e-01, 4.207e-01))
+# ({x7}, (5.437e-01, 4.563e-01))
+# ({x8}, (2.798e-01, 7.202e-01))
+# ({x9}, (7.084e-01, 2.916e-01))
+# ({x10}, (5.776e-01, 4.224e-01))
+# ({x11}, (5.375e-01, 4.625e-01))
+# ({x12}, (3.541e-01, 6.459e-01))
+# ({x13}, (9.038e-01, 9.615e-02))
+# ({x14}, (2.408e-01, 7.592e-01))
+# ({x15}, (6.910e-01, 3.090e-01))
+GLC+_UNICAV_FACTOR                             2.790e-03       6.317e-04       9.520e-03       2.247e-03       N/A             1.000e-09       
+# ({x0}, (3.500e-01, 6.500e-01))
+# ({x1}, (6.444e-01, 3.556e-01))
+# ({x2}, (5.002e-01, 4.998e-01))
+# ({x3}, (3.054e-01, 6.946e-01))
+# ({x4}, (3.699e-01, 6.301e-01))
+# ({x5}, (6.398e-01, 3.601e-01))
+# ({x6}, (5.803e-01, 4.197e-01))
+# ({x7}, (5.435e-01, 4.565e-01))
+# ({x8}, (2.800e-01, 7.199e-01))
+# ({x9}, (7.086e-01, 2.913e-01))
+# ({x10}, (5.748e-01, 4.252e-01))
+# ({x11}, (5.360e-01, 4.640e-01))
+# ({x12}, (3.542e-01, 6.458e-01))
+# ({x13}, (9.037e-01, 9.631e-02))
+# ({x14}, (2.427e-01, 7.573e-01))
+# ({x15}, (6.915e-01, 3.085e-01))
+GLC+_FULLCAV_FACTOR                            3.871e-05       9.116e-06       1.157e-04       1.534e-05       N/A             1.000e-09       
+# ({x0}, (3.500e-01, 6.500e-01))
+# ({x1}, (6.447e-01, 3.553e-01))
+# ({x2}, (4.997e-01, 5.003e-01))
+# ({x3}, (3.049e-01, 6.951e-01))
+# ({x4}, (3.699e-01, 6.301e-01))
+# ({x5}, (6.401e-01, 3.599e-01))
+# ({x6}, (5.793e-01, 4.207e-01))
+# ({x7}, (5.437e-01, 4.563e-01))
+# ({x8}, (2.800e-01, 7.200e-01))
+# ({x9}, (7.084e-01, 2.916e-01))
+# ({x10}, (5.776e-01, 4.224e-01))
+# ({x11}, (5.375e-01, 4.625e-01))
+# ({x12}, (3.542e-01, 6.458e-01))
+# ({x13}, (9.038e-01, 9.617e-02))
+# ({x14}, (2.408e-01, 7.592e-01))
+# ({x15}, (6.910e-01, 3.090e-01))
+GLC+_UNICAV_LOOP3                              2.790e-03       6.317e-04       9.520e-03       2.247e-03       N/A             1.000e-09       
+# ({x0}, (3.500e-01, 6.500e-01))
+# ({x1}, (6.444e-01, 3.556e-01))
+# ({x2}, (5.002e-01, 4.998e-01))
+# ({x3}, (3.054e-01, 6.946e-01))
+# ({x4}, (3.699e-01, 6.301e-01))
+# ({x5}, (6.398e-01, 3.601e-01))
+# ({x6}, (5.803e-01, 4.197e-01))
+# ({x7}, (5.435e-01, 4.565e-01))
+# ({x8}, (2.800e-01, 7.199e-01))
+# ({x9}, (7.086e-01, 2.913e-01))
+# ({x10}, (5.748e-01, 4.252e-01))
+# ({x11}, (5.360e-01, 4.640e-01))
+# ({x12}, (3.542e-01, 6.458e-01))
+# ({x13}, (9.037e-01, 9.631e-02))
+# ({x14}, (2.427e-01, 7.573e-01))
+# ({x15}, (6.915e-01, 3.085e-01))
+GLC+_FULLCAV_LOOP3                             3.871e-05       9.116e-06       1.157e-04       1.534e-05       N/A             1.000e-09       
+# ({x0}, (3.500e-01, 6.500e-01))
+# ({x1}, (6.447e-01, 3.553e-01))
+# ({x2}, (4.997e-01, 5.003e-01))
+# ({x3}, (3.049e-01, 6.951e-01))
+# ({x4}, (3.699e-01, 6.301e-01))
+# ({x5}, (6.401e-01, 3.599e-01))
+# ({x6}, (5.793e-01, 4.207e-01))
+# ({x7}, (5.437e-01, 4.563e-01))
+# ({x8}, (2.800e-01, 7.200e-01))
+# ({x9}, (7.084e-01, 2.916e-01))
+# ({x10}, (5.776e-01, 4.224e-01))
+# ({x11}, (5.375e-01, 4.625e-01))
+# ({x12}, (3.542e-01, 6.458e-01))
+# ({x13}, (9.038e-01, 9.617e-02))
+# ({x14}, (2.408e-01, 7.592e-01))
+# ({x15}, (6.910e-01, 3.090e-01))
+GLC_UNICAV_LOOP4                               5.304e-03       1.690e-03       N/A             N/A             N/A             1.000e-09       
+# ({x0}, (3.512e-01, 6.488e-01))
+# ({x1}, (6.430e-01, 3.570e-01))
+# ({x2}, (4.992e-01, 5.008e-01))
+# ({x3}, (3.049e-01, 6.951e-01))
+# ({x4}, (3.729e-01, 6.271e-01))
+# ({x5}, (6.365e-01, 3.635e-01))
+# ({x6}, (5.773e-01, 4.227e-01))
+# ({x7}, (5.437e-01, 4.563e-01))
+# ({x8}, (2.842e-01, 7.158e-01))
+# ({x9}, (7.030e-01, 2.970e-01))
+# ({x10}, (5.756e-01, 4.244e-01))
+# ({x11}, (5.363e-01, 4.637e-01))
+# ({x12}, (3.555e-01, 6.445e-01))
+# ({x13}, (9.035e-01, 9.653e-02))
+# ({x14}, (2.412e-01, 7.588e-01))
+# ({x15}, (6.913e-01, 3.087e-01))
+GLC_FULLCAV_LOOP4                              5.304e-03       1.690e-03       N/A             N/A             N/A             1.000e-09       
+# ({x0}, (3.512e-01, 6.488e-01))
+# ({x1}, (6.430e-01, 3.570e-01))
+# ({x2}, (4.992e-01, 5.008e-01))
+# ({x3}, (3.049e-01, 6.951e-01))
+# ({x4}, (3.729e-01, 6.271e-01))
+# ({x5}, (6.365e-01, 3.635e-01))
+# ({x6}, (5.773e-01, 4.227e-01))
+# ({x7}, (5.437e-01, 4.563e-01))
+# ({x8}, (2.842e-01, 7.158e-01))
+# ({x9}, (7.030e-01, 2.970e-01))
+# ({x10}, (5.756e-01, 4.244e-01))
+# ({x11}, (5.363e-01, 4.637e-01))
+# ({x12}, (3.555e-01, 6.445e-01))
+# ({x13}, (9.035e-01, 9.653e-02))
+# ({x14}, (2.412e-01, 7.588e-01))
+# ({x15}, (6.913e-01, 3.087e-01))
+GLC+_UNICAV_LOOP4                              5.578e-06       1.980e-06       3.980e-05       5.551e-06       N/A             1.000e-09       
+# ({x0}, (3.500e-01, 6.500e-01))
+# ({x1}, (6.447e-01, 3.553e-01))
+# ({x2}, (4.997e-01, 5.003e-01))
+# ({x3}, (3.049e-01, 6.951e-01))
+# ({x4}, (3.699e-01, 6.301e-01))
+# ({x5}, (6.401e-01, 3.599e-01))
+# ({x6}, (5.793e-01, 4.207e-01))
+# ({x7}, (5.437e-01, 4.563e-01))
+# ({x8}, (2.800e-01, 7.200e-01))
+# ({x9}, (7.083e-01, 2.917e-01))
+# ({x10}, (5.776e-01, 4.224e-01))
+# ({x11}, (5.375e-01, 4.625e-01))
+# ({x12}, (3.542e-01, 6.458e-01))
+# ({x13}, (9.038e-01, 9.616e-02))
+# ({x14}, (2.408e-01, 7.592e-01))
+# ({x15}, (6.910e-01, 3.090e-01))
+GLC+_FULLCAV_LOOP4                             1.512e-07       2.111e-08       3.977e-07       3.838e-08       N/A             1.000e-09       
+# ({x0}, (3.500e-01, 6.500e-01))
+# ({x1}, (6.447e-01, 3.553e-01))
+# ({x2}, (4.997e-01, 5.003e-01))
+# ({x3}, (3.049e-01, 6.951e-01))
+# ({x4}, (3.699e-01, 6.301e-01))
+# ({x5}, (6.401e-01, 3.599e-01))
+# ({x6}, (5.793e-01, 4.207e-01))
+# ({x7}, (5.437e-01, 4.563e-01))
+# ({x8}, (2.800e-01, 7.200e-01))
+# ({x9}, (7.083e-01, 2.917e-01))
+# ({x10}, (5.776e-01, 4.224e-01))
+# ({x11}, (5.375e-01, 4.625e-01))
+# ({x12}, (3.542e-01, 6.458e-01))
+# ({x13}, (9.038e-01, 9.617e-02))
+# ({x14}, (2.408e-01, 7.592e-01))
+# ({x15}, (6.910e-01, 3.090e-01))
+GLC+_UNICAV_LOOP5                              5.578e-06       1.980e-06       3.980e-05       5.551e-06       N/A             1.000e-09       
+# ({x0}, (3.500e-01, 6.500e-01))
+# ({x1}, (6.447e-01, 3.553e-01))
+# ({x2}, (4.997e-01, 5.003e-01))
+# ({x3}, (3.049e-01, 6.951e-01))
+# ({x4}, (3.699e-01, 6.301e-01))
+# ({x5}, (6.401e-01, 3.599e-01))
+# ({x6}, (5.793e-01, 4.207e-01))
+# ({x7}, (5.437e-01, 4.563e-01))
+# ({x8}, (2.800e-01, 7.200e-01))
+# ({x9}, (7.083e-01, 2.917e-01))
+# ({x10}, (5.776e-01, 4.224e-01))
+# ({x11}, (5.375e-01, 4.625e-01))
+# ({x12}, (3.542e-01, 6.458e-01))
+# ({x13}, (9.038e-01, 9.616e-02))
+# ({x14}, (2.408e-01, 7.592e-01))
+# ({x15}, (6.910e-01, 3.090e-01))
+GLC+_FULLCAV_LOOP5                             1.513e-07       2.113e-08       3.977e-07       3.838e-08       N/A             1.000e-09       
+# ({x0}, (3.500e-01, 6.500e-01))
+# ({x1}, (6.447e-01, 3.553e-01))
+# ({x2}, (4.997e-01, 5.003e-01))
+# ({x3}, (3.049e-01, 6.951e-01))
+# ({x4}, (3.699e-01, 6.301e-01))
+# ({x5}, (6.401e-01, 3.599e-01))
+# ({x6}, (5.793e-01, 4.207e-01))
+# ({x7}, (5.437e-01, 4.563e-01))
+# ({x8}, (2.800e-01, 7.200e-01))
+# ({x9}, (7.083e-01, 2.917e-01))
+# ({x10}, (5.776e-01, 4.224e-01))
+# ({x11}, (5.375e-01, 4.625e-01))
+# ({x12}, (3.542e-01, 6.458e-01))
+# ({x13}, (9.038e-01, 9.617e-02))
+# ({x14}, (2.408e-01, 7.592e-01))
+# ({x15}, (6.910e-01, 3.090e-01))
 # testfast.fg
 # METHOD                                       MAX VAR ERR     AVG VAR ERR     MAX FAC ERR     AVG FAC ERR     LOGZ ERROR      MAXDIFF 
 JTREE_MINFILL_HUGIN_MAP