#*******************************************************************
#*   
#*   PCSIM Tutorial Exercises
#*
#*   Exercise 2: Balanced Random Network 
#*                 with Excitatory and Inhibitory Populations
#*   
#*   Part 3: Spatial Networks and  Distance-based random 
#*           Connections
#*      
#*       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 *

###################################################################
# Global parameter values
###################################################################
nNeurons        = 1000;   # number of neurons
ConnP           = 0.1;    # connectivity probability
Frac_EXC        = 0.8;    # fraction of excitatory neurons
nInputNeurons   = 100     # number of neurons which provide initial input (for a time span of Tinp)
inpConnP        = 0.1     # connectivity probability from input neurons to network neurons
inputFiringRate = 10      # firing rate of the input neurons during the initial input [spikes/sec]
Tsim            = 1.0     # simulation time 

# Create an empty pcsim network
net = SingleThreadNetwork()


###################################################################
# CREATING THE NEURONS
###################################################################
# Create the LIF neuron model used to generate
# the neurons in the populations
lifmodel = LifNeuron( Cm=2e-10, 
                      Rm=1e8, 
                      Vthresh=-50e-3, 
                      Vresting=-49e-3, 
                      Vreset=-60e-3, 
                      Trefract=5e-3, 
                      Vinit=-60e-3 ) ;

# create a variation factory based on the model
nrn_factory = SimObjectVariationFactory( lifmodel )

# associate the initial membrane potential parameter with a normal distribution
nrn_factory.set( 'Vinit', BndNormalDistribution(mu = -55e-3, 
                                                cv = 0.1, 
                                                lowerBound = -60e-3, 
                                                upperBound = -50e-3) )

# Create a spatial population composed of two families of neurons (excitatory and inhibitory)
# with ratio 4:1. Neurons are located at integer 3d grid coordinates within a volume 20x10x5.  
all_nrn_popul = SpatialFamilyPopulation(net, [ nrn_factory, nrn_factory], 
                                            RatioBasedFamilies((4,1)),
                                            CuboidIntegerGrid3D(20, 10, 5));

# split the population in two sub-populations, one for each family
exc_nrn_popul, inh_nrn_popul = all_nrn_popul.splitFamilies()

print "Created", exc_nrn_popul.size(), "exc and", inh_nrn_popul.size(), "inh neurons";

###################################################################
# SYNAPTIC CONNECTIONS
###################################################################
# create the excitatory synapse model
exc_syn_model = StaticSpikingSynapse( W=1.62e-11, tau=5e-3, delay=1e-3 )

# create the variation factory based on the model
exc_syn_factory = SimObjectVariationFactory( exc_syn_model )

# associate a gamma distribution with the weight parameter W
exc_syn_factory.set( "W", BndGammaDistribution(mu = 1.62e-11, cv = 0.1, 
                                               upperBound = 2 * 1.62e-11 ) )


print 'Making synaptic connections:'
# randomly connect all the excitatory neurons to all neurons in the network
# with a probability dependent on the distance between the neurons
exc_syn_project = ConnectionsProjection( exc_nrn_popul, all_nrn_popul, 
                                     exc_syn_factory,                                               
                                     EuclideanDistanceRandomConnections( ConnP, 10.0 ))

# randomly connect all the inhibitory neurons to all neurons in the network
# with a probability dependent on the distance between the neurons
inh_syn_project = ConnectionsProjection( inh_nrn_popul, all_nrn_popul,
                                     StaticSpikingSynapse( W=-10e-11, tau=10e-3, delay=1e-3 ), 
                                     EuclideanDistanceRandomConnections( ConnP, 10.0 ))

print " Number of excitatory synapses ", exc_syn_project.size()

print " Number of inhibitory synapses ", inh_syn_project.size()

###################################################################
#  INPUTS
###################################################################
inp_nrn_popul = SimObjectPopulation( net, PoissonInputNeuron( rate = 10, duration = 10), 10) 

inp_syn_project = ConnectionsProjection( inp_nrn_popul, exc_nrn_popul, 
                          StaticSpikingSynapse( W=1.2e-10, tau=5e-3, delay=1e-3 ), 
                          RandomConnections( conn_prob = inpConnP ) );


###################################################################
#  RECORDERS
###################################################################
# population of analog recorders for recording the membrane potential 
# of all the neurons
spike_rec_popul = all_nrn_popul.record(SpikeTimeRecorder())

# population of analog recorders for recording the membrane potential 
# of all the neurons
vm_rec_popul = all_nrn_popul.record( AnalogRecorder(), "Vm" )

###########################################################
# Simulate the network
###########################################################
print "Starting simulation..."
net.simulate( Tsim )
print "Done"


###########################################################
# PRESENTATION OF THE RESULTS
###########################################################

# collect all the spikes from the neurons into a list of numpy arrays
spikes = [ array(spike_rec_popul.object(i).getRecordedValues()) for i in range(spike_rec_popul.size()) ]


# plot the spike raster of the network activity
x, y = create_raster(spikes, 0, Tsim)
plot(x, y, '.')

# Necessary functions to calculate the following figures
def diff(x):
    return [ x[i+1] - x[i] for i in range(len(x)-1) ]


def meanisi(spikes):
    if( len(spikes) > 1):
        return mean(diff(spikes));
    else:
        return None

def cv(spikes):
    if( len(spikes) > 2):
        isi=diff(spikes);
        return std(isi) / mean(isi)
    else:
        return None

    
# Calculate the mean firing rate of the network activity
n = [ len(s) for s in spikes ]
print "spikes total", sum(n)
print "mean rate:", sum(n)/Tsim/len(n)

# Calculate the mean inter-spike interval
isi = [ meanisi(s) for s in spikes ]
isi = filter( lambda x: x != None, isi);
print "mean ISI:", sum(isi)/len(isi)

# Calculate the coefficient of variation
CV = [ cv(s) for s in spikes ]
CV = filter( lambda x: x != None, CV);
print "mean CV:", mean(CV)