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