Improved Gibbs and added FactorGraph::logScore( const std::vector<size_t>& statevec )
[libdai.git] / src / decmap.cpp
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) 2010 Joris Mooij [joris dot mooij at libdai dot org]
8 */
9
10
11 #include <dai/alldai.h>
12
13
14 namespace dai {
15
16
17 using namespace std;
18
19
20 const char *DecMAP::Name = "DECMAP";
21
22
23 void DecMAP::setProperties( const PropertySet &opts ) {
24 DAI_ASSERT( opts.hasKey("ianame") );
25 DAI_ASSERT( opts.hasKey("iaopts") );
26
27 props.ianame = opts.getStringAs<string>("ianame");
28 props.iaopts = opts.getStringAs<PropertySet>("iaopts");
29 if( opts.hasKey("verbose") )
30 props.verbose = opts.getStringAs<size_t>("verbose");
31 else
32 props.verbose = 0;
33 if( opts.hasKey("reinit") )
34 props.reinit = opts.getStringAs<bool>("reinit");
35 else
36 props.reinit = true;
37 }
38
39
40 PropertySet DecMAP::getProperties() const {
41 PropertySet opts;
42 opts.set( "verbose", props.verbose );
43 opts.set( "reinit", props.reinit );
44 opts.set( "ianame", props.ianame );
45 opts.set( "iaopts", props.iaopts );
46 return opts;
47 }
48
49
50 string DecMAP::printProperties() const {
51 stringstream s( stringstream::out );
52 s << "[";
53 s << "verbose=" << props.verbose << ",";
54 s << "reinit=" << props.reinit << ",";
55 s << "ianame=" << props.ianame << ",";
56 s << "iaopts=" << props.iaopts << "]";
57 return s.str();
58 }
59
60
61 DecMAP::DecMAP( const FactorGraph& fg, const PropertySet& opts ) : DAIAlgFG(fg), _state(), _logp(), _maxdiff(), _iters(), props() {
62 setProperties( opts );
63
64 _state = vector<size_t>( nrVars(), 0 );
65 _logp = -INFINITY;
66 }
67
68
69 string DecMAP::identify() const {
70 return string(Name) + printProperties();
71 }
72
73
74 Factor DecMAP::belief( const VarSet& vs ) const {
75 if( vs.size() == 0 )
76 return Factor();
77 else {
78 map<Var, size_t> state;
79 for( VarSet::const_iterator v = vs.begin(); v != vs.end(); v++ )
80 state[*v] = _state[findVar(*v)];
81 return createFactorDelta( vs, calcLinearState( vs, state ) );
82 }
83 }
84
85
86 Factor DecMAP::beliefV( size_t i ) const {
87 return createFactorDelta( var(i), _state[i] );
88 }
89
90
91 vector<Factor> DecMAP::beliefs() const {
92 vector<Factor> result;
93 for( size_t i = 0; i < nrVars(); ++i )
94 result.push_back( beliefV(i) );
95 for( size_t I = 0; I < nrFactors(); ++I )
96 result.push_back( beliefF(I) );
97 return result;
98 }
99
100
101 Real DecMAP::run() {
102 if( props.verbose >= 1 )
103 cerr << "Starting " << identify() << "...";
104 if( props.verbose >= 2 )
105 cerr << endl;
106
107 // the variables which have not been clamped yet
108 SmallSet<size_t> freeVars;
109 for( size_t i = 0; i < nrVars(); i++ )
110 freeVars |= i;
111
112 // prepare the inference algorithm object
113 InfAlg *clamped = newInfAlg( props.ianame, fg(), props.iaopts );
114
115 // decimate until no free variables remain
116 while( freeVars.size() ) {
117 Real md = clamped->run();
118 if( md > _maxdiff )
119 _maxdiff = md;
120 _iters += clamped->Iterations();
121
122 // store the variables that need initialization
123 VarSet varsToInit;
124 SmallSet<size_t> varsToClamp;
125
126 // schedule clamping for the free variables with zero entropy
127 for( SmallSet<size_t>::const_iterator it = freeVars.begin(); it != freeVars.end(); ) {
128 if( clamped->beliefV( *it ).entropy() == 0.0 ) {
129 // this variable should be clamped
130 varsToInit |= var( *it );
131 varsToClamp |= *it;
132 _state[*it] = clamped->beliefV( *it ).p().argmax().first;
133 freeVars.erase( *it );
134 } else
135 it++;
136 }
137
138 // find the free factor with lowest entropy
139 size_t bestI = 0;
140 Real bestEnt = INFINITY;
141 for( size_t I = 0; I < nrFactors(); I++ ) {
142 // check if the factor is still free
143 if( freeVars.intersects( bipGraph().nb2Set(I) ) ) {
144 Real EntI = clamped->beliefF(I).entropy();
145 if( EntI < bestEnt ) {
146 bestI = I;
147 bestEnt = EntI;
148 }
149 }
150 }
151
152 // schedule clamping for the factor with lowest entropy
153 vector<size_t> Istate(1,0);
154 Istate[0] = clamped->beliefF(bestI).p().argmax().first;
155 map<Var, size_t> Istatemap = calcState( factor(bestI).vars(), Istate[0] );
156 foreach( size_t i, bipGraph().nb2Set(bestI) & freeVars ) {
157 varsToInit |= var(i);
158 varsToClamp |= i;
159 _state[i] = Istatemap[var(i)];
160 freeVars.erase(i);
161 }
162
163 // clamp all variables scheduled for clamping
164 foreach( size_t i, varsToClamp )
165 clamped->clamp( i, _state[i], false );
166
167 // initialize clamped for the next run
168 if( props.reinit )
169 clamped->init();
170 else
171 clamped->init( varsToInit );
172 }
173
174 // calculate MAP state
175 map<Var, size_t> state;
176 for( size_t i = 0; i < nrVars(); i++ )
177 state[var(i)] = _state[i];
178 _logp = 0.0;
179 for( size_t I = 0; I < nrFactors(); I++ )
180 _logp += dai::log( factor(I)[calcLinearState( factor(I).vars(), state )] );
181
182 // clean up
183 delete clamped;
184
185 return _maxdiff;
186 }
187
188
189 } // end of namespace dai