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