Miind
NumericalMasterEquationCode.hpp
Go to the documentation of this file.
1 // Copyright (c) 2005 - 2011 Marc de Kamps
2 // All rights reserved.
3 //
4 // Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
5 //
6 // * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
7 // * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation
8 // and/or other materials provided with the distribution.
9 // * Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software
10 // without specific prior written permission.
11 //
12 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
13 // WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY
14 // DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
15 // USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
16 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
17 //
18 
19 // If you use this software in work leading to a scientific publication, you should include a reference there to
20 // the 'currently valid reference', which can be found at http://miind.sourceforge.net
21 #ifndef _CODE_LIBS_GEOMLIB_NUMERICALMASTEREQUATIONCODE_INCLUDE_GUARD
22 #define _CODE_LIBS_GEOMLIB_NUMERICALMASTEREQUATIONCODE_INCLUDE_GUARD
23 
24 #include "BasicDefinitions.hpp"
25 #include "CNZLCacheCode.hpp"
27 #include "GeomLibException.hpp"
28 
29 using namespace GeomLib;
30 using namespace MPILib;
31 
32 namespace
33 {
34 
35  namespace
36  {
38  }
39 
40  template <class Estimator>
42  ( Rate rate_e,
43  const AbstractOdeSystem& system,
44  const vector<typename Estimator::CoverPair>& vec_cover_pair,
45  const double y[],
46  double dydt[]
47  )
48  {
49  Number n_bins = system.NumberOfBins();
50  Index i_reset = system.IndexResetBin();
51 
52  for (Index i = 0; i < n_bins; i++){
53 
54  typename Estimator::CoverPair pair_cover = vec_cover_pair[i];
55  if (pair_cover.first._index == pair_cover.second._index) {
56  if (pair_cover.first._index >= 0 && pair_cover.first._index < static_cast<int>(n_bins)) {
57  Index i_trans = pair_cover.first._index;
58  double lower_bound = (1 - pair_cover.first._alpha);
59  dydt[system.MapPotentialToProbabilityBin(i)] += rate_e * (pair_cover.second._alpha - lower_bound) * y[system.MapPotentialToProbabilityBin(i_trans)];
60  }
61  }
62 
63  if (pair_cover.first._index < pair_cover.second._index) {
64  if (pair_cover.first._index >= 0 && pair_cover.first._index < static_cast<int>(n_bins))
65  dydt[system.MapPotentialToProbabilityBin(i)] += rate_e * pair_cover.first._alpha * y[system.MapPotentialToProbabilityBin(pair_cover.first._index)];
66 
67  if (pair_cover.second._index >= 0 && pair_cover.second._index < static_cast<int>(n_bins))
68  dydt[system.MapPotentialToProbabilityBin(i)] += rate_e * pair_cover.second._alpha * y[system.MapPotentialToProbabilityBin(pair_cover.second._index)];
69 
70  for (int j = pair_cover.first._index + 1; j < pair_cover.second._index; j++)
71  dydt[system.MapPotentialToProbabilityBin(i)] += rate_e * y[system.MapPotentialToProbabilityBin(j)];
72  }
73 
74  dydt[system.MapPotentialToProbabilityBin(i)] -= rate_e * y[system.MapPotentialToProbabilityBin(i)];
75 
76  assert(!(pair_cover.first._index > pair_cover.second._index));
77 
78  if (i == i_reset) {
79  Rate rate = 0;
80  typename Estimator::CoverPair pair = vec_cover_pair[n_bins - 1];
81  rate += rate_e * (1 - pair.second._alpha) * y[system.MapPotentialToProbabilityBin(pair.second._index)];
82 
83  for (Index ir = pair.second._index + 1; ir < n_bins; ir++)
84  rate += rate_e * y[system.MapPotentialToProbabilityBin(ir)];
85  // n_bins, does exist!
86  dydt[n_bins] += rate;
87  }
88  }
89  }
90 
91 
92  template <class Estimator>
94  ( Rate rate_i,
95  const AbstractOdeSystem& system,
96  const vector<typename Estimator::CoverPair>& vec_cover_pair,
97  const double y[],
98  double dydt[]
99  )
100  {
101  Number n_bins = system.NumberOfBins();
102 
103  for (Index i = 0; i < n_bins; i++){
104  typename Estimator::CoverPair pair_cover = vec_cover_pair[i];
105  if (pair_cover.first._index == pair_cover.second._index) {
106  if (pair_cover.first._index >= 0 && pair_cover.first._index < static_cast<int>(n_bins)) {
107  Index i_trans = pair_cover.first._index;
108  double lower_bound = (1 - pair_cover.first._alpha);
109  dydt[system.MapPotentialToProbabilityBin(i)] += rate_i * (pair_cover.second._alpha - lower_bound) * y[system.MapPotentialToProbabilityBin(i_trans)];
110  }
111  }
112 
113  if (pair_cover.first._index < pair_cover.second._index) {
114  if (pair_cover.first._index >= 0 && pair_cover.first._index < static_cast<int>(n_bins))
115  dydt[system.MapPotentialToProbabilityBin(i)] += rate_i * pair_cover.first._alpha * y[system.MapPotentialToProbabilityBin(pair_cover.first._index)];
116 
117  if (pair_cover.second._index >= 0 && pair_cover.second._index < static_cast<int>(n_bins))
118  dydt[system.MapPotentialToProbabilityBin(i)] += rate_i * pair_cover.second._alpha * y[system.MapPotentialToProbabilityBin(pair_cover.second._index)];
119 
120  for (int j = pair_cover.first._index + 1; j < pair_cover.second._index; j++)
121  dydt[system.MapPotentialToProbabilityBin(i)] += rate_i * y[system.MapPotentialToProbabilityBin(j)];
122  }
123 
124  // the first index should never be larger than the second one, unless the second one = -1, which would not require any execution
125  assert( pair_cover.second._index == -1 || (pair_cover.first._index <= pair_cover.second._index) );
126 
127  typename Estimator::CoverPair pair = vec_cover_pair[0];
128  Index i_shift = static_cast<Index>(pair.first._index);
129  assert(i_shift >= 0);
130 
131  if (i == i_shift)
132  dydt[system.MapPotentialToProbabilityBin(i)] -= (rate_i * pair.first._alpha) * y[system.MapPotentialToProbabilityBin(i)];
133 
134  if (i > static_cast<Index>(i_shift))
135  dydt[system.MapPotentialToProbabilityBin(i)] -= rate_i * y[system.MapPotentialToProbabilityBin(i)];
136 
137  }
138  }
139 
140  static inline void SetDyDtZero(double dydt[], Number n_bins)
141  {
142  // mind the equality sign
143  for (Index i = 0; i <=n_bins; i++)
144  dydt[i] = 0.0;
145  }
146 
147  template <class Estimator>
148  int CachedDeriv(double, const double y[], double dydt[], void* params)
149  {
150  MasterParameter<Estimator>* p_par = static_cast<MasterParameter<Estimator>* >(params);
151  const vector<InputParameterSet>& vec_set = *(p_par->_p_vec_set);
152 
153  const AbstractOdeSystem& system = *(p_par->_p_system);
154  SetDyDtZero(dydt,system.NumberOfBins());
155 
156  Index j = 0;
157  for(auto& set: vec_set){
158  if (set._rate_exc > 0)
159  HandleExcitatoryInput<Estimator>(set._rate_exc, system, p_par->_p_cache->List()[j].first, y, dydt);
160  if (set._rate_inh > 0)
161  HandleInhibitoryInput<Estimator>(set._rate_inh, system, p_par->_p_cache->List()[j].second, y, dydt);
162  j++;
163  }
164  return GSL_SUCCESS;
165  }
166 }
167 
168 template <class Estimator>
170 (
171  AbstractOdeSystem& system,
172  const DiffusionParameter& par_diffusion,
173  const CurrentCompensationParameter& par_current
174 ):
175  _system(system),
176  _par_diffusion(par_diffusion),
177  _par_current(par_current),
178  _estimator
179  (
180  _system.InterpretationBuffer(),
181  _system.Par()
182  ),
183  _convertor
184  (
185  _system.Par()._par_pop,
186  par_diffusion,
187  par_current,
188  _system.InterpretationBuffer()
189  ),
190  _integrator
191  (
192  MAXITERATION,
193  &_system.MassBuffer()[0],
194  _system.NumberOfBins() + 1, // TODO: review, rather hacky
195  T_STEP,
196  T_START,
197  PRECISION,
198  CachedDeriv<Estimator>
199  ),
200  _queue(TIME_QUEUE_BATCH),
201  _rate(0.0),
202  _scratch_dense(_system.NumberOfBins()),
203  _scratch_pot(_system.NumberOfBins()+1)
204 {
205 }
206 
207 template <class Estimator>
209 {
210 }
211 
212 template <class Estimator>
214 (
215  const std::vector<MPILib::Rate>& vec_rates,
216  const std::vector<MPILib::DelayedConnection>& vec_cons,
217  const std::vector<MPILib::NodeType>& vec_types
218 
219 )
220 {
221  _convertor.SortConnectionvector(vec_rates,vec_cons,vec_types);
222 }
223 
224 template <class Estimator>
226 {
227  _integrator.Parameter()._nr_bins = _system.NumberOfBins();
228  _integrator.Parameter()._p_system = &_system;
229  _integrator.Parameter()._p_estimator = &_estimator;
230  _integrator.Parameter()._i_reset = _system.IndexResetBin();
231  _integrator.Parameter()._p_vec_set = &_convertor.SolverParameter();
232  _integrator.Parameter()._p_cache = &_cache;
233 }
234 
235 template <class Estimator>
237 {
238  // existence of the extra bin is guaranteed in AbstractOdeSytem
239  MPILib::Probability p = _system.MassBuffer()[_system.NumberOfBins()];
240 
241  _system.MassBuffer()[_system.NumberOfBins()] = 0.0;
242  return p;
243 }
244 
245 template <class Estimator>
247 {
248  InitializeIntegrator();
249  // initialize caching values for the bin estimator
250 
251  _cache.Initialize(_system, _convertor.SolverParameter(), _estimator);
252 
253  MPILib::Time t_integrator = 0;
254 
255  // time is the time until the next bin is added.
256  _time_current = _integrator.CurrentTime();
257  // the integrator requires the total end time
258  t += _time_current;
259 
260  t_integrator = 0.0;
261 
262  double p = 0;
263 
264  while (t_integrator < t){
265  t_integrator = _integrator.Evolve(t);
266  p += RecaptureProbability();
267  }
268 
269 
271  prob._prob = p;
272  prob._time = t + _system.Par()._par_pop._tau_refractive;
273  _queue.push(prob);
274  _rate = p/_system.TStep();
275 
276  p = _queue.CollectAndRemove(t);
277  _system.MassBuffer()[_system.MapPotentialToProbabilityBin(_system.IndexResetBin())] += p;
278 }
279 
280 template <class Estimator>
282 {
283  return _rate;
284 }
285 
286 template <class Estimator>
288 {
289  Potential v_reset = _system.Par()._par_pop._theta;
290  Potential h = _convertor.SolverParameter()[0]._h_exc;
291  Potential v_bound = v_reset - h;
292 
293  Number n_bins = _system.NumberOfBins();
294  vector<Potential>& array_interpretation = _system.InterpretationBuffer();
295  vector<Potential>& array_mass = _system.MassBuffer();
296 
297  int i_bound = 0;
298  for (int i = static_cast<int>(n_bins - 1); i >= 0; i--)
299  if (array_interpretation[i] < v_bound){
300  i_bound = i;
301  break;
302  }
303 
304  if (i_bound == static_cast<int>(n_bins - 1)){
305  return 0.0;
306  }
307 
308  _scratch_pot[i_bound] = v_bound;
309  for (Index i = i_bound + 1; i < n_bins; i++)
310  _scratch_pot[i] = array_interpretation[i];
311 
312  double integ = 0.0;
313  _scratch_pot[n_bins] = v_reset;
314 
315  for (Index i = i_bound; i < n_bins; i++)
316  _scratch_dense[i] = array_mass[_system.MapPotentialToProbabilityBin(i)]/(_scratch_pot[i+1]-_scratch_pot[i]);
317 
318  for (Index i = i_bound; i < n_bins-1; i++)
319  integ += 0.5*(_scratch_dense[i]+_scratch_dense[i+1])*(_scratch_pot[i+1]-_scratch_pot[i]);
320  integ += 0.5*_scratch_dense[n_bins-1]*(_scratch_pot[n_bins]-_scratch_pot[n_bins-1]);
321 
322  Rate rate = _convertor.SolverParameter()[0]._rate_exc;
323 
324  return rate*integ;
325 
326 }
327 
328 template <class Estimator>
330 {
331  return std::accumulate(_system.MassBuffer().begin(),_system.MassBuffer().end(),0.0);
332 }
333 #endif
virtual void apply(MPILib::Time)
Apply the master equations over a period of time.
const NumtoolsLib::Precision PRECISION(1e-7, 0.0)
Default precision for running numerical integrations.
An auxiliary struct to help communicate with GSL C code. This parameter is passed in as void* and rec...
double Potential
int CachedDeriv(double, const double y[], double dydt[], void *params)
const CNZLCache< Estimator > * _p_cache
A geometric grid to represent population densities.
void HandleExcitatoryInput(Rate rate_e, const AbstractOdeSystem &system, const vector< typename Estimator::CoverPair > &vec_cover_pair, const double y[], double dydt[])
const MPILib::Time TIME_QUEUE_BATCH
Batch times for probability queues.
MPILib::Index IndexResetBin() const
Index of the reset bin relative to the interpretation array, i.e. constant during simulation...
When to switch to a two Poisson input approximation, and what input jump to use then.
unsigned int Number
unsigned int Index
double Time
NumericalMasterEquation(AbstractOdeSystem &, const DiffusionParameter &, const CurrentCompensationParameter &)
Parameter for setting current compensation values for the neural models that use it.
double Density
The objective is find a numerical solution for this equation This requires a numerical representation of the density We will work in the state space of a two dimensional and define a mesh there We first give two examples and then define the general procedure and given in Table for given fixed Delta g see Fig and we will denote coordinates in this dimension by a small letter $v The second dimension can be used to represent parameters as varied as and will represented by $w A strip is constructed by choosing two neighbouring points in state e g and integrating the vector field for a time $T that is assumed to be an integer multiple of a period of time Delta t
Definition: 2D.hpp:101
const MPILib::Time T_STEP
Default time step.
const MPILib::Time T_START
Standard start time of simulation.
A time stamped measure of probability.
double Probability
MPILib::Index MapPotentialToProbabilityBin(MPILib::Index i) const
Maintains the current mapping from a probability mass bin to its current poetntial bin in the interpr...
const Number MAXITERATION
Maximum number of iterations before DynamicNetwork reports a failure.
MPILib::Number NumberOfBins() const
Number of bins used in the grid representation.
const vector< InputParameterSet > * _p_vec_set
void HandleInhibitoryInput(Rate rate_i, const AbstractOdeSystem &system, const vector< typename Estimator::CoverPair > &vec_cover_pair, const double y[], double dydt[])
The objective is find a numerical solution for this equation This requires a numerical representation of the density We will work in the state space of a two dimensional and define a mesh there We first give two examples and then define the general procedure and given in Table for given fixed Delta g see Fig and we will denote coordinates in this dimension by a small letter $v The second dimension can be used to represent parameters as varied as and will represented by $w A strip is constructed by choosing two neighbouring points in state e g and integrating the vector field for a time $T that is assumed to be an integer multiple of a period of time Delta which we assume to be a defining characteristic of the grid Let then the set of points the set of points which is quadrilateral in shape The quadrilateral should be but not necessarily as long as they are but it is convenient to number them in order of creation In the we will assume that strip numbers created by the integration procedure start and are so that the numbers i in each identify a unique strip Strip no is reserved for stationary points There may or more cells in strip The number of cells in strip $i denoted by with $i the strip number and $j the cell as the i
Definition: 2D.hpp:145
const AbstractOdeSystem * _p_system
virtual void sortConnectionVector(const std::vector< MPILib::Rate > &, const std::vector< MPILib::DelayedConnection > &, const std::vector< MPILib::NodeType > &)
Interpretation of the input is done by a GeomInputConvertor object. This is a hook to provide that ob...
double Rate