[Frederik Eaton] Added BP_Dual, BBP and CBP algorithms
[libdai.git] / include / dai / bbp.h
1 #ifndef ___defined_libdai_bbp_h
2 #define ___defined_libdai_bbp_h
3
4 #include <vector>
5 #include <utility>
6 #include <ext/hash_map>
7
8 #include <boost/tuple/tuple.hpp>
9
10 #include <dai/prob.h>
11 #include <dai/daialg.h>
12 #include <dai/factorgraph.h>
13 #include <dai/enum.h>
14
15 #include <dai/bp_dual.h>
16
17 namespace dai {
18
19 using namespace std;
20 using namespace __gnu_cxx;
21 using boost::tuple;
22 using boost::get;
23
24 // utility function to subtract 1 from a size_t, or return 0 if the
25 // argument is 0
26 inline size_t oneLess(size_t v) { return v==0?v:v-1; }
27
28 /// function to compute adj_w_unnorm from w, Z_w, adj_w
29 Prob unnormAdjoint(const Prob &w, Real Z_w, const Prob &adj_w);
30
31 vector<size_t> getGibbsState(const BP_dual& fg, size_t iters);
32
33 vector<Prob> get_zero_adj_2(const BP_dual& bp_dual);
34 vector<Prob> get_zero_adj_1(const BP_dual& bp_dual);
35
36 class BBP {
37 // ----------------------------------------------------------------
38 // inputs
39 const BP_dual* _bp_dual;
40 const FactorGraph* _fg;
41 vector<Prob> _adj_b_1, _adj_b_2;
42 vector<Prob> _adj_b_1_unnorm, _adj_b_2_unnorm;
43 // in case the objective function depends on factors as well:
44 vector<Prob> _init_adj_psi_1;
45 vector<Prob> _init_adj_psi_2;
46 // ----------------------------------------------------------------
47 // ouptuts
48 vector<Prob> _adj_psi_1, _adj_psi_2;
49 // the following vectors are length _fg->nr_edges()
50 vector<Prob> _adj_n, _adj_m;
51 vector<Prob> _adj_n_unnorm, _adj_m_unnorm;
52 vector<Prob> _new_adj_n, _new_adj_m;
53 // ----------------------------------------------------------------
54 // auxiliary data
55 typedef vector<hash_map<size_t, vector<size_t> > > _trip_t;
56 typedef vector<size_t> _trip_end_t;
57
58 // members to store indices of VFV and FVF triples
59 // XXX need to have vector going in other direction so that we can
60 // do faster looping on iIj and IiJ
61 _trip_t _VFV;
62 _trip_t _FVF;
63
64 typedef vector<tuple<size_t,size_t,size_t> > _trip_ind_t;
65 _trip_ind_t _VFV_ind;
66 _trip_ind_t _FVF_ind;
67
68 // helper quantities computed from the BP messages
69 vector<Prob> _T,_U,_S,_R;
70 size_t _iters;
71 size_t _max_iter;
72
73 public:
74 void RegenerateInds();
75 void RegenerateT();
76 void RegenerateU();
77 void RegenerateS();
78 void RegenerateR();
79 void RegenerateInputs();
80 /// initialise members for messages and factor adjoints
81 void RegeneratePsiAdjoints();
82 void RegenerateMessageAdjoints();
83 void Regenerate();
84
85 size_t VFV(size_t i, size_t I, size_t j) const { return (*const_cast<_trip_t*>(&_VFV))[i][I][j]; }
86 size_t FVF(size_t I, size_t i, size_t J) const { return (*const_cast<_trip_t*>(&_FVF))[I][i][J]; }
87 DAI_ACCMUT(Prob & T(size_t i, size_t I), { return _T[_fg->VV2E(i,I)]; });
88 DAI_ACCMUT(Prob & U(size_t I, size_t i), { return _U[_fg->VV2E(i,I)]; });
89 DAI_ACCMUT(Prob & S(size_t i, size_t I, size_t j), { return _S[VFV(i,I,j)]; });
90 DAI_ACCMUT(Prob & R(size_t I, size_t i, size_t J), { return _R[FVF(I,i,J)]; });
91
92 DAI_ACCMUT(Prob & T(size_t iI), { return _T[iI]; });
93 DAI_ACCMUT(Prob & U(size_t iI), { return _U[iI]; });
94 DAI_ACCMUT(Prob & S(size_t iIj), { return _S[iIj]; });
95 DAI_ACCMUT(Prob & R(size_t IiJ), { return _R[IiJ]; });
96
97 void calcNewN(size_t iI);
98 void calcNewM(size_t iI);
99 void upMsgM(size_t iI);
100 void upMsgN(size_t iI);
101 void doParUpdate();
102 Real getUnMsgMag();
103 void getMsgMags(Real &s, Real &new_s);
104
105 void zero_adj_b_2() {
106 _adj_b_2.clear();
107 _adj_b_2.reserve(_bp_dual->nrFactors());
108 for(size_t I=0; I<_bp_dual->nrFactors(); I++) {
109 _adj_b_2.push_back(Prob(_bp_dual->factor(I).states(),Real(0.0)));
110 }
111 }
112 public:
113 BBP(const BP_dual &bp_dual) :
114 _bp_dual(&bp_dual), _fg(&bp_dual), _max_iter(50) {}
115
116 DAI_ACCMUT(size_t& maxIter(), { return _max_iter; });
117
118 void init(const vector<Prob> &adj_b_1, const vector<Prob> &adj_b_2,
119 const vector<Prob> &adj_psi_1, const vector<Prob> &adj_psi_2) {
120 _adj_b_1 = adj_b_1;
121 _adj_b_2 = adj_b_2;
122 _init_adj_psi_1 = adj_psi_1;
123 _init_adj_psi_2 = adj_psi_2;
124 Regenerate();
125 }
126 void init(const vector<Prob> &adj_b_1, const vector<Prob> &adj_b_2) {
127 init(adj_b_1, adj_b_2, get_zero_adj_1(*_bp_dual), get_zero_adj_2(*_bp_dual));
128 }
129 void init(const vector<Prob> &adj_b_1) {
130 init(adj_b_1, get_zero_adj_2(*_bp_dual));
131 }
132
133 /// run until change is less than given tolerance
134 void run(Real tol);
135
136 size_t doneIters() { return _iters; }
137
138 DAI_ACCMUT(Prob& adj_psi_1(size_t i), { return _adj_psi_1[i]; });
139 DAI_ACCMUT(Prob& adj_psi_2(size_t I), { return _adj_psi_2[I]; });
140 DAI_ACCMUT(Prob& adj_b_1(size_t i), { return _adj_b_1[i]; });
141 DAI_ACCMUT(Prob& adj_b_2(size_t I), { return _adj_b_2[I]; });
142 DAI_ACCMUT(Prob& adj_n(size_t i, size_t I), { return _adj_n[_fg->VV2E(i,I)]; });
143 DAI_ACCMUT(Prob& adj_m(size_t I, size_t i), { return _adj_m[_fg->VV2E(i,I)]; });
144 };
145
146 /// Cost functions. Not used by BBP class, only used by following functions
147 DAI_ENUM(bbp_cfn_t, cfn_b, cfn_b2, cfn_exp, cfn_var_ent, cfn_factor_ent, cfn_bethe_ent);
148
149 /// initialise BBP using BP_dual and (cfn_type, stateP)
150 void gibbsInitBBPCostFnAdj(BBP& bbp, const BP_dual& fg, bbp_cfn_t cfn_type, const vector<size_t>* stateP);
151
152 /// calculate actual value of cost function (cfn_type, stateP)
153 Real gibbsCostFn(const BP_dual& fg, bbp_cfn_t cfn_type, const vector<size_t> *stateP);
154
155 /// function to test the validity of adjoints computed by BBP
156 double numericBBPTest(const BP_dual& bp_dual, bbp_cfn_t cfn, double bbpTol, double h);
157
158 /// output to "./bbp-data/" a set of files which evaluate the effect
159 /// on the given cost function (for a random Gibbs state) of
160 /// perturbing each graph variable
161 void makeBBPGraph(const BP_dual& bp_dual, bbp_cfn_t cfn);
162
163 }
164
165 #endif