utils/fginfo.cpp

   1 /*  This file is part of libDAI - http://www.libdai.org/
   2  *
   3  *  Copyright (c) 2006-2011, 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 #include <iostream>
  10 #include <cstdlib>
  11 #include <dai/factorgraph.h>
  12 #include <dai/jtree.h>
  13
  14
  15 using namespace std;
  16 using namespace dai;
  17
  18
  19 void findLoopClusters( const FactorGraph & fg, std::set<VarSet> &allcl, VarSet newcl, const Var & root, size_t length, VarSet vars ) {
  20     for( VarSet::const_iterator in = vars.begin(); in != vars.end(); in++ ) {
  21         VarSet ind = fg.delta( *in );
  22         if( (newcl.size()) >= 2 && ind.contains( root ) ) {
  23             allcl.insert( newcl | *in );
  24         }
  25         else if( length > 1 )
  26             findLoopClusters( fg, allcl, newcl | *in, root, length - 1, ind / newcl );
  27     }
  28 }
  29
  30
  31 size_t countLoops( const FactorGraph & fg, size_t loopdepth ) {
  32     set<VarSet> loops;
  33     for( vector<Var>::const_iterator i0 = fg.vars().begin(); i0 != fg.vars().end(); i0++ ) {
  34         VarSet i0d = fg.delta(*i0);
  35         if( loopdepth > 1 )
  36             findLoopClusters( fg, loops, *i0, *i0, loopdepth - 1, i0d );
  37     }
  38     return loops.size();
  39 }
  40
  41
  42 bool hasShortLoops( const std::vector<Factor> &P ) {
  43     bool found = false;
  44     vector<Factor>::const_iterator I, J;
  45     for( I = P.begin(); I != P.end(); I++ ) {
  46         J = I;
  47         J++;
  48         for( ; J != P.end(); J++ )
  49             if( (I->vars() & J->vars()).size() >= 2 ) {
  50                 found = true;
  51                 break;
  52             }
  53         if( found )
  54             break;
  55     }
  56     return found;
  57 }
  58
  59
  60 bool hasNegatives( const std::vector<Factor> &P ) {
  61     bool found = false;
  62     for( size_t I = 0; I < P.size(); I++ )
  63         if( P[I].hasNegatives() ) {
  64             found = true;
  65             break;
  66         }
  67     return found;
  68 }
  69
  70
  71 int main( int argc, char *argv[] ) {
  72     if( argc != 3 ) {
  73         // Display help message if number of command line arguments is incorrect
  74         cout << "This program is part of libDAI - http://www.libdai.org/" << endl << endl;
  75         cout << "Usage: ./fginfo <in.fg> <maxstates>" << endl << endl;
  76         cout << "Reports some detailed information about the factor graph <in.fg>." << endl;
  77         cout << "Also calculates treewidth, with maximum total number of states" << endl;
  78         cout << "given by <maxstates>, where 0 means unlimited." << endl << endl;
  79         return 1;
  80     } else {
  81         // Read factorgraph
  82         FactorGraph fg;
  83         char *infile = argv[1];
  84         size_t maxstates = fromString<size_t>( argv[2] );
  85         fg.ReadFromFile( infile );
  86
  87         // Output various statistics
  88         cout << "Number of variables:   " << fg.nrVars() << endl;
  89         cout << "Number of factors:     " << fg.nrFactors() << endl;
  90         cout << "Connected:             " << fg.isConnected() << endl;
  91         cout << "Tree:                  " << fg.isTree() << endl;
  92         cout << "Has short loops:       " << hasShortLoops(fg.factors()) << endl;
  93         cout << "Has negatives:         " << hasNegatives(fg.factors()) << endl;
  94         cout << "Binary variables?      " << fg.isBinary() << endl;
  95         cout << "Pairwise interactions? " << fg.isPairwise() << endl;
  96
  97         // Calculate treewidth using various heuristics
  98 #ifdef DAI_WITH_JTREE
  99         std::pair<size_t,BigInt> tw;
 100         cout << "Treewidth (MinNeighbors):     ";
 101         try {
 102             tw = boundTreewidth(fg, &eliminationCost_MinNeighbors, maxstates );
 103             cout << tw.first << " (" << tw.second << " states)" << endl;
 104         } catch( Exception &e ) {
 105             if( e.getCode() == Exception::OUT_OF_MEMORY )
 106                 cout << "> " << maxstates << endl;
 107             else
 108                 cout << "an exception occurred" << endl;
 109         }
 110
 111         cout << "Treewidth (MinWeight):        ";
 112         try {
 113             tw = boundTreewidth(fg, &eliminationCost_MinWeight, maxstates );
 114             cout << tw.first << " (" << tw.second << " states)" << endl;
 115         } catch( Exception &e ) {
 116             if( e.getCode() == Exception::OUT_OF_MEMORY )
 117                 cout << "> " << maxstates << endl;
 118             else
 119                 cout << "an exception occurred" << endl;
 120         }
 121
 122         cout << "Treewidth (MinFill):          ";
 123         try {
 124             tw = boundTreewidth(fg, &eliminationCost_MinFill, maxstates );
 125             cout << tw.first << " (" << tw.second << " states)" << endl;
 126         } catch( Exception &e ) {
 127             if( e.getCode() == Exception::OUT_OF_MEMORY )
 128                 cout << "> " << maxstates << endl;
 129             else
 130                 cout << "an exception occurred" << endl;
 131         }
 132
 133         cout << "Treewidth (WeightedMinFill):  ";
 134         try {
 135             tw = boundTreewidth(fg, &eliminationCost_WeightedMinFill, maxstates );
 136             cout << tw.first << " (" << tw.second << " states)" << endl;
 137         } catch( Exception &e ) {
 138             if( e.getCode() == Exception::OUT_OF_MEMORY )
 139                 cout << "> " << maxstates << endl;
 140             else
 141                 cout << "an exception occurred" << endl;
 142         }
 143 #endif
 144
 145         // Calculate total state space
 146         BigInt stsp = 1;
 147         for( size_t i = 0; i < fg.nrVars(); i++ )
 148             stsp *= fg.var(i).states();
 149         cout << "Total state space:   " << stsp << endl;
 150         // Output type of factor graph
 151         cout << "Type: " << (fg.isPairwise() ? "pairwise" : "higher order") << " interactions, " << (fg.isBinary() ? "binary" : "nonbinary") << " variables" << endl;
 152
 153         // Calculate complexity for LCBP
 154         BigInt cavsum_lcbp = 0;
 155         BigInt cavsum_lcbp2 = 0;
 156         BigInt max_Delta_size = 0;
 157         map<size_t,size_t> cavsizes;
 158         for( size_t i = 0; i < fg.nrVars(); i++ ) {
 159             VarSet di = fg.delta(i);
 160             if( cavsizes.count(di.size()) )
 161                 cavsizes[di.size()]++;
 162             else
 163                 cavsizes[di.size()] = 1;
 164             BigInt Ds = fg.Delta(i).nrStates();
 165             if( Ds > max_Delta_size )
 166                 max_Delta_size = Ds;
 167             cavsum_lcbp += di.nrStates();
 168             for( VarSet::const_iterator j = di.begin(); j != di.end(); j++ )
 169                 cavsum_lcbp2 += j->states();
 170         }
 171         cout << "Maximum pancake has " << max_Delta_size << " states" << endl;
 172         cout << "LCBP with full cavities needs " << cavsum_lcbp << " BP runs" << endl;
 173         cout << "LCBP with only pairinteractions needs " << cavsum_lcbp2 << " BP runs" << endl;
 174         cout << "Cavity sizes: ";
 175         for( map<size_t,size_t>::const_iterator it = cavsizes.begin(); it != cavsizes.end(); it++ )
 176             cout << it->first << "(" << it->second << ") ";
 177         cout << endl;
 178
 179         // Calculate girth and length of loops
 180         if( fg.isPairwise() ) {
 181             bool girth_reached = false;
 182             size_t loopdepth;
 183             for( loopdepth = 2; loopdepth <= fg.nrVars() && !girth_reached; loopdepth++ ) {
 184                 size_t nr_loops = countLoops( fg, loopdepth );
 185                 cout << "Loops up to " << loopdepth << " variables: " << nr_loops << endl;
 186                 if( nr_loops > 0 )
 187                     girth_reached = true;
 188             }
 189             if( girth_reached )
 190                 cout << "Girth: " << loopdepth-1 << endl;
 191             else
 192                 cout << "Girth: infinity" << endl;
 193         }
 194
 195         // Output factor state spaces
 196         map<size_t,size_t> facsizes;
 197         for( size_t I = 0; I < fg.nrFactors(); I++ ) {
 198             size_t Isize = fg.factor(I).vars().size();
 199             if( facsizes.count( Isize ) )
 200                 facsizes[Isize]++;
 201             else
 202                 facsizes[Isize] = 1;
 203         }
 204         cout << "Factor sizes: ";
 205         for( map<size_t,size_t>::const_iterator it = facsizes.begin(); it != facsizes.end(); it++ )
 206             cout << it->first << "(" << it->second << ") ";
 207         cout << endl;
 208
 209         return 0;
 210     }
 211 }