76b4ef81264789efd36fde234d15200ca53b768a
[libdai.git] / src / mr.cpp
1 /* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
2 Radboud University Nijmegen, The Netherlands
3
4 This file is part of libDAI.
5
6 libDAI is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 libDAI is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with libDAI; if not, write to the Free Software
18 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
19 */
20
21
22 #include <cstdio>
23 #include <ctime>
24 #include <cmath>
25 #include <cstdlib>
26 #include <dai/mr.h>
27 #include <dai/bp.h>
28 #include <dai/jtree.h>
29 #include <dai/util.h>
30 #include <dai/diffs.h>
31
32
33 namespace dai {
34
35
36 using namespace std;
37
38
39 const char *MR::Name = "MR";
40
41
42 void MR::setProperties( const PropertySet &opts ) {
43 assert( opts.hasKey("tol") );
44 assert( opts.hasKey("verbose") );
45 assert( opts.hasKey("updates") );
46 assert( opts.hasKey("inits") );
47
48 props.tol = opts.getStringAs<double>("tol");
49 props.verbose = opts.getStringAs<size_t>("verbose");
50 props.updates = opts.getStringAs<Properties::UpdateType>("updates");
51 props.inits = opts.getStringAs<Properties::InitType>("inits");
52 }
53
54
55 PropertySet MR::getProperties() const {
56 PropertySet opts;
57 opts.Set( "tol", props.tol );
58 opts.Set( "verbose", props.verbose );
59 opts.Set( "updates", props.updates );
60 opts.Set( "inits", props.inits );
61 return opts;
62 }
63
64
65 string MR::printProperties() const {
66 stringstream s( stringstream::out );
67 s << "[";
68 s << "tol=" << props.tol << ",";
69 s << "verbose=" << props.verbose << ",";
70 s << "updates=" << props.updates << ",";
71 s << "inits=" << props.inits << "]";
72 return s.str();
73 }
74
75
76 // init N, con, nb, tJ, theta
77 void MR::init(size_t Nin, double *_w, double *_th) {
78 size_t i,j;
79
80 N = Nin;
81
82 con.resize(N);
83 nb.resize(N);
84 tJ.resize(N);
85 for(i=0; i<N; i++ ) {
86 nb[i].resize(kmax);
87 tJ[i].resize(kmax);
88 con[i]=0;
89 for(j=0; j<N; j++ )
90 if( _w[i*N+j] != 0.0 ) {
91 nb[i][con[i]] = j;
92 tJ[i][con[i]] = tanh(_w[i*N+j]);
93 con[i]++;
94 }
95 }
96
97 theta.resize(N);
98 for(i=0; i<N; i++)
99 theta[i] = _th[i];
100 }
101
102
103 // calculate cors
104 double MR::init_cor_resp() {
105 size_t j,k,l, runx,i2;
106 double variab1, variab2;
107 double md, maxdev;
108 double thbJsite[kmax];
109 double xinter;
110 double rinter;
111 double res[kmax];
112 size_t s2;
113 size_t flag;
114 size_t concav;
115 size_t runs = 3000;
116 double eps = 0.2;
117 size_t cavity;
118
119 vector<vector<double> > tJ_org;
120 vector<vector<size_t> > nb_org;
121 vector<size_t> con_org;
122 vector<double> theta_org;
123
124 vector<double> xfield(N*kmax,0.0);
125 vector<double> rfield(N*kmax,0.0);
126 vector<double> Hfield(N,0.0);
127 vector<double> devs(N*kmax,0.0);
128 vector<double> devs2(N*kmax,0.0);
129 vector<double> dev(N,0.0);
130 vector<double> avmag(N,0.0);
131
132 // save original tJ, nb
133 nb_org = nb;
134 tJ_org = tJ;
135 con_org = con;
136 theta_org = theta;
137
138 maxdev = 0.0;
139 for(cavity=0; cavity<N; cavity++){ // for each spin to be removed
140 con = con_org;
141 concav=con[cavity];
142
143 nb = nb_org;
144 tJ = tJ_org;
145
146 // Adapt the graph variables nb[], tJ[] and con[]
147 for(size_t i=0; i<con[cavity]; i++) {
148 size_t ij = nb[cavity][i];
149 flag=0;
150 j=0;
151 do{
152 if(nb[ij][j]==cavity){
153 while(j<(con[ij]-1)){
154 nb[ij][j]=nb[ij][j+1];
155 tJ[ij][j] = tJ[ij][j+1];
156 j++;
157 }
158 flag=1;
159 }
160 j++;
161 } while(flag==0);
162 }
163 for(size_t i=0; i<con[cavity]; i++)
164 con[nb[cavity][i]]--;
165 con[cavity] = 0;
166
167
168 // Do everything starting from the new graph********
169
170 makekindex();
171 theta = theta_org;
172
173 for(size_t i=0; i<kmax*N; i++)
174 xfield[i] = 3.0*(2*rnd_uniform()-1.);
175
176 for(i2=0; i2<concav; i2++){ // Subsequently apply a field to each cavity spin ****************
177
178 s2 = nb[cavity][i2]; // identify the index of the cavity spin
179 for(size_t i=0; i<con[s2]; i++)
180 rfield[kmax*s2+i] = 1.;
181
182 runx=0;
183 do { // From here start the response and belief propagation
184 runx++;
185 md=0.0;
186 for(k=0; k<N; k++){
187 if(k!=cavity) {
188 for(size_t i=0; i<con[k]; i++)
189 thbJsite[i] = tJ[k][i];
190 for(l=0; l<con[k]; l++){
191 xinter = 1.;
192 rinter = 0.;
193 if(k==s2) rinter += 1.;
194 for(j=0; j<con[k]; j++)
195 if(j!=l){
196 variab2 = tanh(xfield[kmax*nb[k][j]+kindex[k][j]]);
197 variab1 = thbJsite[j]*variab2;
198 xinter *= (1.+variab1)/(1.-variab1);
199
200 rinter += thbJsite[j]*rfield[kmax*nb[k][j]+kindex[k][j]]*(1-variab2*variab2)/(1-variab1*variab1);
201 }
202
203 variab1 = 0.5*log(xinter);
204 xinter = variab1 + theta[k];
205 devs[kmax*k+l] = xinter-xfield[kmax*k+l];
206 xfield[kmax*k+l] = xfield[kmax*k+l]+devs[kmax*k+l]*eps;
207 if( fabs(devs[kmax*k+l]) > md )
208 md = fabs(devs[kmax*k+l]);
209
210 devs2[kmax*k+l] = rinter-rfield[kmax*k+l];
211 rfield[kmax*k+l] = rfield[kmax*k+l]+devs2[kmax*k+l]*eps;
212 if( fabs(devs2[kmax*k+l]) > md )
213 md = fabs(devs2[kmax*k+l]);
214 }
215 }
216 }
217 } while((md > props.tol)&&(runx<runs)); // Precision condition reached -> BP and RP finished
218 if(runx==runs)
219 if( props.verbose >= 2 )
220 cout << "init_cor_resp: Convergence not reached (md=" << md << ")..." << endl;
221 if(md > maxdev)
222 maxdev = md;
223
224 // compute the observables (i.e. magnetizations and responses)******
225
226 for(size_t i=0; i<concav; i++){
227 rinter = 0.;
228 xinter = 1.;
229 if(i!=i2)
230 for(j=0; j<con[nb[cavity][i]]; j++){
231 variab2 = tanh(xfield[kmax*nb[nb[cavity][i]][j]+kindex[nb[cavity][i]][j]]);
232 variab1 = tJ[nb[cavity][i]][j]*variab2;
233 rinter += tJ[nb[cavity][i]][j]*rfield[kmax*nb[nb[cavity][i]][j]+kindex[nb[cavity][i]][j]]*(1-variab2*variab2)/(1-variab1*variab1);
234 xinter *= (1.+variab1)/(1.-variab1);
235 }
236 xinter = tanh(0.5*log(xinter)+theta[nb[cavity][i]]);
237 res[i] = rinter*(1-xinter*xinter);
238 }
239
240 // *******************
241
242 for(size_t i=0; i<concav; i++)
243 if(nb[cavity][i]!=s2)
244 // if(i!=i2)
245 cors[cavity][i2][i] = res[i];
246 else
247 cors[cavity][i2][i] = 0;
248 } // close for i2 = 0...concav
249 }
250
251 // restore nb, tJ, con
252 tJ = tJ_org;
253 nb = nb_org;
254 con = con_org;
255 theta = theta_org;
256
257 return maxdev;
258 }
259
260
261 double MR::T(size_t i, sub_nb A) {
262 // i is a variable index
263 // A is a subset of nb[i]
264 //
265 // calculate T{(i)}_A as defined in Rizzo&Montanari e-print (2.17)
266
267 sub_nb _nbi_min_A(con[i]);
268 for( size_t __j = 0; __j < A.size(); __j++ )
269 _nbi_min_A -= A[__j];
270
271 double res = theta[i];
272 for( size_t __j = 0; __j < _nbi_min_A.size(); __j++ ) {
273 size_t _j = _nbi_min_A[__j];
274 res += atanh(tJ[i][_j] * M[i][_j]);
275 }
276 return tanh(res);
277 }
278
279
280 double MR::T(size_t i, size_t _j) {
281 sub_nb j(con[i]);
282 j.clear();
283 j += _j;
284 return T(i,j);
285 }
286
287
288 double MR::Omega(size_t i, size_t _j, size_t _l) {
289 sub_nb jl(con[i]);
290 jl.clear();
291 jl += _j;
292 jl += _l;
293 double Tijl = T(i,jl);
294 return Tijl / (1.0 + tJ[i][_l] * M[i][_l] * Tijl);
295 }
296
297
298 double MR::Gamma(size_t i, size_t _j, size_t _l1, size_t _l2) {
299 sub_nb jll(con[i]);
300 jll.clear();
301 jll += _j;
302 double Tij = T(i,jll);
303 jll += _l1;
304 jll += _l2;
305 double Tijll = T(i,jll);
306
307 return (Tijll - Tij) / (1.0 + tJ[i][_l1] * tJ[i][_l2] * M[i][_l1] * M[i][_l2] + tJ[i][_l1] * M[i][_l1] * Tijll + tJ[i][_l2] * M[i][_l2] * Tijll);
308 }
309
310
311 double MR::Gamma(size_t i, size_t _l1, size_t _l2) {
312 sub_nb ll(con[i]);
313 ll.clear();
314 double Ti = T(i,ll);
315 ll += _l1;
316 ll += _l2;
317 double Till = T(i,ll);
318
319 return (Till - Ti) / (1.0 + tJ[i][_l1] * tJ[i][_l2] * M[i][_l1] * M[i][_l2] + tJ[i][_l1] * M[i][_l1] * Till + tJ[i][_l2] * M[i][_l2] * Till);
320 }
321
322
323 double MR::_tJ(size_t i, sub_nb A) {
324 // i is a variable index
325 // A is a subset of nb[i]
326 //
327 // calculate the product of all tJ[i][_j] for _j in A
328
329 size_t Asize = A.size();
330 switch( Asize ) {
331 case 0: return 1.0;
332 // case 1: return tJ[i][A[0]];
333 // case 2: return tJ[i][A[0]] * tJ[i][A[1]];
334 // case 3: return tJ[i][A[0]] * tJ[i][A[1]] * tJ[i][A[2]];
335 default:
336 size_t __j = Asize - 1;
337 size_t _j = A[__j];
338 sub_nb A_j = A - _j;
339 return tJ[i][_j] * _tJ(i, A_j);
340 }
341
342 }
343
344
345 double MR::appM(size_t i, sub_nb A) {
346 // i is a variable index
347 // A is a subset of nb[i]
348 //
349 // calculate the moment of variables in A from M and cors, neglecting higher order cumulants,
350 // defined as the sum over all partitions of A into subsets of cardinality two at most of the
351 // product of the cumulants (either first order, i.e. M, or second order, i.e. cors) of the
352 // entries of the partitions
353
354 size_t Asize = A.size();
355 switch( Asize ) {
356 case 0: return 1.0;
357 // case 1: return M[i][A[0]];
358 // case 2: return M[i][A[0]] * M[i][A[1]] + cors[i][A[0]][A[1]];
359 // case 3: return M[i][A[0]] * M[i][A[1]] * M[i][A[2]] + M[i][A[0]] * cors[i][A[1]][A[2]] + M[i][A[1]] * cors[i][A[0]][A[2]] + M[i][A[2]] * cors[i][A[0]][A[1]];
360 default:
361 size_t __j = Asize - 1;
362 size_t _j = A[__j];
363 sub_nb A_j = A - _j;
364
365 double result = M[i][_j] * appM(i, A_j);
366 for( size_t __k = 0; __k < __j; __k++ ) {
367 size_t _k = A[__k];
368 result += cors[i][_j][_k] * appM(i,A_j - _k);
369 }
370
371 return result;
372 }
373 }
374
375
376 void MR::sum_subs(size_t j, sub_nb A, double *sum_even, double *sum_odd) {
377 // j is a variable index
378 // A is a subset of nb[j]
379
380 // calculate sum over all even/odd subsets B of A of _tJ(j,B) appM(j,B)
381
382 *sum_even = 0.0;
383 *sum_odd = 0.0;
384
385 sub_nb S(A.size());
386 S.clear();
387 do { // for all subsets of A
388
389 // construct subset B of A corresponding to S
390 sub_nb B = A;
391 B.clear();
392 size_t Ssize = S.size();
393 for( size_t bit = 0; bit < Ssize; bit++ )
394 B += A[S[bit]];
395
396 if( S.size() % 2 )
397 *sum_odd += _tJ(j,B) * appM(j,B);
398 else
399 *sum_even += _tJ(j,B) * appM(j,B);
400
401 ++S;
402 } while( !S.empty() );
403 }
404
405
406 void MR::solvemcav() {
407 double sum_even, sum_odd;
408 double maxdev;
409 size_t maxruns = 1000;
410
411 makekindex();
412 for(size_t i=0; i<N; i++)
413 for(size_t _j=0; _j<con[i]; _j++)
414 M[i][_j]=0.1;
415
416 size_t run=0;
417 do {
418 maxdev=0.0;
419 run++;
420 for(size_t i=0; i<N; i++){ // for all i
421 for(size_t _j=0; _j<con[i]; _j++){ // for all j in N_i
422 size_t _i = kindex[i][_j];
423 size_t j = nb[i][_j];
424 assert( nb[j][_i] == i );
425
426 double newM = 0.0;
427 if( props.updates == Properties::UpdateType::FULL ) {
428 // find indices in nb[j] that do not correspond with i
429 sub_nb _nbj_min_i(con[j]);
430 _nbj_min_i -= kindex[i][_j];
431
432 // find indices in nb[i] that do not correspond with j
433 sub_nb _nbi_min_j(con[i]);
434 _nbi_min_j -= _j;
435
436 sum_subs(j, _nbj_min_i, &sum_even, &sum_odd);
437 newM = (tanh(theta[j]) * sum_even + sum_odd) / (sum_even + tanh(theta[j]) * sum_odd);
438
439 sum_subs(i, _nbi_min_j, &sum_even, &sum_odd);
440 double denom = sum_even + tanh(theta[i]) * sum_odd;
441 double numer = 0.0;
442 for(size_t _k=0; _k<con[i]; _k++) if(_k != _j) {
443 sum_subs(i, _nbi_min_j - _k, &sum_even, &sum_odd);
444 numer += tJ[i][_k] * cors[i][_j][_k] * (tanh(theta[i]) * sum_even + sum_odd);
445 }
446 newM -= numer / denom;
447 } else if( props.updates == Properties::UpdateType::LINEAR ) {
448 newM = T(j,_i);
449 for(size_t _l=0; _l<con[i]; _l++) if( _l != _j )
450 newM -= Omega(i,_j,_l) * tJ[i][_l] * cors[i][_j][_l];
451 for(size_t _l1=0; _l1<con[j]; _l1++) if( _l1 != _i )
452 for( size_t _l2=_l1+1; _l2<con[j]; _l2++) if( _l2 != _i)
453 newM += Gamma(j,_i,_l1,_l2) * tJ[j][_l1] * tJ[j][_l2] * cors[j][_l1][_l2];
454 }
455
456 double dev = newM - M[i][_j];
457 // dev *= 0.02;
458 if( fabs(dev) >= maxdev )
459 maxdev = fabs(dev);
460
461 newM = M[i][_j] + dev;
462 if( fabs(newM) > 1.0 )
463 newM = sign(newM);
464 M[i][_j] = newM;
465 }
466 }
467 } while((maxdev>props.tol)&&(run<maxruns));
468
469 _iters = run;
470 if( maxdev > _maxdiff )
471 _maxdiff = maxdev;
472
473 if(run==maxruns){
474 if( props.verbose >= 1 )
475 cout << "solve_mcav: Convergence not reached (maxdev=" << maxdev << ")..." << endl;
476 }
477 }
478
479
480 void MR::solveM() {
481 for(size_t i=0; i<N; i++) {
482 if( props.updates == Properties::UpdateType::FULL ) {
483 // find indices in nb[i]
484 sub_nb _nbi(con[i]);
485
486 // calc numerator1 and denominator1
487 double sum_even, sum_odd;
488 sum_subs(i, _nbi, &sum_even, &sum_odd);
489
490 Mag[i] = (tanh(theta[i]) * sum_even + sum_odd) / (sum_even + tanh(theta[i]) * sum_odd);
491
492 } else if( props.updates == Properties::UpdateType::LINEAR ) {
493 sub_nb empty(con[i]);
494 empty.clear();
495 Mag[i] = T(i,empty);
496
497 for(size_t _l1=0; _l1<con[i]; _l1++)
498 for( size_t _l2=_l1+1; _l2<con[i]; _l2++)
499 Mag[i] += Gamma(i,_l1,_l2) * tJ[i][_l1] * tJ[i][_l2] * cors[i][_l1][_l2];
500 }
501 if(fabs(Mag[i])>1.)
502 Mag[i] = sign(Mag[i]);
503 }
504 }
505
506
507 void MR::init_cor() {
508 for( size_t i = 0; i < nrVars(); i++ ) {
509 vector<Factor> pairq;
510 if( props.inits == Properties::InitType::CLAMPING ) {
511 BP bpcav(*this, PropertySet()("updates", string("SEQMAX"))("tol", string("1e-9"))("maxiter", string("10000UL"))("verbose", string("0UL"))("logdomain", string("0")));
512 bpcav.makeCavity( i );
513 pairq = calcPairBeliefs( bpcav, delta(i), false );
514 } else if( props.inits == Properties::InitType::EXACT ) {
515 JTree jtcav(*this, PropertySet()("updates", string("HUGIN"))("verbose", string("0UL")) );
516 jtcav.makeCavity( i );
517 pairq = calcPairBeliefs( jtcav, delta(i), false );
518 }
519 for( size_t jk = 0; jk < pairq.size(); jk++ ) {
520 VarSet::const_iterator kit = pairq[jk].vars().begin();
521 size_t j = findVar( *(kit) );
522 size_t k = findVar( *(++kit) );
523 pairq[jk].normalize();
524 double cor = (pairq[jk][3] - pairq[jk][2] - pairq[jk][1] + pairq[jk][0]) - (pairq[jk][3] + pairq[jk][2] - pairq[jk][1] - pairq[jk][0]) * (pairq[jk][3] - pairq[jk][2] + pairq[jk][1] - pairq[jk][0]);
525 for( size_t _j = 0; _j < con[i]; _j++ ) if( nb[i][_j] == j )
526 for( size_t _k = 0; _k < con[i]; _k++ ) if( nb[i][_k] == k ) {
527 cors[i][_j][_k] = cor;
528 cors[i][_k][_j] = cor;
529 }
530 }
531 }
532 }
533
534
535 string MR::identify() const {
536 return string(Name) + printProperties();
537 }
538
539
540 double MR::run() {
541 if( supported ) {
542 if( props.verbose >= 1 )
543 cout << "Starting " << identify() << "...";
544
545 double tic = toc();
546 // Diffs diffs(nrVars(), 1.0);
547
548 M.resize(N);
549 for(size_t i=0; i<N; i++)
550 M[i].resize(kmax);
551
552 cors.resize(N);
553 for(size_t i=0; i<N; i++)
554 cors[i].resize(kmax);
555 for(size_t i=0; i<N; i++)
556 for(size_t j=0; j<kmax; j++)
557 cors[i][j].resize(kmax);
558
559 kindex.resize(N);
560 for(size_t i=0; i<N; i++)
561 kindex[i].resize(kmax);
562
563 if( props.inits == Properties::InitType::RESPPROP ) {
564 double md = init_cor_resp();
565 if( md > _maxdiff )
566 _maxdiff = md;
567 } else if( props.inits == Properties::InitType::EXACT )
568 init_cor(); // FIXME no MaxDiff() calculation
569 else if( props.inits == Properties::InitType::CLAMPING )
570 init_cor(); // FIXME no MaxDiff() calculation
571
572 solvemcav();
573
574 Mag.resize(N);
575 solveM();
576
577 if( props.verbose >= 1 )
578 cout << Name << " needed " << toc() - tic << " seconds." << endl;
579
580 return 0.0;
581 } else
582 return 1.0;
583 }
584
585
586 void MR::makekindex() {
587 for(size_t i=0; i<N; i++)
588 for(size_t j=0; j<con[i]; j++) {
589 size_t ij = nb[i][j]; // ij is the j'th neighbour of spin i
590 size_t k=0;
591 while( nb[ij][k] != i )
592 k++;
593 kindex[i][j] = k; // the j'th neighbour of spin i has spin i as its k'th neighbour
594 }
595 }
596
597
598 Factor MR::belief( const Var &n ) const {
599 if( supported ) {
600 size_t i = findVar( n );
601
602 Real x[2];
603 x[0] = 0.5 - Mag[i] / 2.0;
604 x[1] = 0.5 + Mag[i] / 2.0;
605
606 return Factor( n, x );
607 } else
608 return Factor();
609 }
610
611
612 vector<Factor> MR::beliefs() const {
613 vector<Factor> result;
614 for( size_t i = 0; i < nrVars(); i++ )
615 result.push_back( belief( var(i) ) );
616 return result;
617 }
618
619
620
621 MR::MR( const FactorGraph &fg, const PropertySet &opts ) : DAIAlgFG(fg), supported(true), _maxdiff(0.0), _iters(0) {
622 setProperties( opts );
623
624 // check whether all vars in fg are binary
625 // check whether connectivity is <= kmax
626 for( size_t i = 0; i < fg.nrVars(); i++ )
627 if( (fg.var(i).states() > 2) || (fg.delta(i).size() > kmax) ) {
628 supported = false;
629 break;
630 }
631
632 if( !supported )
633 return;
634
635 // check whether all interactions are pairwise or single
636 for( size_t I = 0; I < fg.nrFactors(); I++ )
637 if( fg.factor(I).vars().size() > 2 ) {
638 supported = false;
639 break;
640 }
641
642 if( !supported )
643 return;
644
645 // create w and th
646 size_t Nin = fg.nrVars();
647
648 double *w = new double[Nin*Nin];
649 double *th = new double[Nin];
650
651 for( size_t i = 0; i < Nin; i++ ) {
652 th[i] = 0.0;
653 for( size_t j = 0; j < Nin; j++ )
654 w[i*Nin+j] = 0.0;
655 }
656
657 for( size_t I = 0; I < fg.nrFactors(); I++ ) {
658 const Factor &psi = fg.factor(I);
659 if( psi.vars().size() == 1 ) {
660 size_t i = fg.findVar( *(psi.vars().begin()) );
661 th[i] += 0.5 * log(psi[1] / psi[0]);
662 } else if( psi.vars().size() == 2 ) {
663 size_t i = fg.findVar( *(psi.vars().begin()) );
664 VarSet::const_iterator jit = psi.vars().begin();
665 size_t j = fg.findVar( *(++jit) );
666
667 w[i*Nin+j] += 0.25 * log(psi[3] * psi[0] / (psi[2] * psi[1]));
668 w[j*Nin+i] += 0.25 * log(psi[3] * psi[0] / (psi[2] * psi[1]));
669
670 th[i] += 0.25 * log(psi[3] / psi[2] * psi[1] / psi[0]);
671 th[j] += 0.25 * log(psi[3] / psi[1] * psi[2] / psi[0]);
672 }
673 }
674
675 init(Nin, w, th);
676
677 delete th;
678 delete w;
679 }
680
681
682 } // end of namespace dai