aGrUM  0.14.1
GibbsKL_tpl.h
Go to the documentation of this file.
1 /***************************************************************************
2  * Copyright (C) 2005 by Pierre-Henri WUILLEMIN et Christophe GONZALES *
3  * {prenom.nom}_at_lip6.fr *
4  * *
5  * This program is free software; you can redistribute it and/or modify *
6  * it under the terms of the GNU General Public License as published by *
7  * the Free Software Foundation; either version 2 of the License, or *
8  * (at your option) any later version. *
9  * *
10  * This program is distributed in the hope that it will be useful, *
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13  * GNU General Public License for more details. *
14  * *
15  * You should have received a copy of the GNU General Public License *
16  * along with this program; if not, write to the *
17  * Free Software Foundation, Inc., *
18  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
19  ***************************************************************************/
27 #include <agrum/core/math/math.h>
28 #include <agrum/BN/IBayesNet.h>
32 #include <agrum/core/hashTable.h>
33 
34 #define GIBBSKL_DEFAULT_MAXITER 10000000
35 #define GIBBSKL_DEFAULT_EPSILON 1e-10
36 #define GIBBSKL_DEFAULT_MIN_EPSILON_RATE 1e-10
37 #define GIBBSKL_DEFAULT_PERIOD_SIZE 200
38 #define GIBBSKL_DEFAULT_VERBOSITY false
39 #define GIBBSKL_DEFAULT_BURNIN 2000
40 #define GIBBSKL_DEFAULT_TIMEOUT 6000
41 
42 #define GIBBSKL_POURCENT_DRAWN_SAMPLE 10 // percent drawn
43 #define GIBBSKL_DRAWN_AT_RANDOM false
44 
45 namespace gum {
46 
47 
48  template < typename GUM_SCALAR >
51  BNdistance< GUM_SCALAR >(P, Q),
53  GibbsOperator< GUM_SCALAR >(
54  P,
55  nullptr,
56  1 + (P.size() * GIBBSKL_POURCENT_DRAWN_SAMPLE / 100),
58  GUM_CONSTRUCTOR(GibbsBNdistance);
59 
67  }
68 
69  template < typename GUM_SCALAR >
71  const BNdistance< GUM_SCALAR >& kl) :
72  BNdistance< GUM_SCALAR >(kl),
74  // Gibbs operator with 10% of nodes changes at random between each samples
75  ,
76  GibbsOperator< GUM_SCALAR >(
77  kl.p(),
78  nullptr,
79  1 + (kl.p().size() * GIBBSKL_POURCENT_DRAWN_SAMPLE / 100),
80  true) {
81  GUM_CONSTRUCTOR(GibbsBNdistance);
82 
90  }
91 
92  template < typename GUM_SCALAR >
94  GUM_DESTRUCTOR(GibbsBNdistance);
95  }
96 
97  template < typename GUM_SCALAR >
99  auto Iq = _q.completeInstantiation();
100 
103 
104  // map between particle() variables and _q variables (using name of vars)
106 
107  for (Idx ite = 0; ite < I.nbrDim(); ++ite) {
108  map.insert(&I.variable(ite), &_q.variableFromName(I.variable(ite).name()));
109  }
110 
111  // BURN IN
112  for (Idx i = 0; i < burnIn(); i++)
113  I = this->nextSample(I);
114 
115  // SAMPLING
116  _klPQ = _klQP = _hellinger = _jsd = (GUM_SCALAR)0.0;
117  _errorPQ = _errorQP = 0;
119  GUM_SCALAR delta, ratio, error;
120  delta = ratio = error = (GUM_SCALAR)-1;
121  GUM_SCALAR oldPQ = 0.0;
122  GUM_SCALAR pp, pq, pmid;
123 
124  do {
125  this->disableMinEpsilonRate();
126  I = this->nextSample(I);
128 
129  //_p.synchroInstantiations( Ip,I);
130  Iq.setValsFrom(map, I);
131 
132  pp = _p.jointProbability(I);
133  pq = _q.jointProbability(Iq);
134  pmid = (pp + pq) / 2.0;
135 
136  if (pp != (GUM_SCALAR)0.0) {
137  _hellinger += std::pow(std::sqrt(pp) - std::sqrt(pq), 2) / pp;
138 
139  if (pq != (GUM_SCALAR)0.0) {
140  _bhattacharya += std::sqrt(pq / pp); // std::sqrt(pp*pq)/pp
142  this->enableMinEpsilonRate(); // replace check_rate=true;
143  ratio = pq / pp;
144  delta = (GUM_SCALAR)log2(ratio);
145  _klPQ += delta;
146 
147  // pmid!=0
148  _jsd -= log2(pp / pmid) + ratio * log2(pq / pmid);
149  } else {
150  _errorPQ++;
151  }
152  }
153 
154  if (pq != (GUM_SCALAR)0.0) {
155  if (pp != (GUM_SCALAR)0.0) {
156  // if we are here, it is certain that delta and ratio have been
157  // computed
158  // further lines above. (for now #112-113)
159  _klQP += (GUM_SCALAR)(-delta * ratio);
160  } else {
161  _errorQP++;
162  }
163  }
164 
165  if (this->isEnabledMinEpsilonRate()) { // replace check_rate
166  // delta is used as a temporary variable
167  delta = _klPQ / nbrIterations();
168  error = (GUM_SCALAR)std::abs(delta - oldPQ);
169  oldPQ = delta;
170  }
171  } while (continueApproximationScheme(error)); //
172 
173  _klPQ = -_klPQ / (nbrIterations());
174  _klQP = -_klQP / (nbrIterations());
175  _jsd = -0.5 * _jsd / (nbrIterations());
176  _hellinger = std::sqrt(_hellinger / nbrIterations());
177  _bhattacharya = -std::log(_bhattacharya / (nbrIterations()));
178  }
179 
180  template < typename GUM_SCALAR >
182  this->_burn_in = b;
183  }
184 
185  template < typename GUM_SCALAR >
187  return this->_burn_in;
188  }
189 } // namespace gum
Useful macros for maths.
This file contains Gibbs sampling (for BNs) class definitions.
This file contains general scheme for iteratively convergent algorithms.
#define GIBBSKL_DEFAULT_TIMEOUT
Definition: GibbsKL_tpl.h:40
void disableMinEpsilonRate()
Disable stopping criterion on epsilon rate.
void _computeKL() final
Definition: GibbsKL_tpl.h:98
GUM_SCALAR _klPQ
Definition: BNdistance.h:136
GUM_SCALAR _bhattacharya
Definition: BNdistance.h:139
Idx nbrDim() const final
Returns the number of variables in the Instantiation.
#define GIBBSKL_DRAWN_AT_RANDOM
Definition: GibbsKL_tpl.h:43
Approximation Scheme.
const IBayesNet< GUM_SCALAR > & _p
Definition: BNdistance.h:133
void enableMinEpsilonRate()
Enable stopping criterion on epsilon rate.
#define GIBBSKL_DEFAULT_BURNIN
Definition: GibbsKL_tpl.h:39
void setPeriodSize(Size p)
How many samples between two stopping is enable.
#define GIBBSKL_POURCENT_DRAWN_SAMPLE
Definition: GibbsKL_tpl.h:42
Class representing Bayesian networks.
GUM_SCALAR _jsd
Definition: BNdistance.h:140
const IBayesNet< GUM_SCALAR > & p() const
void initApproximationScheme()
Initialise the scheme.
Class representing the minimal interface for Bayesian Network.
Definition: IBayesNet.h:59
gum is the global namespace for all aGrUM entities
Definition: agrum.h:25
void setMinEpsilonRate(double rate)
Given that we approximate f(t), stopping criterion on d/dt(|f(t+1)-f(t)|).
algorithm for approximated computation KL divergence between BNs using GIBBS sampling ...
void setVerbosity(bool v)
Set the verbosity on (true) or off (false).
The class for generic Hash Tables.
Definition: hashTable.h:676
void setMaxTime(double timeout)
Stopping criterion on timeout.
Size _burn_in
Number of iterations before checking stopping criteria.
Size nbrIterations() const
Returns the number of iterations.
~GibbsBNdistance()
destructor
Definition: GibbsKL_tpl.h:93
const IBayesNet< GUM_SCALAR > & _q
Definition: BNdistance.h:134
bool continueApproximationScheme(double error)
Update the scheme w.r.t the new error.
#define GIBBSKL_DEFAULT_EPSILON
Definition: GibbsKL_tpl.h:35
GUM_SCALAR _hellinger
Definition: BNdistance.h:138
bool isEnabledMinEpsilonRate() const
Returns true if stopping criterion on epsilon rate is enabled, false otherwise.
Instantiation nextSample(Instantiation prev)
draws next sample of Gibbs sampling
GibbsKL computes the KL divergence betweens 2 BNs using an approximation pattern: GIBBS sampling...
#define GIBBSKL_DEFAULT_PERIOD_SIZE
Definition: GibbsKL_tpl.h:37
Size burnIn() const
Returns the number of burn in.
Definition: GibbsKL_tpl.h:186
Class for assigning/browsing values to tuples of discrete variables.
Definition: instantiation.h:80
void setBurnIn(Size b)
Number of burn in for one iteration.
Definition: GibbsKL_tpl.h:181
void setMaxIter(Size max)
Stopping criterion on number of iterations.
GUM_SCALAR _klQP
Definition: BNdistance.h:137
Instantiation monteCarloSample()
draws a Monte Carlo sample
#define GIBBSKL_DEFAULT_MIN_EPSILON_RATE
Definition: GibbsKL_tpl.h:36
void setEpsilon(double eps)
Given that we approximate f(t), stopping criterion on |f(t+1)-f(t)|.
Size Idx
Type for indexes.
Definition: types.h:50
std::size_t Size
In aGrUM, hashed values are unsigned long int.
Definition: types.h:45
value_type & insert(const Key &key, const Val &val)
Adds a new element (actually a copy of this element) into the hash table.
const std::string & name() const
returns the name of the variable
const DiscreteVariable & variable(Idx i) const final
Returns the variable at position i in the tuple.
Class hash tables iterators.
#define GIBBSKL_DEFAULT_VERBOSITY
Definition: GibbsKL_tpl.h:38
class containing all variables and methods required for Gibbssampling
Definition: gibbsOperator.h:47
GibbsBNdistance(const IBayesNet< GUM_SCALAR > &P, const IBayesNet< GUM_SCALAR > &Q)
constructor must give 2 BNs
Definition: GibbsKL_tpl.h:49
#define GIBBSKL_DEFAULT_MAXITER
Definition: GibbsKL_tpl.h:34
void updateApproximationScheme(unsigned int incr=1)
Update the scheme w.r.t the new error and increment steps.