include/dai/jtree.h

   1 /*  This file is part of libDAI - http://www.libdai.org/
   2  *
   3  *  libDAI is licensed under the terms of the GNU General Public License version
   4  *  2, or (at your option) any later version. libDAI is distributed without any
   5  *  warranty. See the file COPYING for more details.
   6  *
   7  *  Copyright (C) 2006-2009  Joris Mooij  [joris dot mooij at libdai dot org]
   8  *  Copyright (C) 2006-2007  Radboud University Nijmegen, The Netherlands
   9  */
  10
  11
  12 /// \file
  13 /// \brief Defines class JTree, which implements the junction tree algorithm
  14
  15
  16 #ifndef __defined_libdai_jtree_h
  17 #define __defined_libdai_jtree_h
  18
  19
  20 #include <vector>
  21 #include <string>
  22 #include <dai/daialg.h>
  23 #include <dai/varset.h>
  24 #include <dai/regiongraph.h>
  25 #include <dai/factorgraph.h>
  26 #include <dai/clustergraph.h>
  27 #include <dai/weightedgraph.h>
  28 #include <dai/enum.h>
  29 #include <dai/properties.h>
  30
  31
  32 namespace dai {
  33
  34
  35 /// Exact inference algorithm using junction tree
  36 /** The junction tree algorithm uses message passing on a junction tree to calculate
  37  *  exact marginal probability distributions ("beliefs") for specified cliques
  38  *  (outer regions) and separators (intersections of pairs of cliques).
  39  *
  40  *  There are two variants, the sum-product algorithm (corresponding to
  41  *  finite temperature) and the max-product algorithm (corresponding to
  42  *  zero temperature).
  43  */
  44 class JTree : public DAIAlgRG {
  45     private:
  46         /// Stores the messages
  47         std::vector<std::vector<Factor> >  _mes;
  48
  49         /// Stores the logarithm of the partition sum
  50         Real _logZ;
  51
  52     public:
  53         /// The junction tree (stored as a rooted tree)
  54         RootedTree RTree;
  55
  56         /// Outer region beliefs
  57         std::vector<Factor> Qa;
  58
  59         /// Inner region beliefs
  60         std::vector<Factor> Qb;
  61
  62         /// Parameters for JTree
  63         struct Properties {
  64             /// Enumeration of possible JTree updates
  65             /** There are two types of updates:
  66              *  - HUGIN similar to those in HUGIN
  67              *  - SHSH Shafer-Shenoy type
  68              */
  69             DAI_ENUM(UpdateType,HUGIN,SHSH);
  70
  71             /// Enumeration of inference variants
  72             /** There are two inference variants:
  73              *  - SUMPROD Sum-Product
  74              *  - MAXPROD Max-Product (equivalent to Min-Sum)
  75              */
  76             DAI_ENUM(InfType,SUMPROD,MAXPROD);
  77
  78             /// Verbosity (amount of output sent to stderr)
  79             size_t verbose;
  80
  81             /// Type of updates
  82             UpdateType updates;
  83
  84             /// Type of inference
  85             InfType inference;
  86         } props;
  87
  88         /// Name of this inference algorithm
  89         static const char *Name;
  90
  91     public:
  92     /// \name Constructors/destructors
  93     //@{
  94         /// Default constructor
  95         JTree() : DAIAlgRG(), _mes(), _logZ(), RTree(), Qa(), Qb(), props() {}
  96
  97         /// Construct from FactorGraph \a fg and PropertySet \a opts
  98         /** \param fg factor graph (which has to be connected);
  99          ** \param opts Parameters @see Properties
 100          *  \param automatic if \c true, construct the junction tree automatically, using the MinFill heuristic.
 101          *  \throw FACTORGRAPH_NOT_CONNECTED if \a fg is not connected
 102          */
 103         JTree( const FactorGraph &fg, const PropertySet &opts, bool automatic=true );
 104     //@}
 105
 106
 107     /// \name General InfAlg interface
 108     //@{
 109         virtual JTree* clone() const { return new JTree(*this); }
 110         virtual std::string identify() const;
 111         virtual Factor belief( const VarSet &vs ) const;
 112         virtual std::vector<Factor> beliefs() const;
 113         virtual Real logZ() const;
 114         virtual void init() {}
 115         virtual void init( const VarSet &/*ns*/ ) {}
 116         virtual Real run();
 117         virtual Real maxDiff() const { return 0.0; }
 118         virtual size_t Iterations() const { return 1UL; }
 119         virtual void setProperties( const PropertySet &opts );
 120         virtual PropertySet getProperties() const;
 121         virtual std::string printProperties() const;
 122     //@}
 123
 124
 125     /// \name Additional interface specific for JTree
 126     //@{
 127         /// Constructs a junction tree based on the cliques \a cl (corresponding to some elimination sequence).
 128         /** First, constructs a weighted graph, where the nodes are the elements of \a cl, and
 129          *  each edge is weighted with the cardinality of the intersection of the state spaces of the nodes.
 130          *  Then, a maximal spanning tree for this weighted graph is calculated.
 131          *  Finally, a corresponding region graph is built:
 132          *    - the outer regions correspond with the cliques and have counting number 1;
 133          *    - the inner regions correspond with the seperators, i.e., the intersections of two
 134          *      cliques that are neighbors in the spanning tree, and have counting number -1;
 135          *    - inner and outer regions are connected by an edge if the inner region is a
 136          *      seperator for the outer region.
 137          */
 138         void GenerateJT( const std::vector<VarSet> &cl );
 139
 140         /// Returns constant reference to the message from outer region \a alpha to its \a _beta 'th neighboring inner region
 141         const Factor & message( size_t alpha, size_t _beta ) const { return _mes[alpha][_beta]; }
 142         /// Returns reference to the message from outer region \a alpha to its \a _beta 'th neighboring inner region
 143         Factor & message( size_t alpha, size_t _beta ) { return _mes[alpha][_beta]; }
 144
 145         /// Runs junction tree algorithm using HUGIN updates
 146         /** \note The initial messages may be arbitrary.
 147          */
 148         void runHUGIN();
 149
 150         /// Runs junction tree algorithm using Shafer-Shenoy updates
 151         /** \note The initial messages may be arbitrary.
 152          */
 153         void runShaferShenoy();
 154
 155         /// Finds an efficient subtree for calculating the marginal of the variables in \a vs
 156         /** First, the current junction tree is reordered such that it gets as root the clique
 157          *  that has maximal state space overlap with \a vs. Then, the minimal subtree
 158          *  (starting from the root) is identified that contains all the variables in \a vs
 159          *  and also the outer region with index \a PreviousRoot (if specified). Finally,
 160          *  the current junction tree is reordered such that this minimal subtree comes
 161          *  before the other edges, and the size of the minimal subtree is returned.
 162          */
 163         size_t findEfficientTree( const VarSet& vs, RootedTree &Tree, size_t PreviousRoot=(size_t)-1 ) const;
 164
 165         /// Calculates the marginal of a set of variables (using cutset conditioning, if necessary)
 166         /** \pre assumes that run() has been called already
 167          */
 168         Factor calcMarginal( const VarSet& vs );
 169
 170         /// Calculates the joint state of all variables that has maximum probability
 171         /** \pre Assumes that run() has been called and that \a props.inference == \c MAXPROD
 172          */
 173         std::vector<std::size_t> findMaximum() const;
 174     //@}
 175 };
 176
 177
 178 /// Calculates upper bound to the treewidth of a FactorGraph, using the MinFill heuristic
 179 /** \relates JTree
 180  *  \return a pair (number of variables in largest clique, number of states in largest clique)
 181  */
 182 std::pair<size_t,size_t> boundTreewidth( const FactorGraph & fg );
 183
 184
 185 } // end of namespace dai
 186
 187
 188 #endif