#*******************************************************************
#*   
#*   PCSIM Tutorial Exercises
#*
#*   Exercise 1: Small Networks Composed of Several Neurons
#*
#*     Part 1: Construct and Simulate the Model
#*
#*
#*       FIAS Theoretical Neuroscience Summer School 2008
#*   
#*       Author: Dejan Pecevski 
#*
#*******************************************************************

# We will need the pypcsim package for simulation,
# the numerical package numpy 
# and the pylab (matplotlib) package for the plotting
from pypcsim import *
from numpy import *
from pylab import *

net = SingleThreadNetwork()

# Create the LIF neuron
nrn = net.create( LifNeuron( Rm = 10e8, # in Ohms
                             Cm = 3e-10, # in Farrads
                             Vinit = -70e-3, # in Volts
                             Vresting = -70e-3, #
                             Vreset = -70e-3,
                             Vthresh = -60e-3,
                             Trefract = 5e-3 ) )

# Create the input neurons that spike Poisson processes
inp_nrns = net.create ( PoissonInputNeuron( rate = 5, duration = 100),10 )

# Connect each of the input neurons to the LIF neuron
for inp in inp_nrns:
    net.connect(inp, nrn, StaticSpikingSynapse( tau = 5e-3,
                                                        delay = 1e-3,
                                                        W = 1e-10))
    
# We want to record the spikes of the LIF Neuron and the membrane potential
# so we setup two recorders for that
spike_rec = net.record(nrn, SpikeTimeRecorder())
vm_rec = net.record(nrn, "Vm", AnalogRecorder())

# simulate the model for 1 second
net.simulate(1.0)

# Retrieve the list of spikes from the spike recorder
recorded_spikes = list(net.object(spike_rec).getRecordedValues())

# Retrieve the recorded values of the membrane potential from the analog recorder
recorded_vm = array(net.object(vm_rec).getRecordedValues())


# plot the spikes
figure(1, figsize = [8,3])
plot(recorded_spikes, zeros(len(recorded_spikes)), '|', markersize=50)
xlim(0,1)   
     
     
# plot the membrane potential
figure(2)
plot(arange(0,1,1e-4), recorded_vm)
xlim(0,1)