
# simulation class
class LIFNeurons:
  """
  Keeps track of membrane potential for multiple realizations of LIF neuron, 
  and performs single step discrete time integration.
  """
  def __init__(self, n, t_ref_mu=0.01, t_ref_sigma=0.002,
               tau=20e-3, el=-60e-3, vr=-70e-3, vth=-50e-3, r=100e6
               ):
    
    # neuron count
    self.n = n

    # neuron parameters
    self.tau = tau        # second
    self.el = el          # milivolt
    self.vr = vr          # milivolt
    self.vth = vth        # milivolt
    self.r = r            # ohm

    # initializes refractory period distribution
    self.t_ref_mu = t_ref_mu
    self.t_ref_sigma = t_ref_sigma
    self.t_ref = self.t_ref_mu + self.t_ref_sigma*np.random.normal(size=self.n)
    self.t_ref[self.t_ref<0] = 0

    # state variables
    self.v = self.el * np.ones(self.n)
    self.spiked = self.v >= self.vth
    self.last_spike = -self.t_ref * np.ones([self.n])
    self.t = 0
    self.steps = 0


  def ode_step(self, dt, i):
    
    # update running time and steps
    self.t += dt
    self.steps += 1

    # one step of discrete time integration of dt
    self.v = self.v + dt/self.tau * (self.el - self.v + self.r*i)

    # spike and clamp
    self.spiked = (self.v >= self.vth)
    self.v[self.spiked] = self.vr
    self.last_spike[self.spiked] = self.t
    clamped = (self.t_ref > self.t-self.last_spike)
    self.v[clamped] = self.vr

    self.last_spike[self.spiked] = self.t
      

# set random number generator
np.random.seed(2020)

# initialize step_end, t_range, n, v_n, syn and raster
t_range = np.arange(0, t_max, dt)
step_end = len(t_range)
n = 500
v_n = el * np.ones([n,step_end])
syn = i_mean * (1 + 0.1*(t_max/dt)**(0.5)*(2*np.random.random([n,step_end])-1))
raster = np.zeros([n,step_end])

# initialize neurons
neurons = LIFNeurons(n)

# loop time steps
for step, t in enumerate(t_range):

  neurons.ode_step(dt, syn[:,step])

  # log v_n and spike history
  v_n[:,step] = neurons.v
  raster[neurons.spiked, step] = 1.

# report running time and steps
print(f'Ran for {neurons.t:.3}s in {neurons.steps} steps.')

# plot multiple realizations of Vm, spikes and mean spike rate
with plt.xkcd():
  plot_all(t_range, v_n, raster)

# simulation class
class LIFNeurons:
  """
  Keeps track of membrane potential for multiple realizations of LIF neuron, 
  and performs single step discrete time integration.
  """
  def __init__(self, n, t_ref_mu=0.01, t_ref_sigma=0.002,
               tau=20e-3, el=-60e-3, vr=-70e-3, vth=-50e-3, r=100e6
               ):
    
    # neuron count
    self.n = n

    # neuron parameters
    self.tau = tau        # second
    self.el = el          # milivolt
    self.vr = vr          # milivolt
    self.vth = vth        # milivolt
    self.r = r            # ohm

    # initializes refractory period distribution
    self.t_ref_mu = t_ref_mu
    self.t_ref_sigma = t_ref_sigma
    self.t_ref = self.t_ref_mu + self.t_ref_sigma*np.random.normal(size=self.n)
    self.t_ref[self.t_ref<0] = 0

    # state variables
    self.v = self.el * np.ones(self.n)
    self.spiked = self.v >= self.vth
    self.last_spike = -self.t_ref * np.ones([self.n])
    self.t = 0
    self.steps = 0


  def ode_step(self, dt, i):
    
    # update running time and steps
    self.t += dt
    self.steps += 1

    # one step of discrete time integration of dt
    self.v = self.v + dt/self.tau * (self.el - self.v + self.r*i)

    # spike and clamp
    self.spiked = (self.v >= self.vth)
    self.v[self.spiked] = self.vr
    self.last_spike[self.spiked] = self.t
    clamped = (self.t_ref > self.t-self.last_spike)
    self.v[clamped] = self.vr

    self.last_spike[self.spiked] = self.t
      

# set random number generator
np.random.seed(2020)

# initialize step_end, t_range, n, v_n, syn and raster
t_range = np.arange(0, t_max, dt)
step_end = len(t_range)
n = 500
v_n = el * np.ones([n,step_end])
syn = i_mean * (1 + 0.1*(t_max/dt)**(0.5)*(2*np.random.random([n,step_end])-1))
raster = np.zeros([n,step_end])

# initialize neurons
neurons = LIFNeurons(n)

# loop time steps
for step, t in enumerate(t_range):

  neurons.ode_step(dt, syn[:,step])

  # log v_n and spike history
  v_n[:,step] = neurons.v
  raster[neurons.spiked, step] = 1.

# report running time and steps
print(f'Ran for {neurons.t:.3}s in {neurons.steps} steps.')

# plot multiple realizations of Vm, spikes and mean spike rate
with plt.xkcd():
  plot_all(t_range, v_n, raster)
# simulation class Keeps track of membrane potential for multiple realizations of LIF neuron,    and performs single step discrete time integration. # neuron count # neuron parameters # second # milivolt # milivolt # milivolt # ohm # initializes refractory period distribution # state variables # update running time and steps # one step of discrete time integration of dt # spike and clamp # set random number generator # initialize step_end, t_range, n, v_n, syn and raster # initialize neurons # loop time steps # log v_n and spike history # report running time and steps # plot multiple realizations of Vm, spikes and mean spike rate