Some small documentation updates
[libdai.git] / src / hak.cpp
1 /* This file is part of libDAI - http://www.libdai.org/
2 *
3 * libDAI is licensed under the terms of the GNU General Public License version
4 * 2, or (at your option) any later version. libDAI is distributed without any
5 * warranty. See the file COPYING for more details.
6 *
7 * Copyright (C) 2006-2010 Joris Mooij [joris dot mooij at libdai dot org]
8 * Copyright (C) 2006-2007 Radboud University Nijmegen, The Netherlands
9 */
10
11
12 #include <map>
13 #include <dai/hak.h>
14 #include <dai/util.h>
15 #include <dai/exceptions.h>
16
17
18 namespace dai {
19
20
21 using namespace std;
22
23
24 const char *HAK::Name = "HAK";
25
26
27 /// Sets factor entries that lie between 0 and \a epsilon to \a epsilon
28 template <class T>
29 TFactor<T>& makePositive( TFactor<T> &f, T epsilon ) {
30 for( size_t t = 0; t < f.states(); t++ )
31 if( (0 < f[t]) && (f[t] < epsilon) )
32 f[t] = epsilon;
33 return f;
34 }
35
36 /// Sets factor entries that are smaller (in absolute value) than \a epsilon to 0
37 template <class T>
38 TFactor<T>& makeZero( TFactor<T> &f, T epsilon ) {
39 for( size_t t = 0; t < f.states(); t++ )
40 if( f[t] < epsilon && f[t] > -epsilon )
41 f[t] = 0;
42 return f;
43 }
44
45
46 void HAK::setProperties( const PropertySet &opts ) {
47 DAI_ASSERT( opts.hasKey("tol") );
48 DAI_ASSERT( opts.hasKey("doubleloop") );
49 DAI_ASSERT( opts.hasKey("clusters") );
50
51 props.tol = opts.getStringAs<Real>("tol");
52 props.doubleloop = opts.getStringAs<bool>("doubleloop");
53 props.clusters = opts.getStringAs<Properties::ClustersType>("clusters");
54
55 if( opts.hasKey("maxiter") )
56 props.maxiter = opts.getStringAs<size_t>("maxiter");
57 else
58 props.maxiter = 10000;
59 if( opts.hasKey("maxtime") )
60 props.maxtime = opts.getStringAs<Real>("maxtime");
61 else
62 props.maxtime = INFINITY;
63 if( opts.hasKey("verbose") )
64 props.verbose = opts.getStringAs<size_t>("verbose");
65 else
66 props.verbose = 0;
67 if( opts.hasKey("loopdepth") )
68 props.loopdepth = opts.getStringAs<size_t>("loopdepth");
69 else
70 DAI_ASSERT( props.clusters != Properties::ClustersType::LOOP );
71 if( opts.hasKey("damping") )
72 props.damping = opts.getStringAs<Real>("damping");
73 else
74 props.damping = 0.0;
75 if( opts.hasKey("init") )
76 props.init = opts.getStringAs<Properties::InitType>("init");
77 else
78 props.init = Properties::InitType::UNIFORM;
79 }
80
81
82 PropertySet HAK::getProperties() const {
83 PropertySet opts;
84 opts.set( "tol", props.tol );
85 opts.set( "maxiter", props.maxiter );
86 opts.set( "maxtime", props.maxtime );
87 opts.set( "verbose", props.verbose );
88 opts.set( "doubleloop", props.doubleloop );
89 opts.set( "clusters", props.clusters );
90 opts.set( "init", props.init );
91 opts.set( "loopdepth", props.loopdepth );
92 opts.set( "damping", props.damping );
93 return opts;
94 }
95
96
97 string HAK::printProperties() const {
98 stringstream s( stringstream::out );
99 s << "[";
100 s << "tol=" << props.tol << ",";
101 s << "maxiter=" << props.maxiter << ",";
102 s << "maxtime=" << props.maxtime << ",";
103 s << "verbose=" << props.verbose << ",";
104 s << "doubleloop=" << props.doubleloop << ",";
105 s << "clusters=" << props.clusters << ",";
106 s << "init=" << props.init << ",";
107 s << "loopdepth=" << props.loopdepth << ",";
108 s << "damping=" << props.damping << "]";
109 return s.str();
110 }
111
112
113 void HAK::construct() {
114 // Create outer beliefs
115 if( props.verbose >= 3 )
116 cerr << "Constructing outer beliefs" << endl;
117 _Qa.clear();
118 _Qa.reserve(nrORs());
119 for( size_t alpha = 0; alpha < nrORs(); alpha++ )
120 _Qa.push_back( Factor( OR(alpha) ) );
121
122 // Create inner beliefs
123 if( props.verbose >= 3 )
124 cerr << "Constructing inner beliefs" << endl;
125 _Qb.clear();
126 _Qb.reserve(nrIRs());
127 for( size_t beta = 0; beta < nrIRs(); beta++ )
128 _Qb.push_back( Factor( IR(beta) ) );
129
130 // Create messages
131 if( props.verbose >= 3 )
132 cerr << "Constructing messages" << endl;
133 _muab.clear();
134 _muab.reserve( nrORs() );
135 _muba.clear();
136 _muba.reserve( nrORs() );
137 for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
138 _muab.push_back( vector<Factor>() );
139 _muba.push_back( vector<Factor>() );
140 _muab[alpha].reserve( nbOR(alpha).size() );
141 _muba[alpha].reserve( nbOR(alpha).size() );
142 foreach( const Neighbor &beta, nbOR(alpha) ) {
143 _muab[alpha].push_back( Factor( IR(beta) ) );
144 _muba[alpha].push_back( Factor( IR(beta) ) );
145 }
146 }
147 }
148
149
150 HAK::HAK( const RegionGraph &rg, const PropertySet &opts ) : DAIAlgRG(rg), _Qa(), _Qb(), _muab(), _muba(), _maxdiff(0.0), _iters(0U), props() {
151 setProperties( opts );
152
153 construct();
154 }
155
156
157 void HAK::findLoopClusters( const FactorGraph & fg, std::set<VarSet> &allcl, VarSet newcl, const Var & root, size_t length, VarSet vars ) {
158 for( VarSet::const_iterator in = vars.begin(); in != vars.end(); in++ ) {
159 VarSet ind = fg.delta( fg.findVar( *in ) );
160 if( (newcl.size()) >= 2 && ind.contains( root ) )
161 allcl.insert( newcl | *in );
162 else if( length > 1 )
163 findLoopClusters( fg, allcl, newcl | *in, root, length - 1, ind / newcl );
164 }
165 }
166
167
168 HAK::HAK(const FactorGraph & fg, const PropertySet &opts) : DAIAlgRG(), _Qa(), _Qb(), _muab(), _muba(), _maxdiff(0.0), _iters(0U), props() {
169 setProperties( opts );
170
171 if( props.verbose >= 3 )
172 cerr << "Constructing clusters" << endl;
173
174 vector<VarSet> cl;
175 if( props.clusters == Properties::ClustersType::MIN ) {
176 cl = fg.maximalFactorDomains();
177 constructCVM( fg, cl );
178 } else if( props.clusters == Properties::ClustersType::DELTA ) {
179 cl.reserve( fg.nrVars() );
180 for( size_t i = 0; i < fg.nrVars(); i++ )
181 cl.push_back( fg.Delta(i) );
182 constructCVM( fg, cl );
183 } else if( props.clusters == Properties::ClustersType::LOOP ) {
184 cl = fg.maximalFactorDomains();
185 set<VarSet> scl;
186 if( props.verbose >= 2 )
187 cerr << "Searching loops...";
188 for( size_t i0 = 0; i0 < fg.nrVars(); i0++ ) {
189 VarSet i0d = fg.delta(i0);
190 if( props.loopdepth > 1 )
191 findLoopClusters( fg, scl, fg.var(i0), fg.var(i0), props.loopdepth - 1, fg.delta(i0) );
192 }
193 if( props.verbose >= 2 )
194 cerr << "done" << endl;
195 for( set<VarSet>::const_iterator c = scl.begin(); c != scl.end(); c++ )
196 cl.push_back(*c);
197 if( props.verbose >= 3 ) {
198 cerr << Name << " uses the following clusters: " << endl;
199 for( vector<VarSet>::const_iterator cli = cl.begin(); cli != cl.end(); cli++ )
200 cerr << *cli << endl;
201 }
202 constructCVM( fg, cl );
203 } else if( props.clusters == Properties::ClustersType::BETHE ) {
204 // Copy factor graph structure
205 if( props.verbose >= 3 )
206 cerr << "Copying factor graph" << endl;
207 FactorGraph::operator=( fg );
208
209 // Construct inner regions (single variables)
210 if( props.verbose >= 3 )
211 cerr << "Constructing inner regions" << endl;
212 _IRs.reserve( fg.nrVars() );
213 for( size_t i = 0; i < fg.nrVars(); i++ )
214 _IRs.push_back( Region( fg.var(i), 1.0 ) );
215
216 // Construct graph
217 if( props.verbose >= 3 )
218 cerr << "Constructing graph" << endl;
219 _G = BipartiteGraph( 0, nrIRs() );
220
221 // Construct outer regions:
222 // maximal factors become new outer regions
223 // non-maximal factors are assigned an outer region that contains them
224 if( props.verbose >= 3 )
225 cerr << "Construct outer regions" << endl;
226 _fac2OR.reserve( nrFactors() );
227 queue<pair<size_t, size_t> > todo;
228 for( size_t I = 0; I < fg.nrFactors(); I++ ) {
229 size_t J = fg.maximalFactor( I );
230 if( J == I ) {
231 // I is maximal; add it to the outer regions
232 _fac2OR.push_back( nrORs() );
233 // Construct outer region (with counting number 1.0)
234 _ORs.push_back( FRegion( fg.factor(I), 1.0 ) );
235 // Add node and edges to graph
236 SmallSet<size_t> irs = fg.bipGraph().nb2Set( I );
237 _G.addNode1( irs.begin(), irs.end(), irs.size() );
238 } else if( J < I ) {
239 // J is larger and has already been assigned to an outer region
240 // so I should belong to the same outer region as J
241 _fac2OR.push_back( _fac2OR[J] );
242 _ORs[_fac2OR[J]] *= fg.factor(I);
243 } else {
244 // J is larger but has not yet been assigned to an outer region
245 // we handle this case later
246 _fac2OR.push_back( -1 );
247 todo.push( make_pair( I, J ) );
248 }
249 }
250 // finish the construction
251 while( !todo.empty() ) {
252 size_t I = todo.front().first;
253 size_t J = todo.front().second;
254 todo.pop();
255 _fac2OR[I] = _fac2OR[J];
256 _ORs[_fac2OR[J]] *= fg.factor(I);
257 }
258
259 // Calculate inner regions' counting numbers
260 for( size_t beta = 0; beta < nrIRs(); beta++ )
261 _IRs[beta].c() = 1.0 - _G.nb2(beta).size();
262 } else
263 DAI_THROW(UNKNOWN_ENUM_VALUE);
264
265 construct();
266
267 if( props.verbose >= 3 )
268 cerr << Name << " regiongraph: " << *this << endl;
269 }
270
271
272 string HAK::identify() const {
273 return string(Name) + printProperties();
274 }
275
276
277 void HAK::init( const VarSet &ns ) {
278 for( size_t alpha = 0; alpha < nrORs(); alpha++ )
279 if( _Qa[alpha].vars().intersects( ns ) ) {
280 if( props.init == Properties::InitType::UNIFORM )
281 _Qa[alpha].setUniform();
282 else
283 _Qa[alpha].randomize();
284 _Qa[alpha] *= OR(alpha);
285 _Qa[alpha].normalize();
286 }
287
288 for( size_t beta = 0; beta < nrIRs(); beta++ )
289 if( IR(beta).intersects( ns ) ) {
290 if( props.init == Properties::InitType::UNIFORM )
291 _Qb[beta].fill( 1.0 );
292 else
293 _Qb[beta].randomize();
294 foreach( const Neighbor &alpha, nbIR(beta) ) {
295 size_t _beta = alpha.dual;
296 if( props.init == Properties::InitType::UNIFORM ) {
297 muab( alpha, _beta ).fill( 1.0 );
298 muba( alpha, _beta ).fill( 1.0 );
299 } else {
300 muab( alpha, _beta ).randomize();
301 muba( alpha, _beta ).randomize();
302 }
303 }
304 }
305 }
306
307
308 void HAK::init() {
309 for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
310 if( props.init == Properties::InitType::UNIFORM )
311 _Qa[alpha].setUniform();
312 else
313 _Qa[alpha].randomize();
314 _Qa[alpha] *= OR(alpha);
315 _Qa[alpha].normalize();
316 }
317
318 for( size_t beta = 0; beta < nrIRs(); beta++ )
319 if( props.init == Properties::InitType::UNIFORM )
320 _Qb[beta].setUniform();
321 else
322 _Qb[beta].randomize();
323
324 for( size_t alpha = 0; alpha < nrORs(); alpha++ )
325 foreach( const Neighbor &beta, nbOR(alpha) ) {
326 size_t _beta = beta.iter;
327 if( props.init == Properties::InitType::UNIFORM ) {
328 muab( alpha, _beta ).setUniform();
329 muba( alpha, _beta ).setUniform();
330 } else {
331 muab( alpha, _beta ).randomize();
332 muba( alpha, _beta ).randomize();
333 }
334 }
335 }
336
337
338 Real HAK::doGBP() {
339 if( props.verbose >= 1 )
340 cerr << "Starting " << identify() << "...";
341 if( props.verbose >= 3)
342 cerr << endl;
343
344 double tic = toc();
345
346 // Check whether counting numbers won't lead to problems
347 for( size_t beta = 0; beta < nrIRs(); beta++ )
348 DAI_ASSERT( nbIR(beta).size() + IR(beta).c() != 0.0 );
349
350 // Keep old beliefs to check convergence
351 vector<Factor> oldBeliefsV;
352 oldBeliefsV.reserve( nrVars() );
353 for( size_t i = 0; i < nrVars(); i++ )
354 oldBeliefsV.push_back( beliefV(i) );
355 vector<Factor> oldBeliefsF;
356 oldBeliefsF.reserve( nrFactors() );
357 for( size_t I = 0; I < nrFactors(); I++ )
358 oldBeliefsF.push_back( beliefF(I) );
359
360 // do several passes over the network until maximum number of iterations has
361 // been reached or until the maximum belief difference is smaller than tolerance
362 Real maxDiff = INFINITY;
363 for( _iters = 0; _iters < props.maxiter && maxDiff > props.tol; _iters++ ) {
364 for( size_t beta = 0; beta < nrIRs(); beta++ ) {
365 foreach( const Neighbor &alpha, nbIR(beta) ) {
366 size_t _beta = alpha.dual;
367 muab( alpha, _beta ) = _Qa[alpha].marginal(IR(beta)) / muba(alpha,_beta);
368 /* TODO: INVESTIGATE THIS PROBLEM
369 *
370 * In some cases, the muab's can have very large entries because the muba's have very
371 * small entries. This may cause NANs later on (e.g., multiplying large quantities may
372 * result in +inf; normalization then tries to calculate inf / inf which is NAN).
373 * A fix of this problem would consist in normalizing the messages muab.
374 * However, it is not obvious whether this is a real solution, because it has a
375 * negative performance impact and the NAN's seem to be a symptom of a fundamental
376 * numerical unstability.
377 */
378 muab(alpha,_beta).normalize();
379 }
380
381 Factor Qb_new;
382 foreach( const Neighbor &alpha, nbIR(beta) ) {
383 size_t _beta = alpha.dual;
384 Qb_new *= muab(alpha,_beta) ^ (1 / (nbIR(beta).size() + IR(beta).c()));
385 }
386
387 Qb_new.normalize();
388 if( Qb_new.hasNaNs() ) {
389 // TODO: WHAT TO DO IN THIS CASE?
390 cerr << Name << "::doGBP: Qb_new has NaNs!" << endl;
391 return 1.0;
392 }
393 /* TODO: WHAT IS THE PURPOSE OF THE FOLLOWING CODE?
394 *
395 * _Qb[beta] = Qb_new.makeZero(1e-100);
396 */
397
398 if( props.doubleloop || props.damping == 0.0 )
399 _Qb[beta] = Qb_new; // no damping for double loop
400 else
401 _Qb[beta] = (Qb_new^(1.0 - props.damping)) * (_Qb[beta]^props.damping);
402
403 foreach( const Neighbor &alpha, nbIR(beta) ) {
404 size_t _beta = alpha.dual;
405 muba(alpha,_beta) = _Qb[beta] / muab(alpha,_beta);
406
407 /* TODO: INVESTIGATE WHETHER THIS HACK (INVENTED BY KEES) TO PREVENT NANS MAKES SENSE
408 *
409 * muba(beta,*alpha).makePositive(1e-100);
410 *
411 */
412
413 Factor Qa_new = OR(alpha);
414 foreach( const Neighbor &gamma, nbOR(alpha) )
415 Qa_new *= muba(alpha,gamma.iter);
416 Qa_new ^= (1.0 / OR(alpha).c());
417 Qa_new.normalize();
418 if( Qa_new.hasNaNs() ) {
419 cerr << Name << "::doGBP: Qa_new has NaNs!" << endl;
420 return 1.0;
421 }
422 /* TODO: WHAT IS THE PURPOSE OF THE FOLLOWING CODE?
423 *
424 * _Qb[beta] = Qb_new.makeZero(1e-100);
425 */
426
427 if( props.doubleloop || props.damping == 0.0 )
428 _Qa[alpha] = Qa_new; // no damping for double loop
429 else
430 // FIXME: GEOMETRIC DAMPING IS SLOW!
431 _Qa[alpha] = (Qa_new^(1.0 - props.damping)) * (_Qa[alpha]^props.damping);
432 }
433 }
434
435 // Calculate new single variable beliefs and compare with old ones
436 maxDiff = -INFINITY;
437 for( size_t i = 0; i < nrVars(); ++i ) {
438 Factor b = beliefV(i);
439 maxDiff = std::max( maxDiff, dist( b, oldBeliefsV[i], DISTLINF ) );
440 oldBeliefsV[i] = b;
441 }
442 for( size_t I = 0; I < nrFactors(); ++I ) {
443 Factor b = beliefF(I);
444 maxDiff = std::max( maxDiff, dist( b, oldBeliefsF[I], DISTLINF ) );
445 oldBeliefsF[I] = b;
446 }
447
448 if( props.verbose >= 3 )
449 cerr << Name << "::doGBP: maxdiff " << maxDiff << " after " << _iters+1 << " passes" << endl;
450 }
451
452 if( maxDiff > _maxdiff )
453 _maxdiff = maxDiff;
454
455 if( props.verbose >= 1 ) {
456 if( maxDiff > props.tol ) {
457 if( props.verbose == 1 )
458 cerr << endl;
459 cerr << Name << "::doGBP: WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << maxDiff << endl;
460 } else {
461 if( props.verbose >= 2 )
462 cerr << Name << "::doGBP: ";
463 cerr << "converged in " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
464 }
465 }
466
467 return maxDiff;
468 }
469
470
471 Real HAK::doDoubleLoop() {
472 if( props.verbose >= 1 )
473 cerr << "Starting " << identify() << "...";
474 if( props.verbose >= 3)
475 cerr << endl;
476
477 double tic = toc();
478
479 // Save original outer regions
480 vector<FRegion> org_ORs = _ORs;
481
482 // Save original inner counting numbers and set negative counting numbers to zero
483 vector<Real> org_IR_cs( nrIRs(), 0.0 );
484 for( size_t beta = 0; beta < nrIRs(); beta++ ) {
485 org_IR_cs[beta] = IR(beta).c();
486 if( IR(beta).c() < 0.0 )
487 IR(beta).c() = 0.0;
488 }
489
490 // Keep old beliefs to check convergence
491 vector<Factor> oldBeliefsV;
492 oldBeliefsV.reserve( nrVars() );
493 for( size_t i = 0; i < nrVars(); i++ )
494 oldBeliefsV.push_back( beliefV(i) );
495 vector<Factor> oldBeliefsF;
496 oldBeliefsF.reserve( nrFactors() );
497 for( size_t I = 0; I < nrFactors(); I++ )
498 oldBeliefsF.push_back( beliefF(I) );
499
500 size_t outer_maxiter = props.maxiter;
501 Real outer_tol = props.tol;
502 size_t outer_verbose = props.verbose;
503 Real org_maxdiff = _maxdiff;
504
505 // Set parameters for inner loop
506 props.maxiter = 5;
507 props.verbose = outer_verbose ? outer_verbose - 1 : 0;
508
509 size_t outer_iter = 0;
510 size_t total_iter = 0;
511 Real maxDiff = INFINITY;
512 for( outer_iter = 0; outer_iter < outer_maxiter && maxDiff > outer_tol && (toc() - tic) < props.maxtime; outer_iter++ ) {
513 // Calculate new outer regions
514 for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
515 OR(alpha) = org_ORs[alpha];
516 foreach( const Neighbor &beta, nbOR(alpha) )
517 OR(alpha) *= _Qb[beta] ^ ((IR(beta).c() - org_IR_cs[beta]) / nbIR(beta).size());
518 }
519
520 // Inner loop
521 if( isnan( doGBP() ) )
522 return 1.0;
523
524 // Calculate new single variable beliefs and compare with old ones
525 maxDiff = -INFINITY;
526 for( size_t i = 0; i < nrVars(); ++i ) {
527 Factor b = beliefV(i);
528 maxDiff = std::max( maxDiff, dist( b, oldBeliefsV[i], DISTLINF ) );
529 oldBeliefsV[i] = b;
530 }
531 for( size_t I = 0; I < nrFactors(); ++I ) {
532 Factor b = beliefF(I);
533 maxDiff = std::max( maxDiff, dist( b, oldBeliefsF[I], DISTLINF ) );
534 oldBeliefsF[I] = b;
535 }
536
537 total_iter += Iterations();
538
539 if( props.verbose >= 3 )
540 cerr << Name << "::doDoubleLoop: maxdiff " << maxDiff << " after " << total_iter << " passes" << endl;
541 }
542
543 // restore _maxiter, _verbose and _maxdiff
544 props.maxiter = outer_maxiter;
545 props.verbose = outer_verbose;
546 _maxdiff = org_maxdiff;
547
548 _iters = total_iter;
549 if( maxDiff > _maxdiff )
550 _maxdiff = maxDiff;
551
552 // Restore original outer regions
553 _ORs = org_ORs;
554
555 // Restore original inner counting numbers
556 for( size_t beta = 0; beta < nrIRs(); ++beta )
557 IR(beta).c() = org_IR_cs[beta];
558
559 if( props.verbose >= 1 ) {
560 if( maxDiff > props.tol ) {
561 if( props.verbose == 1 )
562 cerr << endl;
563 cerr << Name << "::doDoubleLoop: WARNING: not converged after " << total_iter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << maxDiff << endl;
564 } else {
565 if( props.verbose >= 3 )
566 cerr << Name << "::doDoubleLoop: ";
567 cerr << "converged in " << total_iter << " passes (" << toc() - tic << " seconds)." << endl;
568 }
569 }
570
571 return maxDiff;
572 }
573
574
575 Real HAK::run() {
576 if( props.doubleloop )
577 return doDoubleLoop();
578 else
579 return doGBP();
580 }
581
582
583 Factor HAK::belief( const VarSet &ns ) const {
584 vector<Factor>::const_iterator beta;
585 for( beta = _Qb.begin(); beta != _Qb.end(); beta++ )
586 if( beta->vars() >> ns )
587 break;
588 if( beta != _Qb.end() )
589 return( beta->marginal(ns) );
590 else {
591 vector<Factor>::const_iterator alpha;
592 for( alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
593 if( alpha->vars() >> ns )
594 break;
595 if( alpha == _Qa.end() )
596 DAI_THROW(BELIEF_NOT_AVAILABLE);
597 return( alpha->marginal(ns) );
598 }
599 }
600
601
602 vector<Factor> HAK::beliefs() const {
603 vector<Factor> result;
604 for( size_t beta = 0; beta < nrIRs(); beta++ )
605 result.push_back( Qb(beta) );
606 for( size_t alpha = 0; alpha < nrORs(); alpha++ )
607 result.push_back( Qa(alpha) );
608 return result;
609 }
610
611
612 Real HAK::logZ() const {
613 Real s = 0.0;
614 for( size_t beta = 0; beta < nrIRs(); beta++ )
615 s += IR(beta).c() * Qb(beta).entropy();
616 for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
617 s += OR(alpha).c() * Qa(alpha).entropy();
618 s += (OR(alpha).log(true) * Qa(alpha)).sum();
619 }
620 return s;
621 }
622
623
624 } // end of namespace dai