include/dai/cbp.h

   1 /*  Copyright (C) 2009  Frederik Eaton [frederik at ofb dot net]
   2
   3     This file is part of libDAI.
   4
   5     libDAI is free software; you can redistribute it and/or modify
   6     it under the terms of the GNU General Public License as published by
   7     the Free Software Foundation; either version 2 of the License, or
   8     (at your option) any later version.
   9
  10     libDAI is distributed in the hope that it will be useful,
  11     but WITHOUT ANY WARRANTY; without even the implied warranty of
  12     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13     GNU General Public License for more details.
  14
  15     You should have received a copy of the GNU General Public License
  16     along with libDAI; if not, write to the Free Software
  17     Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  18 */
  19
  20
  21 /// \file
  22 /// \brief Defines class CBP [\ref EaG09]
  23 /// \todo Improve documentation
  24 /// \todo Clean up
  25
  26
  27 #ifndef __defined_libdai_cbp_h
  28 #define __defined_libdai_cbp_h
  29
  30
  31 #include <fstream>
  32
  33 #include <boost/shared_ptr.hpp>
  34
  35 #include <dai/daialg.h>
  36
  37 #include <dai/cbp.h>
  38 #include <dai/bbp.h>
  39
  40
  41 namespace dai {
  42
  43
  44 /// Find a variable to clamp using BBP (goes with maximum adjoint)
  45 /// \see BBP
  46 std::pair<size_t, size_t> bbpFindClampVar(const InfAlg &in_bp, bool clampingVar,
  47             const PropertySet &bbp_props, bbp_cfn_t cfn,
  48             Real *maxVarOut);
  49
  50
  51 /// Class for CBP (Clamped Belief Propagation)
  52 /** This algorithm uses configurable heuristics to choose a variable
  53  *  x_i and a state x_i*. Inference is done with x_i "clamped" to x_i*
  54  *  (e.g. conditional on x_i==x_i*), and also with the negation of this
  55  *  condition. Clamping is done recursively up to a fixed number of
  56  *  levels (other stopping criteria are also implemented, see
  57  *  "recursion" property). The resulting approximate marginals are
  58  *  combined using logZ estimates.
  59  */
  60 class CBP : public DAIAlgFG {
  61     std::vector<Factor> _beliefsV;
  62     std::vector<Factor> _beliefsF;
  63     double _logZ;
  64     double _est_logZ;
  65     double _old_est_logZ;
  66
  67     double _sum_level;
  68     size_t _num_leaves;
  69
  70     boost::shared_ptr<std::ofstream> _clamp_ofstream;
  71
  72     bbp_cfn_t BBP_cost_function() {
  73       return props.bbp_cfn;
  74     }
  75
  76     void printDebugInfo();
  77
  78     /// Called by 'run', and by itself. Implements the main algorithm.
  79     /** Chooses a variable to clamp, recurses, combines the logZ and
  80      *  beliefs estimates of the children, and returns the improved
  81      *  estimates in \a lz_out and \a beliefs_out to its parent
  82      */
  83     void runRecurse(InfAlg *bp,
  84                     double orig_logZ,
  85                     std::vector<size_t> clamped_vars_list,
  86                     size_t &num_leaves,
  87                     size_t &choose_count,
  88                     double &sum_level,
  89                     Real &lz_out,
  90                     std::vector<Factor>& beliefs_out);
  91
  92     /// Choose the next variable to clamp
  93     /** Choose the next variable to clamp, given a converged InfAlg (\a bp),
  94      *  and a vector of variables that are already clamped (\a
  95      *  clamped_vars_list). Returns the chosen variable in \a i, and
  96      *  the set of states in \a xis. If \a maxVarOut is non-NULL and
  97      *  props.choose==CHOOSE_BBP then it is used to store the
  98      *  adjoint of the chosen variable
  99      */
 100     virtual bool chooseNextClampVar(InfAlg* bp,
 101                                     std::vector<size_t> &clamped_vars_list,
 102                                     size_t &i, std::vector<size_t> &xis, Real *maxVarOut);
 103
 104     /// Return the InfAlg to use at each step of the recursion.
 105     /// \todo At present, only returns a BP instance
 106     InfAlg* getInfAlg();
 107
 108     size_t _iters;
 109     double _maxdiff;
 110
 111     void setBeliefs(const std::vector<Factor>& bs, double logZ) {
 112         size_t i=0;
 113         _beliefsV.clear(); _beliefsV.reserve(nrVars());
 114         _beliefsF.clear(); _beliefsF.reserve(nrFactors());
 115         for(i=0; i<nrVars(); i++) _beliefsV.push_back(bs[i]);
 116         for(; i<nrVars()+nrFactors(); i++) _beliefsF.push_back(bs[i]);
 117         _logZ = logZ;
 118     }
 119
 120     void construct();
 121
 122   public:
 123
 124     //----------------------------------------------------------------
 125
 126     /// construct CBP object from FactorGraph
 127     CBP(const FactorGraph &fg, const PropertySet &opts) : DAIAlgFG(fg) {
 128         props.set(opts);
 129         construct();
 130     }
 131
 132     static const char *Name;
 133
 134     /// @name General InfAlg interface
 135     //@{
 136     virtual CBP* clone() const { return new CBP(*this); }
 137     virtual CBP* create() const { DAI_THROW(NOT_IMPLEMENTED); }
 138     virtual std::string identify() const { return std::string(Name) + props.toString(); }
 139     virtual Factor belief (const Var &n) const { return _beliefsV[findVar(n)]; }
 140     virtual Factor belief (const VarSet &) const { DAI_THROW(NOT_IMPLEMENTED); }
 141     virtual std::vector<Factor> beliefs() const { return concat(_beliefsV, _beliefsF); }
 142     virtual Real logZ() const { return _logZ; }
 143     virtual void init() {};
 144     virtual void init( const VarSet & ) {};
 145     virtual double run();
 146     virtual double maxDiff() const { return _maxdiff; }
 147     virtual size_t Iterations() const { return _iters; }
 148     //@}
 149
 150     Factor beliefV (size_t i) const { return _beliefsV[i]; }
 151     Factor beliefF (size_t I) const { return _beliefsF[I]; }
 152
 153     //----------------------------------------------------------------
 154
 155  public:
 156 /* PROPERTIES(props,CBP) {
 157         typedef BP::Properties::UpdateType UpdateType;
 158         DAI_ENUM(RecurseType,REC_FIXED,REC_LOGZ,REC_BDIFF);
 159         DAI_ENUM(ChooseMethodType,CHOOSE_RANDOM,CHOOSE_MAXENT,CHOOSE_BBP,CHOOSE_BP_L1,CHOOSE_BP_CFN);
 160         DAI_ENUM(ClampType,CLAMP_VAR,CLAMP_FACTOR);
 161
 162         size_t verbose = 0;
 163
 164         /// Tolerance to use in BP
 165         double tol;
 166         /// Update style for BP
 167         UpdateType updates;
 168         /// Maximum number of iterations for BP
 169         size_t maxiter;
 170
 171         /// Tolerance to use for controlling recursion depth (\a recurse
 172         /// is REC_LOGZ or REC_BDIFF)
 173         double rec_tol;
 174         /// Maximum number of levels of recursion (\a recurse is REC_FIXED)
 175         size_t max_levels = 10;
 176         /// If choose=CHOOSE_BBP and maximum adjoint is less than this value, don't recurse
 177         double min_max_adj;
 178         /// Heuristic for choosing clamping variable
 179         ChooseMethodType choose;
 180         /// Method for deciding when to stop recursing
 181         RecurseType recursion;
 182         /// Whether to clamp variables or factors
 183         ClampType clamp;
 184         /// Properties to pass to BBP
 185         PropertySet bbp_props;
 186         /// Cost function to use for BBP
 187         bbp_cfn_t bbp_cfn;
 188         /// Random seed
 189         size_t rand_seed = 0;
 190
 191         /// If non-empty, write clamping choices to this file
 192         std::string clamp_outfile = "";
 193    }
 194 */
 195 /* {{{ GENERATED CODE: DO NOT EDIT. Created by
 196     ./scripts/regenerate-properties include/dai/cbp.h src/cbp.cpp
 197 */
 198     struct Properties {
 199         typedef BP::Properties::UpdateType UpdateType;
 200         DAI_ENUM(RecurseType,REC_FIXED,REC_LOGZ,REC_BDIFF);
 201         DAI_ENUM(ChooseMethodType,CHOOSE_RANDOM,CHOOSE_MAXENT,CHOOSE_BBP,CHOOSE_BP_L1,CHOOSE_BP_CFN);
 202         DAI_ENUM(ClampType,CLAMP_VAR,CLAMP_FACTOR);
 203         size_t verbose;
 204         /// Tolerance to use in BP
 205         double tol;
 206         /// Update style for BP
 207         UpdateType updates;
 208         /// Maximum number of iterations for BP
 209         size_t maxiter;
 210         /// Tolerance to use for controlling recursion depth (\a recurse
 211         /// is REC_LOGZ or REC_BDIFF)
 212         double rec_tol;
 213         /// Maximum number of levels of recursion (\a recurse is REC_FIXED)
 214         size_t max_levels;
 215         /// If choose=CHOOSE_BBP and maximum adjoint is less than this value, don't recurse
 216         double min_max_adj;
 217         /// Heuristic for choosing clamping variable
 218         ChooseMethodType choose;
 219         /// Method for deciding when to stop recursing
 220         RecurseType recursion;
 221         /// Whether to clamp variables or factors
 222         ClampType clamp;
 223         /// Properties to pass to BBP
 224         PropertySet bbp_props;
 225         /// Cost function to use for BBP
 226         bbp_cfn_t bbp_cfn;
 227         /// Random seed
 228         size_t rand_seed;
 229         /// If non-empty, write clamping choices to this file
 230         std::string clamp_outfile;
 231
 232         /// Set members from PropertySet
 233         void set(const PropertySet &opts);
 234         /// Get members into PropertySet
 235         PropertySet get() const;
 236         /// Convert to a string which can be parsed as a PropertySet
 237         std::string toString() const;
 238     } props;
 239 /* }}} END OF GENERATED CODE */
 240
 241     Properties::ChooseMethodType ChooseMethod() { return props.choose; }
 242     Properties::RecurseType Recursion() { return props.recursion; }
 243     Properties::ClampType Clamping() { return props.clamp; }
 244     size_t maxClampLevel() { return props.max_levels; }
 245     double minMaxAdj() { return props.min_max_adj; }
 246     double recTol() { return props.rec_tol; }
 247 };
 248
 249 /// Given a sorted vector of states \a xis and total state count \a n_states, return a vector of states not in \a xis
 250 std::vector<size_t> complement( std::vector<size_t>& xis, size_t n_states );
 251
 252 /// Computes \f$\frac{\exp(a)}{\exp(a)+\exp(b)}\f$
 253 Real unSoftMax( Real a, Real b );
 254
 255 /// Computes log of sum of exponents, i.e., \f$\log\left(\exp(a) + \exp(b)\right)\f$
 256 Real logSumExp( Real a, Real b );
 257
 258 /// Compute sum of pairwise L-infinity distances of the first \a nv factors in each vector
 259 Real dist( const std::vector<Factor>& b1, const std::vector<Factor>& b2, size_t nv );
 260
 261
 262 } // end of namespace dai
 263
 264
 265 #endif