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