aGrUM 3.1.1
a C++ library for (probabilistic) graphical models
loopyBeliefPropagation_tpl.h
Go to the documentation of this file.
1/****************************************************************************
2 * This file is part of the aGrUM/pyAgrum library. *
3 * *
4 * Copyright (c) 2005-2026 by *
5 * - Pierre-Henri WUILLEMIN(_at_LIP6) *
6 * - Christophe GONZALES(_at_AMU) *
7 * *
8 * The aGrUM/pyAgrum library is free software; you can redistribute it *
9 * and/or modify it under the terms of either : *
10 * *
11 * - the GNU Lesser General Public License as published by *
12 * the Free Software Foundation, either version 3 of the License, *
13 * or (at your option) any later version, *
14 * - the MIT license (MIT), *
15 * - or both in dual license, as here. *
16 * *
17 * (see https://agrum.gitlab.io/articles/dual-licenses-lgplv3mit.html) *
18 * *
19 * This aGrUM/pyAgrum library is distributed in the hope that it will be *
20 * useful, but WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, *
21 * INCLUDING BUT NOT LIMITED TO THE WARRANTIES MERCHANTABILITY or FITNESS *
22 * FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE *
23 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
24 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
25 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR *
26 * OTHER DEALINGS IN THE SOFTWARE. *
27 * *
28 * See LICENCES for more details. *
29 * *
30 * SPDX-FileCopyrightText: Copyright 2005-2026 *
31 * - Pierre-Henri WUILLEMIN(_at_LIP6) *
32 * - Christophe GONZALES(_at_AMU) *
33 * SPDX-License-Identifier: LGPL-3.0-or-later OR MIT *
34 * *
35 * Contact : info_at_agrum_dot_org *
36 * homepage : http://agrum.gitlab.io *
37 * gitlab : https://gitlab.com/agrumery/agrum *
38 * *
39 ****************************************************************************/
40
41#pragma once
42
43
48#include <agrum/BN/inference/loopyBeliefPropagation.h> // to ease IDE parser
49#ifndef DOXYGEN_SHOULD_SKIP_THIS
50
51# include <algorithm>
52# include <sstream>
53# include <string>
54
55constexpr auto LBP_DEFAULT_MAXITER = 100;
56constexpr auto LBP_DEFAULT_EPSILON = 1e-8;
57constexpr auto LBP_DEFAULT_MIN_EPSILON_RATE = 1e-10;
58constexpr auto LBP_DEFAULT_PERIOD_SIZE = 1;
59constexpr auto LBP_DEFAULT_VERBOSITY = false;
60
61
62// to ease parsing for IDE
64
65namespace gum {
66
68 template < GUM_Numeric GUM_SCALAR >
70 ApproximateInference< GUM_SCALAR >(bn) {
71 // for debugging purposes
72 GUM_CONSTRUCTOR(LoopyBeliefPropagation)
73
74 this->setEpsilon(LBP_DEFAULT_EPSILON);
75 this->setMinEpsilonRate(LBP_DEFAULT_MIN_EPSILON_RATE);
76 this->setMaxIter(LBP_DEFAULT_MAXITER);
77 this->setVerbosity(LBP_DEFAULT_VERBOSITY);
78 this->setPeriodSize(LBP_DEFAULT_PERIOD_SIZE);
79
80 _init_messages_();
81 }
82
84 template < GUM_Numeric GUM_SCALAR >
85 LoopyBeliefPropagation< GUM_SCALAR >::~LoopyBeliefPropagation() {
86 GUM_DESTRUCTOR(LoopyBeliefPropagation)
87 }
88
89 template < GUM_Numeric GUM_SCALAR >
90 void LoopyBeliefPropagation< GUM_SCALAR >::_init_messages_() {
91 _messages_.clear();
92 for (const auto& tail: this->BN().nodes()) {
93 Tensor< GUM_SCALAR > p;
94 p.add(this->BN().variable(tail));
95 p.fill(static_cast< GUM_SCALAR >(1));
96
97 for (const auto& head: this->BN().children(tail)) {
98 _messages_.insert(Arc(head, tail), p);
99 _messages_.insert(Arc(tail, head), p);
100 }
101 }
102 }
103
104 template < GUM_Numeric GUM_SCALAR >
105 void LoopyBeliefPropagation< GUM_SCALAR >::updateOutdatedStructure_() {
106 _init_messages_();
107 }
108
109 template < GUM_Numeric GUM_SCALAR >
110 Tensor< GUM_SCALAR > LoopyBeliefPropagation< GUM_SCALAR >::_computeProdPi_(NodeId X) {
111 const auto& varX = this->BN().variable(X);
112
113 auto piX = this->BN().cpt(X);
114 for (const auto& U: this->BN().parents(X)) {
115 piX *= _messages_[Arc(U, X)];
116 }
117 piX = piX.sumIn({&varX});
118
119 return piX;
120 }
121
122 template < GUM_Numeric GUM_SCALAR >
123 Tensor< GUM_SCALAR > LoopyBeliefPropagation< GUM_SCALAR >::_computeProdPi_(NodeId X,
124 NodeId except) {
125 const auto& varX = this->BN().variable(X);
126 const auto& varExcept = this->BN().variable(except);
127 auto piXexcept = this->BN().cpt(X);
128 for (const auto& U: this->BN().parents(X)) {
129 if (U != except) { piXexcept *= _messages_[Arc(U, X)]; }
130 }
131 piXexcept = piXexcept.sumIn({&varX, &varExcept});
132 return piXexcept;
133 }
134
135 template < GUM_Numeric GUM_SCALAR >
136 Tensor< GUM_SCALAR > LoopyBeliefPropagation< GUM_SCALAR >::_computeProdLambda_(NodeId X) {
137 Tensor< GUM_SCALAR > lamX;
138 if (this->hasEvidence(X)) {
139 lamX = *(this->evidence()[X]);
140 } else {
141 lamX.add(this->BN().variable(X));
142 lamX.fill(1);
143 }
144 for (const auto& Y: this->BN().children(X)) {
145 lamX *= _messages_[Arc(Y, X)];
146 }
147
148 return lamX;
149 }
150
151 template < GUM_Numeric GUM_SCALAR >
152 Tensor< GUM_SCALAR > LoopyBeliefPropagation< GUM_SCALAR >::_computeProdLambda_(NodeId X,
153 NodeId except) {
154 Tensor< GUM_SCALAR > lamXexcept;
155 if (this->hasEvidence(X)) { //
156 lamXexcept = *this->evidence()[X];
157 } else {
158 lamXexcept.add(this->BN().variable(X));
159 lamXexcept.fill(1);
160 }
161 for (const auto& Y: this->BN().children(X)) {
162 if (Y != except) { lamXexcept *= _messages_[Arc(Y, X)]; }
163 }
164
165 return lamXexcept;
166 }
167
168 template < GUM_Numeric GUM_SCALAR >
169 GUM_SCALAR LoopyBeliefPropagation< GUM_SCALAR >::_updateNodeMessage_(NodeId X) {
170 auto piX = _computeProdPi_(X);
171 auto lamX = _computeProdLambda_(X);
172
173 GUM_SCALAR KL = 0;
174 Arc argKL(0, 0);
175
176 // update lambda_par (for arc U->x)
177 for (const auto& U: this->BN().parents(X)) {
178 auto newLambda = (_computeProdPi_(X, U) * lamX).sumIn({&this->BN().variable(U)});
179 newLambda.normalize();
180 auto ekl = static_cast< GUM_SCALAR >(0);
181 try {
182 ekl = _messages_[Arc(X, U)].KL(newLambda);
183 } catch (InvalidArgument const&) {
184 GUM_ERROR(InvalidArgument, "Not compatible pi during computation")
185 } catch (FatalError const&) { // 0 misplaced
186 ekl = std::numeric_limits< GUM_SCALAR >::infinity();
187 }
188 if (ekl > KL) {
189 KL = ekl;
190 argKL = Arc(X, U);
191 }
192 _messages_.set(Arc(X, U), newLambda);
193 }
194
195 // update pi_child (for arc x->child)
196 for (const auto& Y: this->BN().children(X)) {
197 auto newPi = (piX * _computeProdLambda_(X, Y));
198 newPi.normalize();
199 GUM_SCALAR ekl = KL;
200 try {
201 ekl = _messages_[Arc(X, Y)].KL(newPi);
202 } catch (InvalidArgument const&) {
203 GUM_ERROR(InvalidArgument, "Not compatible pi during computation")
204 } catch (FatalError const&) { // 0 misplaced
205 ekl = std::numeric_limits< GUM_SCALAR >::infinity();
206 }
207 if (ekl > KL) {
208 KL = ekl;
209 argKL = Arc(X, Y);
210 }
211 _messages_.set(Arc(X, Y), newPi);
212 }
213
214 return KL;
215 }
216
217 template < GUM_Numeric GUM_SCALAR >
218 void LoopyBeliefPropagation< GUM_SCALAR >::_initStats_() {
219 _init_messages_();
220 for (const auto& node: this->BN().topologicalOrder()) {
221 _updateNodeMessage_(node);
222 }
223 }
224
226 template < GUM_Numeric GUM_SCALAR >
227 void LoopyBeliefPropagation< GUM_SCALAR >::makeInference_() {
228 _initStats_();
229 this->initApproximationScheme();
230
231 std::vector< NodeId > shuffleIds;
232 for (const auto& node: this->BN().nodes())
233 shuffleIds.push_back(node);
234
235 auto engine = std::default_random_engine{};
236
237 GUM_SCALAR error = 0.0;
238 do {
239 std::shuffle(std::begin(shuffleIds), std::end(shuffleIds), engine);
240 this->updateApproximationScheme();
241 for (const auto& node: shuffleIds) {
242 GUM_SCALAR e = _updateNodeMessage_(node);
243 if (e > error) error = e;
244 }
245 } while (this->continueApproximationScheme(error));
246 }
247
249 template < GUM_Numeric GUM_SCALAR >
250 const Tensor< GUM_SCALAR >& LoopyBeliefPropagation< GUM_SCALAR >::posterior_(NodeId id) {
251 auto p = _computeProdPi_(id) * _computeProdLambda_(id);
252 p.normalize();
253 _posteriors_.set(id, p);
254
255 return _posteriors_[id];
256 }
257
258 template < GUM_Numeric GUM_SCALAR >
259 void LoopyBeliefPropagation< GUM_SCALAR >::onStateChanged_() {}
260
261 template < GUM_Numeric GUM_SCALAR >
262 void LoopyBeliefPropagation< GUM_SCALAR >::onEvidenceAdded_(const NodeId id,
263 bool isHardEvidence) {}
264
265 template < GUM_Numeric GUM_SCALAR >
266 void LoopyBeliefPropagation< GUM_SCALAR >::onEvidenceErased_(const NodeId id,
267 bool isHardEvidence) {}
268
269 template < GUM_Numeric GUM_SCALAR >
270 void LoopyBeliefPropagation< GUM_SCALAR >::onAllEvidenceErased_(bool contains_hard_evidence) {}
271
272 template < GUM_Numeric GUM_SCALAR >
273 void LoopyBeliefPropagation< GUM_SCALAR >::onEvidenceChanged_(const NodeId id,
274 bool hasChangedSoftHard) {}
275
276 template < GUM_Numeric GUM_SCALAR >
277 void LoopyBeliefPropagation< GUM_SCALAR >::onModelChanged_(const GraphicalModel* bn) {}
278
279 template < GUM_Numeric GUM_SCALAR >
280 void LoopyBeliefPropagation< GUM_SCALAR >::updateOutdatedTensors_() {}
281
282 template < GUM_Numeric GUM_SCALAR >
283 void LoopyBeliefPropagation< GUM_SCALAR >::onMarginalTargetAdded_(const NodeId id) {}
284
285 template < GUM_Numeric GUM_SCALAR >
286 void LoopyBeliefPropagation< GUM_SCALAR >::onMarginalTargetErased_(const NodeId id) {}
287
288 template < GUM_Numeric GUM_SCALAR >
289 void LoopyBeliefPropagation< GUM_SCALAR >::onAllMarginalTargetsAdded_() {}
290
291 template < GUM_Numeric GUM_SCALAR >
292 void LoopyBeliefPropagation< GUM_SCALAR >::onAllMarginalTargetsErased_() {}
293
294} /* namespace gum */
295
296#endif // DOXYGEN_SHOULD_SKIP_THIS
KL is the base class for KL computation betweens 2 BNs.
Exception : fatal (unknown ?) error.
Class representing the minimal interface for Bayesian network with no numerical data.
Definition IBayesNet.h:75
Exception: at least one argument passed to a function is not what was expected.
LoopyBeliefPropagation(const IBayesNet< GUM_SCALAR > *bn)
Default constructor.
#define GUM_ERROR(type, msg)
Definition exceptions.h:76
This file contains gibbs sampling (for BNs) class definitions.
Sequence< NodeId > topologicalOrder(const G &g)
Returns a topological ordering of the nodes of g (Kahn's algorithm).
gum is the global namespace for all aGrUM entities
Definition agrum.h:46