aGrUM  0.13.2
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 <cmath>
28 
29 #include <agrum/BN/IBayesNet.h>
33 #include <agrum/core/hashTable.h>
34 
35 #define GIBBSKL_DEFAULT_MAXITER 10000000
36 #define GIBBSKL_DEFAULT_EPSILON 1e-10
37 #define GIBBSKL_DEFAULT_MIN_EPSILON_RATE 1e-10
38 #define GIBBSKL_DEFAULT_PERIOD_SIZE 200
39 #define GIBBSKL_DEFAULT_VERBOSITY false
40 #define GIBBSKL_DEFAULT_BURNIN 2000
41 #define GIBBSKL_DEFAULT_TIMEOUT 6000
42 
43 #define GIBBSKL_POURCENT_DRAWN_SAMPLE 10 // percent drawn
44 #define GIBBSKL_DRAWN_AT_RANDOM false
45 
46 namespace gum {
47 
48 
49  template < typename GUM_SCALAR >
51  const IBayesNet< GUM_SCALAR >& Q) :
52  KL< GUM_SCALAR >(P, Q),
54  GibbsOperator< GUM_SCALAR >(
55  P,
56  nullptr,
57  1 + (P.size() * GIBBSKL_POURCENT_DRAWN_SAMPLE / 100),
59  GUM_CONSTRUCTOR(GibbsKL);
60 
68  }
69 
70  template < typename GUM_SCALAR >
72  KL< GUM_SCALAR >(kl), ApproximationScheme()
73  // Gibbs operator with 10% of nodes changes at random between each samples
74  ,
75  GibbsOperator< GUM_SCALAR >(
76  kl.p(),
77  nullptr,
78  1 + (kl.p().size() * GIBBSKL_POURCENT_DRAWN_SAMPLE / 100),
79  true) {
80  GUM_CONSTRUCTOR(GibbsKL);
81 
89  }
90 
91  template < typename GUM_SCALAR >
93  GUM_DESTRUCTOR(GibbsKL);
94  }
95 
96  template < typename GUM_SCALAR >
98  auto Iq = _q.completeInstantiation();
99 
102 
103  // map between particle() variables and _q variables (using name of vars)
105 
106  for (Idx ite = 0; ite < I.nbrDim(); ++ite) {
107  map.insert(&I.variable(ite), &_q.variableFromName(I.variable(ite).name()));
108  }
109 
110  // BURN IN
111  for (Idx i = 0; i < burnIn(); i++)
112  I = this->nextSample(I);
113 
114  // SAMPLING
115  _klPQ = _klQP = _hellinger = (GUM_SCALAR)0.0;
116  _errorPQ = _errorQP = 0;
118  GUM_SCALAR delta, ratio, error;
119  delta = ratio = error = (GUM_SCALAR)-1;
120  GUM_SCALAR oldPQ = 0.0;
121  GUM_SCALAR pp, pq;
122 
123  do {
124  this->disableMinEpsilonRate();
125  I = this->nextSample(I);
127 
128  //_p.synchroInstantiations( Ip,I);
129  Iq.setValsFrom(map, I);
130 
131  pp = _p.jointProbability(I);
132  pq = _q.jointProbability(Iq);
133 
134  if (pp != (GUM_SCALAR)0.0) {
135  _hellinger += std::pow(std::sqrt(pp) - std::sqrt(pq), 2) / pp;
136 
137  if (pq != (GUM_SCALAR)0.0) {
138  _bhattacharya += std::sqrt(pq / pp); // std::sqrt(pp*pq)/pp
140  this->enableMinEpsilonRate(); // replace check_rate=true;
141  ratio = pq / pp;
142  delta = (GUM_SCALAR)log2(ratio);
143  _klPQ += delta;
144  } else {
145  _errorPQ++;
146  }
147  }
148 
149  if (pq != (GUM_SCALAR)0.0) {
150  if (pp != (GUM_SCALAR)0.0) {
151  // if we are here, it is certain that delta and ratio have been
152  // computed
153  // further lines above. (for now #112-113)
154  _klQP += (GUM_SCALAR)(-delta * ratio);
155  } else {
156  _errorQP++;
157  }
158  }
159 
160  if (this->isEnabledMinEpsilonRate()) { // replace check_rate
161  // delta is used as a temporary variable
162  delta = _klPQ / nbrIterations();
163  error = (GUM_SCALAR)std::abs(delta - oldPQ);
164  oldPQ = delta;
165  }
166  } while (continueApproximationScheme(error)); //
167 
168  _klPQ = -_klPQ / (nbrIterations());
169  _klQP = -_klQP / (nbrIterations());
170  _hellinger = std::sqrt(_hellinger / nbrIterations());
171  _bhattacharya = -std::log(_bhattacharya);
172  }
173 
174  template < typename GUM_SCALAR >
176  this->_burn_in = b;
177  }
178 
179  template < typename GUM_SCALAR >
181  return this->_burn_in;
182  }
183 } // namespace gum
GUM_SCALAR _klPQ
Definition: KL.h:132
const IBayesNet< GUM_SCALAR > & p() const
Definition: KL_tpl.h:108
unsigned long Size
In aGrUM, hashed values are unsigned long int.
Definition: types.h:50
GibbsKL(const IBayesNet< GUM_SCALAR > &P, const IBayesNet< GUM_SCALAR > &Q)
constructor must give 2 BNs
Definition: GibbsKL_tpl.h:50
This file contains Gibbs sampling (for BNs) class definitions.
bool isEnabledMinEpsilonRate() const
Returns true if stopping criterion on epsilon rate is enabled, false otherwise.
This file contains general scheme for iteratively convergent algorithms.
#define GIBBSKL_DEFAULT_TIMEOUT
Definition: GibbsKL_tpl.h:41
void disableMinEpsilonRate()
Disable stopping criterion on epsilon rate.
Idx nbrDim() const final
Returns the number of variables in the Instantiation.
#define GIBBSKL_DRAWN_AT_RANDOM
Definition: GibbsKL_tpl.h:44
Approximation Scheme.
const IBayesNet< GUM_SCALAR > & _p
Definition: KL.h:129
GibbsKL computes the KL divergence betweens 2 BNs using an approximation pattern: GIBBS sampling...
Definition: GibbsKL.h:76
void enableMinEpsilonRate()
Enable stopping criterion on epsilon rate.
#define GIBBSKL_DEFAULT_BURNIN
Definition: GibbsKL_tpl.h:40
Size burnIn() const
Returns the number of burn in.
Definition: GibbsKL_tpl.h:180
void _computeKL() final
Definition: GibbsKL_tpl.h:97
GUM_SCALAR _bhattacharya
Definition: KL.h:135
void setPeriodSize(Size p)
How many samples between two stopping is enable.
#define GIBBSKL_POURCENT_DRAWN_SAMPLE
Definition: GibbsKL_tpl.h:43
Class representing Bayesian networks.
void initApproximationScheme()
Initialise the scheme.
~GibbsKL()
destructor
Definition: GibbsKL_tpl.h:92
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)|).
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.
void setBurnIn(Size b)
Number of burn in for one iteration.
Definition: GibbsKL_tpl.h:175
bool continueApproximationScheme(double error)
Update the scheme w.r.t the new error.
const IBayesNet< GUM_SCALAR > & _q
Definition: KL.h:130
#define GIBBSKL_DEFAULT_EPSILON
Definition: GibbsKL_tpl.h:36
Instantiation nextSample(Instantiation prev)
draws next sample of Gibbs sampling
#define GIBBSKL_DEFAULT_PERIOD_SIZE
Definition: GibbsKL_tpl.h:38
KL is the base class for KL computation betweens 2 BNs.
Definition: KL.h:65
Size _errorPQ
Definition: KL.h:137
Class for assigning/browsing values to tuples of discrete variables.
Definition: instantiation.h:80
void setMaxIter(Size max)
Stopping criterion on number of iterations.
Instantiation monteCarloSample()
draws a Monte Carlo sample
algorithm for approximated computation KL divergence between BNs using GIBBS sampling ...
#define GIBBSKL_DEFAULT_MIN_EPSILON_RATE
Definition: GibbsKL_tpl.h:37
void setEpsilon(double eps)
Given that we approximate f(t), stopping criterion on |f(t+1)-f(t)|.
Size nbrIterations() const
Returns the number of iterations.
value_type & insert(const Key &key, const Val &val)
Adds a new element (actually a copy of this element) into the hash table.
unsigned long Idx
Type for indexes.
Definition: types.h:43
const std::string & name() const
returns the name of the variable
GUM_SCALAR _klQP
Definition: KL.h:133
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:39
Size _errorQP
Definition: KL.h:138
class containing all variables and methods required for Gibbssampling
Definition: gibbsOperator.h:47
#define GIBBSKL_DEFAULT_MAXITER
Definition: GibbsKL_tpl.h:35
void updateApproximationScheme(unsigned int incr=1)
Update the scheme w.r.t the new error and increment steps.
GUM_SCALAR _hellinger
Definition: KL.h:134