From dc680d679e5fa129d3df5697227f58e460eca28e Mon Sep 17 00:00:00 2001 From: Joris Mooij Date: Fri, 19 Sep 2008 13:23:20 +0200 Subject: [PATCH] Replaced complex numbers by real numbers --- include/dai/bp.h | 2 +- include/dai/daialg.h | 2 +- include/dai/factor.h | 18 +++++------------- include/dai/hak.h | 2 +- include/dai/jtree.h | 2 +- include/dai/lc.h | 2 +- include/dai/mf.h | 2 +- include/dai/mr.h | 2 +- include/dai/prob.h | 32 ++++++++++---------------------- include/dai/treeep.h | 2 +- matlab/dai.cpp | 2 +- src/bp.cpp | 6 +++--- src/daialg.cpp | 22 +++++++--------------- src/hak.cpp | 8 ++++---- src/jtree.cpp | 8 ++++---- src/mf.cpp | 6 +++--- src/treeep.cpp | 8 ++++---- tests/test.cpp | 2 +- 18 files changed, 50 insertions(+), 78 deletions(-) diff --git a/include/dai/bp.h b/include/dai/bp.h index 9bb1b51..c972949 100644 --- a/include/dai/bp.h +++ b/include/dai/bp.h @@ -91,7 +91,7 @@ class BP : public DAIAlgFG { Factor belief (const Var &n) const; Factor belief (const VarSet &n) const; std::vector beliefs() const; - Complex logZ() const; + Real logZ() const; void init( const VarSet &ns ); void undoProbs( const VarSet &ns ) { FactorGraph::undoProbs(ns); init(ns); } diff --git a/include/dai/daialg.h b/include/dai/daialg.h index 32c1ece..049d8a4 100644 --- a/include/dai/daialg.h +++ b/include/dai/daialg.h @@ -132,7 +132,7 @@ class InfAlg { virtual std::vector beliefs() const = 0; /// Get log partition sum - virtual Complex logZ() const = 0; + virtual Real logZ() const = 0; /// Clear messages and beliefs virtual void init() = 0; diff --git a/include/dai/factor.h b/include/dai/factor.h index 41b5eaf..6f2e5b0 100644 --- a/include/dai/factor.h +++ b/include/dai/factor.h @@ -36,16 +36,15 @@ namespace dai { template class TFactor; typedef TFactor Factor; -typedef TFactor CFactor; // predefine friends template Real dist( const TFactor & x, const TFactor & y, Prob::DistType dt ); -template Complex KL_dist( const TFactor & p, const TFactor & q ); +template Real KL_dist( const TFactor & p, const TFactor & q ); template std::ostream& operator<< (std::ostream& os, const TFactor& P); -// T should be castable from and to double and to complex +// T should be castable from and to double template class TFactor { protected: VarSet _vs; @@ -206,13 +205,6 @@ template class TFactor { return l0; } - CFactor clog0() const { - CFactor l0; - l0._vs = _vs; - l0._p = _p.clog0(); - return l0; - } - T normalize( typename Prob::NormType norm ) { return _p.normalize( norm ); } TFactor normalized( typename Prob::NormType norm ) const { TFactor result; @@ -247,7 +239,7 @@ template class TFactor { T totalSum() const { return _p.totalSum(); } T maxAbs() const { return _p.maxAbs(); } T maxVal() const { return _p.maxVal(); } - Complex entropy() const { return _p.entropy(); } + Real entropy() const { return _p.entropy(); } T strength( const Var &i, const Var &j ) const; friend Real dist( const TFactor & x, const TFactor & y, Prob::DistType dt ) { @@ -260,7 +252,7 @@ template class TFactor { return dist( x._p, y._p, dt ); } } - friend Complex KL_dist <> (const TFactor & p, const TFactor & q); + friend Real KL_dist <> (const TFactor & p, const TFactor & q); template friend std::ostream& operator<< (std::ostream& os, const TFactor& P); }; @@ -302,7 +294,7 @@ template TFactor TFactor::operator* (const TFactor& Q) cons } -template Complex KL_dist(const TFactor & P, const TFactor & Q) { +template Real KL_dist(const TFactor & P, const TFactor & Q) { if( P._vs.empty() || Q._vs.empty() ) return -1; else { diff --git a/include/dai/hak.h b/include/dai/hak.h index f1adfc4..9420298 100644 --- a/include/dai/hak.h +++ b/include/dai/hak.h @@ -88,7 +88,7 @@ class HAK : public DAIAlgRG { Factor belief( const Var &n ) const; Factor belief( const VarSet &ns ) const; std::vector beliefs() const; - Complex logZ () const; + Real logZ () const; void init( const VarSet &ns ); void undoProbs( const VarSet &ns ) { RegionGraph::undoProbs( ns ); init( ns ); } diff --git a/include/dai/jtree.h b/include/dai/jtree.h index e6bf12f..3244dd9 100644 --- a/include/dai/jtree.h +++ b/include/dai/jtree.h @@ -79,7 +79,7 @@ class JTree : public DAIAlgRG { Factor belief( const Var &n ) const; Factor belief( const VarSet &ns ) const; std::vector beliefs() const; - Complex logZ() const; + Real logZ() const; void init( const VarSet &/*ns*/ ) {} void undoProbs( const VarSet &ns ) { RegionGraph::undoProbs( ns ); init( ns ); } diff --git a/include/dai/lc.h b/include/dai/lc.h index fbfaf9a..4af42dc 100644 --- a/include/dai/lc.h +++ b/include/dai/lc.h @@ -87,7 +87,7 @@ class LC : public DAIAlgFG { return Factor(); } std::vector beliefs() const { return _beliefs; } - Complex logZ() const { + Real logZ() const { DAI_THROW(NOT_IMPLEMENTED); return 0.0; } diff --git a/include/dai/mf.h b/include/dai/mf.h index 7f35546..56e0272 100644 --- a/include/dai/mf.h +++ b/include/dai/mf.h @@ -64,7 +64,7 @@ class MF : public DAIAlgFG { Factor belief (const Var &n) const; Factor belief (const VarSet &ns) const; std::vector beliefs() const; - Complex logZ() const; + Real logZ() const; void init( const VarSet &ns ); void undoProbs( const VarSet &ns ) { FactorGraph::undoProbs(ns); init(ns); } diff --git a/include/dai/mr.h b/include/dai/mr.h index ef7b6b0..2dfb2b9 100644 --- a/include/dai/mr.h +++ b/include/dai/mr.h @@ -77,7 +77,7 @@ class MR : public DAIAlgFG { return Factor(); } std::vector beliefs() const; - Complex logZ() const { + Real logZ() const { DAI_THROW(NOT_IMPLEMENTED); return 0.0; } diff --git a/include/dai/prob.h b/include/dai/prob.h index d11cda1..6037339 100644 --- a/include/dai/prob.h +++ b/include/dai/prob.h @@ -23,7 +23,6 @@ #define __defined_libdai_prob_h -#include #include #include #include @@ -38,15 +37,13 @@ namespace dai { typedef double Real; -typedef std::complex Complex; template class TProb; typedef TProb Prob; -typedef TProb CProb; /// TProb implements a probability vector of type T. -/// T should be castable from and to double and to complex. +/// T should be castable from and to double. template class TProb { protected: /// The entries @@ -54,7 +51,7 @@ template class TProb { private: /// Calculate x times log(x), or 0 if x == 0 - Complex xlogx( Real x ) const { return( x == 0.0 ? 0.0 : Complex(x) * std::log(Complex(x))); } + Real xlogx( Real x ) const { return( x == 0.0 ? 0.0 : x * std::log(x)); } public: /// NORMPROB means that the sum of all entries should be 1 @@ -331,15 +328,6 @@ template class TProb { return l0; } - /// Pointwise (complex) log (or 0 if == 0) -/* CProb clog0() const { - CProb l0; - l0._p.reserve( size() ); - for( size_t i = 0; i < size(); i++ ) - l0._p.push_back( (_p[i] == 0.0) ? 0.0 : std::log( Complex( _p[i] ) ) ); - return l0; - }*/ - /// Return distance of p and q friend Real dist( const TProb & p, const TProb & q, DistType dt ) { #ifdef DAI_DEBUG @@ -369,16 +357,16 @@ template class TProb { return result; } - /// Return (complex) Kullback-Leibler distance with q - friend Complex KL_dist( const TProb & p, const TProb & q ) { + /// Return Kullback-Leibler distance with q + friend Real KL_dist( const TProb & p, const TProb & q ) { #ifdef DAI_DEBUG assert( p.size() == q.size() ); #endif - Complex result = 0.0; + Real result = 0.0; for( size_t i = 0; i < p.size(); i++ ) { if( (Real) p[i] != 0.0 ) { - Complex p_i = p[i]; - Complex q_i = q[i]; + Real p_i = p[i]; + Real q_i = q[i]; result += p_i * (std::log(p_i) - std::log(q_i)); } } @@ -444,9 +432,9 @@ template class TProb { return (std::find_if( _p.begin(), _p.end(), std::bind2nd( std::less_equal(), 0.0 ) ) != _p.end()); } - /// Returns (complex) entropy - Complex entropy() const { - Complex S = 0.0; + /// Returns entropy + Real entropy() const { + Real S = 0.0; for( size_t i = 0; i < size(); i++ ) S -= xlogx(_p[i]); return S; diff --git a/include/dai/treeep.h b/include/dai/treeep.h index 4d3e417..edef1ce 100644 --- a/include/dai/treeep.h +++ b/include/dai/treeep.h @@ -111,7 +111,7 @@ class TreeEP : public JTree { std::string identify() const; void init(); double run(); - Complex logZ() const; + Real logZ() const; bool offtree( size_t I ) const { return (fac2OR[I] == -1U); } diff --git a/matlab/dai.cpp b/matlab/dai.cpp index c4e766b..7ad11b8 100644 --- a/matlab/dai.cpp +++ b/matlab/dai.cpp @@ -103,7 +103,7 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[] ) // Save logZ - double logZ = real( obj->logZ() ); + double logZ = obj->logZ(); // Save maxdiff double maxdiff = obj->MaxDiff(); diff --git a/src/bp.cpp b/src/bp.cpp index 79d751f..580a4c6 100644 --- a/src/bp.cpp +++ b/src/bp.cpp @@ -390,10 +390,10 @@ Factor BP::beliefF (size_t I) const { } -Complex BP::logZ() const { - Complex sum = 0.0; +Real BP::logZ() const { + Real sum = 0.0; for(size_t i = 0; i < nrVars(); ++i ) - sum += Complex(1.0 - nbV(i).size()) * beliefV(i).entropy(); + sum += (1.0 - nbV(i).size()) * beliefV(i).entropy(); for( size_t I = 0; I < nrFactors(); ++I ) sum -= KL_dist( beliefF(I), factor(I) ); return sum; diff --git a/src/daialg.cpp b/src/daialg.cpp index a688f36..062de76 100644 --- a/src/daialg.cpp +++ b/src/daialg.cpp @@ -38,7 +38,7 @@ Factor calcMarginal( const InfAlg & obj, const VarSet & ns, bool reInit ) { if( !reInit ) clamped->init(); - Complex logZ0; + Real logZ0 = 0.0; for( State s(ns); s.valid(); s++ ) { // save unclamped factors connected to ns clamped->saveProbs( ns ); @@ -54,18 +54,16 @@ Factor calcMarginal( const InfAlg & obj, const VarSet & ns, bool reInit ) { clamped->init(); clamped->run(); - Complex Z; + Real Z; if( s == 0 ) { logZ0 = clamped->logZ(); Z = 1.0; } else { // subtract logZ0 to avoid very large numbers Z = exp(clamped->logZ() - logZ0); - if( fabs(imag(Z)) > 1e-5 ) - cout << "Marginal:: WARNING: complex Z (" << Z << ")" << endl; } - Pns[s] = real(Z); + Pns[s] = Z; // restore clamped factors clamped->undoProbs( ns ); @@ -98,7 +96,7 @@ vector calcPairBeliefs( const InfAlg & obj, const VarSet& ns, bool reIni if( !reInit ) clamped->init(); - Complex logZ0; + Real logZ0 = 0.0; for( size_t j = 0; j < N; j++ ) { // clamp Var j to its possible values for( size_t j_val = 0; j_val < vns[j].states(); j_val++ ) { @@ -120,10 +118,7 @@ vector calcPairBeliefs( const InfAlg & obj, const VarSet& ns, bool reIni logZ0 = clamped->logZ(); } else { // subtract logZ0 to avoid very large numbers - Complex Z = exp(clamped->logZ() - logZ0); - if( fabs(imag(Z)) > 1e-5 ) - cout << "calcPairBelief:: Warning: complex Z: " << Z << endl; - Z_xj = real(Z); + Z_xj = exp(clamped->logZ() - logZ0); } for( size_t k = 0; k < N; k++ ) @@ -177,7 +172,7 @@ vector calcPairBeliefsNew( const InfAlg & obj, const VarSet& ns, bool re if( !reInit ) clamped->init(); - Complex logZ0; + Real logZ0 = 0.0; VarSet::const_iterator nj = ns.begin(); for( long j = 0; j < (long)ns.size() - 1; j++, nj++ ) { size_t k = 0; @@ -201,10 +196,7 @@ vector calcPairBeliefsNew( const InfAlg & obj, const VarSet& ns, bool re logZ0 = clamped->logZ(); } else { // subtract logZ0 to avoid very large numbers - Complex Z = exp(clamped->logZ() - logZ0); - if( fabs(imag(Z)) > 1e-5 ) - cout << "calcPairBelief:: Warning: complex Z: " << Z << endl; - Z_xj = real(Z); + Z_xj = exp(clamped->logZ() - logZ0); } // we assume that j.label() < k.label() diff --git a/src/hak.cpp b/src/hak.cpp index 5a5a19b..5570dee 100644 --- a/src/hak.cpp +++ b/src/hak.cpp @@ -415,12 +415,12 @@ vector HAK::beliefs() const { } -Complex HAK::logZ() const { - Complex sum = 0.0; +Real HAK::logZ() const { + Real sum = 0.0; for( size_t beta = 0; beta < nrIRs(); beta++ ) - sum += Complex(IR(beta).c()) * Qb(beta).entropy(); + sum += IR(beta).c() * Qb(beta).entropy(); for( size_t alpha = 0; alpha < nrORs(); alpha++ ) { - sum += Complex(OR(alpha).c()) * Qa(alpha).entropy(); + sum += OR(alpha).c() * Qa(alpha).entropy(); sum += (OR(alpha).log0() * Qa(alpha)).totalSum(); } return sum; diff --git a/src/jtree.cpp b/src/jtree.cpp index c7652e1..ce619f7 100644 --- a/src/jtree.cpp +++ b/src/jtree.cpp @@ -292,12 +292,12 @@ double JTree::run() { } -Complex JTree::logZ() const { - Complex sum = 0.0; +Real JTree::logZ() const { + Real sum = 0.0; for( size_t beta = 0; beta < nrIRs(); beta++ ) - sum += Complex(IR(beta).c()) * _Qb[beta].entropy(); + sum += IR(beta).c() * _Qb[beta].entropy(); for( size_t alpha = 0; alpha < nrORs(); alpha++ ) { - sum += Complex(OR(alpha).c()) * _Qa[alpha].entropy(); + sum += OR(alpha).c() * _Qa[alpha].entropy(); sum += (OR(alpha).log0() * _Qa[alpha]).totalSum(); } return sum; diff --git a/src/mf.cpp b/src/mf.cpp index 5587a3e..921692b 100644 --- a/src/mf.cpp +++ b/src/mf.cpp @@ -168,8 +168,8 @@ vector MF::beliefs() const { } -Complex MF::logZ() const { - Complex sum = 0.0; +Real MF::logZ() const { + Real sum = 0.0; for(size_t i=0; i < nrVars(); i++ ) sum -= beliefV(i).entropy(); @@ -181,7 +181,7 @@ Complex MF::logZ() const { Factor piet; piet = factor(I).log0(); piet *= henk; - sum -= Complex( piet.totalSum() ); + sum -= piet.totalSum(); } return -sum; diff --git a/src/treeep.cpp b/src/treeep.cpp index 1fee924..2a60738 100644 --- a/src/treeep.cpp +++ b/src/treeep.cpp @@ -210,7 +210,7 @@ TreeEP::TreeEP( const FactorGraph &fg, const Properties &opts ) : JTree(fg, opts if( piet.vars() >> (v_i | *j) ) { piet = piet.marginal( v_i | *j ); Factor pietf = piet.marginal(v_i) * piet.marginal(*j); - wg[UEdge(i,findVar(*j))] = real( KL_dist( piet, pietf ) ); + wg[UEdge(i,findVar(*j))] = KL_dist( piet, pietf ); } else wg[UEdge(i,findVar(*j))] = 0; } @@ -455,14 +455,14 @@ double TreeEP::run() { } -Complex TreeEP::logZ() const { +Real TreeEP::logZ() const { double sum = 0.0; // entropy of the tree for( size_t beta = 0; beta < nrIRs(); beta++ ) - sum -= real(_Qb[beta].entropy()); + sum -= _Qb[beta].entropy(); for( size_t alpha = 0; alpha < nrORs(); alpha++ ) - sum += real(_Qa[alpha].entropy()); + sum += _Qa[alpha].entropy(); // energy of the on-tree factors for( size_t alpha = 0; alpha < nrORs(); alpha++ ) diff --git a/tests/test.cpp b/tests/test.cpp index 54294ed..fe42f43 100644 --- a/tests/test.cpp +++ b/tests/test.cpp @@ -107,7 +107,7 @@ class TestAI { obj->run(); time += toc() - tic; try { - logZ = real(obj->logZ()); + logZ = obj->logZ(); has_logZ = true; } catch( Exception &e ) { has_logZ = false; -- 2.20.1