include/dai/mr.h

   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 #ifndef __defined_libdai_mr_h
  23 #define __defined_libdai_mr_h
  24
  25
  26 #include <vector>
  27 #include <string>
  28 #include <dai/factorgraph.h>
  29 #include <dai/daialg.h>
  30 #include <dai/enum.h>
  31
  32
  33 namespace dai {
  34
  35
  36 class sub_nb;
  37
  38
  39 class MR : public DAIAlgFG {
  40     private:
  41         bool supported;                                            // is the underlying factor graph supported?
  42
  43         std::vector<size_t>                             con;       // con[i] = connectivity of spin i
  44         std::vector<std::vector<size_t> >               nb;        // nb[i] are the neighbours of spin i
  45         std::vector<std::vector<double> >               tJ;        // tJ[i][_j] is the tanh of the interaction between spin i and its neighbour nb[i][_j]
  46         std::vector<double>                             theta;     // theta[i] is the local field on spin i
  47         std::vector<std::vector<double> >               M;         // M[i][_j] is M^{(i)}_j
  48         std::vector<std::vector<size_t> >               kindex;    // the _j'th neighbour of spin i has spin i as its kindex[i][_j]'th neighbour
  49         std::vector<std::vector<std::vector<double> > > cors;
  50
  51         static const size_t kmax = 31;
  52
  53         size_t N;
  54
  55         std::vector<double> Mag;
  56
  57     public:
  58         ENUM2(UpdateType,FULL,LINEAR)
  59         ENUM3(InitType,RESPPROP,CLAMPING,EXACT)
  60
  61         UpdateType Updates() const { return GetPropertyAs<UpdateType>("updates"); }
  62         InitType Inits() const { return GetPropertyAs<InitType>("inits"); }
  63
  64         MR( const FactorGraph & fg, const Properties &opts );
  65         void init(size_t _N, double *_w, double *_th);
  66         void makekindex();
  67         void read_files();
  68         void init_cor();
  69         double init_cor_resp();
  70         void solvemcav();
  71         void solveM();
  72         double run();
  73         Factor belief( const Var &n ) const;
  74         Factor belief( const VarSet &/*ns*/ ) const { assert( 0 == 1 ); }
  75         std::vector<Factor> beliefs() const;
  76         Complex logZ() const { return NAN; }
  77         void init() { assert( checkProperties() ); }
  78         static const char *Name;
  79         std::string identify() const;
  80         double _tJ(size_t i, sub_nb A);
  81
  82         double Omega(size_t i, size_t _j, size_t _l);
  83         double T(size_t i, sub_nb A);
  84         double T(size_t i, size_t _j);
  85         double Gamma(size_t i, size_t _j, size_t _l1, size_t _l2);
  86         double Gamma(size_t i, size_t _l1, size_t _l2);
  87
  88         double appM(size_t i, sub_nb A);
  89         void sum_subs(size_t j, sub_nb A, double *sum_even, double *sum_odd);
  90
  91         double sign(double a) { return (a >= 0) ? 1.0 : -1.0; }
  92         MR* clone() const { assert( 0 == 1 ); }
  93
  94         bool checkProperties();
  95
  96 };
  97
  98
  99 // represents a subset of nb[i] as a binary number
 100 // the elements of a subset should be thought of as indices in nb[i]
 101 class sub_nb {
 102     private:
 103         size_t s;
 104         size_t bits;
 105
 106     public:
 107         // construct full subset containing nr_elmt elements
 108         sub_nb(size_t nr_elmt) {
 109 #ifdef DEBUG
 110             assert( nr_elmt < sizeof(size_t) / sizeof(char) * 8 );
 111 #endif
 112             bits = nr_elmt;
 113             s = (1U << bits) - 1;
 114         }
 115
 116         // copy constructor
 117         sub_nb( const sub_nb & x ) : s(x.s), bits(x.bits) {}
 118
 119         // assignment operator
 120         sub_nb & operator=( const sub_nb &x ) {
 121             if( this != &x ) {
 122                 s = x.s;
 123                 bits = x.bits;
 124             }
 125             return *this;
 126         }
 127
 128         // returns number of elements
 129         size_t size() {
 130             size_t size = 0;
 131             for(size_t bit = 0; bit < bits; bit++)
 132                 if( s & (1U << bit) )
 133                     size++;
 134             return size;
 135         }
 136
 137         // increases s by one (for enumeration in lexicographical order)
 138         sub_nb operator++() {
 139             s++;
 140             if( s >= (1U << bits) )
 141                 s = 0;
 142             return *this;
 143         }
 144
 145         // return i'th element of this subset
 146         size_t operator[](size_t i) {
 147             size_t bit;
 148             for(bit = 0; bit < bits; bit++ )
 149                 if( s & (1U << bit) ) {
 150                     if( i == 0 )
 151                         break;
 152                     else
 153                         i--;
 154                 }
 155 #ifdef DEBUG
 156             assert( bit < bits );
 157 #endif
 158             return bit;
 159         }
 160
 161         // add index _j to this subset
 162         sub_nb &operator +=(size_t _j) {
 163             s |= (1U << _j);
 164             return *this;
 165         }
 166
 167         // return copy with index _j
 168         sub_nb operator+(size_t _j) {
 169             sub_nb x = *this;
 170             x += _j;
 171             return x;
 172         }
 173
 174         // delete index _j from this subset
 175         sub_nb &operator -=(size_t _j) {
 176             s &= ~(1U << _j);
 177             return *this;
 178         }
 179
 180         // return copy without index _j
 181         sub_nb operator-(size_t _j) {
 182             sub_nb x = *this;
 183             x -= _j;
 184             return x;
 185         }
 186
 187         // empty this subset
 188         sub_nb & clear() {
 189             s = 0;
 190             return *this;
 191         }
 192
 193         // returns true if subset is empty
 194         bool empty() { return (s == 0); }
 195
 196         // return 1 if _j is contained, 0 otherwise ("is element of")
 197         size_t operator&(size_t _j) { return s & (1U << _j); }
 198
 199         friend std::ostream& operator<< (std::ostream& os, const sub_nb x) {
 200             if( x.bits == 0 )
 201                 os << "empty";
 202             else {
 203                 for(size_t bit = x.bits; bit > 0; bit-- )
 204                     if( x.s & (1U << (bit-1)) )
 205                         os << "1";
 206                     else
 207                         os << "0";
 208             }
 209             return os;
 210         }
 211 };
 212
 213
 214 } // end of namespace dai
 215
 216
 217 #endif