00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00015
00016
00017 #ifndef __defined_libdai_factor_h
00018 #define __defined_libdai_factor_h
00019
00020
00021 #include <iostream>
00022 #include <functional>
00023 #include <cmath>
00024 #include <dai/prob.h>
00025 #include <dai/varset.h>
00026 #include <dai/index.h>
00027 #include <dai/util.h>
00028
00029
00030 namespace dai {
00031
00032
00034
00058 template <typename T> class TFactor {
00059 private:
00061 VarSet _vs;
00063 TProb<T> _p;
00064
00065 public:
00067
00068
00069 TFactor ( T p = 1 ) : _vs(), _p(1,p) {}
00070
00072 TFactor( const Var &v ) : _vs(v), _p(v.states()) {}
00073
00075 TFactor( const VarSet& vars ) : _vs(vars), _p(_vs.nrStates()) {}
00076
00078 TFactor( const VarSet& vars, T p ) : _vs(vars), _p(_vs.nrStates(),p) {}
00079
00081
00085 template<typename S>
00086 TFactor( const VarSet& vars, const std::vector<S> &x ) : _vs(vars), _p(x.begin(), x.begin() + _vs.nrStates(), _vs.nrStates()) {
00087 DAI_ASSERT( x.size() == vars.nrStates() );
00088 }
00089
00091
00094 TFactor( const VarSet& vars, const T* p ) : _vs(vars), _p(p, p + _vs.nrStates(), _vs.nrStates()) {}
00095
00097 TFactor( const VarSet& vars, const TProb<T> &p ) : _vs(vars), _p(p) {
00098 DAI_DEBASSERT( _vs.nrStates() == _p.size() );
00099 }
00100
00102 TFactor( const std::vector<Var> &vars, const std::vector<T> &p ) : _vs(vars.begin(), vars.end(), vars.size()), _p(p.size()) {
00103 Permute permindex(vars);
00104 for( size_t li = 0; li < p.size(); ++li )
00105 _p[permindex.convertLinearIndex(li)] = p[li];
00106 }
00108
00110
00111
00112 const TProb<T>& p() const { return _p; }
00113
00115 TProb<T>& p() { return _p; }
00116
00118 T operator[] (size_t i) const { return _p[i]; }
00119
00121 T& operator[] (size_t i) { return _p[i]; }
00122
00124 const VarSet& vars() const { return _vs; }
00125
00127 VarSet& vars() { return _vs; }
00128
00130
00132 size_t states() const { return _p.size(); }
00133
00135 T entropy() const { return _p.entropy(); }
00136
00138 T max() const { return _p.max(); }
00139
00141 T min() const { return _p.min(); }
00142
00144 T sum() const { return _p.sum(); }
00145
00147 T maxAbs() const { return _p.maxAbs(); }
00148
00150 bool hasNaNs() const { return _p.hasNaNs(); }
00151
00153 bool hasNegatives() const { return _p.hasNegatives(); }
00154
00156 T strength( const Var &i, const Var &j ) const;
00158
00160
00161
00162 TFactor<T> abs() const {
00163
00164
00165
00166 TFactor<T> x;
00167 x._vs = _vs;
00168 x._p = _p.abs();
00169 return x;
00170 }
00171
00173 TFactor<T> exp() const {
00174 TFactor<T> x;
00175 x._vs = _vs;
00176 x._p = _p.exp();
00177 return x;
00178 }
00179
00181
00183 TFactor<T> log(bool zero=false) const {
00184 TFactor<T> x;
00185 x._vs = _vs;
00186 x._p = _p.log(zero);
00187 return x;
00188 }
00189
00191
00193 TFactor<T> inverse(bool zero=true) const {
00194 TFactor<T> x;
00195 x._vs = _vs;
00196 x._p = _p.inverse(zero);
00197 return x;
00198 }
00199
00201
00203 TFactor<T> normalized( typename TProb<T>::NormType norm=TProb<T>::NORMPROB ) const {
00204 TFactor<T> x;
00205 x._vs = _vs;
00206 x._p = _p.normalized( norm );
00207 return x;
00208 }
00210
00212
00213
00214 TFactor<T> & randomize () { _p.randomize(); return *this; }
00215
00217 TFactor<T>& setUniform () { _p.setUniform(); return *this; }
00218
00220
00222 T normalize( typename TProb<T>::NormType norm=TProb<T>::NORMPROB ) { return _p.normalize( norm ); }
00224
00226
00227
00228 TFactor<T> & fill (T x) { _p.fill( x ); return *this; }
00229
00231 TFactor<T>& operator+= (T x) { _p += x; return *this; }
00232
00234 TFactor<T>& operator-= (T x) { _p -= x; return *this; }
00235
00237 TFactor<T>& operator*= (T x) { _p *= x; return *this; }
00238
00240 TFactor<T>& operator/= (T x) { _p /= x; return *this; }
00241
00243 TFactor<T>& operator^= (T x) { _p ^= x; return *this; }
00245
00247
00248
00249 TFactor<T> operator+ (T x) const {
00250
00251
00252
00253
00254 TFactor<T> result;
00255 result._vs = _vs;
00256 result._p = p() + x;
00257 return result;
00258 }
00259
00261 TFactor<T> operator- (T x) const {
00262 TFactor<T> result;
00263 result._vs = _vs;
00264 result._p = p() - x;
00265 return result;
00266 }
00267
00269 TFactor<T> operator* (T x) const {
00270 TFactor<T> result;
00271 result._vs = _vs;
00272 result._p = p() * x;
00273 return result;
00274 }
00275
00277 TFactor<T> operator/ (T x) const {
00278 TFactor<T> result;
00279 result._vs = _vs;
00280 result._p = p() / x;
00281 return result;
00282 }
00283
00285 TFactor<T> operator^ (T x) const {
00286 TFactor<T> result;
00287 result._vs = _vs;
00288 result._p = p() ^ x;
00289 return result;
00290 }
00292
00294
00295
00296
00300 template<typename binOp> TFactor<T>& binaryOp( const TFactor<T> &g, binOp op ) {
00301 if( _vs == g._vs )
00302 _p.pwBinaryOp( g._p, op );
00303 else {
00304 TFactor<T> f(*this);
00305 _vs |= g._vs;
00306 size_t N = _vs.nrStates();
00307
00308 IndexFor i_f( f._vs, _vs );
00309 IndexFor i_g( g._vs, _vs );
00310
00311 _p.p().clear();
00312 _p.p().reserve( N );
00313 for( size_t i = 0; i < N; i++, ++i_f, ++i_g )
00314 _p.p().push_back( op( f._p[i_f], g._p[i_g] ) );
00315 }
00316 return *this;
00317 }
00318
00320
00324 TFactor<T>& operator+= (const TFactor<T>& g) { return binaryOp( g, std::plus<T>() ); }
00325
00327
00331 TFactor<T>& operator-= (const TFactor<T>& g) { return binaryOp( g, std::minus<T>() ); }
00332
00334
00338 TFactor<T>& operator*= (const TFactor<T>& g) { return binaryOp( g, std::multiplies<T>() ); }
00339
00341
00345 TFactor<T>& operator/= (const TFactor<T>& g) { return binaryOp( g, fo_divides0<T>() ); }
00347
00349
00350
00351
00355 template<typename binOp> TFactor<T> binaryTr( const TFactor<T> &g, binOp op ) const {
00356
00357
00358 TFactor<T> result;
00359 if( _vs == g._vs ) {
00360 result._vs = _vs;
00361 result._p = _p.pwBinaryTr( g._p, op );
00362 } else {
00363 result._vs = _vs | g._vs;
00364 size_t N = result._vs.nrStates();
00365
00366 IndexFor i_f( _vs, result.vars() );
00367 IndexFor i_g( g._vs, result.vars() );
00368
00369 result._p.p().clear();
00370 result._p.p().reserve( N );
00371 for( size_t i = 0; i < N; i++, ++i_f, ++i_g )
00372 result._p.p().push_back( op( _p[i_f], g._p[i_g] ) );
00373 }
00374 return result;
00375 }
00376
00378
00382 TFactor<T> operator+ (const TFactor<T>& g) const {
00383 return binaryTr(g,std::plus<T>());
00384 }
00385
00387
00391 TFactor<T> operator- (const TFactor<T>& g) const {
00392 return binaryTr(g,std::minus<T>());
00393 }
00394
00396
00400 TFactor<T> operator* (const TFactor<T>& g) const {
00401 return binaryTr(g,std::multiplies<T>());
00402 }
00403
00405
00409 TFactor<T> operator/ (const TFactor<T>& g) const {
00410 return binaryTr(g,fo_divides0<T>());
00411 }
00413
00415
00416
00417
00428 TFactor<T> slice( const VarSet& vars, size_t varsState ) const;
00429
00431
00436 TFactor<T> embed(const VarSet & vars) const {
00437 DAI_ASSERT( vars >> _vs );
00438 if( _vs == vars )
00439 return *this;
00440 else
00441 return (*this) * TFactor<T>(vars / _vs, (T)1);
00442 }
00443
00445 TFactor<T> marginal(const VarSet &vars, bool normed=true) const;
00446
00448 TFactor<T> maxMarginal(const VarSet &vars, bool normed=true) const;
00450 };
00451
00452
00453 template<typename T> TFactor<T> TFactor<T>::slice( const VarSet& vars, size_t varsState ) const {
00454 DAI_ASSERT( vars << _vs );
00455 VarSet varsrem = _vs / vars;
00456 TFactor<T> result( varsrem, T(0) );
00457
00458
00459 IndexFor i_vars (vars, _vs);
00460 IndexFor i_varsrem (varsrem, _vs);
00461 for( size_t i = 0; i < states(); i++, ++i_vars, ++i_varsrem )
00462 if( (size_t)i_vars == varsState )
00463 result._p[i_varsrem] = _p[i];
00464
00465 return result;
00466 }
00467
00468
00469 template<typename T> TFactor<T> TFactor<T>::marginal(const VarSet &vars, bool normed) const {
00470 VarSet res_vars = vars & _vs;
00471
00472 TFactor<T> res( res_vars, 0.0 );
00473
00474 IndexFor i_res( res_vars, _vs );
00475 for( size_t i = 0; i < _p.size(); i++, ++i_res )
00476 res._p[i_res] += _p[i];
00477
00478 if( normed )
00479 res.normalize( TProb<T>::NORMPROB );
00480
00481 return res;
00482 }
00483
00484
00485 template<typename T> TFactor<T> TFactor<T>::maxMarginal(const VarSet &vars, bool normed) const {
00486 VarSet res_vars = vars & _vs;
00487
00488 TFactor<T> res( res_vars, 0.0 );
00489
00490 IndexFor i_res( res_vars, _vs );
00491 for( size_t i = 0; i < _p.size(); i++, ++i_res )
00492 if( _p[i] > res._p[i_res] )
00493 res._p[i_res] = _p[i];
00494
00495 if( normed )
00496 res.normalize( TProb<T>::NORMPROB );
00497
00498 return res;
00499 }
00500
00501
00502 template<typename T> T TFactor<T>::strength( const Var &i, const Var &j ) const {
00503 DAI_DEBASSERT( _vs.contains( i ) );
00504 DAI_DEBASSERT( _vs.contains( j ) );
00505 DAI_DEBASSERT( i != j );
00506 VarSet ij(i, j);
00507
00508 T max = 0.0;
00509 for( size_t alpha1 = 0; alpha1 < i.states(); alpha1++ )
00510 for( size_t alpha2 = 0; alpha2 < i.states(); alpha2++ )
00511 if( alpha2 != alpha1 )
00512 for( size_t beta1 = 0; beta1 < j.states(); beta1++ )
00513 for( size_t beta2 = 0; beta2 < j.states(); beta2++ )
00514 if( beta2 != beta1 ) {
00515 size_t as = 1, bs = 1;
00516 if( i < j )
00517 bs = i.states();
00518 else
00519 as = j.states();
00520 T f1 = slice( ij, alpha1 * as + beta1 * bs ).p().divide( slice( ij, alpha2 * as + beta1 * bs ).p() ).max();
00521 T f2 = slice( ij, alpha2 * as + beta2 * bs ).p().divide( slice( ij, alpha1 * as + beta2 * bs ).p() ).max();
00522 T f = f1 * f2;
00523 if( f > max )
00524 max = f;
00525 }
00526
00527 return std::tanh( 0.25 * std::log( max ) );
00528 }
00529
00530
00532
00534 template<typename T> std::ostream& operator<< (std::ostream& os, const TFactor<T>& f) {
00535 os << "(" << f.vars() << ", (";
00536 for( size_t i = 0; i < f.states(); i++ )
00537 os << (i == 0 ? "" : ", ") << f[i];
00538 os << "))";
00539 return os;
00540 }
00541
00542
00544
00547 template<typename T> T dist( const TFactor<T> &f, const TFactor<T> &g, typename TProb<T>::DistType dt ) {
00548 if( f.vars().empty() || g.vars().empty() )
00549 return -1;
00550 else {
00551 DAI_DEBASSERT( f.vars() == g.vars() );
00552 return dist( f.p(), g.p(), dt );
00553 }
00554 }
00555
00556
00558
00561 template<typename T> TFactor<T> max( const TFactor<T> &f, const TFactor<T> &g ) {
00562 DAI_ASSERT( f._vs == g._vs );
00563 return TFactor<T>( f._vs, max( f.p(), g.p() ) );
00564 }
00565
00566
00568
00571 template<typename T> TFactor<T> min( const TFactor<T> &f, const TFactor<T> &g ) {
00572 DAI_ASSERT( f._vs == g._vs );
00573 return TFactor<T>( f._vs, min( f.p(), g.p() ) );
00574 }
00575
00576
00578
00581 template<typename T> T MutualInfo(const TFactor<T> &f) {
00582 DAI_ASSERT( f.vars().size() == 2 );
00583 VarSet::const_iterator it = f.vars().begin();
00584 Var i = *it; it++; Var j = *it;
00585 TFactor<T> projection = f.marginal(i) * f.marginal(j);
00586 return dist( f.normalized(), projection, TProb<T>::DISTKL );
00587 }
00588
00589
00591 typedef TFactor<Real> Factor;
00592
00593
00595
00598 Factor createFactorIsing( const Var &x, Real h );
00599
00600
00602
00606 Factor createFactorIsing( const Var &x1, const Var &x2, Real J );
00607
00608
00610
00615 Factor createFactorExpGauss( const VarSet &vs, Real beta );
00616
00617
00619
00623 Factor createFactorPotts( const Var &x1, const Var &x2, Real J );
00624
00625
00627
00630 Factor createFactorDelta( const Var &v, size_t state );
00631
00632
00633 }
00634
00635
00636 #endif