21 #ifndef _CODE_LIBS_GEOMLIB_NUMERICALMASTEREQUATIONCODE_INCLUDE_GUARD
22 #define _CODE_LIBS_GEOMLIB_NUMERICALMASTEREQUATIONCODE_INCLUDE_GUARD
40 template <
class Estimator>
44 const vector<typename Estimator::CoverPair>& vec_cover_pair,
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);
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))
67 if (pair_cover.second._index >= 0 && pair_cover.second._index < static_cast<int>(n_bins))
70 for (
int j = pair_cover.first._index + 1; j < pair_cover.second._index; j++)
76 assert(!(pair_cover.first._index > pair_cover.second._index));
80 typename Estimator::CoverPair pair = vec_cover_pair[n_bins - 1];
83 for (
Index ir = pair.second._index + 1; ir < n_bins; ir++)
92 template <
class Estimator>
96 const vector<typename Estimator::CoverPair>& vec_cover_pair,
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);
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))
117 if (pair_cover.second._index >= 0 && pair_cover.second._index < static_cast<int>(n_bins))
120 for (
int j = pair_cover.first._index + 1; j < pair_cover.second._index; j++)
125 assert( pair_cover.second._index == -1 || (pair_cover.first._index <= pair_cover.second._index) );
127 typename Estimator::CoverPair pair = vec_cover_pair[0];
128 Index i_shift =
static_cast<Index>(pair.first._index);
129 assert(i_shift >= 0);
134 if (
i > static_cast<Index>(i_shift))
140 static inline void SetDyDtZero(
double dydt[],
Number n_bins)
147 template <
class Estimator>
148 int CachedDeriv(
double,
const double y[],
double dydt[],
void* params)
151 const vector<InputParameterSet>& vec_set = *(p_par->
_p_vec_set);
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);
168 template <
class Estimator>
176 _par_diffusion(par_diffusion),
177 _par_current(par_current),
180 _system.InterpretationBuffer(),
185 _system.Par()._par_pop,
188 _system.InterpretationBuffer()
193 &_system.MassBuffer()[0],
194 _system.NumberOfBins() + 1,
198 CachedDeriv<Estimator>
202 _scratch_dense(_system.NumberOfBins()),
203 _scratch_pot(_system.NumberOfBins()+1)
207 template <
class Estimator>
212 template <
class Estimator>
215 const std::vector<MPILib::Rate>& vec_rates,
216 const std::vector<MPILib::DelayedConnection>& vec_cons,
217 const std::vector<MPILib::NodeType>& vec_types
221 _convertor.SortConnectionvector(vec_rates,vec_cons,vec_types);
224 template <
class Estimator>
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;
235 template <
class Estimator>
241 _system.MassBuffer()[_system.NumberOfBins()] = 0.0;
245 template <
class Estimator>
248 InitializeIntegrator();
251 _cache.Initialize(_system, _convertor.SolverParameter(), _estimator);
256 _time_current = _integrator.CurrentTime();
264 while (t_integrator < t){
265 t_integrator = _integrator.Evolve(t);
266 p += RecaptureProbability();
272 prob.
_time = t + _system.Par()._par_pop._tau_refractive;
274 _rate = p/_system.TStep();
276 p = _queue.CollectAndRemove(t);
277 _system.MassBuffer()[_system.MapPotentialToProbabilityBin(_system.IndexResetBin())] += p;
280 template <
class Estimator>
286 template <
class Estimator>
289 Potential v_reset = _system.Par()._par_pop._theta;
290 Potential h = _convertor.SolverParameter()[0]._h_exc;
293 Number n_bins = _system.NumberOfBins();
294 vector<Potential>& array_interpretation = _system.InterpretationBuffer();
295 vector<Potential>& array_mass = _system.MassBuffer();
298 for (
int i = static_cast<int>(n_bins - 1);
i >= 0;
i--)
299 if (array_interpretation[
i] < v_bound){
304 if (i_bound == static_cast<int>(n_bins - 1)){
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];
313 _scratch_pot[n_bins] = v_reset;
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]);
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]);
322 Rate rate = _convertor.SolverParameter()[0]._rate_exc;
328 template <
class Estimator>
331 return std::accumulate(_system.MassBuffer().begin(),_system.MassBuffer().end(),0.0);
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...
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.
NumericalMasterEquation(AbstractOdeSystem &, const DiffusionParameter &, const CurrentCompensationParameter &)
Parameter for setting current compensation values for the neural models that use it.
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
virtual Rate getTransitionRate() const
const MPILib::Time T_STEP
Default time step.
Rate IntegralRate() const
const MPILib::Time T_START
Standard start time of simulation.
A time stamped measure of 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
Density RecaptureProbability()
virtual ~NumericalMasterEquation()
Destructor.
void InitializeIntegrator()
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
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...