Several changes by Giuseppe Passino
[libdai.git] / include / dai / factor.h
1 /* Copyright (C) 2006-2008 Joris Mooij [joris dot mooij at tuebingen dot mpg dot de]
2 Radboud University Nijmegen, The Netherlands /
3 Max Planck Institute for Biological Cybernetics, Germany
4
5 Copyright (C) 2002 Martijn Leisink [martijn@mbfys.kun.nl]
6 Radboud University Nijmegen, The Netherlands
7
8 This file is part of libDAI.
9
10 libDAI is free software; you can redistribute it and/or modify
11 it under the terms of the GNU General Public License as published by
12 the Free Software Foundation; either version 2 of the License, or
13 (at your option) any later version.
14
15 libDAI is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 GNU General Public License for more details.
19
20 You should have received a copy of the GNU General Public License
21 along with libDAI; if not, write to the Free Software
22 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
23 */
24
25
26 /// \file
27 /// \brief Defines TFactor<T> and Factor classes
28
29
30 #ifndef __defined_libdai_factor_h
31 #define __defined_libdai_factor_h
32
33
34 #include <iostream>
35 #include <cmath>
36 #include <dai/prob.h>
37 #include <dai/varset.h>
38 #include <dai/index.h>
39
40
41 namespace dai {
42
43
44 // predefine TFactor<T> class
45 template<typename T> class TFactor;
46
47
48 /// Represents a factor with probability entries represented as Real
49 typedef TFactor<Real> Factor;
50
51
52 /// Represents a probability factor.
53 /** A \e factor is a function of the Cartesian product of the state
54 * spaces of some set of variables to the nonnegative real numbers.
55 * More formally, if \f$x_i \in X_i\f$ for all \f$i\f$, then a factor
56 * depending on the variables \f$\{x_i\}\f$ is a function defined
57 * on \f$\prod_i X_i\f$ with values in \f$[0,\infty)\f$.
58 *
59 * A Factor has two components: a VarSet, defining the set of variables
60 * that the factor depends on, and a TProb<T>, containing the values of
61 * the factor for all possible joint states of the variables.
62 *
63 * \tparam T Should be castable from and to double.
64 */
65 template <typename T> class TFactor {
66 private:
67 VarSet _vs;
68 TProb<T> _p;
69
70 public:
71 /// Iterator over factor entries
72 typedef typename TProb<T>::iterator iterator;
73
74 /// Const iterator over factor entries
75 typedef typename TProb<T>::const_iterator const_iterator;
76
77 /// Construct Factor with empty VarSet
78 TFactor ( Real p = 1.0 ) : _vs(), _p(1,p) {}
79
80 /// Construct Factor from VarSet
81 TFactor( const VarSet& ns ) : _vs(ns), _p(_vs.nrStates()) {}
82
83 /// Construct Factor from VarSet and initial value
84 TFactor( const VarSet& ns, Real p ) : _vs(ns), _p(_vs.nrStates(),p) {}
85
86 /// Construct Factor from VarSet and initial array
87 TFactor( const VarSet& ns, const Real *p ) : _vs(ns), _p(_vs.nrStates(),p) {}
88
89 /// Construct Factor from VarSet and TProb<T>
90 TFactor( const VarSet& ns, const TProb<T>& p ) : _vs(ns), _p(p) {
91 #ifdef DAI_DEBUG
92 assert( _vs.nrStates() == _p.size() );
93 #endif
94 }
95
96 /// Construct Factor from Var
97 TFactor( const Var& n ) : _vs(n), _p(n.states()) {}
98
99 /// Copy constructor
100 TFactor( const TFactor<T> &x ) : _vs(x._vs), _p(x._p) {}
101
102 /// Assignment operator
103 TFactor<T> & operator= (const TFactor<T> &x) {
104 if( this != &x ) {
105 _vs = x._vs;
106 _p = x._p;
107 }
108 return *this;
109 }
110
111 /// Returns const reference to probability entries
112 const TProb<T> & p() const { return _p; }
113 /// Returns reference to probability entries
114 TProb<T> & p() { return _p; }
115
116 /// Returns const reference to variables
117 const VarSet & vars() const { return _vs; }
118
119 /// Returns the number of possible joint states of the variables
120 size_t states() const { return _p.size(); }
121
122 /// Returns a copy of the i'th probability value
123 T operator[] (size_t i) const { return _p[i]; }
124
125 /// Returns a reference to the i'th probability value
126 T& operator[] (size_t i) { return _p[i]; }
127
128 /// Returns iterator pointing to first entry
129 iterator begin() { return _p.begin(); }
130 /// Returns const iterator pointing to first entry
131 const_iterator begin() const { return _p.begin(); }
132 /// Returns iterator pointing beyond last entry
133 iterator end() { return _p.end(); }
134 /// Returns const iterator pointing beyond last entry
135 const_iterator end() const { return _p.end(); }
136
137 /// Sets all probability entries to p
138 TFactor<T> & fill (T p) { _p.fill( p ); return(*this); }
139
140 /// Fills all probability entries with random values
141 TFactor<T> & randomize () { _p.randomize(); return(*this); }
142
143 /// Returns product of *this with x
144 TFactor<T> operator* (T x) const {
145 Factor result = *this;
146 result.p() *= x;
147 return result;
148 }
149
150 /// Multiplies each probability entry with x
151 TFactor<T>& operator*= (T x) {
152 _p *= x;
153 return *this;
154 }
155
156 /// Returns quotient of *this with x
157 TFactor<T> operator/ (T x) const {
158 Factor result = *this;
159 result.p() /= x;
160 return result;
161 }
162
163 /// Divides each probability entry by x
164 TFactor<T>& operator/= (T x) {
165 _p /= x;
166 return *this;
167 }
168
169 /// Returns product of *this with another Factor
170 TFactor<T> operator* (const TFactor<T>& Q) const;
171
172 /// Returns quotient of *this with another Factor
173 TFactor<T> operator/ (const TFactor<T>& Q) const;
174
175 /// Multiplies *this with another Factor
176 TFactor<T>& operator*= (const TFactor<T>& Q) { return( *this = (*this * Q) ); }
177
178 /// Divides *this by another Factor
179 TFactor<T>& operator/= (const TFactor<T>& Q) { return( *this = (*this / Q) ); }
180
181 /// Returns sum of *this and another Factor (their vars() should be identical)
182 TFactor<T> operator+ (const TFactor<T>& Q) const {
183 #ifdef DAI_DEBUG
184 assert( Q._vs == _vs );
185 #endif
186 TFactor<T> sum(*this);
187 sum._p += Q._p;
188 return sum;
189 }
190
191 /// Returns difference of *this and another Factor (their vars() should be identical)
192 TFactor<T> operator- (const TFactor<T>& Q) const {
193 #ifdef DAI_DEBUG
194 assert( Q._vs == _vs );
195 #endif
196 TFactor<T> sum(*this);
197 sum._p -= Q._p;
198 return sum;
199 }
200
201 /// Adds another Factor to *this (their vars() should be identical)
202 TFactor<T>& operator+= (const TFactor<T>& Q) {
203 #ifdef DAI_DEBUG
204 assert( Q._vs == _vs );
205 #endif
206 _p += Q._p;
207 return *this;
208 }
209
210 /// Subtracts another Factor from *this (their vars() should be identical)
211 TFactor<T>& operator-= (const TFactor<T>& Q) {
212 #ifdef DAI_DEBUG
213 assert( Q._vs == _vs );
214 #endif
215 _p -= Q._p;
216 return *this;
217 }
218
219 /// Adds scalar to *this
220 TFactor<T>& operator+= (T q) {
221 _p += q;
222 return *this;
223 }
224
225 /// Subtracts scalar from *this
226 TFactor<T>& operator-= (T q) {
227 _p -= q;
228 return *this;
229 }
230
231 /// Returns sum of *this and a scalar
232 TFactor<T> operator+ (T q) const {
233 TFactor<T> result(*this);
234 result._p += q;
235 return result;
236 }
237
238 /// Returns difference of *this with a scalar
239 TFactor<T> operator- (T q) const {
240 TFactor<T> result(*this);
241 result._p -= q;
242 return result;
243 }
244
245 /// Returns *this raised to some power
246 TFactor<T> operator^ (Real a) const { TFactor<T> x; x._vs = _vs; x._p = _p^a; return x; }
247
248 /// Raises *this to some power
249 TFactor<T>& operator^= (Real a) { _p ^= a; return *this; }
250
251 /// Sets all entries that are smaller than epsilon to zero
252 TFactor<T>& makeZero( Real epsilon ) {
253 _p.makeZero( epsilon );
254 return *this;
255 }
256
257 /// Sets all entries that are smaller than epsilon to epsilon
258 TFactor<T>& makePositive( Real epsilon ) {
259 _p.makePositive( epsilon );
260 return *this;
261 }
262
263 /// Returns inverse of *this
264 TFactor<T> inverse() const {
265 TFactor<T> inv;
266 inv._vs = _vs;
267 inv._p = _p.inverse(true); // FIXME
268 return inv;
269 }
270
271 /// Returns *this divided by another Factor
272 TFactor<T> divided_by( const TFactor<T>& denom ) const {
273 #ifdef DAI_DEBUG
274 assert( denom._vs == _vs );
275 #endif
276 TFactor<T> quot(*this);
277 quot._p /= denom._p;
278 return quot;
279 }
280
281 /// Divides *this by another Factor
282 TFactor<T>& divide( const TFactor<T>& denom ) {
283 #ifdef DAI_DEBUG
284 assert( denom._vs == _vs );
285 #endif
286 _p /= denom._p;
287 return *this;
288 }
289
290 /// Returns exp of *this
291 TFactor<T> exp() const {
292 TFactor<T> e;
293 e._vs = _vs;
294 e._p = _p.exp();
295 return e;
296 }
297
298 /// Returns absolute value of *this
299 TFactor<T> abs() const {
300 TFactor<T> e;
301 e._vs = _vs;
302 e._p = _p.abs();
303 return e;
304 }
305
306 /// Returns logarithm of *this
307 TFactor<T> log() const {
308 TFactor<T> l;
309 l._vs = _vs;
310 l._p = _p.log();
311 return l;
312 }
313
314 /// Returns logarithm of *this (defining log(0)=0)
315 TFactor<T> log0() const {
316 TFactor<T> l0;
317 l0._vs = _vs;
318 l0._p = _p.log0();
319 return l0;
320 }
321
322 /// Normalizes *this Factor
323 T normalize( typename Prob::NormType norm = Prob::NORMPROB ) { return _p.normalize( norm ); }
324
325 /// Returns a normalized copy of *this
326 TFactor<T> normalized( typename Prob::NormType norm = Prob::NORMPROB ) const {
327 TFactor<T> result;
328 result._vs = _vs;
329 result._p = _p.normalized( norm );
330 return result;
331 }
332
333 /// Returns a slice of this factor, where the subset ns is in state ns_state
334 Factor slice( const VarSet & ns, size_t ns_state ) const {
335 assert( ns << _vs );
336 VarSet nsrem = _vs / ns;
337 Factor result( nsrem, 0.0 );
338
339 // OPTIMIZE ME
340 IndexFor i_ns (ns, _vs);
341 IndexFor i_nsrem (nsrem, _vs);
342 for( size_t i = 0; i < states(); i++, ++i_ns, ++i_nsrem )
343 if( (size_t)i_ns == ns_state )
344 result._p[i_nsrem] = _p[i];
345
346 return result;
347 }
348
349 /// Returns unnormalized marginal; ns should be a subset of vars()
350 TFactor<T> partSum(const VarSet & ns) const;
351
352 /// Returns (normalized by default) marginal; ns should be a subset of vars()
353 TFactor<T> marginal(const VarSet & ns, bool normed = true) const { if(normed) return partSum(ns).normalized(); else return partSum(ns); }
354
355 /// Sums out all variables except those in ns
356 TFactor<T> notSum(const VarSet & ns) const { return partSum(vars() ^ ns); }
357
358 /// Embeds this factor in a larger VarSet
359 TFactor<T> embed(const VarSet & ns) const {
360 VarSet vs = vars();
361 assert( ns >> vs );
362 if( vs == ns )
363 return *this;
364 else
365 return (*this) * Factor(ns / vs, 1.0);
366 }
367
368 /// Returns true if *this has NANs
369 bool hasNaNs() const { return _p.hasNaNs(); }
370
371 /// Returns true if *this has negative entries
372 bool hasNegatives() const { return _p.hasNegatives(); }
373
374 /// Returns total sum of probability entries
375 T totalSum() const { return _p.totalSum(); }
376
377 /// Returns maximum absolute value of probability entries
378 T maxAbs() const { return _p.maxAbs(); }
379
380 /// Returns maximum value of probability entries
381 T maxVal() const { return _p.maxVal(); }
382
383 /// Returns minimum value of probability entries
384 T minVal() const { return _p.minVal(); }
385
386 /// Returns entropy of *this
387 Real entropy() const { return _p.entropy(); }
388
389 /// Returns strength of *this, between variables i and j, using (52) of [\ref MoK07b]
390 T strength( const Var &i, const Var &j ) const;
391 };
392
393
394 template<typename T> TFactor<T> TFactor<T>::partSum(const VarSet & ns) const {
395 #ifdef DAI_DEBUG
396 assert( ns << _vs );
397 #endif
398
399 TFactor<T> res( ns, 0.0 );
400
401 IndexFor i_res( ns, _vs );
402 for( size_t i = 0; i < _p.size(); i++, ++i_res )
403 res._p[i_res] += _p[i];
404
405 return res;
406 }
407
408
409 template<typename T> TFactor<T> TFactor<T>::operator* (const TFactor<T>& Q) const {
410 TFactor<T> prod( _vs | Q._vs, 0.0 );
411
412 IndexFor i1(_vs, prod._vs);
413 IndexFor i2(Q._vs, prod._vs);
414
415 for( size_t i = 0; i < prod._p.size(); i++, ++i1, ++i2 )
416 prod._p[i] += _p[i1] * Q._p[i2];
417
418 return prod;
419 }
420
421
422 template<typename T> TFactor<T> TFactor<T>::operator/ (const TFactor<T>& Q) const {
423 TFactor<T> quot( _vs + Q._vs, 0.0 );
424
425 IndexFor i1(_vs, quot._vs);
426 IndexFor i2(Q._vs, quot._vs);
427
428 for( size_t i = 0; i < quot._p.size(); i++, ++i1, ++i2 )
429 quot._p[i] += _p[i1] / Q._p[i2];
430
431 return quot;
432 }
433
434
435 template<typename T> T TFactor<T>::strength( const Var &i, const Var &j ) const {
436 #ifdef DAI_DEBUG
437 assert( _vs.contains( i ) );
438 assert( _vs.contains( j ) );
439 assert( i != j );
440 #endif
441 VarSet ij(i, j);
442
443 T max = 0.0;
444 for( size_t alpha1 = 0; alpha1 < i.states(); alpha1++ )
445 for( size_t alpha2 = 0; alpha2 < i.states(); alpha2++ )
446 if( alpha2 != alpha1 )
447 for( size_t beta1 = 0; beta1 < j.states(); beta1++ )
448 for( size_t beta2 = 0; beta2 < j.states(); beta2++ )
449 if( beta2 != beta1 ) {
450 size_t as = 1, bs = 1;
451 if( i < j )
452 bs = i.states();
453 else
454 as = j.states();
455 T f1 = slice( ij, alpha1 * as + beta1 * bs ).p().divide( slice( ij, alpha2 * as + beta1 * bs ).p() ).maxVal();
456 T f2 = slice( ij, alpha2 * as + beta2 * bs ).p().divide( slice( ij, alpha1 * as + beta2 * bs ).p() ).maxVal();
457 T f = f1 * f2;
458 if( f > max )
459 max = f;
460 }
461
462 return std::tanh( 0.25 * std::log( max ) );
463 }
464
465
466 /// Writes a Factor to an output stream
467 template<typename T> std::ostream& operator<< (std::ostream& os, const TFactor<T>& P) {
468 os << "(" << P.vars() << " <";
469 for( size_t i = 0; i < P.states(); i++ )
470 os << P[i] << " ";
471 os << ">)";
472 return os;
473 }
474
475
476 /// Returns distance between two Factors (with identical vars())
477 template<typename T> Real dist( const TFactor<T> & x, const TFactor<T> & y, Prob::DistType dt ) {
478 if( x.vars().empty() || y.vars().empty() )
479 return -1;
480 else {
481 #ifdef DAI_DEBUG
482 assert( x.vars() == y.vars() );
483 #endif
484 return dist( x.p(), y.p(), dt );
485 }
486 }
487
488
489 /// Returns the pointwise maximum of two Factors
490 template<typename T> TFactor<T> max( const TFactor<T> & P, const TFactor<T> & Q ) {
491 assert( P._vs == Q._vs );
492 return TFactor<T>( P._vs, min( P.p(), Q.p() ) );
493 }
494
495
496 /// Returns the pointwise minimum of two Factors
497 template<typename T> TFactor<T> min( const TFactor<T> & P, const TFactor<T> & Q ) {
498 assert( P._vs == Q._vs );
499 return TFactor<T>( P._vs, max( P.p(), Q.p() ) );
500 }
501
502
503 /// Calculates the mutual information between the two variables in P
504 template<typename T> Real MutualInfo(const TFactor<T> & P) {
505 assert( P.vars().size() == 2 );
506 VarSet::const_iterator it = P.vars().begin();
507 Var i = *it; it++; Var j = *it;
508 TFactor<T> projection = P.marginal(i) * P.marginal(j);
509 return real( dist( P.normalized(), projection, Prob::DISTKL ) );
510 }
511
512
513 } // end of namespace dai
514
515
516 #endif