// predefine friends
template<typename T> Real dist( const TFactor<T> & x, const TFactor<T> & y, Prob::DistType dt );
template<typename T> Real KL_dist( const TFactor<T> & p, const TFactor<T> & q );
+template<typename T> Real MutualInfo( const TFactor<T> & p );
+template<typename T> TFactor<T> max( const TFactor<T> & P, const TFactor<T> & Q );
+template<typename T> TFactor<T> min( const TFactor<T> & P, const TFactor<T> & Q );
template<typename T> std::ostream& operator<< (std::ostream& os, const TFactor<T>& P);
TProb<T> _p;
public:
- // Default constructor
- TFactor () : _vs(), _p(1,1.0) {}
-
+ // Construct Factor with empty VarSet but nonempty _p
+ TFactor ( Real p = 1.0 ) : _vs(), _p(1,p) {}
+
// Construct Factor from VarSet
TFactor( const VarSet& ns ) : _vs(ns), _p(_vs.states()) {}
TFactor( const VarSet& ns, const Real* p ) : _vs(ns), _p(_vs.states(),p) {}
// Construct Factor from VarSet and TProb<T>
- TFactor( const VarSet& ns, const TProb<T> p ) : _vs(ns), _p(p) {
+ TFactor( const VarSet& ns, const TProb<T>& p ) : _vs(ns), _p(p) {
#ifdef DAI_DEBUG
assert( _vs.states() == _p.size() );
#endif
return *this;
}
TFactor<T> operator* (const TFactor<T>& Q) const;
+ TFactor<T> operator/ (const TFactor<T>& Q) const;
TFactor<T>& operator*= (const TFactor<T>& Q) { return( *this = (*this * Q) ); }
+ TFactor<T>& operator/= (const TFactor<T>& Q) { return( *this = (*this / Q) ); }
TFactor<T> operator+ (const TFactor<T>& Q) const {
#ifdef DAI_DEBUG
assert( Q._vs == _vs );
_p -= Q._p;
return *this;
}
+ TFactor<T>& operator+= (T q) {
+ _p += q;
+ return *this;
+ }
+ TFactor<T>& operator-= (T q) {
+ _p -= q;
+ return *this;
+ }
+ TFactor<T> operator+ (T q) const {
+ TFactor<T> result(*this);
+ result._p += q;
+ return result;
+ }
+ TFactor<T> operator- (T q) const {
+ TFactor<T> result(*this);
+ result._p -= q;
+ return result;
+ }
TFactor<T> operator^ (Real a) const { TFactor<T> x; x._vs = _vs; x._p = _p^a; return x; }
TFactor<T>& operator^= (Real a) { _p ^= a; return *this; }
_p.makeZero( epsilon );
return *this;
}
+
+ TFactor<T>& makePositive( Real epsilon ) {
+ _p.makePositive( epsilon );
+ return *this;
+ }
TFactor<T> inverse() const {
TFactor<T> inv;
return e;
}
+ TFactor<T> abs() const {
+ TFactor<T> e;
+ e._vs = _vs;
+ e._p = _p.abs();
+ return e;
+ }
+
TFactor<T> log() const {
TFactor<T> l;
l._vs = _vs;
return l0;
}
- T normalize( typename Prob::NormType norm ) { return _p.normalize( norm ); }
- TFactor<T> normalized( typename Prob::NormType norm ) const {
+ T normalize( typename Prob::NormType norm = Prob::NORMPROB ) { return _p.normalize( norm ); }
+ TFactor<T> normalized( typename Prob::NormType norm = Prob::NORMPROB ) const {
TFactor<T> result;
result._vs = _vs;
result._p = _p.normalized( norm );
return result;
}
- // returns unnormalized marginal
- TFactor<T> part_sum(const VarSet & ns) const;
- // returns normalized marginal
- TFactor<T> marginal(const VarSet & ns) const { return part_sum(ns).normalized( Prob::NORMPROB ); }
+ // returns unnormalized marginal; ns should be a subset of vars()
+ TFactor<T> partSum(const VarSet & ns) const;
+ // returns (normalized by default) marginal; ns should be a subset of vars()
+ TFactor<T> marginal(const VarSet & ns, bool normed = true) const { if(normed) return partSum(ns).normalized(); else return partSum(ns); }
+ // sums out all variables except those in ns
+ TFactor<T> notSum(const VarSet & ns) const { return partSum(vars() ^ ns); }
+
+ // embeds this factor in larger varset ns
+ TFactor<T> embed(const VarSet & ns) const {
+ VarSet vs = vars();
+ assert( ns >> vs );
+ if( vs == ns )
+ return *this;
+ else
+ return (*this) * Factor(ns / vs, 1.0);
+ }
bool hasNaNs() const { return _p.hasNaNs(); }
bool hasNegatives() const { return _p.hasNegatives(); }
T totalSum() const { return _p.totalSum(); }
T maxAbs() const { return _p.maxAbs(); }
T maxVal() const { return _p.maxVal(); }
+ T minVal() const { return _p.minVal(); }
Real entropy() const { return _p.entropy(); }
T strength( const Var &i, const Var &j ) const;
}
}
friend Real KL_dist <> (const TFactor<T> & p, const TFactor<T> & q);
+ friend Real MutualInfo <> ( const TFactor<T> & P );
template<class U> friend std::ostream& operator<< (std::ostream& os, const TFactor<U>& P);
};
-template<typename T> TFactor<T> TFactor<T>::part_sum(const VarSet & ns) const {
+template<typename T> TFactor<T> TFactor<T>::partSum(const VarSet & ns) const {
#ifdef DAI_DEBUG
assert( ns << _vs );
#endif
}
+template<typename T> TFactor<T> TFactor<T>::operator/ (const TFactor<T>& Q) const {
+ TFactor<T> quot( _vs + Q._vs, 0.0 );
+
+ IndexFor i1(_vs, quot._vs);
+ IndexFor i2(Q._vs, quot._vs);
+
+ for( size_t i = 0; i < quot._p.size(); i++, ++i1, ++i2 )
+ quot._p[i] += _p[i1] / Q._p[i2];
+
+ return quot;
+}
+
+
template<typename T> Real KL_dist(const TFactor<T> & P, const TFactor<T> & Q) {
if( P._vs.empty() || Q._vs.empty() )
return -1;
}
+// calculate mutual information of x_i and x_j where P.vars() = \{x_i,x_j\}
+template<typename T> Real MutualInfo(const TFactor<T> & P) {
+ assert( P._vs.size() == 2 );
+ VarSet::const_iterator it = P._vs.begin();
+ Var i = *it; it++; Var j = *it;
+ TFactor<T> projection = P.marginal(i) * P.marginal(j);
+ return real( KL_dist( P.normalized(), projection ) );
+}
+
+
+template<typename T> TFactor<T> max( const TFactor<T> & P, const TFactor<T> & Q ) {
+ assert( P._vs == Q._vs );
+ return TFactor<T>( P._vs, min( P.p(), Q.p() ) );
+}
+
+template<typename T> TFactor<T> min( const TFactor<T> & P, const TFactor<T> & Q ) {
+ assert( P._vs == Q._vs );
+ return TFactor<T>( P._vs, max( P.p(), Q.p() ) );
+}
+
// calculate N(psi, i, j)
template<typename T> T TFactor<T>::strength( const Var &i, const Var &j ) const {
#ifdef DAI_DEBUG
assert( _vs.contains( j ) );
assert( i != j );
#endif
- VarSet ij = i | j;
+ VarSet ij(i, j);
T max = 0.0;
for( size_t alpha1 = 0; alpha1 < i.states(); alpha1++ )
VarSet vars = psi.vars();
for( size_t iter = 0; iter < 100; iter++ ) {
for( VarSet::const_iterator n = vars.begin(); n != vars.end(); n++ )
- result = result * result.part_sum(*n).inverse();
- result.normalize( Prob::NORMPROB );
+ result = result * result.partSum(*n).inverse();
+ result.normalize();
}
return result;
for( size_t i = _RTree.size(); (i--) != 0; ) {
// Make outer region _RTree[i].n1 consistent with outer region _RTree[i].n2
// IR(i) = seperator OR(_RTree[i].n1) && OR(_RTree[i].n2)
- Factor new_Qb = _Qa[_RTree[i].n2].part_sum( IR( i ) );
+ Factor new_Qb = _Qa[_RTree[i].n2].partSum( IR( i ) );
_logZ += log(new_Qb.normalize( Prob::NORMPROB ));
_Qa[_RTree[i].n1] *= new_Qb.divided_by( _Qb[i] );
_Qb[i] = new_Qb;
foreach( const Neighbor &k, nbOR(i) )
if( k != e )
piet *= message( i, k.iter );
- message( j, _e ) = piet.part_sum( IR(e) );
+ message( j, _e ) = piet.partSum( IR(e) );
_logZ += log( message(j,_e).normalize( Prob::NORMPROB ) );
}
_Qa[T[i].n2] *= piet;
}
- Factor new_Qb = _Qa[T[i].n2].part_sum( IR( b[i] ) );
+ Factor new_Qb = _Qa[T[i].n2].partSum( IR( b[i] ) );
logZ += log(new_Qb.normalize( Prob::NORMPROB ));
_Qa[T[i].n1] *= new_Qb.divided_by( _Qb[b[i]] );
_Qb[b[i]] = new_Qb;
Factor piet( nsrem, 0.0 );
piet[s] = exp(logZ);
- Pns += piet * _Qa[T[0].n1].part_sum( ns / nsrem ); // OPTIMIZE ME
+ Pns += piet * _Qa[T[0].n1].partSum( ns / nsrem ); // OPTIMIZE ME
// Restore clamped beliefs
for( map<size_t,Factor>::const_iterator alpha = _Qa_old.begin(); alpha != _Qa_old.end(); alpha++ )