#*******************************************************************
#*   
#*   PCSIM Tutorial Exercises
#*
#*   Exercise 2: Balanced Random Network 
#*                 with Excitatory and Inhibitory Populations
#*   
#*      
#*       FIAS Theoretical Neuroscience Summer School, August 2008
#*   
#*       Author: Dejan Pecevski 
#*
#*******************************************************************

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

Tsim = 100.0

# Create an empty network
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 = Tsim), 100 )

# Connect each of the input neurons to the LIF neuron
syns = []
for inp in inp_nrns:
    syn_id = net.connect(inp, nrn, StaticStdpSynapse( tau = 5e-3,
                                             delay = 1e-3,
                                             Winit = 1e-11,
                                             useFroemkeDanSTDP = False, 
                                             Apos = 2e-13,
                                             Aneg = -2.2e-13,
                                             taupos = 30e-3,
                                             tauneg = 30e-3,
                                             mupos = 0, 
                                             muneg = 0, 
                                             Wex = 2 * 1e-11))
    syns.append( syn_id )
    
# We want to record the spikes of the LIF Neuron 
spike_rec = net.record(nrn, SpikeTimeRecorder())

# Create a population from the list of synapses IDs
syn_popul = SimObjectPopulation( net, syns )

# use the synapse population to create population of recorders
# of the synaptic weights of all synapses
syn_rec_popul = syn_popul.record(AnalogRecorder(samplingTime = 1000), "W")

# simulate the model for 100 seconds
net.simulate(Tsim)

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

# collect all weight recordings in a two dimensional array
weights = vstack([ array(syn_rec_popul.object(i).getRecordedValues()) for i in range(syn_rec_popul.size()) ])


# plot the weight change of the first synapse
figure(1)
plot(arange(0,Tsim,1e-1), weights[0])
xlim(0,Tsim)   

# plot a histogram of the weights at the beginning of the simulation
figure(2)
subplot(2,1,1)
hist(weights[:,0])
   
# plot a histogram of the weights at the end of the simulation
subplot(2,1,2)
hist(weights[:,-1])