include/dai/glc.h

   1 /*  This file is part of libDAI - http://www.libdai.org/
   2  *
   3  *  Copyright (c) 2006-2012, The libDAI authors. All rights reserved.
   4  *
   5  *  Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
   6  */
   7
   8
   9 /// \file
  10 /// \brief Defines classes GLC and Cobweb, which implement the "Generalized Loop Correction method"
  11 /// \todo Fix the init of GLC
  12
  13
  14 #ifndef __defined_libdai_glc_h
  15 #define __defined_libdai_glc_h
  16
  17
  18 #include <dai/dai_config.h>
  19 #ifdef DAI_WITH_GLC
  20
  21
  22 #include <algorithm>
  23 #include <set>
  24 #include <string>
  25 #include <dai/util.h>
  26 #include <dai/daialg.h>
  27 #include <dai/enum.h>
  28 #include <dai/cobwebgraph.h>
  29 #include <dai/properties.h>
  30 #include <dai/exceptions.h>
  31
  32
  33 namespace dai {
  34
  35
  36 /// Represents a region for the GLC algorithm.
  37 /** This region contains all the proper factors from the factor graph. It also contains a factor for incoming messages and the
  38  *  cavity distribution over this region. The outgoing messages also receive a factor mainly so that we can easily calculate the marginal
  39  *  of the region over the outgoing message variables.
  40  *
  41  *  \author Siamak Ravanbakhsh
  42  */
  43 class Cobweb {
  44     public:
  45         /// Pointer to the approximate inference algorithm and factor graph over region R (provided in parameters of GLC)
  46         DAIAlgFG *cav;
  47
  48         /// Global to local factor index conversion.
  49         /** Factor(I) in original graph corresponds to cav->factor(_g2l[I]) in this region
  50          *  the j'th incoming message to this region M(R,j) corresponds to cav->factor(_g2l[-j-1])
  51          *  the j'th cavity factor (for pairwise case),_cavitydists[R][j] corresponds to cav->factor(_g2l[-m-j-1])
  52          *      when m = M(R).size()
  53          *      note that for full cavity the _cavitydists[R] has a single factor and for uniform cavity it is empty
  54          *  the j'th outgoing message from this region outM(R,j) corresponds to cav->factor(_g2l[-c-m-j-1])
  55          *      when m  = M(R).size() and c = _cavitydists[R].size()
  56          */
  57         std::map<int, size_t> _g2l;
  58
  59         /// Default constructor
  60         Cobweb(){}
  61
  62         /// Construct from global-to-local factor index conversion structure
  63         Cobweb( const std::map<int, size_t>& g2l ): cav(0), _g2l(g2l) {}
  64
  65         /// Default destructor
  66         virtual ~Cobweb(){ delete cav; }
  67
  68         /// Sets the factor graph and inference algorithm for this region
  69         void setInfAlg( DAIAlgFG* alg ) {
  70             cav = alg;
  71         }
  72
  73         /// Returns the factor of this region globaly indexed by \a I (mainly used to access messages to and from this region)
  74         Factor factor( size_t I ) {
  75             return cav->factor( _g2l[I] );
  76         }
  77
  78         /// Initializes the inference algorithm over this region
  79         void initialize() {
  80             cav->init();
  81         }
  82
  83         /// Calculates the marginal over \a ns after temporarily removing the factors indexed by \a rmgind
  84         Factor marginal( const VarSet& ns, const std::vector<size_t>& rmgind ) {
  85             // Set the factors in rmgind to one
  86             VarSet vs;
  87             std::vector<size_t> rmlinds; // local index
  88             rmlinds.reserve( rmgind.size() );
  89             for( std::vector<size_t>::const_iterator it = rmgind.begin(); it != rmgind.end(); it++ ) {
  90                 vs |= cav->factor( _g2l[*it] ).vars();
  91                 rmlinds.push_back( _g2l[*it] );
  92             }
  93             cav->makeRegionCavity( rmlinds, true );
  94             // Initialize if necessary
  95             cav->init( vs );
  96             cav->run();
  97             Factor result;
  98             result = cav->belief( ns );
  99             cav->restoreFactors();
 100             return result;
 101         }
 102
 103         /// Calculates the marginal over the variables in outgoing message indexed by \a outmsgind after temporarily removing the factors indexed by \a rmgind
 104         /** This function is used to calculate outgoing messages from this region
 105          */
 106         Factor marginal( const int& outmsgind, const std::vector<size_t>& rmgind ) {
 107             // set the factors in rmgind to one
 108             VarSet vs;
 109             std::vector<size_t> rmlinds; // local index
 110             rmlinds.reserve( rmgind.size() );
 111             for( std::vector<size_t>::const_iterator it = rmgind.begin(); it != rmgind.end(); it++ ) {
 112                 vs |= cav->factor( _g2l[*it] ).vars();
 113                 rmlinds.push_back( _g2l[*it] );
 114             }
 115             cav->makeRegionCavity( rmlinds, true );
 116             // Initialize if necessary
 117             cav->init( vs );
 118             cav->run();
 119             Factor result;
 120             result = cav->beliefF( _g2l[outmsgind] );
 121             cav->restoreFactors();
 122             return result;
 123         }
 124
 125         /// Sets (or multiplies, if \a multiply == \c true) the factor indexed by \a gind by factor \a msg
 126         void updateFactor( int gind, const Factor& msg, bool multiply=false ) {
 127             if( !multiply )
 128                 cav->setFactor( _g2l[gind], msg, false );
 129             else
 130                 cav->setFactor( _g2l[gind], (msg*(cav->factor( _g2l[gind] ))).normalized(), false );
 131         }
 132
 133         /// Runs inference and returns belief of variable \a v
 134         Factor belief( const Var v ) {
 135             cav->run();
 136             return cav->belief( v );
 137         }
 138
 139         /// Runs inference and returns belief of variables \a vs
 140         Factor belief( const VarSet& vs ) {
 141             cav->run();
 142             return cav->belief( vs );
 143         }
 144 };
 145
 146
 147 /// Approximate inference algorithm "Generalized Loop Correction" [\ref RYG12]
 148 /** Inherits from a CobwebGraph which in turn inherits from a FactorGraph.
 149  *  Although a CobwebGraph represents the graph with cobweb clusters (for design reasons?),
 150  *  the actual vector of cobwebs is contained in GLC class. All the other functionalities
 151  *  needed to represent the cobwebgraph is in its own class. The cobweb provides the function
 152  *  of marginalization using the provided inference method. This can be beneficial when using
 153  *  uniform cavities (or a subset of pairwise cavities, but this is not implemented yet).
 154  *  However, in general the performance does not improve using Junction Tree instead of exact marginalization...
 155  *
 156  *  \author Siamak Ravanbakhsh
 157  */
 158 class GLC : public DAIAlgCG {
 159     private:
 160         /// Stores for each variable the approximate cavity distribution
 161         std::vector<Cobweb> _CWs;
 162
 163         /// Gives the starting global index for the factors of outgoing messages of each region
 164         std::vector<int> _outOffset;
 165
 166         /// Cavity distributions
 167         /** For each region it contains:
 168          *  - Nothing if using UniCAV
 169          *  - One Factor if using FullCAV
 170          *  - All pairwise caivity factors if using PairCAV or Pair2CAV
 171          */
 172         std::vector<std::vector<Factor> > _cavitydists;
 173
 174         /// Single variable beliefs
 175         std::vector<Factor> _beliefs;
 176
 177         /// Beliefs over factors
 178         std::map<VarSet, Factor> _factorBeliefs;
 179
 180         /// Maximum difference encountered so far
 181         Real _maxdiff;
 182
 183         /// Number of iterations needed
 184         size_t _iters;
 185
 186     public:
 187         /// Parameters for GLC
 188         struct Properties {
 189             /// Enumeration of possible ways to initialize the cavities
 190             /** The following initialization methods are defined:
 191              *  - FULL calculates the marginal using calcMarginal()
 192              *  - PAIR calculates only second order interactions using calcPairBeliefs() with \a accurate == \c false
 193              *  - PAIR2 calculates only second order interactions using calcPairBeliefs() with \a accurate == \c true
 194              *  - UNIFORM uses a uniform distribution
 195              */
 196             DAI_ENUM(CavityType,FULL,PAIR,PAIR2,UNIFORM);
 197
 198             /// Enumeration of different possible types of region (\ominus r in the paper)
 199             /** The following initialization methods are defined:
 200              *  - SINGLE using single variable
 201              *  - FACTOR partitioning regions of that form corresponding to GLC (in the paper)
 202              *  - LOOP partitioning regions of that form corresponding to GLC (in the paper)
 203              *  - DELTA partitioning regions of that form corresponding to GLC (in the paper)
 204              *  - OVFACTOR, OVLOOP, OVDELTA the corresponding form of region in the overlapping form corresponding to GLC+ (in the paper)
 205              *  For example, using OVFACTOR each factor becomes a region (\ominus r), while using FACTOR the regions contain only
 206              *  non-overlapping factors and if some variables can't be covered by any remaining factors single variable regions
 207              *  are constructed for them.
 208              *  With this convention OVDELTA here corresponds to DELTA using CVM (HAK).
 209              */
 210             DAI_ENUM(RegionType,SINGLE,FACTOR,OVFACTOR,LOOP,OVLOOP,DELTA,OVDELTA)
 211
 212             /// Enumeration of different update schedules
 213             /** The following update schedules are defined:
 214              *  - SEQFIX sequential fixed schedule
 215              *  - SEQRND sequential random schedule
 216              */
 217             DAI_ENUM(UpdateType,SEQFIX,SEQRND);
 218
 219             /// Verbosity (amount of output sent to stderr)
 220             size_t verbose;
 221
 222             // Maximum length of the loops to use in construction of loopy regions
 223             size_t loopdepth;
 224
 225             /// Maximum number of iterations
 226             size_t maxiter;
 227
 228             /// The algorithm will quit and report the current result at reaching this runtime
 229             Real maxtime;
 230
 231             /// Tolerance for convergence test
 232             Real tol;
 233
 234             /// Complete or partial reinitialization of the messages in calculating the cavity distribution?
 235             bool reinit;
 236
 237             /// Whether or not to include auxiliary factors in a region.
 238             /** The auxilary factors are the ones that only depend on variables on the boundary of a region (\oplus r - \ominus r).
 239              *  This sometimes produces orders of magnitude improvement (e.g., in grids), but it is not discussed in the paper.
 240              */
 241             bool auxfactors;
 242
 243             /// How to initialize the cavities
 244             CavityType cavity;
 245
 246             // Type of region
 247             RegionType rgntype;
 248
 249             /// What update schedule to use
 250             UpdateType updates;
 251
 252             /** Indicates what happens if a subset of variables in the boundary of a region (\ominus r_p) is shared by
 253              *  some neighbors such that one (\ominus r_{p,q1}) is a subset of another (\ominus r_{p,q2}).
 254              *  - ALL     all such neighbors are included in the the updates.
 255              *  - TOP     only (\ominus r_{p,q2}) is included unless (\ominus r_{p,q2} = \ominus r_{p,q1}) in which case both are included
 256              *  - 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
 257              *  \note Not important in perfomance!
 258              */
 259             CobwebGraph::NeighborType neighbors;
 260
 261             /// Name of the algorithm used to initialize the cavity distributions
 262             std::string cavainame;
 263
 264             /// Parameters for the algorithm used to initialize the cavity distributions
 265             PropertySet cavaiopts;
 266
 267             // Name of the algorithm used for marginalization over a region (EXACT is fine)
 268             std::string inainame;
 269
 270             // Parameters for algorithm used for marginalization over a region
 271             PropertySet inaiopts;
 272         } props;
 273
 274
 275     public:
 276         /// Default constructor
 277         GLC() : DAIAlgCG(), _CWs(),_outOffset(), _cavitydists(), _beliefs(), _factorBeliefs(), _maxdiff(), _iters(), props() {}
 278
 279         /// Construct from FactorGraph \a fg and PropertySet \a opts
 280         /** \param fg Factor graph.
 281          *  \param opts Parameters @see Properties
 282          */
 283         GLC( const FactorGraph &fg, const PropertySet &opts );
 284
 285     /// \name General InfAlg interface
 286     //@{
 287         virtual GLC* clone() const { return new GLC(*this); }
 288         virtual GLC* construct( const FactorGraph &fg, const PropertySet &opts ) const { return new GLC( fg, opts ); }
 289         virtual std::string name() const { return "GLC"; }
 290         virtual Factor belief( const Var &v ) const { return beliefV( findVar( v ) ); }
 291         virtual Factor belief( const VarSet &vs ) const;
 292         virtual Factor beliefV( size_t i ) const { return _beliefs[i]; }
 293         virtual std::vector<Factor> beliefs() const;
 294         virtual Real logZ() const { DAI_THROW(NOT_IMPLEMENTED); return 0.0; }
 295         virtual void init(){ initCWs(); }
 296         virtual void init( const VarSet &/*ns*/ ) { init(); }
 297         virtual Real run();
 298         virtual Real maxDiff() const { return _maxdiff; }
 299         virtual size_t Iterations() const { return _iters; }
 300         virtual void setMaxIter( size_t maxiter ) { props.maxiter = maxiter; }
 301         virtual void setProperties( const PropertySet &opts );
 302         virtual PropertySet getProperties() const;
 303         virtual std::string printProperties() const;
 304     //@}
 305
 306     /// \name Additional interface specific for GLC
 307     //@{
 308         /// Returns constant reference to outer region \a alpha
 309         const Cobweb& CW( size_t alpha ) const {
 310             DAI_DEBASSERT( alpha < nrCWs() );
 311             return _CWs[alpha];
 312         }
 313
 314         /// Returns reference to outer region \a alpha
 315         Cobweb& CW( size_t alpha ) {
 316             DAI_DEBASSERT( alpha < nrCWs() );
 317             return _CWs[alpha];
 318         }
 319
 320         /// Approximates the cavity distribution of variable \a i, using the inference algorithm \a name with parameters \a opts
 321         Real CalcCavityDist( size_t i, const std::string &name, const PropertySet &opts );
 322
 323         /// Approximates all cavity distributions using inference algorithm \a name with parameters \a opts
 324         Real InitCavityDists( const std::string &name, const PropertySet &opts );
 325
 326         /// Updates the belief of the Markov blanket of region \a R using its intersection with \a _R2'th neighbor (partitioning case).
 327         void NewPancake( size_t R, size_t _R2 );
 328
 329         /// Updates the belief of the Markov blanket of region \a R using all its overlapping neighbors (non-partitioning case).
 330         /** \note The implementation uses a more compact form than the message passing over neighborhood region-graph explained in the paper.
 331          */
 332         void OVNewPancake( size_t R );
 333
 334         /// 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
 335         void CalcBelief( size_t i, bool isFinal = false );
 336
 337         /// Calculates the belief of factor \a I
 338         void CalcFactorBelief( size_t I );
 339
 340         /// Designates the variables in each region of the graph based on the rgntype
 341         std::vector<SmallSet<size_t> > calcRegions();
 342
 343         /// Initializes the cobweb regions and all the messages between regions.
 344         void initCWs();
 345
 346         /// Initializes all the cobweb regions.
 347         /** Defines their factors using the cavity distribution (as factor(s))
 348          *  and all the anticipated messages to be exchanged with other regions (treating a message as a factor),
 349          *  and also it builds a mapping from cobweb factors to factors of the original factor graph and messages.
 350          */
 351         void setCWs( const std::string &name, const PropertySet &opts );
 352
 353         /// Helper function for partitioning case to build the regions by finding loops of maximum given size.
 354         /** If a variable is contained in a loop it is not considered for inclusion in any other loop.
 355          */
 356         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 );
 357
 358         /// Helper function for general case to build the regions by finding loops of maximum given size.
 359         /** 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.
 360          */
 361         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 );
 362     //@}
 363 };
 364
 365
 366 } // end of namespace dai
 367
 368
 369 #endif
 370
 371
 372 #endif