libDAI
factor.h
Go to the documentation of this file.
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 
11 
12 
13 #ifndef __defined_libdai_factor_h
14 #define __defined_libdai_factor_h
15 
16 
17 #include <iostream>
18 #include <functional>
19 #include <cmath>
20 #include <dai/prob.h>
21 #include <dai/varset.h>
22 #include <dai/index.h>
23 #include <dai/util.h>
24 
25 
26 namespace dai {
27 
28 
30 
54 template <typename T>
55 class TFactor {
56  private:
61 
62  public:
64 
65 
66  TFactor ( T p = 1 ) : _vs(), _p(1,p) {}
67 
69  TFactor( const Var &v ) : _vs(v), _p(v.states()) {}
70 
72  TFactor( const VarSet& vars ) : _vs(vars), _p() {
74  }
75 
77  TFactor( const VarSet& vars, T p ) : _vs(vars), _p() {
78  _p = TProb<T>( BigInt_size_t( _vs.nrStates() ), p );
79  }
80 
82 
86  template<typename S>
87  TFactor( const VarSet& vars, const std::vector<S> &x ) : _vs(vars), _p() {
88  DAI_ASSERT( (BigInt)x.size() == vars.nrStates() );
89  _p = TProb<T>( x.begin(), x.end(), x.size() );
90  }
91 
93 
96  TFactor( const VarSet& vars, const T* p ) : _vs(vars), _p() {
97  size_t N = BigInt_size_t( _vs.nrStates() );
98  _p = TProb<T>( p, p + N, N );
99  }
100 
102  TFactor( const VarSet& vars, const TProb<T> &p ) : _vs(vars), _p(p) {
103  DAI_ASSERT( _vs.nrStates() == (BigInt)_p.size() );
104  }
105 
107  TFactor( const std::vector<Var> &vars, const std::vector<T> &p ) : _vs(vars.begin(), vars.end(), vars.size()), _p(p.size()) {
108  BigInt nrStates = 1;
109  for( size_t i = 0; i < vars.size(); i++ )
110  nrStates *= vars[i].states();
111  DAI_ASSERT( nrStates == p.size() );
112  Permute permindex(vars);
113  for( size_t li = 0; li < p.size(); ++li )
114  _p.set( permindex.convertLinearIndex(li), p[li] );
115  }
117 
119 
120 
121  void set( size_t i, T val ) { _p.set( i, val ); }
122 
124  T get( size_t i ) const { return _p[i]; }
126 
128 
129 
130  const TProb<T>& p() const { return _p; }
131 
133  TProb<T>& p() { return _p; }
134 
136  T operator[] (size_t i) const { return _p[i]; }
137 
139  const VarSet& vars() const { return _vs; }
140 
142  VarSet& vars() { return _vs; }
143 
145 
147  size_t nrStates() const { return _p.size(); }
148 
150  T entropy() const { return _p.entropy(); }
151 
153  T max() const { return _p.max(); }
154 
156  T min() const { return _p.min(); }
157 
159  T sum() const { return _p.sum(); }
160 
162  T sumAbs() const { return _p.sumAbs(); }
163 
165  T maxAbs() const { return _p.maxAbs(); }
166 
168  bool hasNaNs() const { return _p.hasNaNs(); }
169 
171  bool hasNegatives() const { return _p.hasNegatives(); }
172 
174  T strength( const Var &i, const Var &j ) const;
175 
177  bool operator==( const TFactor<T>& y ) const {
178  return (_vs == y._vs) && (_p == y._p);
179  }
181 
183 
184 
186  // Note: the alternative (shorter) way of implementing this,
187  // return TFactor<T>( _vs, _p.abs() );
188  // is slower because it invokes the copy constructor of TProb<T>
189  TFactor<T> x;
190  x._vs = _vs;
191  x._p = -_p;
192  return x;
193  }
194 
196  TFactor<T> abs() const {
197  TFactor<T> x;
198  x._vs = _vs;
199  x._p = _p.abs();
200  return x;
201  }
202 
204  TFactor<T> exp() const {
205  TFactor<T> x;
206  x._vs = _vs;
207  x._p = _p.exp();
208  return x;
209  }
210 
212 
214  TFactor<T> log(bool zero=false) const {
215  TFactor<T> x;
216  x._vs = _vs;
217  x._p = _p.log(zero);
218  return x;
219  }
220 
222 
224  TFactor<T> inverse(bool zero=true) const {
225  TFactor<T> x;
226  x._vs = _vs;
227  x._p = _p.inverse(zero);
228  return x;
229  }
230 
232 
234  TFactor<T> normalized( ProbNormType norm=NORMPROB ) const {
235  TFactor<T> x;
236  x._vs = _vs;
237  x._p = _p.normalized( norm );
238  return x;
239  }
241 
243 
244 
245  TFactor<T>& randomize() { _p.randomize(); return *this; }
246 
248  TFactor<T>& setUniform() { _p.setUniform(); return *this; }
249 
251  TFactor<T>& takeAbs() { _p.takeAbs(); return *this; }
252 
254  TFactor<T>& takeExp() { _p.takeExp(); return *this; }
255 
257 
259  TFactor<T>& takeLog( bool zero = false ) { _p.takeLog(zero); return *this; }
260 
262 
264  T normalize( ProbNormType norm=NORMPROB ) { return _p.normalize( norm ); }
266 
268 
269 
270  TFactor<T>& fill (T x) { _p.fill( x ); return *this; }
271 
273  TFactor<T>& operator+= (T x) { _p += x; return *this; }
274 
276  TFactor<T>& operator-= (T x) { _p -= x; return *this; }
277 
279  TFactor<T>& operator*= (T x) { _p *= x; return *this; }
280 
282  TFactor<T>& operator/= (T x) { _p /= x; return *this; }
283 
285  TFactor<T>& operator^= (T x) { _p ^= x; return *this; }
287 
289 
290 
291  TFactor<T> operator+ (T x) const {
292  // Note: the alternative (shorter) way of implementing this,
293  // TFactor<T> result(*this);
294  // result._p += x;
295  // is slower because it invokes the copy constructor of TFactor<T>
296  TFactor<T> result;
297  result._vs = _vs;
298  result._p = p() + x;
299  return result;
300  }
301 
303  TFactor<T> operator- (T x) const {
304  TFactor<T> result;
305  result._vs = _vs;
306  result._p = p() - x;
307  return result;
308  }
309 
311  TFactor<T> operator* (T x) const {
312  TFactor<T> result;
313  result._vs = _vs;
314  result._p = p() * x;
315  return result;
316  }
317 
319  TFactor<T> operator/ (T x) const {
320  TFactor<T> result;
321  result._vs = _vs;
322  result._p = p() / x;
323  return result;
324  }
325 
327  TFactor<T> operator^ (T x) const {
328  TFactor<T> result;
329  result._vs = _vs;
330  result._p = p() ^ x;
331  return result;
332  }
334 
336 
337 
338 
342  template<typename binOp> TFactor<T>& binaryOp( const TFactor<T> &g, binOp op ) {
343  if( _vs == g._vs ) // optimize special case
344  _p.pwBinaryOp( g._p, op );
345  else {
346  TFactor<T> f(*this); // make a copy
347  _vs |= g._vs;
348  size_t N = BigInt_size_t( _vs.nrStates() );
349 
350  IndexFor i_f( f._vs, _vs );
351  IndexFor i_g( g._vs, _vs );
352 
353  _p.p().clear();
354  _p.p().reserve( N );
355  for( size_t i = 0; i < N; i++, ++i_f, ++i_g )
356  _p.p().push_back( op( f._p[i_f], g._p[i_g] ) );
357  }
358  return *this;
359  }
360 
362 
366  TFactor<T>& operator+= (const TFactor<T>& g) { return binaryOp( g, std::plus<T>() ); }
367 
369 
373  TFactor<T>& operator-= (const TFactor<T>& g) { return binaryOp( g, std::minus<T>() ); }
374 
376 
380  TFactor<T>& operator*= (const TFactor<T>& g) { return binaryOp( g, std::multiplies<T>() ); }
381 
383 
387  TFactor<T>& operator/= (const TFactor<T>& g) { return binaryOp( g, fo_divides0<T>() ); }
389 
391 
392 
393 
397  template<typename binOp> TFactor<T> binaryTr( const TFactor<T> &g, binOp op ) const {
398  // Note that to prevent a copy to be made, it is crucial
399  // that the result is declared outside the if-else construct.
400  TFactor<T> result;
401  if( _vs == g._vs ) { // optimize special case
402  result._vs = _vs;
403  result._p = _p.pwBinaryTr( g._p, op );
404  } else {
405  result._vs = _vs | g._vs;
406  size_t N = BigInt_size_t( result._vs.nrStates() );
407 
408  IndexFor i_f( _vs, result.vars() );
409  IndexFor i_g( g._vs, result.vars() );
410 
411  result._p.p().clear();
412  result._p.p().reserve( N );
413  for( size_t i = 0; i < N; i++, ++i_f, ++i_g )
414  result._p.p().push_back( op( _p[i_f], g[i_g] ) );
415  }
416  return result;
417  }
418 
420 
424  TFactor<T> operator+ (const TFactor<T>& g) const {
425  return binaryTr(g,std::plus<T>());
426  }
427 
429 
433  TFactor<T> operator- (const TFactor<T>& g) const {
434  return binaryTr(g,std::minus<T>());
435  }
436 
438 
442  TFactor<T> operator* (const TFactor<T>& g) const {
443  return binaryTr(g,std::multiplies<T>());
444  }
445 
447 
451  TFactor<T> operator/ (const TFactor<T>& g) const {
452  return binaryTr(g,fo_divides0<T>());
453  }
455 
457 
458 
459 
470  TFactor<T> slice( const VarSet& vars, size_t varsState ) const;
471 
473 
478  TFactor<T> embed(const VarSet & vars) const {
479  DAI_ASSERT( vars >> _vs );
480  if( _vs == vars )
481  return *this;
482  else
483  return (*this) * TFactor<T>(vars / _vs, (T)1);
484  }
485 
487  TFactor<T> marginal(const VarSet &vars, bool normed=true) const;
488 
490  TFactor<T> maxMarginal(const VarSet &vars, bool normed=true) const;
492 };
493 
494 
495 template<typename T> TFactor<T> TFactor<T>::slice( const VarSet& vars, size_t varsState ) const {
496  DAI_ASSERT( vars << _vs );
497  VarSet varsrem = _vs / vars;
498  TFactor<T> result( varsrem, T(0) );
499 
500  // OPTIMIZE ME
501  IndexFor i_vars (vars, _vs);
502  IndexFor i_varsrem (varsrem, _vs);
503  for( size_t i = 0; i < nrStates(); i++, ++i_vars, ++i_varsrem )
504  if( (size_t)i_vars == varsState )
505  result.set( i_varsrem, _p[i] );
506 
507  return result;
508 }
509 
510 
511 template<typename T> TFactor<T> TFactor<T>::marginal(const VarSet &vars, bool normed) const {
512  VarSet res_vars = vars & _vs;
513 
514  TFactor<T> res( res_vars, 0.0 );
515 
516  IndexFor i_res( res_vars, _vs );
517  for( size_t i = 0; i < _p.size(); i++, ++i_res )
518  res.set( i_res, res[i_res] + _p[i] );
519 
520  if( normed )
521  res.normalize( NORMPROB );
522 
523  return res;
524 }
525 
526 
527 template<typename T> TFactor<T> TFactor<T>::maxMarginal(const VarSet &vars, bool normed) const {
528  VarSet res_vars = vars & _vs;
529 
530  TFactor<T> res( res_vars, 0.0 );
531 
532  IndexFor i_res( res_vars, _vs );
533  for( size_t i = 0; i < _p.size(); i++, ++i_res )
534  if( _p[i] > res._p[i_res] )
535  res.set( i_res, _p[i] );
536 
537  if( normed )
538  res.normalize( NORMPROB );
539 
540  return res;
541 }
542 
543 
544 template<typename T> T TFactor<T>::strength( const Var &i, const Var &j ) const {
545  DAI_DEBASSERT( _vs.contains( i ) );
546  DAI_DEBASSERT( _vs.contains( j ) );
547  DAI_DEBASSERT( i != j );
548  VarSet ij(i, j);
549 
550  T max = 0.0;
551  for( size_t alpha1 = 0; alpha1 < i.states(); alpha1++ )
552  for( size_t alpha2 = 0; alpha2 < i.states(); alpha2++ )
553  if( alpha2 != alpha1 )
554  for( size_t beta1 = 0; beta1 < j.states(); beta1++ )
555  for( size_t beta2 = 0; beta2 < j.states(); beta2++ )
556  if( beta2 != beta1 ) {
557  size_t as = 1, bs = 1;
558  if( i < j )
559  bs = i.states();
560  else
561  as = j.states();
562  T f1 = slice( ij, alpha1 * as + beta1 * bs ).p().divide( slice( ij, alpha2 * as + beta1 * bs ).p() ).max();
563  T f2 = slice( ij, alpha2 * as + beta2 * bs ).p().divide( slice( ij, alpha1 * as + beta2 * bs ).p() ).max();
564  T f = f1 * f2;
565  if( f > max )
566  max = f;
567  }
568 
569  return std::tanh( 0.25 * std::log( max ) );
570 }
571 
572 
574 
576 template<typename T> std::ostream& operator<< (std::ostream& os, const TFactor<T>& f) {
577  os << "(" << f.vars() << ", (";
578  for( size_t i = 0; i < f.nrStates(); i++ )
579  os << (i == 0 ? "" : ", ") << f[i];
580  os << "))";
581  return os;
582 }
583 
584 
586 
589 template<typename T> T dist( const TFactor<T> &f, const TFactor<T> &g, ProbDistType dt ) {
590  if( f.vars().empty() || g.vars().empty() )
591  return -1;
592  else {
593  DAI_DEBASSERT( f.vars() == g.vars() );
594  return dist( f.p(), g.p(), dt );
595  }
596 }
597 
598 
600 
603 template<typename T> TFactor<T> max( const TFactor<T> &f, const TFactor<T> &g ) {
604  DAI_ASSERT( f.vars() == g.vars() );
605  return TFactor<T>( f.vars(), max( f.p(), g.p() ) );
606 }
607 
608 
610 
613 template<typename T> TFactor<T> min( const TFactor<T> &f, const TFactor<T> &g ) {
614  DAI_ASSERT( f.vars() == g.vars() );
615  return TFactor<T>( f.vars(), min( f.p(), g.p() ) );
616 }
617 
618 
620 
623 template<typename T> T MutualInfo(const TFactor<T> &f) {
624  DAI_ASSERT( f.vars().size() == 2 );
625  VarSet::const_iterator it = f.vars().begin();
626  Var i = *it; it++; Var j = *it;
627  TFactor<T> projection = f.marginal(i) * f.marginal(j);
628  return dist( f.normalized(), projection, DISTKL );
629 }
630 
631 
634 
635 
637 
640 Factor createFactorIsing( const Var &x, Real h );
641 
642 
644 
648 Factor createFactorIsing( const Var &x1, const Var &x2, Real J );
649 
650 
652 
657 Factor createFactorExpGauss( const VarSet &vs, Real beta );
658 
659 
661 
665 Factor createFactorPotts( const Var &x1, const Var &x2, Real J );
666 
667 
669 
672 Factor createFactorDelta( const Var &v, size_t state );
673 
674 
676 
679 Factor createFactorDelta( const VarSet& vs, size_t state );
680 
681 
682 } // end of namespace dai
683 
684 
685 #endif