File: HistogramsAverageAps.h

package info (click to toggle)
cain 1.10+dfsg-2
  • links: PTS, VCS
  • area: main
  • in suites: stretch
  • size: 29,856 kB
  • sloc: cpp: 49,612; python: 14,988; xml: 11,654; ansic: 3,644; makefile: 133; sh: 2
file content (231 lines) | stat: -rw-r--r-- 7,708 bytes parent folder | download | duplicates (4)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
// -*- C++ -*-

/*!
  \file stochastic/HistogramsAverageAps.h
  \brief Collect the average values of species populations. Use the all possible states method.
*/

#if !defined(__stochastic_HistogramsAverageAps_h__)
#define __stochastic_HistogramsAverageAps_h__

#include "HistogramsAverageBase.h"

namespace stochastic {

//! Accumulate a histogram of equilibrium state using Gillespie's direct method.
/*!
  \param _DiscreteGenerator Random deviate generator for the discrete,
  finite distribution with reaction propensities as scaled probabilities.
  \param _ExponentialGenerator Random deviate generator for the exponential
  distribution. By default the ziggurat algorithm is used.
  \param _PropensitiesFunctor Can calculate propensities as a function of the
  reaction index and the populations.
*/
template < class _DiscreteGenerator,
         class _ExponentialGenerator =
         numerical::ExponentialGeneratorZiggurat<>,
         class _PropensitiesFunctor = PropensitiesSingle<true> >
class HistogramsAverageAps :
   public HistogramsAverageBase < _DiscreteGenerator, _ExponentialGenerator,
      _PropensitiesFunctor > {
   //
   // Private types.
   //
private:

   typedef HistogramsAverageBase < _DiscreteGenerator, _ExponentialGenerator,
           _PropensitiesFunctor > Base;

   //
   // Public types.
   //
public:

   //! The propensities functor.
   typedef typename Base::PropensitiesFunctor PropensitiesFunctor;
   //! The exponential generator.
   typedef typename Base::ExponentialGenerator ExponentialGenerator;
   //! The discrete, finite generator.
   typedef typename Base::DiscreteGenerator DiscreteGenerator;
   //! The discrete, uniform generator.
   typedef typename Base::DiscreteUniformGenerator DiscreteUniformGenerator;

   //
   // Member data.
   //
protected:

   using Base::_state;
   using Base::_time;
   using Base::_tau;
   using Base::_exponentialGenerator;
   using Base::_discreteGenerator;
   using Base::_recordedSpecies;
   using Base::_histograms;
   using Base::_reactionInfluence;

   //
   // Not implemented.
   //
private:

   //! Default constructor not implemented.
   HistogramsAverageAps();
   //! Copy constructor not implemented.
   HistogramsAverageAps(const HistogramsAverageAps&);
   //! Assignment operator not implemented.
   HistogramsAverageAps&
   operator=(const HistogramsAverageAps&);

   //--------------------------------------------------------------------------
   //! \name Constructors etc.
   //@{
public:

   //! Construct.
   HistogramsAverageAps
   (const State& state,
    const PropensitiesFunctor& propensitiesFunctor,
    const array::StaticArrayOfArrays<std::size_t>& reactionInfluence,
    const std::vector<std::size_t>& recordedSpecies,
    const std::size_t numberOfBins,
    const std::size_t multiplicity,
    const double maxSteps) :
      Base(state, propensitiesFunctor, reactionInfluence, recordedSpecies,
           numberOfBins, multiplicity, maxSteps) {
   }

   // Use the default destructor.

   //@}
   //--------------------------------------------------------------------------
   //! \name Simulation.
   //@{
public:

   //! Initialize the state with the initial populations and time.
   using Base::initialize;

   //! Generate a trajectory and record the state in the histograms.
   void
   simulate(const double equilibrationTime, const double recordingTime) {
      // Let the process equilibrate.
      Base::equilibrate(equilibrationTime);
      // Step until we have reached the end time.
      const double endTime = _time + recordingTime;
      while (true) {
         // Check that we have not exceeded the allowed number of steps.
         if (! Base::incrementStepCount()) {
            setStepCountError();
            break;
         }

         // Compute an exponential deviate and the time step. Note that we don't
         // use computeTau() because we need the exponential deviate and not just
         // the time step.
         const double exponentialDeviate = _exponentialGenerator();
         if (_discreteGenerator.isValid()) {
            const double mean = 1.0 / _discreteGenerator.sum();
            _time.updateEpoch(mean);
            _tau = mean * exponentialDeviate;
         }
         else {
            _tau = std::numeric_limits<double>::max();
         }

         // If we have not passed the end time.
         if (_time + _tau < endTime) {
            // Increment the time.
            _time += _tau;
         }
         else {
            // Reduce the time step and indicate that this is the last step.
            _tau = endTime - _time;
            _time = endTime;
         }

         // For each reaction.
         for (std::size_t reactionIndex = 0;
               reactionIndex != _state.getNumberOfReactions(); ++reactionIndex) {
            // Skip reactions that have zero probability.
            if (_discreteGenerator[reactionIndex] == 0) {
               continue;
            }
            // This will hold the sum of the propensities after firing the reaction.
            double sum = _discreteGenerator.sum();
            // Fire the reaction.
            _state.fireReaction(reactionIndex);
            // Compute the sum of the propensities after firing the reaction.
            for (typename array::StaticArrayOfArrays<std::size_t>::const_iterator
                  i = _reactionInfluence.begin(reactionIndex);
                  i != _reactionInfluence.end(reactionIndex); ++i) {
               sum += this->_propensitiesFunctor(*i, _state.getPopulations()) -
                      _discreteGenerator[*i];
            }
            // Record the probabilities for this state.
            // P(reaction) == _discreteGenerator[reactionIndex] /
            //                _discreteGenerator.sum()
            // Time to hold state == exponentialDeviate / sum
            // The weight is the product of these two.
            const double weight = _discreteGenerator[reactionIndex] *
                                  exponentialDeviate / (_discreteGenerator.sum() * sum);
            for (std::size_t i = 0; i != _recordedSpecies.size(); ++i) {
               _histograms(i).accumulate(_state.getPopulation(_recordedSpecies[i]),
                                         weight);
            }
            // Un-fire the reaction.
            _state.unFireReaction(reactionIndex);
         }

         // If we have reached the end time.
         if (_time >= endTime) {
            // End the simulation.
            return;
         }

         // Determine the reaction to fire.
         const std::size_t reactionIndex = _discreteGenerator();
#ifdef DEBUG_stlib
         assert(_discreteGenerator[reactionIndex] > 0);
#endif
         // Fire the reaction.
         _state.fireReaction(reactionIndex);
         // Recompute the propensities and update the discrete, finite generator.
         Base::updatePropensities(reactionIndex);
      }
   }

   //! Synchronize the two sets of histograms so that corresponding historams have the same lower bounds and widths.
   using Base::synchronize;

protected:

   //! Record a step count error message. Record the current time.
   using Base::setStepCountError;

   //@}
   //--------------------------------------------------------------------------
   //! \name Accessors.
   //@{
public:

   //! Return a const reference to the state.
   using Base::getState;

   //! Return a const reference to the discrete, uniform generator.
   using Base::getDiscreteUniformGenerator;

   //! Return the vector of recorded species.
   using Base::getRecordedSpecies;

   //! Return the set of histograms.
   using Base::getHistograms;

   //@}
};

//@}

} // namespace stochastic

#endif