TemperedMarkovChain.h
Go to the documentation of this file.
1 /*!
2  *
3  *
4  * \brief -
5  *
6  * \author -
7  * \date -
8  *
9  *
10  * \par Copyright 1995-2015 Shark Development Team
11  *
12  * <BR><HR>
13  * This file is part of Shark.
14  * <http://image.diku.dk/shark/>
15  *
16  * Shark is free software: you can redistribute it and/or modify
17  * it under the terms of the GNU Lesser General Public License as published
18  * by the Free Software Foundation, either version 3 of the License, or
19  * (at your option) any later version.
20  *
21  * Shark is distributed in the hope that it will be useful,
22  * but WITHOUT ANY WARRANTY; without even the implied warranty of
23  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24  * GNU Lesser General Public License for more details.
25  *
26  * You should have received a copy of the GNU Lesser General Public License
27  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
28  *
29  */
30 #ifndef SHARK_UNSUPERVISED_RBM_SAMPLING_TEMPEREDMARKOVCHAIN_H
31 #define SHARK_UNSUPERVISED_RBM_SAMPLING_TEMPEREDMARKOVCHAIN_H
32 
33 #include <shark/Data/Dataset.h>
36 #include <vector>
37 #include "Impl/SampleTypes.h"
38 namespace shark{
39 
40 
41 //\brief models a set of tempered Markov chains given a TransitionOperator.
42 // e.g. TemperedMarkovChain<GibbsOperator<RBM> > chain, leads to the set of chains
43 // used for parallel tempering.
44 template<class Operator>
46 private:
47  typedef typename Operator::HiddenSample HiddenSample;
48  typedef typename Operator::VisibleSample VisibleSample;
49 public:
50 
51  ///\brief The MarkovChain can't be used to compute several samples at once.
52  ///
53  /// The tempered markov chain ues it's batch capabilities allready to compute the samples for all temperatures
54  /// At the same time. Also it is much more powerfull when all samples are drawn one after another for a higher mixing rate.
55  static const bool computesBatch = false;
56 
57  ///\brief The type of the RBM the operator is working with.
58  typedef typename Operator::RBM RBM;
59 
60  ///\brief A batch of samples containing hidden and visible samples as well as the energies.
62 
63  ///\brief Mutable reference to an element of the batch.
64  typedef typename SampleBatch::reference reference;
65 
66  ///\brief Immutable reference to an element of the batch.
67  typedef typename SampleBatch::const_reference const_reference;
68 
69 private:
70  SampleBatch m_temperedChains;
71  RealVector m_betas;
72  Operator m_operator;
73 
74  void metropolisSwap(reference low, double betaLow, reference high, double betaHigh){
75  RealVector const& baseRate = transitionOperator().rbm()->visibleNeurons().baseRate();
76  double betaDiff = betaLow - betaHigh;
77  double energyDiff = low.energy - high.energy;
78  double baseRateDiff = inner_prod(low.visible.state,baseRate) - inner_prod(high.visible.state,baseRate);
79  double r = betaDiff * energyDiff + betaDiff*baseRateDiff;
80 
81 
82  Uniform<typename RBM::RngType> uni(m_operator.rbm()->rng(),0,1);
83  double z = uni();
84  if( r >= 0 || (z > 0 && std::log(z) < r) ){
85  swap(high,low);
86  }
87  }
88 
89 public:
90  TemperedMarkovChain(RBM* rbm):m_operator(rbm){}
91 
92  const Operator& transitionOperator()const{
93  return m_operator;
94  }
95  Operator& transitionOperator(){
96  return m_operator;
97  }
98 
99 
100  /// \brief Sets the number of temperatures and initializes the tempered chains accordingly.
101  ///
102  /// @param temperatures number of temperatures
103  void setNumberOfTemperatures(std::size_t temperatures){
104  std::size_t visibles=m_operator.rbm()->numberOfVN();
105  std::size_t hiddens=m_operator.rbm()->numberOfHN();
106  m_temperedChains = SampleBatch(temperatures,visibles,hiddens);
107  m_betas.resize(temperatures);
108  }
109 
110  /// \brief Sets the number of temperatures and initializes them in a uniform spacing
111  ///
112  /// Temperatures are spaced equally between 0 and 1.
113  /// @param temperatures number of temperatures
114  void setUniformTemperatureSpacing(std::size_t temperatures){
115  setNumberOfTemperatures(temperatures);
116  for(std::size_t i = 0; i != temperatures; ++i){
117  double factor = temperatures - 1.0;
118  setBeta(i,1.0 - i/factor);
119  }
120  }
121 
122 
123  /// \brief Returns the number Of temperatures.
124  std::size_t numberOfTemperatures()const{
125  return m_betas.size();
126  }
127 
128  void setBatchSize(std::size_t batchSize){
129  SHARK_CHECK(batchSize == 1, "[TemperedMarkovChain::setBatchSize] markov chain can only compute batches of size 1.");
130  }
131  std::size_t batchSize(){
132  return 1;
133  }
134 
135  void setBeta(std::size_t i, double beta){
136  SIZE_CHECK(i < m_betas.size());
137  m_betas(i) = beta;
138  }
139 
140  double beta(std::size_t i)const{
141  SIZE_CHECK(i < m_betas.size());
142  return m_betas(i);
143  }
144 
145  RealVector const& beta()const{
146  return m_betas;
147  }
148 
149  ///\brief Returns the current state of the chain for beta = 1.
150  const_reference sample()const{
151  return const_reference(m_temperedChains,0);
152  }
153  ///\brief Returns the current state of the chain for all beta values.
154  SampleBatch const& samples()const{
155  return m_temperedChains;
156  }
157 
158  /// \brief Returns the current batch of samples of the Markov chain.
159  SampleBatch& samples(){
160  return m_temperedChains;
161  }
162 
163  ///\brief Initializes the markov chain using samples drawn uniformly from the set.
164  ///
165  /// Be aware that the number of chains and the temperatures need to bee specified previously.
166  /// @param dataSet the data set
167  void initializeChain(Data<RealVector> const& dataSet){
168  if(m_temperedChains.size()==0)
169  throw SHARKEXCEPTION("you did not initialize the number of temperatures bevor initializing the chain!");
170  DiscreteUniform<typename RBM::RngType> uni(m_operator.rbm()->rng(),0,dataSet.numberOfElements()-1);
171  std::size_t visibles = m_operator.rbm()->numberOfVN();
172  RealMatrix sampleData(m_temperedChains.size(),visibles);
173 
174  for(std::size_t i = 0; i != m_temperedChains.size(); ++i){
175  noalias(row(sampleData,i)) = dataSet.element(uni());
176  }
177  initializeChain(sampleData);
178  }
179 
180  /// \brief Initializes with data points from a batch of points
181  ///
182  /// Be aware that the number of chains and the temperatures need to bee specified previously.
183  /// @param sampleData the data set
184  void initializeChain(RealMatrix const& sampleData){
185  if(m_temperedChains.size()==0)
186  throw SHARKEXCEPTION("you did not initialize the number of temperatures bevor initializing the chain!");
187 
188  m_operator.createSample(m_temperedChains.hidden,m_temperedChains.visible,sampleData,m_betas);
189 
190  m_temperedChains.energy = m_operator.calculateEnergy(
191  m_temperedChains.hidden,
192  m_temperedChains.visible
193  );
194  }
195  //updates the chain using the current sample
196  void step(unsigned int k){
197  for(std::size_t i = 0; i != k; ++i){
198  //do one step of the tempered the Markov chains at the same time
199  m_operator.stepVH(m_temperedChains.hidden, m_temperedChains.visible,1,m_betas);
200 
201  //calculate energy for samples at all temperatures
202  m_temperedChains.energy = m_operator.calculateEnergy(
203  m_temperedChains.hidden,
204  m_temperedChains.visible
205  );
206 
207  //EVEN phase
208  std::size_t elems = m_temperedChains.size();
209  for(std::size_t i = 0; i < elems-1; i+=2){
210  metropolisSwap(
211  reference(m_temperedChains,i),m_betas(i),
212  reference(m_temperedChains,i+1),m_betas(i+1)
213  );
214  }
215  //ODD phase
216  for(std::size_t i = 1; i < elems-1; i+=2){
217  metropolisSwap(
218  reference(m_temperedChains,i),m_betas(i),
219  reference(m_temperedChains,i+1),m_betas(i+1)
220  );
221  }
222  m_operator.rbm()->hiddenNeurons().sufficientStatistics(
223  m_temperedChains.hidden.input,m_temperedChains.hidden.statistics, m_betas
224  );
225  }
226  }
227 };
228 
229 }
230 #endif