// gen_events.hoc
// The gen_events function returns a vector of spike times that are
// based on a poisson firing rate created from periodic repetition of
// a gaussian probability distribution function see below for
// documentation on the arguments passed to the gen_events function:

load_file("gauss.hoc")

objref _poisson_rate
objref store_gauss_centers, time_vec, r_
objref r_
renormalization=1 // set to 0 if no renomalization
seed_ = 1 // make random stream reproducible
r_ = new Random(seed_)
r_.uniform(0,1) // random numbers between 0 and 1 will be generated by repick


obfunc poisson() { local i localobj probability, spike_times // receives 
// a poisson rate vector and returns a vector of spike times
  probability = $o1.mul(dt)
  spike_times = new Vector()
  for i=0, probability.size()-1 {
    if (probability.x[i]>r_.repick()) {
      spike_times.append(dt*i) // create a spike where randomly appropriate
    }
  }
  return spike_times
}

obfunc gen_poisson() { local total_time, period, gauss_center, gauss_half_width, gauss_peak  localobj spike_times_vec, probability
// args are tstop, T (period), gauss_center, gauss_half_width, gauss_peak
//          $1,    $2,         $3,           $4,               $5
// Here the gauss_center, gauss_half_width are in milliseconds (ms) and
// the gauss_peak is in spikes/second (Hz).
// The spike_times_vec returns the spike times in ms,
// _poisson_rate probability time series is returned and also is a global
// vector that stores the poisson firing rate at each time step

total_time = $1
// just use dt instead of delta_t = $2
period = $2
gauss_center = $3
gauss_half_width = $4
gauss_peak = $5
print "Creating stimuli's poissons"
spike_times_vec = new Vector()

num_of_stimuli = int(total_time/period+1) // each period is considered a stimuli
//  - actually a period contains a poisson rate and the rate is kind of like a stimuli

time_vec = new Vector()
time_vec.indgen(0, total_time, dt)
_poisson_rate = new Vector(time_vec.size())

for i=-1, num_of_stimuli+1 {
  mu = i*period+gauss_center
  sigma = gauss_half_width
  _poisson_rate.add(time_vec.c.apply("uni_gauss_x"))
}

if (renormalization) {
  renorm = _poisson_rate.max()

  if (renorm != 1) {
    _poisson_rate.div(renorm) // renormalize to a peak of 1 because
    print "gaussians overlapped significantly so renormalized from ",renorm," to 1"
  }                                    // if the gaussians from different 
  // periods overlap significantly their combined peaks will be higher than gauss_peak
}

_poisson_rate.mul(gauss_peak)

// _poisson_rate is now ready to create spikes
//
// loop over the time steps and generate a spike or not at each time

return _poisson_rate
}


obfunc gen_events() { local total_time, period, gauss_center, gauss_half_width, gauss_peak  localobj spike_times_vec, probability
// args are tstop, T (period), gauss_center, gauss_half_width, gauss_peak
//          $1,    $2,         $3,           $4,               $5
// Here the gauss_center, gauss_half_width are in milliseconds (ms) and
// the gauss_peak is in spikes/second (Hz).
// The spike_times_vec returns the spike times in ms,
// _poisson_rate probability time series is a temporary
// vector that stores the poisson firing rate at each time step

return poisson(gen_poisson($1, $2, $3, $4, $5).c.div(1000)) // divide by 1000 because rate is Hz not spikes/ms
}