Information transmission in cerebellar granule cell models (Rossert et al. 2014)

 Download zip file 
Help downloading and running models
Accession:156733
" ... In this modeling study we analyse how electrophysiological granule cell properties and spike sampling influence information coded by firing rate modulation, assuming no signal-related, i.e., uncorrelated inhibitory feedback (open-loop mode). A detailed one-compartment granule cell model was excited in simulation by either direct current or mossy-fiber synaptic inputs. Vestibular signals were represented as tonic inputs to the flocculus modulated at frequencies up to 20 Hz (approximate upper frequency limit of vestibular-ocular reflex, VOR). Model outputs were assessed using estimates of both the transfer function, and the fidelity of input-signal reconstruction measured as variance-accounted-for. The detailed granule cell model with realistic mossy-fiber synaptic inputs could transmit infoarmation faithfully and linearly in the frequency range of the vestibular-ocular reflex. ... "
Reference:
1 . Rössert C, Solinas S, D'Angelo E, Dean P, Porrill J (2014) Model cerebellar granule cells can faithfully transmit modulated firing rate signals. Front Cell Neurosci 8:304 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell; Synapse;
Brain Region(s)/Organism: Cerebellum;
Cell Type(s): Cerebellum interneuron granule GLU cell;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON; Python;
Model Concept(s): Action Potentials; Markov-type model;
Implementer(s): Solinas, Sergio [solinas at unipv.it]; Roessert, Christian [christian.a at roessert.de];
Search NeuronDB for information about:  Cerebellum interneuron granule GLU cell;
# -*- coding: utf-8 -*-
"""
@author: chris

Modified from THOMAS MCTAVISH (2010-11-04).

mpiexec -f ~/machinefile -enable-x -n 96 python Population.py --noplot
"""

from __future__ import with_statement
from __future__ import division

import sys
sys.path.append('../NET/sheff/weasel/')
sys.path.append('../NET/sheffprk/template/')

import os

#use_pc = True

import sys
argv = sys.argv

if "-python" in argv:
    use_pc = True
else:
    use_pc = False  

if use_pc == True:
    
    from neuron import h
    pc = h.ParallelContext()
    rank = int(pc.id())
    nhost = pc.nhost()
        
else:
    
    from mpi4py import MPI
    from neuron import h
    rank = MPI.COMM_WORLD.rank

#print sys.version

if __name__ == '__main__':
    
    import argparse
    parser = argparse.ArgumentParser()
    parser.add_argument('-o', action='store', dest='opt')
    parser.add_argument('--noplot', action='store_true')
    parser.add_argument('--norun', action='store_true')
    parser.add_argument('--noconst', action='store_true')
    parser.add_argument('--noqual', action='store_true')
    pars, unknown = parser.parse_known_args(['-o','--noplot','--norun','--noconst','--noqual'])

if __name__ == '__main__':
    
    import matplotlib
    if rank == 0: 
        matplotlib.use('Tkagg', warn=True)    
    else:    
        matplotlib.use('Agg', warn=True) 

if __name__ == '__main__':
    
    do_plot = 1
    if results.noplot:  # do not plot to windows
        matplotlib.use('Agg', warn=True)
        if rank == 0: print "- No plotting"
        do_plot = 0

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab

import random as rnd
import neuronpy.util.spiketrain

#set_printoptions(threshold='nan')

from Stimulation import *
from Stimhelp import *
from units import *

from cells.PassiveCell import *

from itertools import izip

try:
    import cPickle as pickle
except:
    import pickle
    
import gzip
import h5py

from templates.synapse.synapse import Synapse
from synapsepfpurk import Synapse as Synapse2
if use_pc is False: import mdp
                                       
import time as ttime
from scipy.optimize import fmin, leastsq

from NeuroTools import stgen, signals

import md5

#from guppy import hpy
#hpy = hpy()
    

class Population:
    """
    A population of N cells
    """
    
    def __init__(self, cellimport = [], celltype = None, N = [10], temperature = 6.3, cell_exe = 0, ihold = [0*nA], ihold_sigma = [0*nA], amp = [0*nA], amod = [None], anoise = [None], give_freq = False, do_run = 1, pickle_prefix = "default", istart = 0, istop = 0.07, di = 0.001, dt = 0.025*ms, use_mpi = True, use_pc = False):
        """
        :param N: Number of cells.
        :param fluct_m: 
        :param fluct_s: 
        :param fluct_tau: 
        """

        self.use_pc = use_pc
        
        if type(celltype) is not list: celltype = [celltype] #convert to list if it is not given as one
        self.celltype = celltype
        
        if type(cell_exe) is not list: cell_exe = [cell_exe] #convert to list if it is not given as one
        self.cell_exe = cell_exe  
        
        if cellimport is not None:
            if cellimport == []: 
                for n in range(len(celltype)):
                    cellimport.append("from cells." + self.celltype[n] + " import *")
        self.cellimport = cellimport
        
        if type(N) is not list: N = [N]
        self.N = N # Total number of cells in the net
        
        self.n_celltypes = len(self.N)
        self.a_celltype = [0]     # celltype to analyse
        
        self.factor_celltype = [1]*self.n_celltypes
        
        self.set_init(ihold, ihold_sigma, amp, amod)
        
        self.CF_var = False

        self.inh_hold_sigma = [0]
        self.intr_hold_sigma = [0]
        
        #self.sigma_inh_hold = 0
        #self.sigma_ihold = 0
        
 
        if type(anoise) is not list: anoise = [anoise]*self.n_celltypes
        if len(anoise) < self.n_celltypes: anoise = [anoise[0]]*self.n_celltypes
        self.anoise = anoise # RUN self.set_i()
        
        self.give_freq = give_freq # RUN self.set_i()
       
        self.temperature = temperature
        
        self.gid_count = 0
        self.gidlist = []           # List of global identifiers on this host
        self.global_gidlist = []    # List of global identifiers
        self.cells = []             # Cells on this host
                
        self.t_vec = []
        self.id_vec = []
        self.rec_v = []
                
        for n in range(self.n_celltypes):
            if use_mpi:
                self.t_vec.append(h.Vector()) # np.array([0])
                self.id_vec.append(h.Vector()) # np.array([-1], dtype=int)
            else:
                self.t_vec.append([])
                
            self.rec_v.append(h.Vector())
            
        #self.t_vec = h.Vector(np.array([0]))     # Spike time of all cells on this host
        #self.id_vec = h.Vector(np.array([-1]))    # Ids of spike times on this host
        
        self.flucts = []            # Fluctuating inputs on this host
        self.fluct_m = 0            # [nA]
        self.fluct_s = [0]          # [nA]
        self.fluct_tau = 0*ms          # [ms]
        
        self.noises = []            # Random number generators on this host
        self.plays = []             # Play inputs on this host
        self.rec_is = []
        
        self.trains = []
        self.vecstim = []
        self.nc_vecstim = []
        self.spike_vec = []       
        
        self.syn_tau1 = 5*ms        # Synapse of virtual target neuron
        self.syn_tau2 = 5*ms        # Synapse of virtual target neuron
        self.tmax = 10*sec          # maximum length of plot that should be plotted!!
        
        self.nc_delay = 0 #500*ms  # only important if syn_output is used, not used currently
        self.dt = dt
        self.bin_width = dt
        self.jitter = 0*ms
        self.delta_t = 0*ms
        
        self.istart = istart
        self.istop = istop
        self.di = di
        
        self.ic_holds = []
        self.i_holdrs = []
        self.i_holds = []
        self.ic_starts = [] 
        self.vc_starts = []
        self.ic_steps = []
        
        self.rec_step = []
        
        self.tvecs = []
        self.ivecs = []  

        self.noises = []
        
        self.record_syn = []
        self.id_all_vec_input = []
        self.t_all_vec_input = []
        
        if len(self.N) == len(self.cell_exe) == len(self.celltype):
            pass
        else:
            raise ValueError('N, cell_exe, celltype do NOT have equal length!')

        self.use_mpi = use_mpi
        self.use_pc = use_pc
        
        if self.use_mpi:
            
            #### Make a new ParallelContext object
            self.pc = h.ParallelContext()
            self.id = self.pc.id()
            self.nhost = int(self.pc.nhost())
            
            if self.use_pc == False:

                s = "mpi4py thinks I am %d of %d on %s, NEURON thinks I am %d of %d\n"
                processorname = MPI.Get_processor_name()
                self.comm = MPI.COMM_WORLD
                
                if self.id == 0:
                    print s % (self.comm.rank, self.comm.size, processorname, self.id, self.nhost)
                    
            else:
               
                s = "NEURON thinks I am %d of %d\n"
                if self.id == 0:
                    print s % (self.id, self.nhost)
               
            self.barrier()
            
        else:
            self.id = 0
            self.nhost = 1
            
        self.do_run = do_run

        self.first_run = True
                
        self.set_numcells()  # Build the portion of cells on this host.  
        
        self.pickle_prefix = pickle_prefix
        
        # plot options
        self.ymax = 0 
        self.ax = None 
        self.linewidth = 1.5
        self.color_vec = None 
        self.alpha = 0.8   
        self.method_interpol = np.array(['bin','syn'])  
        self.dumpsave = 1  
        self.called_syn_out_all = False
        self.no_fmean=False
        
        self.tau1_ex=[0*ms]*self.n_celltypes
        self.tau2_ex=[10*ms]*self.n_celltypes
        self.tau1_inh=[0*ms]*self.n_celltypes
        self.tau2_inh=[100*ms]*self.n_celltypes
        
        self.n_syn_ex = [0]*self.n_celltypes 
        self.g_syn_ex = [1]*self.n_celltypes
        self.g_syn_ex_s = [0]*self.n_celltypes
        self.mglufac_ex = [1,0]       
        
        self.noise_syn = [0]*self.n_celltypes 
        self.noise_syn_tau = [0*ms]*self.n_celltypes
        self.noise_syn_inh = [0]*self.n_celltypes
        self.noise_syn_tau_inh = [0*ms]*self.n_celltypes
        
        self.noise_a = [1e9]*self.n_celltypes
        self.noise_a_inh = [1e9]*self.n_celltypes
        
        self.inh_hold = [0]*self.n_celltypes
        self.n_syn_inh = [0]*self.n_celltypes
        self.g_syn_inh = [1]*self.n_celltypes
        self.g_syn_inh_s = [0]*self.n_celltypes
        self.intr_hold = [0]*self.n_celltypes
        self.n_syn_intr = [0]*self.n_celltypes
        self.g_syn_intr = [0]*self.n_celltypes
        self.syn_max_mf = [1]*self.n_celltypes # possible mossy fibres per synapse
        self.syn_max_inh = [1]*self.n_celltypes # possible Golgi cells per synapse
        self.syn_max_intr = [1]*self.n_celltypes # possible Intruding cells per synapse
        
        
        self.seed = 50
        
        self.force_run = False
        self.give_psd = False
        self.do_if = True
        
        self.fluct_g_e0 = []
        self.fluct_g_i0 = []
        self.fluct_std_e = []  
        self.fluct_std_i = []  
        self.fluct_tau_e = [] 
        self.fluct_tau_i = []   
        
        self.adjinh = True # adjust inhibition to get CFo instead of g_ex
        self.adjfinh = True # adjust frequnecy of inhibition to get CFo instead of g_ex
        
        self.syn_ex_dist = []
        self.syn_inh_dist = []
        
        self.stdp_used = False
        self.xmax = 20
        self.use_multisplit = False
        self.use_local_dt = False
        self.simstep = 0
        self.plot_train = True
        self.inh_delay = 0 # in ms
        self.plot_input = True
        self.delay_baseline = 8
        
        self.tstop_if = 1
        self.gsyn_in_fac = []
        
        self.netcons = [] # keeping track of!
        self.nclist = []
        
        self.ST_stims = []
        self.PF_stims = []
        
        self.data_dir = "./data"
        self.minimal_dir = False
        

    def set_init(self, ihold,  ihold_sigma, amp, amod):
        # important for all methods:
        if type(ihold) is not list: ihold = [ihold] #convert to list if it is not given as one
        self.ihold = ihold
        self.ihold_orig = ihold
        
        if type(amp) is not list: amp = [amp]
        if len(amp) < self.n_celltypes: amp = [amp[0]]*self.n_celltypes
        self.amp = amp  
        
        if type(amod) is not list: amod = [amod]*self.n_celltypes
        self.amod = amod # RUN self.set_i()
        
        self.ihold_sigma = ihold_sigma
        
    def barrier(self):
        if self.use_mpi:
            if self.use_pc == True:
                self.pc.barrier()
            else:
                self.comm.Barrier()
    
    def broadcast(self, vec, root = 0, fast = False):
        if self.use_mpi: 
            if self.use_pc:
            
                if fast:
                    hvec = h.Vector(vec)
                    v = self.pc.broadcast(hvec,root)
                    vec = np.array(hvec)
                else:
                
                    sendlist = [None]*self.nhost
                    if self.id == root:
                        for i in range(self.nhost):
                            sendlist[i] = vec                             
                    getlist = self.pc.py_alltoall(sendlist)
                    vec = getlist[root]
                
            else:
                #vec = np.array(vec, dtype=np.float64)
                #self.comm.Bcast([vec, MPI.DOUBLE])
                vec = self.comm.bcast(vec, root=0)

        return vec        
            
    def set_numcells(self, N = []):
        """
        Create, layout, and connect N cells.
        """
        self.set_gids(N)
        self.create_cells()

        #self.syn_output() # generate synaptic "output" in neuron
        #self.connect_cells()
       

    def set_gids(self, N = []):
        """Set the gidlist on this host.
        Round-robin counting. Each host as an id from 0 to pc.nhost()-1.
        Example:
        if N = 5 cells and nhost() = 3
        node id() = 0 will get cells [0, 3]
        node id() = 1 will get cells [1, 4]
        node id() = 2 will get cells [2]        
        """
        
        self.gidlist = [] 
        
        if N == []:
            N = self.N
        
        # borders where another celltype begins
        self.global_gidlist = []
        self.n_borders = [0]
        for l in range(1,self.n_celltypes+1):
            self.n_borders.append(sum(N[0:l]))
            self.global_gidlist.append(range(self.n_borders[-2], self.n_borders[-1]))

        for n in range(self.n_celltypes): # create list in list  
            self.gidlist.append([])        
            
        for i in range(int(self.id), sum(N), int(self.nhost)): # loop over all cells
            
            n = np.where((np.array(self.n_borders)-i)>0)[0][0]-1 # find out what cell type this is
            self.gidlist[n].append(i) # put in specific gidlist for that celltype
        
        self.gid_count = self.gid_count + sum(N)
        
        if self.id == 0: print "nodeid:" , self.id , ", gidlist:" , self.gidlist , ", total gids:" , len(self.global_gidlist) , ", sum(N):" , sum(N)    # check gids of node
        
        
    def del_cells(self):
        if self.cells != []:  
            for n in range(self.n_celltypes): 
                for m in self.cells[n]:
                    print "deleting cell", m
                    del m 
            del self.cells
            self.cells = []   
        if self.use_mpi: self.pc.gid_clear() 


    def create_cells(self):
        """
        Create cell objects on this host.
        """
        if self.do_run:
            
            self.del_cells()
            
            if self.id == 0: print "creating cells"
            
            for n in range(self.n_celltypes): 
                self.cells.append([]) # create list in list  
                
                #print self.cellimport[n]
                exec self.cellimport[n]
                
                #print self.gidlist
                for i in self.gidlist[n]:
                    
                    #if "sigma" not in self.cell_exe[n]:
                    #    exec self.cell_exe[n]
                    #    cell.gid = i # tell cell it's gid!
                    #    print i
                    #else:
                    
                    if (self.celltype[n] == "IfCell") or (self.celltype[n] == "Grc"):
                        
                        # add gid to cell and execute!
                        if self.cell_exe[n][-2] == "(":
                            exec self.cell_exe[n][0:-1] + "gid=" + str(i) + ")"
                        else:
                            exec self.cell_exe[n][0:-1] + ", gid=" + str(i) + ")"
                            
                    else:
                        exec self.cell_exe[n]    
                        cell.gid = i
                    
                    self.cells[n].append(cell)  # add to (local) list
    
                    if self.use_mpi:
                        #### Tell this host it has this gid
                        #### gids can be any integer, they just need to be unique.
                        #### In this simple case, we set the gid to i.
                        self.pc.set_gid2node(i, int(self.id))
                        self.pc.cell(i, cell.nc_spike) # Associate the cell with this host and gid
                        
                        ## NOT NECESSARY ANYMORE ##
                        #### Means to tell the ParallelContext that this cell is a source.
                        #nc = cell.connect_target(None)
                        #self.ncs[n].append(nc)  
                                        
                        #### Record spikes of this cell
                        self.pc.spike_record(i, self.t_vec[n], self.id_vec[n])
                        
                        #print n, self.cells[n][-1].nc_spike.thresh
                    else:
                        
                        self.t_vec[n].append(h.Vector())
                        cell.nc_spike.record(self.t_vec[n][-1])                    
                    


    def connect_cells(self, conntype=[], stdp=[], tend=1e9):
        """
        Connect cells as specified.
        """
        
        if self.do_run:
            
            stdp = stdp[:]
            conntype = conntype[:]
            
            if len(stdp) == 0:
                for i in conntype:
                    stdp.append({'wmax':0, 'taupre':0, 'taupost':0, 'apre':0, 'apost':0})    
            else:
                self.stdp_used = True
    
            for i, conn in enumerate(conntype): 
                
                typ = conn['type']
                conv = conn['conv']
                src = conn['src']
                tgt = conn['tgt']
                w0 = conn['w']
                var = conn['var']
                tau1 = conn['tau1']
                tau2 = conn['tau2']
                
                if 'mgr2' in conn.keys():
                    mgr2 = conn['mgr2']
                    mgr2_var = conn['mgr2_var']
                else:
                    mgr2 = 0
                    mgr2_var = 0
                    
                if 'e_inh' in conn.keys():    
                    e_inh = conn['e_inh']
                else:
                    e_inh = -65
                    
                if 'e_ex' in conn.keys():    
                    e_ex = conn['e_ex']
                else:
                    e_ex = 0
                    
                wmax = stdp[i]['wmax']
                taupre = stdp[i]['taupre']
                taupost = stdp[i]['taupost']
                apre = stdp[i]['apre']
                apost = stdp[i]['apost']
                    
                # Connect conv cells of celltype src to every cell of celltype tgt
                for ni, i in enumerate(self.cells[tgt]):
                    
                    rnd.seed(i.gid*10*self.seed)
                    
                    if conv >= len(self.global_gidlist[src]):
                        gids = self.global_gidlist[src]
                        if self.id == 0: print "more or equal conv to len(self.global_gidlist[src])"
                    else:
                        gids = rnd.sample(self.global_gidlist[src],conv) 
    
                    if self.id == 0: print conn['type'], ":", ni, ":", gids[0], "\n"
                    
                    for ng, g in enumerate(gids):
                        
                        np.random.seed(g*12) 
                        #np.random.seed(int(g%10+1)*12)  
                        
                        if len(shape(w0))>0: # array is given
                            print "w array is given"
                            
                            if len(w0[ng]) == self.N[0]:
                                w = w0[ng][ni]
                       
                        elif (var > 0) and (w0>0):
                            w = np.random.normal(w0, w0*var, 1).clip(min=0)
                        else:
                            w = w0
                            
                        if (mgr2_var > 0) and (mgr2>0):
                            mg = np.random.normal(mgr2, mgr2*mgr2_var, 1).clip(min=0)
                        else:
                            mg = mgr2
                
                            
                        #print conn['type'], ":", i.gid, ":", g, ", w:", w, "\n"
                        
                        if self.celltype[tgt] == 'IfCell':
                            
                            if typ == 'gogr':
                                
                                i.whatami = "grc"
                                i.synlist_inh.append(Synapse('goc', i, i.soma, nrel=0, record_all=0, weight_gmax=w))
                                i0 = int(len(i.synlist_inh)-1)
                                
                                i.nc_inh.append(self.pc.gid_connect(g, i.synlist_inh[i0].input))
                                i.nc_inh[-1].delay = 1
                                i.nc_inh[-1].weight[0] = 1
                            
                            if typ == 'grgo':
                                
                                i.whatami = "goc"
                                i.synlist.append(Synapse('grc', i, i.soma, syntype = 'D', nrel=0, record_all=0, weight_gmax=w))
                                e0 = int(len(i.synlist)-1)
                                
                                i.nc.append(self.pc.gid_connect(g, i.synlist[e0].input))
                                i.nc[-1].delay = 1
                                i.nc[-1].weight[0] = 1
                                
                            if typ == 'grgom':
                                
                                i.whatami = "goc"
                                i.synlist.append(Synapse('grc', i, i.soma, syntype = 'DM', nrel=0, record_all=0, weight_gmax=w, mglufac = mg))
                                e0 = int(len(i.synlist)-1)
                                
                                i.nc.append(self.pc.gid_connect(g, i.synlist[e0].input))
                                i.nc[-1].delay = 1
                                i.nc[-1].weight[0] = 1
                              
                                
                            if typ == 'e2inh':
                                
                                i.create_synapses(n_inh=1, tau1_inh=tau1, tau2_inh=tau2, e_inh=e_inh, w = w, wmax = wmax, taupre = taupre, taupost = taupost, apre = apre, apost = apost, tend=tend)
                                i0 = len(i.synlist_inh)-1
                                
                                if self.use_mpi:
                                    if wmax == 0:
                                        i.pconnect_target(self.pc, source=g, target=i0, syntype='inh', weight=w, delay=1)
                                    else:
                                        i.pconnect_target(self.pc, source=g, target=i0, syntype='inh', weight=1, delay=1)
                                        
                                else:  
                                    if wmax == 0:
                                        i.nc_inh.append(self.cells[1][g-self.N[0]].connect_target(target=i.synlist_inh[i0], weight=w, delay=1))
                                    else:
                                        i.nc_inh.append(self.cells[1][g-self.N[0]].connect_target(target=i.synlist_inh[i0], weight=1, delay=1))
                            
                            if typ == 'e2ex':
                                
                                i.create_synapses(n_ex = 1, tau1 = tau1, tau2 = tau2, e_ex=e_ex, w = w, wmax = wmax, taupre = taupre, taupost = taupost, apre = apre, apost = apost, tend=tend)
                                e0 = len(i.synlist)-1
    
                                if self.use_mpi:
                                    if wmax == 0:
                                        i.pconnect_target(self.pc, source=g, target=e0, syntype='ex', weight=w, delay=1) 
                                    else:
                                        i.pconnect_target(self.pc, source=g, target=e0, syntype='ex', weight=1, delay=1)  
                                        
                                else:  
                                    if wmax == 0:
                                        i.nc.append(self.cells[0][g].connect_target(target=i.synlist[e0], weight=w, delay=1))
                                    else:
                                        i.nc.append(self.cells[0][g].connect_target(target=i.synlist[e0], weight=1, delay=1))
    
                        else: # No IfCell
                            
                            if typ == 'gogr':
                                i.createsyn(ngoc = 1, weight_gmax=w) # multiplication factor
                                i0 = len(i.GOC_L)-1 # get number of current synapse!
                                i.pconnect(self.pc,g,i0,'goc')
                                
                            if typ == 'grgo':
                                i.createsyn(ngrc = 1, weight_gmax=w) # multiplication factor
                                i0 = len(i.GRC_L)-1 # get number of current synapse!
                                i.pconnect(self.pc,g,i0,'grc',conduction_speed=0,grc_positions=[1])
                                
                            if typ == 'grgom':
                                #print w, mg
                                i.createsyn(ngrcm = 1, weight_gmax=w, mglufac = mg) # multiplication factor
                                i0 = len(i.GRC_L)-1 # get number of current synapse!
                                i.pconnect(self.pc,g,i0,'grc',conduction_speed=0,grc_positions=[1])
                                
                            if typ == 'grstl':
                                i.createsyn(ngrc = 1, weight_gmax=w) # multiplication factor
                                i0 = len(i.GRC_L)-1 # get number of current synapse!
                                i.pconnect(self.pc,g,i0,'grc',conduction_speed=0,grc_positions=[1])
                                   
                                    
                            if 'e2' in typ:
                                
                                if 'inh' in typ:
                                    Erev = -65
                                elif 'ex' in typ:
                                    Erev = 0
                                
                                if tau1 == 0:
                                    syn = h.ExpSyn(i.soma(0.5))
                                    syn.tau = tau2/ms
                                else:                
                                    if wmax == 0:
                                        syn = h.Exp2Syn(i.soma(0.5))
                                        syn.tau1 = tau1/ms
                                        syn.tau2 = tau2/ms
                
                                    else: # STDP
                                        syn = h.stdpE2S(i.soma(0.5))
                                        syn.tau1 = tau1/ms
                                        syn.tau2 = tau2/ms
                                        
                                        syn.on = 1
                                        syn.thresh = -20
                                        
                                        syn.wmax = wmax
                                        syn.w = w
                                        
                                        syn.taupre  = taupre/ms	
                                        syn.taupost  = taupost/ms
                                        syn.apre    = apre
                                        syn.apost   = apost
    
                                syn.e = Erev/mV
                                
                                if self.celltype[tgt] == 'Grc':
                                    
                                    i.GOC_L.append(syn)
                                    i0 = int(len(i.GOC_L)-1) # get number of current synapse!
        
                                    i.gocncpc.append(self.pc.gid_connect(g, i.GOC_L[i0]))
                                    i.gocncpc[-1].delay = 1
                                    
                                    if wmax == 0:
                                        i.gocncpc[-1].weight[0] = w
                                    else:
                                        i.gocncpc[-1].weight[0] = 1
                                        
                                elif self.celltype[tgt] == 'Goc':
                                    
                                    i.GRC_L.append(syn)
                                    e0 = int(len(i.GRC_L)-1) # get number of current synapse!
        
                                    i.pfncpc.append(self.pc.gid_connect(g, i.GRC_L[e0]))
                                    i.pfncpc[-1].delay = 1
                                    i.pfncpc[-1].weight[0] = w
                                    
                                    if wmax == 0:
                                        i.pfncpc[-1].weight[0] = w
                                    else:
                                        i.pfncpc[-1].weight[0] = 1
                                                
            #self.rec_s1 = h.Vector()
            #self.rec_s1.record(self.cells[0][0].synlist_inh[0]._ref_g)    
            #self.rec_s2 = h.Vector()
            #self.rec_s2.record(self.cells[1][0].synlist_inh[0]._ref_g) 
                    
                
    def syn_output(self):
        """
        Connect cell n to target cell sum(self.N) + 100.
        """
        
        if self.id == 0:  # create target cell

            tgt_gid = self.gid_count
            self.gid_count = self.gid_count + 1 
            
            # Synaptic integrated response
            self.rec_g = h.Vector()
            self.passive_target = PassiveCell()
            if self.use_mpi: self.pc.set_gid2node(tgt_gid, 0)  # Tell this host it has this gid
            
            syn = self.passive_target.create_synapses(tau1 = self.syn_tau1, tau2 = self.syn_tau2)  # if tau1=tau2: alpha synapse!
            
            for i in range(self.n_borders[self.a_celltype[0]],self.n_borders[self.a_celltype[0]+1]): # take all cells, corresponding to self.a_celltype, not just the ones in self.gidlist:
            
                src_gid = i
                
                if self.use_mpi:
                    nc = self.pc.gid_connect(src_gid, syn)
                    nc.weight[0] = 1
                    nc.delay = self.nc_delay/ms #0.05  # MUST be larger than dt!!!
                
                else:
                    nc = self.cells[self.a_celltype[0]][src_gid].connect_target(target=syn, weight=1, delay=self.nc_delay/ms)
                
                self.nclist.append(nc) 
                
            self.rec_g.record(syn._ref_g)
            
            
    def syn_out_all(self, tau1 = 1*ms, tau2 = 30*ms):
        
        if self.do_run:
        
            for n in range(self.n_celltypes): 
                for i, gid in enumerate(self.gidlist[n]):
                    
                    self.cells[n][i].start_record(tau1 = tau1/ms, tau2 = tau2/ms)
                    
            self.called_syn_out_all = True
            
            
    def get_i(self, a, n, do_plot = True):
            
        import md5
        m = md5.new()
        
        if ", sigma" in self.cell_exe[n]: 
            cell_exe_new = self.cell_exe[n].split(", sigma")[0] + ")"
        else:
            cell_exe_new = self.cell_exe[n]
        
        m.update(cell_exe_new)
        filename = self.data_dir + '/if_' + self.celltype[n] + '_' + m.hexdigest() + '.p'
        
        #print filename
        
        if self.id == 0:
            is_there = os.path.isfile(filename)
        else:
            is_there = None
        
        is_there = self.broadcast(is_there)
        
        if (is_there is not True) or (self.force_run is True): # run i/f estimation
        
            if self.id == 0: print '- running i/f estimation for ', self.celltype[n], ' id: ' , m.hexdigest() 
            exec self.cellimport[n]
            exec cell_exe_new
            sim = Stimulation(cell, temperature = self.temperature, use_multisplit = self.use_multisplit)
            sim.spikes_from_neuron = False
            sim.celltype =  self.celltype[n]
            current_vector, freq_vector, freq_onset_vector = sim.get_if(istart = self.istart, istop = self.istop, di = self.di, tstop = self.tstop_if) 
            
            sim = None
            cell = None
            
            if self.id == 0:
                if do_plot:
                    plt.figure(99)
                    plt.plot(current_vector, freq_vector, 'r*-')
                    plt.plot(current_vector, freq_onset_vector, 'b*-')
                    plt.savefig("./figs/dump/latest_if_" + self.celltype[n]  + ".pdf", dpi = 300)  # save it 
                    plt.clf()
                    #plt.show()
            
                ifv = {'i':current_vector,'f':freq_vector}
                print ifv
            
                pickle.dump(ifv, gzip.GzipFile(filename, "wb" ))
            
            self.barrier()
            
        else:
            
            if self.id == 0:            
                ifv = pickle.load(gzip.GzipFile(filename, "rb" ))
                #print ifv
        
        self.barrier()
        
        if self.id == 0:
            
            f =  ifv.get('f') 
            i =  ifv.get('i')
        
            i = i[~isnan(f)]
            f = f[~isnan(f)]
        
            iin = if_extrap(a, f, i)
        
        else:
            
            iin = [0]
            
        iin = self.broadcast(iin, root=0, fast = True)
        self.barrier()
        
        return iin


    def set_i(self, ihold = [0]):
        
        ihold = list(ihold)
        self.ihold_orig = list(ihold)
        
        self.barrier()   # wait for other nodes
           
        # Ihold given as frequency, convert to current
        
        if ((self.give_freq)): 
            
            ihold0 = [[] for _ in range(self.n_celltypes)]
            
            for n in range(self.n_celltypes):
                a = np.array([ihold[n]])
                #print "a:", a
                iin = self.get_i(a, n)
                #print "iin:", iin
                ihold0[n] = iin[0]
                
            if self.id == 0: print '- ihold: ', ihold, 'Hz, => ihold: ', ihold0, 'nA' 
   
        # Modulation depth given, not always applied to current!
        for n in range(self.n_celltypes):
            
            if self.amod[n] is not None:
            
                if self.give_freq:
                    
                    # Apply to amplitude:
                    a = np.array([ihold[n]]) + self.amod[n]*np.array([ihold[n]])
                    self.amp[n] = self.get_i(a, n) - ihold0[n]
                    
                    if self.id == 0:
                        print '- amp: ihold: ', ihold[n], 'Hz , amod: ', self.amod[n], ', => amp: ', self.amp[n], 'nA (' #, self.get_i(a, n), ')'
                        
                elif self.n_syn_ex[n] > 0:
                
                    if self.id == 0:
                        print '- amp: ihold: ', ihold[n], 'Hz , amod: ', self.amod[n], ', => amp will be set for each spike generator'

                else:
                    
                    self.amp[n] = self.amod[n] * ihold[n]     
                
                    if self.id == 0:
                        print '- amp: ihold: ', ihold[n], 'nA , amod: ', self.amod[n], ', => amp: ', self.amp[n], 'nA'
                      
            # Noise depth given, not always applied to current!
            if self.anoise[n] is not None:
                
                if (self.give_freq is True) or (self.n_syn_ex[n] > 0):
                    
                    # Apply to amplitude:
                    a = np.array([ihold[n]]) + self.anoise[n]*np.array([ihold[n]])
                    self.fluct_s[n] = ((self.get_i(a, n) - ihold0[n]))/2. # adjust with /2 so that noise = +-2*std
                    
                    if self.id == 0:
                        print '- noise: ihold: ', ihold[n], 'Hz , anoise: ', self.anoise[n], ', => fluct_s: ', self.fluct_s[n], 'nA'
    
                else:
                    
                    self.fluct_s[n] = self.anoise[n] * ihold[n] 
                
                    if self.id == 0:
                        print '- noise: ihold: ', ihold[n], 'nA , anoise: ', self.anoise[n], ', => fluct_s: ', self.fluct_s[n], 'nA'
                        
        
        if self.give_freq is True:        
            ihold = ihold0
                
        return ihold
        
        
    def calc_fmean(self, t_vec, t_startstop):
        
        #t_startstop[0] = 1
        #t_startstop[1] = 5
        
        f_cells_mean = 0
        f_cells_cv = np.nan
        f_cells_std = np.nan
                      
        if len(t_vec) > 0: 
            
            f_start_in = mlab.find(t_vec >= t_startstop[0]) # 1
            f_stop_in = mlab.find(t_vec <= t_startstop[1]) # 5
            
            if (len(f_start_in) > 0) & (len(f_stop_in) > 0):
                
                f_start = f_start_in[0]  
                f_stop = f_stop_in[-1]+1      
                use_spikes = t_vec[f_start:f_stop]*1e3
                
                if len(use_spikes) > 1:
                    s1 = signals.SpikeTrain(use_spikes)
                    isi = s1.isi()
                    f_cells_mean = s1.mean_rate() # use mean of single cells
                    f_cells_cv = np.std(isi)/np.mean(isi)
                    f_cells_std = np.std(isi)
                    
            #f_start_in = mlab.find(t_vec >= 1) 
            #f_stop_in = mlab.find(t_vec <= 2) 
            
            #if (len(f_start_in) > 0) & (len(f_stop_in) > 0):
                
            #    f_start = f_start_in[0]  
            #    f_stop = f_stop_in[-1]+1      
            #    use_spikes = t_vec[f_start:f_stop]*1e3
                
            #    if len(use_spikes) > 1:
            #        s1 = signals.SpikeTrain(use_spikes)
            #        isi = s1.isi()
            #        f_cells_cv = np.std(isi)/np.mean(isi)
        
        return f_cells_mean, f_cells_cv, f_cells_std 
                    
                    
    def get_fmean(self, t_all_vec_vecn, id_all_vec_vecn, t_startstop, gidlist, facborder = 3): # 1e9
        
        f_cells_mean = zeros(len(gidlist))
        f_cells_base = zeros(len(gidlist))
        f_cells_std = nans(len(gidlist))
        f_cells_cv = nans(len(gidlist))
        f_cells_gid = nans(len(gidlist))
        
        fbase = np.nan
        fmean = np.nan
        fmax = np.nan
        fmstd = np.nan
        fcvm = np.nan
        fstdm = np.nan
        
        f_cells_mean_all = []
        f_cells_base_all = []
        f_cells_cv_all = []
        f_cells_std_all = []
        
        gid_del = np.array([])
        
        if self.no_fmean == False:
            
            if self.id == 0: print "- sorting for fmean"

            for i, l in enumerate(gidlist):
                
                t_0_vec = t_all_vec_vecn[where(id_all_vec_vecn==l)]
                f_cells_mean[i], f_cells_cv[i], f_cells_std[i] = self.calc_fmean(t_0_vec, t_startstop)
                f_cells_base[i], _, _ = self.calc_fmean(t_0_vec, [self.delay_baseline-4,self.delay_baseline])
                f_cells_gid[i] = l
                        
            if self.id == 0:  print "- gather fmean"            
            f_cells_mean_all = self.do_gather(f_cells_mean)
            f_cells_base_all = self.do_gather(f_cells_base)
            f_cells_std_all = self.do_gather(f_cells_std)
            f_cells_cv_all = self.do_gather(f_cells_cv)
            f_cells_gid_all = self.do_gather(f_cells_gid)

            if self.id == 0:
                
                #print f_cells_mean_all
                
                f_cells_mean_all = np.nan_to_num(f_cells_mean_all)
                fmean = mean(f_cells_mean_all)  # compute mean of mean rate for all cells
                fmstd = std(f_cells_mean_all) 
                fmax = max(f_cells_mean_all)
                
                f_cells_base_all = np.nan_to_num(f_cells_base_all)
                fbase = mean(f_cells_base_all)  # compute mean of mean rate for all cells
                
                f_cells_cv_all = f_cells_cv_all[~np.isnan(f_cells_cv_all)]
                f_cells_std_all = f_cells_std_all[~np.isnan(f_cells_std_all)]
                fcvm = mean(f_cells_cv_all)
                fstdm = mean(f_cells_std_all)
            
                print "- get_fmean, fmean: ",fmean, "fmax: ",fmax, "Hz", "fmstd: ",fmstd, "Hz", "fcvm: ",fcvm, "fstdm: ",fstdm, "Hz" ,"fbase: ", fbase, "Hz"
                
                if facborder < 1e9:
                
                    fborder = fmean + facborder*fmstd
                    i = mlab.find(f_cells_mean_all > fborder)
                    gid_del = f_cells_gid_all[i]
                    
                #    f_cells_mean_all[i] = 0
                #    f_cells_cv_all[i] = np.nan
                #    f_cells_std_all[i] = np.nan
                    
                #    fmean2 = mean(np.nan_to_num(f_cells_mean_all))  # compute mean of mean rate for all cells
                #    fmstd2 = std(np.nan_to_num(f_cells_mean_all)) 
                #    fmax2 = max(np.nan_to_num(f_cells_mean_all))
                
                #    fcvm2 = mean(f_cells_cv_all[~np.isnan(f_cells_cv_all)])
                #    fstdm2 = mean(f_cells_std_all[~np.isnan(f_cells_std_all)])
            
                #    print "- after facborder: get_fmean, fmean: ",fmean2, "fmax: ",fmax2, "Hz", "fmstd: ",fmstd2, "Hz", "fcvm: ",fcvm2, "fstdm: ",fstdm2, "Hz, gid_del: ", gid_del
 

        return fmean, fmax, fmstd, fcvm, fstdm, gid_del, f_cells_mean_all, f_cells_cv_all, f_cells_std_all, fbase, f_cells_base_all 

                
    def connect_fluct(self):
        """
        Create fluctuating input onto every cell.
        """
        
        if self.do_run:
        
            for m in self.flucts:
                del m 
            del self.flucts
            
            for m in self.noises:
                del m 
            del self.noises
            
            self.flucts = []
            self.noises = []
            
            for n in range(self.n_celltypes):
                
                for i, gid in enumerate(self.gidlist[n]):  # for every cell in the gidlist 
                                
                    #h.mcell_ran4_init(gid)
                
                    noiseRandObj = h.Random()  # provides NOISE with random stream
                    self.noises.append(noiseRandObj)  # has to be set here not inside the nmodl function!!   
                    
                    # print str(gid) + ": " + str(noiseRandObj.normal(0,1))
                                    
                    fluct = h.Ifluct2(self.cells[n][i].soma(0.5))
                    fluct.m = self.fluct_m/nA      # [nA]
                    fluct.s = self.fluct_s[n]/nA      # [nA]
                    fluct.tau = self.fluct_tau/ms    # [ms]
                    self.flucts.append(fluct)   # add to list     
                    self.flucts[-1].noiseFromRandom(self.noises[-1])  # connect random generator!
                    
                    self.noises[-1].MCellRan4(1, gid+1)  # set lowindex to gid+1, set highindex to > 0 
                    self.noises[-1].normal(0,1)
                    
                    
    def connect_gfluct(self, E_e=0, E_i=-65):
        """
        Create fluctuating conductance input onto every cell.
        """
        if self.do_run:
        
            for m in self.flucts:
                del m 
            del self.flucts
            
            for m in self.noises:
                del m 
            del self.noises
            
            self.flucts = []
            self.noises = []
            
            for n in range(self.n_celltypes):
                
                fluct_g_i0_n = self.fluct_g_i0[n]
            
                if type(fluct_g_i0_n) is not ndarray: fluct_g_i0_n = np.array([fluct_g_i0_n])
            
                if len(fluct_g_i0_n) == len(self.global_gidlist[n]):
                    pass
                else:
                    fluct_g_i0_n = np.ones(int(len(self.global_gidlist[n])))*fluct_g_i0_n[0]
                    if self.id == 0: print "- single value in fluct_g_i0_n"
                
                #print fluct_g_i0_n
                
                for i, gid in enumerate(self.gidlist[n]):  # for every cell in the gidlist 
                                
                    #h.mcell_ran4_init(gid)
                
                    noiseRandObj = h.Random()  # provides NOISE with random stream
                    self.noises.append(noiseRandObj)  # has to be set here not inside the nmodl function!!   
                    
                    # print str(gid) + ": " + str(noiseRandObj.normal(0,1))
                                    
                    fluct = h.Gfluct3(self.cells[n][i].soma(0.5))
                    fluct.E_e = E_e/mV  # [mV]
                    fluct.E_i = E_i/mV  # [mV]
                    fluct.g_e0 = self.fluct_g_e0[n]/uS  # [uS]
                    fluct.g_i0 = fluct_g_i0_n[i]/uS  # [uS]
                    fluct.std_e = self.fluct_std_e[n]/uS  # [uS] 
                    fluct.std_i = self.fluct_std_i[n]/uS  # [uS] 
                    fluct.tau_e = self.fluct_tau_e/ms #tau_e/ms  # [ms] 
                    fluct.tau_i = self.fluct_tau_i/ms #tau_i/ms  # [ms]
            
                    self.flucts.append(fluct)   # add to list     
                    self.flucts[-1].noiseFromRandom(self.noises[-1])  # connect random generator!
                    
                    self.noises[-1].MCellRan4(1, gid+1)  # set lowindex to gid+1, set highindex to > 0 
                    self.noises[-1].normal(0,1)
                    
                    
    def connect_synfluct(self, PF_BG_rate=6, PF_BG_cv=1, STL_BG_rate=20, STL_BG_cv=1):
        """
        Create fluctuating synaptic input onto every cell.
        """
        
        if self.do_run:
            
            for m in self.ST_stims:
                del m 
            del self.ST_stims
            
            for m in self.PF_stims:
                del m 
            del self.PF_stims
            
            self.ST_stims = []
            self.PF_stims = []
            
            
            for n in range(self.n_celltypes):
                
                for i, gid in enumerate(self.gidlist[n]):  # for every cell in the gidlist 
            
                    PF_syn_list = self.cells[n][i].createsyn_PF()
            
                    for d in PF_syn_list:
                        d.input.newnetstim.number = 1e9
                        d.input.newnetstim.noise = PF_BG_cv
                        d.input.newnetstim.interval = 1000.0 / PF_BG_rate
                        d.input.newnetstim.start = 0
                    
                    self.PF_stims.append(PF_syn_list)
                    
                    ST_stim_list = self.cells[n][i].createsyn_ST(record_all=0)

                    for d in ST_stim_list:
                        d.newnetstim.number = 1e9
                        d.newnetstim.noise = STL_BG_cv
                        d.newnetstim.interval =  1000.0 / STL_BG_rate
                        d.newnetstim.start = 0
                        
                    self.ST_stims.append(ST_stim_list)
                                
            if self.id == 0: print "- PF and ST stimulation added."
            
                

    def set_IStim(self, ihold = None, ihold_sigma = None, random_start = True, tstart_offset = 0):
        """
        Add (random) ihold for each cell and offset!
        """
        if self.do_run:
            
            # if not given, use the one in self
            if ihold == None:
                ihold = self.ihold
            if ihold_sigma == None:
                ihold_sigma = self.ihold_sigma
            
            if ihold[self.a_celltype[0]] != 0:
                ihold = self.set_i(ihold)        
            
            for m in self.ic_holds:
                #m.destroy()
                del m 
            del self.ic_holds
            
            for m in self.ic_starts:
                #m.destroy()
                del m 
            del self.ic_starts
            
            for m in self.vc_starts:
                #m.destroy()
                del m 
            del self.vc_starts
            
            self.ic_holds = []
            self.ic_starts = [] 
            self.vc_starts = []
            self.i_holdrs = []
            self.i_holds = ihold
                            
            for n in range(self.n_celltypes):
                self.i_holdrs.append([])
                            
                for i, gid in enumerate(self.gidlist[n]):  # for every cell in the gidlist 
                            
                    np.random.seed(gid*20)
                    
                    tis = 1
                    
                    if random_start == True:
                                                
                        # random start time
                        tstart = np.random.uniform(tstart_offset+0, tstart_offset+0.5)
                        #if self.id == 0: print "tstart:", tstart
                        vc_start = h.SEClamp(self.cells[n][i].soma(0.5))
                        vc_start.dur1 = tstart/ms
                        vc_start.amp1 = -80
                        self.vc_starts.append(vc_start)
                        tis = 0
                                                
                    else:
                        
                        tis = 0    
                    
                    
                    if ihold_sigma[n] != 0:
                        #print ihold_sigma[n], ihold[n]
                        ihold_r = np.random.normal(ihold[n], ihold[n]*ihold_sigma[n], 1).clip(min=0)
                        #ihold_r = np.random.uniform(ihold[n]*ihold_sigma[n], ihold[n])
                        
                    elif self.CF_var is not False:  # CF gets not adapted to current but final frequnecy!
                        
                        r_ok = False
                        while r_ok == False:
                            r_temp = np.random.normal(self.ihold_orig[n], self.CF_var[n][1], 1)    
                            if (r_temp <= self.CF_var[n][2]) and (r_temp >= self.CF_var[n][0]): # check borders!
                                r_ok = True
                                
                        #print r_temp 
                        ihold_r = self.get_i(r_temp, n)
                        #print ihold_r
                        #if self.id == 0: 
                        print "set self.CF_var", r_temp, ihold_r
                    
                    else:  # same ihold for all cells!
                        ihold_r = ihold[n]
                        
                    self.i_holdrs[n].append(ihold_r)
                    
                    if ihold_r != 0:
                        
                        if hasattr(self.cells[n][i], 'input_vec'):
            
                            ic_hold = []
                            for vec in self.cells[n][i].input_vec:
                                for inv in vec:
                                    #print ihold_r
                                    ic_hold.append(h.IClamp(inv(0.5))) 
                                    ic_hold[-1].amp = self.cells[n][i].ifac * ihold_r / self.cells[n][i].n_input_spiny / nA
                                    ic_hold[-1].delay = tis/ms
                                    ic_hold[-1].dur = 1e9
                                
                        else:                        

                            # holding current
                            ic_hold = h.IClamp(self.cells[n][i].soma(0.5))
                            ic_hold.delay = tis/ms
                            ic_hold.dur = 1e9
                            ic_hold.amp = ihold_r/nA
                        
                        self.ic_holds.append(ic_hold)
    
            if self.id == 0: print "set_IStim finished. ihold: ", ihold, ", ihold_sigma: ", ihold_sigma
        
        
    def set_IStep(self, istep = [0], istep_sigma = [0], tstep = 5, tdur = 1e6, give_freq = True):
        """
        Add istep for each cell and offset!
        """
        if self.do_run:
            #for m in self.ic_steps:
            #    m.destroy()
            #    del m 
            #del self.ic_steps
            
            #self.ic_steps = []
            
            istep = list(istep)
            neg = False
            
            for n in range(self.n_celltypes):
                
                if istep[n] < 0: 
                    neg = True
                    istep[n] = abs(istep[n]) # make positive again
                                    
                if istep[n] != 0:
                    if give_freq is True:
                        a = np.array([istep[n]])
                        iin = self.get_i(a, n)[0]
                        if self.id == 0: print "celltype: ", n, " istep: ", istep[n], "Hz => ", iin, " nA"
                        istep[n] = iin  
                        
            for n in range(self.n_celltypes):
                
                for i, gid in enumerate(self.gidlist[n]):  # for every cell in the gidlist 
                            
                    np.random.seed(gid*30)
                    
                    if self.i_holdrs == []:
        
                        if istep_sigma[n] != 0:
                            istep_r = np.random.normal(istep[n], istep[n]*istep_sigma[n], 1).clip(min=0)
                        else:  # same ihold for all cells!
                            istep_r = istep[n]
                        
                    else: # ihold has been set!
                        
                        if istep_sigma[n] != 0:
                            istep_r = np.random.normal(istep[n]-self.i_holds[n], (istep[n]-self.i_holds[n])*istep_sigma[n], 1).clip(min=0) # delta now! put on top of hold!
                        else:  # same ihold for all cells!
                            istep_r = istep[n]-self.i_holds[n] # delta now! put on top of hold!
                    
                    if neg:
                        istep_r = -1*istep_r
                        
                    if istep[n] == 0:
                        istep_r = -1*self.i_holdrs[n][i]  
                        
                    #print 'is:' + str(istep_r) + 'was:' + str(self.i_holdrs[n][i])
                    
                    if istep_r != 0:     
                        # step current
                        ic_step = h.IClamp(self.cells[n][i].soma(0.5))
                        ic_step.delay = tstep/ms
                        ic_step.dur = tdur/ms
                        ic_step.amp = istep_r/nA
                        self.ic_steps.append(ic_step)
                        
                    
            if self.id == 0: print "set_IStep finished. istep: ", istep, ", istep_sigma: ", istep_sigma
                        

    def set_IPlay(self, stimulus, t):
        """
        Initializes values for current clamp to play a signal. 
        """
        
        if self.do_run:
            
            for m in self.tvecs:
                #m.destroy()
                del m 
            del self.tvecs
            
            for m in self.ivecs:
                #m.destroy()
                del m 
            del self.ivecs
            
            for m in self.plays:
                #m.destroy()
                del m 
            del self.plays
            
            self.tvecs = []
            self.ivecs = []
            self.plays = []
            
            for i, gid in enumerate(self.gidlist[self.a_celltype[0]]):  # for every cell in the gidlist 
                        
                tvec = h.Vector(t/ms)
                ivec = h.Vector(stimulus/nA)
            
                play = h.IClamp(self.cells[self.a_celltype[0]][i].soma(0.5))
                play.delay = 0
                play.dur = 1e9
                
                ivec.play(play._ref_amp, tvec, 1)
                
                self.plays.append(play)   # add to list
                self.tvecs.append(tvec)   # add to list
                self.ivecs.append(ivec)   # add to list 
                    
            if self.id == 0: print "set_IPlay finished."
                
        
    def set_IPlay2(self, stimulus, t):
        """
        Initializes values for current clamp to play a signal. 
        """
        
        if self.do_run:
            
            for m in self.tvecs:
                #m.destroy()
                del m 
            del self.tvecs
            
            for m in self.ivecs:
                #m.destroy()
                del m 
            del self.ivecs
            
            for m in self.plays:
                #m.destroy()
                del m 
            del self.plays
            
            self.tvecs = []
            self.ivecs = []
            self.plays = []
                
            for j in self.a_celltype:
                
                tvec = h.Vector(t/ms)
                ivec = []
                for s in stimulus:
                    if hasattr(self.cells[j][0], 'input_vec'):
                        ivec.append(h.Vector(self.factor_celltype[j] * self.cells[j][0].ifac * s / self.cells[j][0].n_input_spiny / nA))
                    else:
                        ivec.append(h.Vector(self.factor_celltype[j]*s/nA))

                self.tvecs.append(tvec)   # add to list
                self.ivecs.append(ivec)   # add to list  
    
                for i, gid in enumerate(self.gidlist[j]):  # for every cell in the gidlist 

                    if hasattr(self.cells[j][i], 'input_vec'):
                        
                        play = []
                        for iloc, vec in enumerate(self.cells[j][i].input_vec):
                            isig = self.syn_ex_dist[j][iloc]-1
                            #print isig
                            for inv in vec:
                                play.append(h.IClamp(inv(0.5))) 
                                play[-1].delay = 0
                                play[-1].dur = 1e9
                                ivec[isig].play(play[-1]._ref_amp, tvec, 1)
                                
                    else:                        
                        #fluctuating current
                        play = h.IClamp(self.cells[j][i].soma(0.5))
                        play.delay = 0
                        play.dur = 1e9
                        ivec[0].play(play._ref_amp, tvec, 1)
                    
                    self.plays.append(play)   # add to list

                    
            if self.id == 0: print "set_IPlay2 finished."
        
        
    def set_IPlay3(self, stimulus, t, amp = None):
        """
        Initializes values for current clamp to play a signal. 
        """
        
        if self.do_run:
            
            for m in self.tvecs:
                #m.destroy()
                del m 
            del self.tvecs
            
            for m in self.ivecs:
                #m.destroy()
                del m 
            del self.ivecs
            
            for m in self.plays:
                #m.destroy()
                del m 
            del self.plays
            
            self.tvecs = []
            self.ivecs = []
            self.plays = []
                
            for j in self.a_celltype:
                
                if amp is None:
                    amp0 = 0
                else:
                    amp0 = amp[j]
                
                tvec = h.Vector(t/ms)
                self.tvecs.append(tvec)   # add to list
    
                for i, gid in enumerate(self.gidlist[j]):  # for every cell in the gidlist 
    
                    if isinstance(self.factor_celltype[j], ( int, long ) ):                  
                        ivec = h.Vector(self.factor_celltype[j]*(stimulus*amp0)/nA)    
                    else:
                        np.random.seed(gid*40)
                        rnd.seed(gid*40)
                        if self.factor_celltype[j][1] > 0:
                            f = np.random.normal(self.factor_celltype[j][0], self.factor_celltype[j][1], 1).clip(min=0)
                        else:
                            f = self.factor_celltype[j][0]    
                        if self.factor_celltype[j][2] > 0: # add inverted input with 50% probability, in future versions this will indicate the propability for -1 and 1
                            f = rnd.sample([-1,1],1)[0] * f
                            if self.id == 0: print "- inverted input with 50% probability:", f 
                        if self.id == 0: print "- randomize play stimulus height"   
                        ivec = h.Vector(f*(stimulus*amp0)/nA)
    
                    self.ivecs.append(ivec)   # add to list  
                                    
                    #fluctuating current
                    play = h.IClamp(self.cells[j][i].soma(0.5))
                    play.delay = 0
                    play.dur = 1e9
                    ivec.play(play._ref_amp, tvec, 1)
                    
                    self.plays.append(play)   # add to list
                    
            if self.id == 0: print "set_IPlay3 finished."
        
        
    def set_PulseStim(self, start_time=[100*ms], dur=[1500*ms], steadyf=[100*Hz], pulsef=[150*Hz], pulse_start=[500*ms], pulse_len=[500*ms], weight0=1, tau01=[1*ms], tau02=[20*ms], weight1=1, tau11=[0*ms], tau12=[1*ms], noise = 1):
        
        if self.do_run:
            
            modulation_vec = []
            
            for n in range(self.n_celltypes):
                
                t_input = np.arange(0, dur[n], self.dt) # create stimulus time vector has to be in ms!!     
                mod = np.concatenate(([np.zeros(round(start_time[n]/self.dt)), steadyf[n]*np.ones(round((pulse_start[n]-start_time[n])/self.dt)), pulsef[n]*np.ones(round(pulse_len[n]/self.dt)),steadyf[n]*np.ones(round((dur[n]-pulse_start[n]-pulse_len[n])/self.dt)) ])) 
                modulation = (t_input, mod)
                
                #print shape(t_input), shape(mod), shape(modulation)
    
                for i, gid in enumerate(self.gidlist[n]):  # for every cell in the gidlist 
                
                    if dur[n] > 0:
                    
                        if self.celltype[n] == 'Grc':
                            
                            nmf = 4
    
                            for j in range(nmf):
                                
                                self.cells[n][i].createsyn(nmf = 1, ngoc = 0, weight = weight0) 
                                e0 = len(self.cells[n][i].MF_L)-1 # get number of current synapse!
                                
                                pulse_gid = int(self.gid_count + gid*1000 + j)
                                
                                train = mod_spike_train(modulation, noise = noise, seed = pulse_gid)
                                            
                                self.setup_Play_train(train = train, input_gid = pulse_gid)
                        
                                self.cells[n][i].pconnect(self.pc,pulse_gid,int(e0),'mf')   
                                
                        elif self.celltype[n] == 'Goc':
                            
                            nmf = 53
     
                            for j in range(nmf):
                                
                                self.cells[n][i].createsyn(nmf = 1, weight = weight1)
                                e0 = len(self.cells[n][i].MF_L)-1 # get number of current synapse!
                                
                                pulse_gid = int(self.gid_count + gid*1000 + j)
                                
                                train = mod_spike_train(modulation, noise = noise, seed = pulse_gid)
                                
                                self.setup_Play_train(train = train, input_gid = pulse_gid)
                
                                self.cells[n][i].pconnect(self.pc,pulse_gid,int(e0),'mf')  
                                
                        
                        elif self.celltype[n] == 'Goc_noloop':
                            
                            ngrc = 100
     
                            for j in range(ngrc):
                                
                                self.cells[n][i].createsyn(ngrc = 1, weight = weight0)
                                e0 = len(self.cells[n][i].GRC_L)-1 # get number of current synapse!
                                
                                pulse_gid = int(self.gid_count + gid*1000 + j)
                                
                                train = mod_spike_train(modulation, noise = noise, seed=pulse_gid)
                                
                                self.setup_Play_train(train = train, input_gid = pulse_gid)
                                
                                self.cells[n][i].pconnect(self.pc,pulse_gid,int(e0),'grc')  
    
                        else:
                            
                            pulse_gid = int(self.gid_count + gid*1000 + 100)
                        
                            train = mod_spike_train(modulation, noise = noise, seed = pulse_gid)
                            self.trains.append(train)
                
                            setup_Play_train(train = train, input_gid = pulse_gid)
                            
                            # NMDA
                            self.cells[n][i].create_synapses(n_ex=1, tau1=tau01[n], tau2=tau02[n])
                            e0 = len(self.cells[n][i].synlist)-1
                    
                            weight=weight0[n]
                            np.random.seed(gid*60)
                            #weight = np.random.normal(weight, weight*0.5, 1).clip(min=0)
                            self.cells[n][i].pconnect_target(self.pc, source=pulse_gid, target=e0, syntype='ex', weight=weight, delay=1)
                            
                            # AMPA
                            self.cells[n][i].create_synapses(n_ex=1, tau1=tau11[n], tau2=tau12[n])
                            e0 = len(self.cells[n][i].synlist)-1
                    
                            weight=weight1[n]
                            np.random.seed(gid*60)
                            #weight = np.random.normal(weight, weight*0.5, 1).clip(min=0)
                            self.cells[n][i].pconnect_target(self.pc, source=pulse_gid, target=e0, syntype='ex', weight=weight, delay=1)
                            
                
                modulation = (t_input, mod) # mack to s!
                modulation_vec.append(modulation)            
                        
            return modulation_vec
        
        
    def connect_Synapse(self, pulse_gid, nt, i, n, gid, j, syntype = "ex", nsyn=0):   
        
        if self.do_run:
            
            if 'gsyn_in' in self.method_interpol:        
                if isinstance(self.factor_celltype[nt], ( int, long ) ):
                    f = self.factor_celltype[nt] 
                else:
                    f = self.factor_celltype[nt][0]       
 
            if syntype == "ex":
            
                # each cell can receive different g_syn_ex !  
                if type(self.g_syn_ex[nt]) is ndarray:
                    if len(self.g_syn_ex[nt]) == len(self.global_gidlist[nt]):
                        w = self.g_syn_ex[nt][n]
                    else:
                        w = self.g_syn_ex[nt] 
                else:
                    w = self.g_syn_ex[nt] 
    
                seed = int(10000 + 10*gid + j)
                np.random.seed(seed*41)
                
                if self.g_syn_ex_s[nt] > 0:
                    w = np.random.normal(w, w*self.g_syn_ex_s[nt], 1).clip(min=0)  #  self.g_syn_ex_s[nt]                          
            
                if self.celltype[nt] == 'Grc':
                
                    # delete old
                    if j == 0: 
                        self.cells[nt][i].MF_L = []
                        self.cells[nt][i].mfncpc = []
                    
                    if "gr" not in str(self.tau1_ex[nt]):
    
                        if "amfit" in str(self.tau1_ex[nt]):
                            syn = h.ExpZSyn(self.cells[nt][i].soma(0.5)) 
                            
                            syn.tau1_ampa = 0.254
                            syn.tau2_ampa = 0.254
                            syn.tau3_ampa = 0.363
                            syn.tau4_ampa = 6.523
                            syn.f1_ampa = 8.8376e-05
                            syn.f2_ampa = 5.5257e-05
                        
                            syn.f1_nmda = 0
                            
                        elif "nmfit" in str(self.tau1_ex[nt]):
                            syn = h.ExpYSyn(self.cells[nt][i].soma(0.5))
                            
                            syn.f1_ampa = 0
                            syn.f2_ampa = 0
                        
                            syn.tau1_nmda = 1.902
                            syn.tau2_nmda = 82.032
                            syn.f1_nmda = 7.853857483005277e-05
                            
                        elif "fit" in str(self.tau1_ex[nt]):    
                            syn = h.ExpGrcSyn(self.cells[nt][i].soma(0.5))
       
                            syn.tau1_ampa = 0.254
                            syn.tau2_ampa = 0.254
                            syn.tau3_ampa = 0.363
                            syn.tau4_ampa = 6.523
                            syn.f1_ampa = 8.8376e-05
                            syn.f2_ampa = 5.5257e-05
                        
                            syn.tau1_nmda = 1.902
                            syn.tau2_nmda = 82.032
                            syn.f1_nmda = 7.853857483005277e-05
                            
                        else:
                            tau1 = self.tau1_ex[nt]
                            tau2 = self.tau2_ex[nt]
                            
                            if tau1 == 0:
                                syn = h.ExpSyn(self.cells[nt][i].soma(0.5))
                                syn.tau = tau2/ms
                    
                            else:                
                                syn = h.Exp2Syn(self.cells[nt][i].soma(0.5))
                                syn.tau1 = tau1/ms
                                syn.tau2 = tau2/ms
                    
                        syn.e = 0/mV
                        
                        self.cells[nt][i].MF_L.append(syn)
                        
                        e0 = len(self.cells[nt][i].MF_L)-1 # get number of current synapse!
                        
                        syn_idx = int(e0)
                        
                        source = int(pulse_gid)
                        self.cells[nt][i].mfncpc.append(self.pc.gid_connect(source, self.cells[nt][i].MF_L[syn_idx]))
                        self.cells[nt][i].mfncpc[-1].delay = 1
                        self.cells[nt][i].mfncpc[-1].weight[0] = w
                        
                        if 'gsyn_in' in self.method_interpol:
                            self.record_syn.append(h.Vector())
                            self.record_syn[-1].record(self.cells[nt][i].MF_L[-1]._ref_g)
                            self.gsyn_in_fac.append(f)
                          
                    else:
                        
                        nrel = 0
                        
                        if "stoch" in str(self.tau1_ex[nt]):
                            nrel = 4
                        
                        self.cells[nt][i].createsyn(nmf = 1, ngoc = 0, weight_gmax = w, nrel=nrel) 
                        
                        if "ampa" in str(self.tau1_ex[nt]):
                            self.cells[nt][i].MF_L[-1].postsyns['NMDA'][0].gmax_factor = 0
                            if "nopre" in str(self.tau1_ex[nt]):
                                print "- no pre"
                                self.cells[nt][i].MF_L[-1].postsyns['AMPA'][0].tau_rec = 1e-9
                                self.cells[nt][i].MF_L[-1].postsyns['AMPA'][0].tau_facil  = 1e-9
                                self.cells[nt][i].MF_L[-1].postsyns['AMPA'][0].tau_1 = 0
                            
                        if "nostdampa" in str(self.tau1_ex[nt]):
                            self.cells[nt][i].MF_L[-1].postsyns['NMDA'][0].gmax_factor = 0
                            self.cells[nt][i].MF_L[-1].postsyns['AMPA'][0].tau_rec = 1e-9
                            self.cells[nt][i].MF_L[-1].postsyns['AMPA'][0].tau_facil  = 1e-9
                            self.cells[nt][i].MF_L[-1].postsyns['AMPA'][0].tau_1 = 0
                            self.cells[nt][i].MF_L[-1].postsyns['AMPA'][0].r6FIX = 0
                            
                        if "nostdnmda" in str(self.tau1_ex[nt]):
                            self.cells[nt][i].MF_L[-1].postsyns['AMPA'][0].gmax_factor = 0
                            self.cells[nt][i].MF_L[-1].postsyns['NMDA'][0].tau_rec = 1e-9
                            self.cells[nt][i].MF_L[-1].postsyns['NMDA'][0].tau_facil  = 1e-9
                            self.cells[nt][i].MF_L[-1].postsyns['NMDA'][0].tau_1 = 0
                            self.cells[nt][i].MF_L[-1].postsyns['NMDA'][0].RdRate = 0	
                            
                        if "nmda" in str(self.tau1_ex[nt]):
                            self.cells[nt][i].MF_L[-1].postsyns['AMPA'][0].gmax_factor = 0
                            if "nopre" in str(self.tau1_ex[nt]):
                                self.cells[nt][i].MF_L[-1].postsyns['NMDA'][0].tau_rec = 1e-9
                                self.cells[nt][i].MF_L[-1].postsyns['NMDA'][0].tau_facil  = 1e-9
                                self.cells[nt][i].MF_L[-1].postsyns['NMDA'][0].tau_1 = 0
                        
                        if "nostdgr" in str(self.tau1_ex[nt]):
                            self.cells[nt][i].MF_L[-1].postsyns['AMPA'][0].r6FIX	= 0 #1.12	
                            self.cells[nt][i].MF_L[-1].postsyns['NMDA'][0].RdRate = 0 #12e-3
                            print "- no std"
                        
                        if "nomggr" in str(self.tau1_ex[nt]):    
                            self.cells[nt][i].MF_L[-1].postsyns['NMDA'][0].v0_block = -1e9
                            print "- no mg block"
                        
                        e0 = len(self.cells[nt][i].MF_L)-1 # get number of current synapse!
                        
                        self.cells[nt][i].pconnect(self.pc,pulse_gid,int(e0),'mf') 
                        
                        if 'gsyn_in' in self.method_interpol:
                            self.record_syn.append(h.Vector())
                            self.record_syn[-1].record(self.cells[nt][i].MF_L[-1].postsyns['AMPA'][0]._ref_g)
                            self.record_syn.append(h.Vector())
                            self.record_syn[-1].record(self.cells[nt][i].MF_L[-1].postsyns['NMDA'][0]._ref_g)
                            self.gsyn_in_fac.append(f)
                            self.gsyn_in_fac.append(f)
    
                
                elif self.celltype[nt] == 'Goc':
                
                    # delete old
                    if j == 0: 
                        self.cells[nt][i].MF_L = []
                        self.cells[nt][i].mfncpc = []
                    
                    if "go" not in str(self.tau1_ex[nt]):
    
                        tau1 = self.tau1_ex[nt]
                        tau2 = self.tau2_ex[nt]
                            
                        if tau1 == 0:
                            syn = h.ExpSyn(self.cells[nt][i].soma(0.5))
                            syn.tau = tau2/ms
                
                        else:                
                            syn = h.Exp2Syn(self.cells[nt][i].soma(0.5))
                            syn.tau1 = tau1/ms
                            syn.tau2 = tau2/ms
                    
                        syn.e = 0/mV
                        
                        self.cells[nt][i].MF_L.append(syn)
                        
                        e0 = len(self.cells[nt][i].MF_L)-1 # get number of current synapse!
                        
                        syn_idx = int(e0)
                        
                        source = int(pulse_gid)
                        self.cells[nt][i].mfncpc.append(self.pc.gid_connect(source, self.cells[nt][i].MF_L[syn_idx]))
                        self.cells[nt][i].mfncpc[-1].delay = 1
                        self.cells[nt][i].mfncpc[-1].weight[0] = w
                        
                        if 'gsyn_in' in self.method_interpol:
                            self.record_syn.append(h.Vector())
                            self.record_syn[-1].record(self.cells[nt][i].MF_L[-1]._ref_g)
                            self.gsyn_in_fac.append(f)
                    else:
                        
                        nrel = 0
                        
                        mg = self.mglufac_ex[0]
                        if self.mglufac_ex[1] > 0:
                            mg = np.random.normal(self.mglufac_ex[0], self.mglufac_ex[1]*self.mglufac_ex[0], 1).clip(min=0)  #  self.g_syn_ex_s[nt]              
                        
                        if "stoch" in str(self.tau1_ex[nt]):
                            nrel = 4
                        
                        self.cells[nt][i].createsyn(nmf = 1, weight_gmax = w, nrel=nrel, mglufac = mg) 
                        
                        e0 = len(self.cells[nt][i].MF_L)-1 # get number of current synapse!
                        
                        self.cells[nt][i].pconnect(self.pc,pulse_gid,int(e0),'mf') 
                        
                        if 'gsyn_in' in self.method_interpol:
                            self.record_syn.append(h.Vector())
                            self.record_syn[-1].record(self.cells[nt][i].MF_L[-1].postsyns['AMPA'][0]._ref_g)
                            self.record_syn.append(h.Vector())
                            self.record_syn[-1].record(self.cells[nt][i].MF_L[-1].postsyns['NMDA'][0]._ref_g)
                            self.gsyn_in_fac.append(f)
                            self.gsyn_in_fac.append(f)
                            
                elif self.celltype[nt] == 'IfCell':  
                    
                    # delete old
                    if j == 0: 
                        self.cells[nt][i].synlist = []
                        self.cells[nt][i].nc = []
                    
                    if "gr" in str(self.tau1_ex[nt]):
                        
                        self.cells[nt][i].whatami = "grc"
                        
                        nrel = 0
                        if "stoch" in str(self.tau1_ex[nt]):
                            nrel = 4
                            
                        self.cells[nt][i].MF_L = self.cells[nt][i].synlist
                        self.cells[nt][i].synlist.append(Synapse('glom', self.cells[nt][i], self.cells[nt][i].soma, nrel=nrel, record_all=0, weight_gmax = w))
                        
                        if "ampa" in str(self.tau1_ex[nt]):
                            self.cells[nt][i].synlist[-1].postsyns['NMDA'][0].gmax_factor = 0
                            if "nopre" in str(self.tau1_ex[nt]):
                                print "- no pre"
                                self.cells[nt][i].synlist[-1].postsyns['AMPA'][0].tau_rec = 1e-9
                                self.cells[nt][i].synlist[-1].postsyns['AMPA'][0].tau_facil  = 1e-9
                                self.cells[nt][i].synlist[-1].postsyns['AMPA'][0].tau_1 = 0
                            
                        if "nmda" in str(self.tau1_ex[nt]):
                            self.cells[nt][i].synlist[-1].postsyns['AMPA'][0].gmax_factor = 0
                            if "nopre" in str(self.tau1_ex[nt]):
                                self.cells[nt][i].synlist[-1].postsyns['NMDA'][0].tau_rec = 1e-9
                                self.cells[nt][i].synlist[-1].postsyns['NMDA'][0].tau_facil  = 1e-9
                                self.cells[nt][i].synlist[-1].postsyns['NMDA'][0].tau_1 = 0
                            
                        if "nostdampa" in str(self.tau1_ex[nt]):
                            self.cells[nt][i].synlist[-1].postsyns['AMPA'][0].tau_rec = 1e-9
                            self.cells[nt][i].synlist[-1].postsyns['AMPA'][0].tau_facil  = 1e-9
                            self.cells[nt][i].synlist[-1].postsyns['AMPA'][0].tau_1 = 0
                            self.cells[nt][i].synlist[-1].postsyns['AMPA'][0].r6FIX	= 0 #1.12	
                            
                        if "nostdnmda" in str(self.tau1_ex[nt]):
                            self.cells[nt][i].synlist[-1].postsyns['NMDA'][0].tau_rec = 1e-9
                            self.cells[nt][i].synlist[-1].postsyns['NMDA'][0].tau_facil  = 1e-9
                            self.cells[nt][i].synlist[-1].postsyns['NMDA'][0].tau_1 = 0
                            self.cells[nt][i].synlist[-1].postsyns['NMDA'][0].RdRate = 0	
                        
                        if "nostdgr" in str(self.tau1_ex[nt]):
                            self.cells[nt][i].synlist[-1].postsyns['AMPA'][0].r6FIX	= 0 #1.12	
                            self.cells[nt][i].synlist[-1].postsyns['NMDA'][0].RdRate = 0 #12e-3
                            print "- no std"
                            
                        if "nomggr" in str(self.tau1_ex[nt]):    
                            self.cells[nt][i].synlist[-1].postsyns['NMDA'][0].v0_block = -1e9 #.k_block  = 1e-9
                            print "- no mg block"
                        
                        e0 = len(self.cells[nt][i].synlist)-1
                        syn_idx = int(e0)
                    
                        source = int(pulse_gid)
                        self.cells[nt][i].nc.append(self.pc.gid_connect(source, self.cells[nt][i].synlist[syn_idx].input))
                        self.cells[nt][i].nc[-1].delay = 1
                        self.cells[nt][i].nc[-1].weight[0] = 1
                        
                        if 'gsyn_in' in self.method_interpol:
                            self.record_syn.append(h.Vector())
                            self.record_syn[-1].record(self.cells[nt][i].synlist[syn_idx].postsyns['AMPA'][0]._ref_g)
                            self.record_syn.append(h.Vector())
                            self.record_syn[-1].record(self.cells[nt][i].synlist[syn_idx].postsyns['NMDA'][0]._ref_g) 
                            self.gsyn_in_fac.append(f)
                            self.gsyn_in_fac.append(f)
                    else:
                        
                        if "amfit" in str(self.tau1_ex):
                            
                            syn = h.ExpGrcSyn(self.cells[nt][i].soma(0.5)) 
                            
                            syn.tau1_ampa = 0.254
                            syn.tau2_ampa = 0.254
                            syn.tau3_ampa = 0.363
                            syn.tau4_ampa = 6.523
                            syn.f1_ampa = 8.8376e-05
                            syn.f2_ampa = 5.5257e-05
                        
                            syn.f1_nmda = 0
                            
                            self.cells[nt][i].synlist.append(syn) # synlist is defined in Cell 
                            
                        elif "nmfit" in str(self.tau1_ex):
                            
                            syn = h.ExpGrcSyn(self.cells[nt][i].soma(0.5))
                            
                            syn.f1_ampa = 0
                            syn.f2_ampa = 0
                        
                            syn.tau1_nmda = 1.902
                            syn.tau2_nmda = 82.032
                            syn.f1_nmda = 7.853857483005277e-05
                            
                            self.cells[nt][i].synlist.append(syn) # synlist is defined in Cell 
                            
                        elif "fit" in str(self.tau1_ex):    
                            
                            syn = h.ExpGrcSyn(self.cells[nt][i].soma(0.5))
       
                            syn.tau1_ampa = 0.254
                            syn.tau2_ampa = 0.254
                            syn.tau3_ampa = 0.363
                            syn.tau4_ampa = 6.523
                            syn.f1_ampa = 8.8376e-05
                            syn.f2_ampa = 5.5257e-05
                        
                            syn.tau1_nmda = 1.902
                            syn.tau2_nmda = 82.032
                            syn.f1_nmda = 7.853857483005277e-05 
                            
                            self.cells[nt][i].synlist.append(syn) # synlist is defined in Cell 
                        
                        else:
                            
                            self.cells[nt][i].create_synapses(n_ex=1, tau1=self.tau1_ex[nt], tau2=self.tau2_ex[nt])
                            
                        
                        e0 = len(self.cells[nt][i].synlist)-1
                        syn_idx = int(e0)
                        
                        self.cells[nt][i].pconnect_target(self.pc, source=pulse_gid, target=int(e0), syntype='ex', weight=w, delay=1)
                        
                        if 'gsyn_in' in self.method_interpol:
                            self.record_syn.append(h.Vector())
                            self.record_syn[-1].record(self.cells[nt][i].synlist[syn_idx]._ref_g)
                            self.gsyn_in_fac.append(f)
                
                elif self.celltype[nt] == 'Prk':
                
                    # delete old
                    if j == 0: 
                        self.cells[nt][i].PF_Lsync = []
                        self.cells[nt][i].spk_nc_pfsync = []
                        self.cells[nt][i].pfrand = []
                        
                        m = len(self.cells[nt][i].dendrange)
                        
                        seed = int(4*gid)
                        np.random.seed(seed)
                
                        for k in xrange(nsyn):
                            m -= 1
                	    mi = np.random.randint(0, m)	    
                	    self.cells[nt][i].dendrange[mi], self.cells[nt][i].dendrange[m] = self.cells[nt][i].dendrange[m], self.cells[nt][i].dendrange[mi]
                	    self.cells[nt][i].pfrand.append(self.cells[nt][i].dendrange[m])
                     
                        #print self.cells[nt][i].pfrand

                    if "prk" not in str(self.tau1_ex[nt]):
                        pass
                    else:
                        self.cells[nt][i].PF_Lsync.append(Synapse2('pf',self.cells[nt][i],self.cells[nt][i].pfrand[j],record_all=0))

                        e0 = len(self.cells[nt][i].PF_Lsync)-1 # get number of current synapse!
                        syn_idx = int(e0)

                        self.cells[nt][i].spk_nc_pfsync.append(self.pc.gid_connect(pulse_gid, self.cells[nt][i].PF_Lsync[syn_idx].input.newnetstim))
                        self.cells[nt][i].spk_nc_pfsync[-1].delay = 1
                        self.cells[nt][i].spk_nc_pfsync[-1].weight[0] = 1
                        
                        if 'gsyn_in' in self.method_interpol:
                            self.record_syn.append(h.Vector())
                            self.record_syn[-1].record(self.cells[nt][i].PF_Lsync[-1].postsyns['AMPA'][0]._ref_g)
                            self.gsyn_in_fac.append(f)  
            
            elif syntype == "inh":
                
                w = self.g_syn_inh[nt]
                
                seed = int(10000 + 10*gid + j)
                np.random.seed(seed*42)
                
                if self.g_syn_inh_s[nt] > 0:
                    w = np.random.normal(w, w*self.g_syn_inh_s[nt], 1).clip(min=w*0.1)  # self.g_syn_inh_s[nt]               
    
                if self.celltype[nt] == 'Grc':
        
                    if j == 0: 
                        self.cells[nt][i].GOC_L = []
                        self.cells[nt][i].gocncpc = []
                    
                    if "gr" not in str(self.tau1_inh[nt]):
    
                        tau1 = self.tau1_inh[nt]
                        tau2 = self.tau2_inh[nt]
    
                        if tau1 == 0:
                            syn = h.ExpSyn(self.cells[nt][i].soma(0.5))
                            syn.tau = tau2/ms
                
                        else:                
                            syn = h.Exp2Syn(self.cells[nt][i].soma(0.5))
                            syn.tau1 = tau1/ms
                            syn.tau2 = tau2/ms
                
                        syn.e = -65
                        
                        self.cells[nt][i].GOC_L.append(syn)
                        
                        i0 = len(self.cells[nt][i].GOC_L)-1 # get number of current synapse!
                        
                        syn_idx = int(i0)
                        source = int(pulse_gid)
                        self.cells[nt][i].gocncpc.append(self.pc.gid_connect(source, self.cells[nt][i].GOC_L[syn_idx]))
                        self.cells[nt][i].gocncpc[-1].delay = 1
                        self.cells[nt][i].gocncpc[-1].weight[0] = w
                          
                    else:
                        
                        self.cells[nt][i].createsyn(nmf = 0, ngoc = 1, weight_gmax = w) 
                        i0 = len(self.cells[nt][i].GOC_L)-1 # get number of current synapse!
                        self.cells[nt][i].pconnect(self.pc,pulse_gid,int(i0),'goc')
                        
                
                if self.celltype[nt] == 'IfCell': 
                    
                    if j == 0: 
                        self.cells[nt][i].synlist_inh = []
                        self.cells[nt][i].nc_inh = []
                    
                    if "gr" in str(self.tau1_inh[nt]):
                        
                        nrel = 0
                        if "stoch" in str(self.tau1_ex[nt]):
                            nrel = 4
                        
                        self.cells[nt][i].GOC_L = self.cells[nt][i].synlist
                        self.cells[nt][i].whatami = "grc"
                        self.cells[nt][i].synlist_inh.append(Synapse('goc', self.cells[nt][i], self.cells[nt][i].soma, nrel=nrel, record_all=0, weight_gmax = w))
                        
                        i0 = len(self.cells[nt][i].synlist_inh)-1
                        syn_idx = int(i0)
                        
                        source = int(pulse_gid)
                        self.cells[nt][i].nc_inh.append(self.pc.gid_connect(source, self.cells[nt][i].synlist_inh[syn_idx].input))
                        self.cells[nt][i].nc_inh[-1].delay = 1
                        self.cells[nt][i].nc_inh[-1].weight[0] = 1
                        
                        if "gaba" in str(self.tau1_ex[nt]):
                            
                            if 'gsyn_in' in self.method_interpol:
                                
                                if "nostdgaba" in str(self.tau1_ex[nt]):
                                    
                                    self.cells[nt][i].synlist_inh[syn_idx].postsyns['GABA'][0].tau_rec = 1e-9 
                                    self.cells[nt][i].synlist_inh[syn_idx].postsyns['GABA'][0].tau_facil = 1e-9       
                                    self.cells[nt][i].synlist_inh[syn_idx].postsyns['GABA'][0].tau_1 = 0      
                                    self.cells[nt][i].synlist_inh[syn_idx].postsyns['GABA'][0].d3 = 0      
                                    self.cells[nt][i].synlist_inh[syn_idx].postsyns['GABA'][0].d1d2 = 0      
                                    self.cells[nt][i].synlist_inh[syn_idx].postsyns['GABA'][0].d1 = 0      
                                    self.cells[nt][i].synlist_inh[syn_idx].postsyns['GABA'][0].d2 = 0      
                                    self.cells[nt][i].synlist_inh[syn_idx].postsyns['GABA'][0].d3_a6 = 0      
                                    self.cells[nt][i].synlist_inh[syn_idx].postsyns['GABA'][0].d1d2_a6 = 0      
                                    self.cells[nt][i].synlist_inh[syn_idx].postsyns['GABA'][0].d1_a6 = 0      
                                    self.cells[nt][i].synlist_inh[syn_idx].postsyns['GABA'][0].d2_a6 = 0        
                                    
                                self.record_syn.append(h.Vector())
                                self.record_syn[-1].record(self.cells[nt][i].synlist_inh[syn_idx].postsyns['GABA'][0]._ref_g)
                                self.gsyn_in_fac.append(f)
                                                    
                    else:
                        
                        self.cells[nt][i].create_synapses(n_inh=1, tau1_inh=self.tau1_inh[nt], tau2_inh=self.tau2_inh[nt], e_inh=-65)
                        i0 = len(self.cells[nt][i].synlist_inh)-1
                        syn_idx = int(i0)
                        self.cells[nt][i].pconnect_target(self.pc, source=pulse_gid, target=int(i0), syntype='inh', weight=w, delay=1)
                        
                        
            elif syntype == "intr":
                
                if self.celltype[nt] == 'Prk':
                
                    pass

            
    def set_SynPlay(self, farray, tarray, N = [], t_startstop = [], amode = 1):
        
        if self.do_run:
            
            delay = 1
            if (self.use_pc is False):
                delay = 0.1
                                    
            if N == []:
                N = self.N
                            
            self.pulse_list = [] 
            self.global_pulse_list = []
            self.global_pulse_list_inh = []
            self.global_pulse_list_intr = []
    
            f_cells_mean_local = []
            f_cells_cv_local = []
            f_cells_std_local = []
            
            for nt in range(self.n_celltypes): # loop over all cells
            
                if (self.n_syn_ex[nt] > 0) or (self.n_syn_inh[nt] > 0) or (self.n_syn_intr[nt] > 0):

                    local_gid_count = 0
                    local_gid_count_type = []
                
                
                    # EXCITATION
                    if str(type(self.g_syn_ex[nt] )) is not ndarray: self.g_syn_ex[nt]  = np.array([self.g_syn_ex[nt] ]) # each cell can receive different g_syn_ex !
                    
                    if len(self.g_syn_ex[nt]) == len(self.global_gidlist[nt]):
                        pass
                    else:
                        self.g_syn_ex[nt] = np.ones(len(self.global_gidlist[nt]))*self.g_syn_ex[nt][0]
                        #print "- single value in g_syn_ex, cells:", len(self.global_gidlist[nt])
                    
                    self.global_pulse_list.append([])
                    for ns in range(self.n_syn_ex[nt]): # loop over all excitatory synapses!
                        self.global_pulse_list[-1].append([])
                        for n in range(self.syn_max_mf[nt]): # number of cells of this celltype
                            self.global_pulse_list[-1][-1].append(local_gid_count+self.gid_count)
                            local_gid_count += 1
                            local_gid_count_type.append([])
                            local_gid_count_type[-1].append('ex')
                            local_gid_count_type[-1].append(n) # number of cell within their population 0..N[nt]
                            local_gid_count_type[-1].append(ns) # number of synapse 
                            
                            
                    # INHIBITION                        
                    if np.array(self.inh_hold[nt]).size <= 1:
                        self.inh_hold[nt] = np.ones(len(self.global_gidlist[nt]))*self.inh_hold[nt]
                        #print "- single value in inh_hold", self.inh_hold[nt] 
        
                        
                    self.global_pulse_list_inh.append([])
                    for ns in range(self.n_syn_inh[nt]): # loop over all inhibitory synapses!
                        self.global_pulse_list_inh[-1].append([])
                        for n in range(self.syn_max_inh[nt]): # number of cells of this celltype
                            self.global_pulse_list_inh[-1][-1].append(local_gid_count+self.gid_count)
                            local_gid_count += 1
                            local_gid_count_type.append([])
                            local_gid_count_type[-1].append('inh')
                            local_gid_count_type[-1].append(n) # number of cell within their population 0..N[nt]
                            local_gid_count_type[-1].append(ns) # number of synapse  

                    
                    # INTRUDER SYNAPSE
                    if str(type(self.g_syn_intr[nt] )) is not ndarray: self.g_syn_intr[nt]  = np.array([self.g_syn_intr[nt] ]) # each cell can receive different g_syn_intr !
                    
                    if len(self.g_syn_intr[nt]) == len(self.global_gidlist[nt]):
                        pass 
                    else:
                        self.g_syn_intr[nt] = np.ones(len(self.global_gidlist[nt]))*self.g_syn_intr[nt][0]
                        #print "- single value in g_syn_intr, cells:", len(self.global_gidlist[nt])
                    
                    self.global_pulse_list_intr.append([])
                    for ns in range(self.n_syn_intr[nt]): # loop over all intruding synapses!
                        self.global_pulse_list_intr[-1].append([])
                        for n in range(self.syn_max_intr[nt]): # number of generators for this celltype
                            self.global_pulse_list_intr[-1][-1].append(local_gid_count+self.gid_count)
                            local_gid_count += 1
                            local_gid_count_type.append([])
                            local_gid_count_type[-1].append('intr')
                            local_gid_count_type[-1].append(n) # number of cell within their population 0..N[nt]
                            local_gid_count_type[-1].append(ns) # number of synapse 
                            
                    
                    t_vec_input = np.array([]) # input trains 
                    id_vec_input = np.array([]) # input trains id
                    fs = 1 / self.dt
                    ih_use_v = []
                    
                    for i in range(int(self.id), local_gid_count, int(self.nhost)): # loop over all train generators and generate them
                    
                        self.pulse_list.append(i+self.gid_count)
                        pulse_gid = self.pulse_list[-1] 
                        gid = local_gid_count_type[i][1] # should correspond to this gid when multiple values inserted
                        
                        if local_gid_count_type[i][0] == 'ex':
                            
                            seed = int(10001 + pulse_gid) #  unique gid for generators! 
                            np.random.seed(seed*423)
                            
                            if self.ihold_sigma[nt] > 0:
                                ih_use = np.random.normal(self.ihold[nt], self.ihold[nt]*self.ihold_sigma[nt], 1).clip(min=0) # self.ihold[nt]*self.ihold_sigma[nt] 
                            
                            elif self.ihold_sigma[nt] < 0:
                                ih_use = np.random.uniform(0.1, self.ihold[nt])
                              
                            else:
                                ih_use = self.ihold[nt]                    
                            
                            ih_use_v.append(ih_use)
                            
                            if ih_use > 0:
                                # train has to be contructed here, to insert different train into each "dendrite"
                                ## different ihold has to be implemented here!!
                                iholdvec = concatenate((zeros(round(fs)), ones(round(len(tarray) - 1 * fs)) * ih_use))
    
                                if isinstance(self.syn_ex_dist[nt], ( tuple ) ): # distribution of amplitude, only one noise source!
                                    
                                    np.random.seed(pulse_gid*40)
                                    if self.syn_ex_dist[nt][1] > 0:
                                        f = np.random.normal(self.syn_ex_dist[nt][0], self.syn_ex_dist[nt][1], 1).clip(min=0)
                                    else:
                                        f = self.syn_ex_dist[nt][0]
                                        
                                    f2 = f
                                    rnd.seed(pulse_gid*40) # use gid so type 1, 2 is identical for each cell
                                    #rnd.seed(gid*40) # use gid so type 1, 2 is identical for each cell
                                    if self.syn_ex_dist[nt][2] > 0: # add inverted input with 50% probability, in future versions this will indicate the propability for -1 and 1                                    
                                        f2 = rnd.sample([-1,1],1)[0] * f
                                        #f2 = f
    
                                    if amode == 1:
                                        inamp = (f2 * self.amod[nt] * ih_use)
                                    elif amode == 2:
                                        inamp = (f2 * self.amod[nt] * self.ihold[nt])  
                                         
                                    modulation = (tarray, inamp * farray[0] + iholdvec)
                                    
                                    #if self.id == 0: print "- randomize play stimulus height, pulse_gid=", pulse_gid, " gid=", gid ," f=", f   
                                    if (gid==0): print "- randomize play stimulus height, pulse_gid=", pulse_gid, " gid=", gid ," f2=", f2,"inamp=",inamp   
                                                          
                                    #rnd.seed(local_gid_count_type[i][1]*300) # pick seed based on number of cell
                                    #nj = rnd.sample(range(len(farray)),1)[0] 
                                    nj = 1
                                                        
                                else: # different noise sources can be used at different synapses, linear combination test in openloop
                              
                                    nj = self.syn_ex_dist[nt][local_gid_count_type[i][2]]
                                    
                                    if nj == 0:
                                        modulation = (tarray, iholdvec)
                                    else:
                                        if amode == 1:
                                            inamp = (self.factor_celltype[nt] * self.amod[nt] * ih_use)
                                        elif amode == 2:
                                            inamp = (self.factor_celltype[nt] * self.amod[nt] * self.ihold[nt])  

                                        modulation = (tarray, inamp * farray[nj-1] + iholdvec)
                                        if self.id == 0: print "ex farray number:", nj-1, "ih_use:", ih_use, "self.amod[nt]:", self.amod[nt], "inamp: ", inamp
    
                                
                                # will be done n_syn_ex * number of cells!
                                if self.noise_syn_tau[nt] < 0: # variable threshold
                                    no = self.noise_syn[nt]
                                else:                   
                                    no = self.noise_syn[nt]*ih_use

                                train, self.n_train_ex = mod_spike_train(modulation, noise = no, seed = seed, noise_tau = self.noise_syn_tau[nt], noise_a = self.noise_a[nt]) 
                                
                                #plt.figure("input")
                                #plt.plot(train, train*0, '|')
                                #plt.show()
                                
                                t_vec_input = np.append(t_vec_input, train*ms).flatten() # use ms to save!!
                                id_vec_input = np.append(id_vec_input, np.ones(len(train))*pulse_gid).flatten()
                                
                                f_cells_mean_local0, f_cells_cv_local0, f_cells_std_local0 = self.calc_fmean(train*ms, t_startstop)
                                f_cells_mean_local.append(f_cells_mean_local0); f_cells_cv_local.append(f_cells_cv_local0); f_cells_std_local.append(f_cells_std_local0)
                                
                                if self.id == 0: print "TRAIN: requ. mean:", ih_use ,"eff. mean:", f_cells_mean_local0, "cv: " , f_cells_cv_local0, "std:" , f_cells_std_local0
                            
                            else:
                                train = []
                                self.n_train_ex = []
                                


                        elif local_gid_count_type[i][0] == 'intr':
                                    
                            # train has to be contructed here, to insert different train into each "dendrite"
                            nj = 0
                            
                            seed = int(10001 + pulse_gid)
                            np.random.seed(seed*4411)
                            
                            if self.intr_hold_sigma[nt] > 0:                    
                                ih_use = np.random.normal(self.intr_hold[nt], self.intr_hold[nt]*self.intr_hold_sigma[nt], 1).clip(min=0) 
                            else:
                                ih_use = self.intr_hold[nt]
                                
                            ih_use_v.append(ih_use)
                            
                            if ih_use > 0:     
                                                
                                iholdvec = concatenate((zeros(round(fs)), ones(round(len(tarray) - 1 * fs)) * ih_use))
                                modulation = (tarray, iholdvec)
                                                                    
                                # will be done n_syn_in * number of cells!                
                                if self.noise_syn_tau_intr[nt] < 0: # variable threshold
                                    no = self.noise_syn_intr[nt]
                                else:                   
                                    no = self.noise_syn_intr[nt]*ih_use
                                
                                if self.noise_syn_tau_intr[nt] >= -1:
                                    train, _ = mod_spike_train(modulation, noise = no, seed = seed, noise_tau = self.noise_syn_tau_intr[nt], noise_a = self.noise_a_intr[nt]) # train in ms
                                else:
                                    train = oscill_spike_train(sor = 4, spike_prob = 1/4, noise_fraction = 4, end_time = tarray[-1]/ms, seed = seed) 
                                
                                
                        elif local_gid_count_type[i][0] == 'inh':
                                    
                            # train has to be contructed here, to insert different train into each "dendrite"
            
                            seed = int(10001 + pulse_gid)
                                
                            np.random.seed(seed*44)
                            
                            if self.inh_hold_sigma[nt] > 0:                    
                                ih_use = np.random.normal(self.inh_hold[nt][gid], self.inh_hold[nt][gid]*self.inh_hold_sigma[nt], 1).clip(min=0) 
                            else:
                                ih_use = self.inh_hold[nt][gid]
                                
                                                
                            iholdvec = concatenate((zeros(round(fs)), ones(round(len(tarray) - 1 * fs)) * ih_use))
                            
                            nj = self.syn_inh_dist[nt][local_gid_count_type[i][2]]
                            if nj == 0:
                                modulation = (tarray, iholdvec)
                            else:
                                inamp = (self.amod[nt] * ih_use)
                                modulation = (tarray, inamp * farray[nj-1] + iholdvec)
                                #print "inh farray number:", nj-1, "ih_use:", ih_use, "amp: ", inamp #old: nj-1+nemax
                                
                            # will be done n_syn_in * number of cells!                
                            if self.noise_syn_tau_inh[nt] < 0: # variable threshold
                                no = self.noise_syn_inh[nt]
                            else:                   
                                no = self.noise_syn_inh[nt]*ih_use
                            
                            train, _ = mod_spike_train(modulation, noise = no, seed = seed, noise_tau = self.noise_syn_tau_inh[nt], noise_a = self.noise_a_inh[nt]) # train in ms
                            #print train
                        
                        #print train
                        if len(train) > 0:
                            if self.id == 0: 
                                print "-", pulse_gid, local_gid_count_type[i], "seed: ", seed, "ih_use:", ih_use, no, nj #, "first spike: ", train[0]             
                            self.setup_Play_train(train = train+self.inh_delay, input_gid = pulse_gid, delay = delay) # train in ms
                                    
                        
                    self.gid_count += local_gid_count # increase gid count
                    
                    self.barrier()
                    
                    for i, gid in enumerate(self.gidlist[nt]):  # for all input cells
                    
                        rnd.seed(gid*200)
                        n = self.global_gidlist[nt].index(gid) # index of cell within their population 0..N[nt]
                        # i is index on this node only!
                        
                        self.record_syn = []
                        for j in range(self.n_syn_ex[nt]):
                            if N[nt] == len(self.global_pulse_list[nt][j]):
                                pulse_gid = self.global_pulse_list[nt][j][n] #every cell of this type receives one pulse gid 
                                if self.id == 0: print "- gid:", gid ," n:", n ," one ex train for each synapse:", pulse_gid, "self.g_syn_ex[nt][n]:", self.g_syn_ex[nt][n] 
                            else:
                                pulse_gid = rnd.sample(self.global_pulse_list[nt][j],1)[0] # not enough, just pick one at random, for inh/f search only one synapse available!
                                if self.id == 0: print "- gid:", gid ," n:", n ," one ex train from", len(self.global_pulse_list[nt][j]), ":", pulse_gid, "self.g_syn_ex[nt][n]:", self.g_syn_ex[nt][n] 
                            
                            if "gaba" in str(self.tau1_ex[nt]):
                                self.connect_Synapse(pulse_gid, nt, i, n, gid, j, syntype = "inh") 
                            else:
                                self.connect_Synapse(pulse_gid, nt, i, n, gid, j, syntype = "ex", nsyn = self.n_syn_ex[nt]) 
                            
                        
                        if self.n_syn_inh[nt] > 0:
                            for j in range(self.n_syn_inh[nt]):
                                
                                if N[nt] == len(self.global_pulse_list_inh[nt][j]):
                                    pulse_gid = self.global_pulse_list_inh[nt][j][n] #every cell of this type receives one pulse gid   
                                    if self.id == 0: print "- one inh train for each synapse:", pulse_gid
                                else:
                                    pulse_gid = rnd.sample(self.global_pulse_list_inh[nt][j],1)[0] # not enough, just pick one at random 
                                    if self.id == 0: print "- one inh train from", len(self.global_pulse_list_inh[nt][j]), ":", pulse_gid
                                
                                self.connect_Synapse(pulse_gid, nt, i, n, gid, j, syntype = "inh") 
                                
                                
                        if self.n_syn_intr[nt] > 0:
                            for j in range(self.n_syn_intr[nt]):
                                
                                if N[nt] == len(self.global_pulse_list_intr[nt][j]):
                                    pulse_gid = self.global_pulse_list_intr[nt][j][n] #every cell of this type receives one pulse gid   
                                    if self.id == 0: print "- one intruding train for each synapse:", pulse_gid
                                else:
                                    pulse_gid = rnd.sample(self.global_pulse_list_intr[nt][j],1)[0] # not enough, just pick one at random 
                                    if self.id == 0: print "- one intruding train from", len(self.global_pulse_list_intr[nt][j]), ":", pulse_gid
                                
                                if (self.use_pc is False):
                                    
                                    if self.celltype[nt] == 'Prk': self.cells[nt][i].delrerun() 
                                    
                                    (msg,CF_input) = self.cells[nt][i].createsyn_CF(record_all=0,factor=self.g_syn_intr[nt][0],cf_setup_select='old')
                                    CF_input.number = 3 # three bursts
                                    CF_input.start = -0.3 # See synapsepfpurk.py
                                    CF_input.interval = 3 # 3 ms interval between bursts

                                    self.cells[nt][i].input_to_CF_nc.append(h.NetCon(self.vecstim[j], CF_input, 0, 0.1, 1))
                                    self.netcons.append(self.cells[nt][i].input_to_CF_nc[-1])
                                        
                                else:
                                    print "NOT IMPLEMENTED"
                                    
                                    
                    if self.id == 0: print "trains connected"
                    
                    if local_gid_count_type[i][0] == 'intr':
                        pass
                    else:
                        self.id_all_vec_input.append(self.do_gather(id_vec_input, dtype = 'i')) 
                        self.t_all_vec_input.append(self.do_gather(t_vec_input)) 
                    
                        f_cells_mean = self.do_gather(f_cells_mean_local) 
                        f_cells_cv = self.do_gather(f_cells_cv_local)  
                        f_cells_std = self.do_gather(f_cells_std_local) 
                    
                    self.fmean_input = np.nan
                    self.fmax_input = np.nan
                    self.fmstd_input = np.nan
                    self.fcvm_input = np.nan
                    self.fstdm_input = np.nan
                    
                    ih_use_v_all = self.do_gather(ih_use_v)
                                
                    if self.id == 0 and local_gid_count_type[i][0] != 'intr':
                        
                        self.fmean_input = mean(np.nan_to_num(f_cells_mean))  # compute mean of mean rate for all cells
                        self.fmstd_input = std(np.nan_to_num(f_cells_mean)) 
                        self.fmax_input = max(np.nan_to_num(f_cells_mean))
                        
                        self.fcvm_input = mean(f_cells_cv[~np.isnan(f_cells_cv)])
                        self.fstdm_input = mean(f_cells_std[~np.isnan(f_cells_std)])
                        
                        self.ih_use_max = max(ih_use_v_all)
                    
                        print "- trains, fmean: ",self.fmean_input, "fmax: ",self.fmax_input, "Hz", "fmstd: ",self.fmstd_input, "Hz", "fcvm: ",self.fcvm_input, "fstdm: ",self.fstdm_input, "Hz, ih_use_max:", self.ih_use_max        
                
                else:
                    self.global_pulse_list.append([])
                    self.global_pulse_list_inh.append([])
                  


    def do_gather(self, v_local, dtype = 'd'):
    
        if self.use_mpi:
            
            self.barrier()
            
            #v_local = v_local.astype(dtype).flatten()
            v_local = np.array(v_local, dtype=dtype).flatten() 
            
            if self.use_pc == False:

                v_global = None
                counts_local = np.array(len(v_local), dtype='i')
                
                counts = 0
                if self.id == 0:
                    counts = np.empty(self.nhost, dtype='i')
                    
                self.comm.Gather(sendbuf=[counts_local, MPI.INT], recvbuf=[counts, MPI.INT], root=0)
                
                if self.id == 0:
                    v_global = np.empty(sum(counts), dtype=dtype)
                    
    
                if dtype == 'd':
                    self.comm.Gatherv(sendbuf=[v_local, MPI.DOUBLE], recvbuf=[v_global, (counts, None), MPI.DOUBLE], root=0)
                elif dtype == 'i':
                    self.comm.Gatherv(sendbuf=[v_local, MPI.INT], recvbuf=[v_global, (counts, None), MPI.INT], root=0) 
                    
                #v_global = np.hstack(v_global)
                    
            else:
                sendlist = [None]*self.nhost                
                sendlist[0] =  v_local
                getlist = self.pc.py_alltoall(sendlist)
                    
                v_global = np.hstack(getlist)  
                
        else:
            
            v_global = np.hstack(v_local)
                
        return v_global
        

    def setup_Play_train(self, train = [], input_gid = 0, delay = 1):
    
        self.trains.append(train)

        # possibility to play spikes into the cells!
        self.vecstim.append(h.VecStim(.5))
        self.nc_vecstim.append(h.NetCon(self.vecstim[-1],None))
        self.nc_vecstim[-1].delay = delay

        self.spike_vec.append(h.Vector(self.trains[-1]))
        self.vecstim[-1].play(self.spike_vec[-1]) 

        if (self.use_mpi):
            self.pc.set_gid2node(input_gid, self.id)  # associate gid with this host
            self.pc.cell(input_gid,self.nc_vecstim[-1])  # associate gid with spike detector
 

    def record(self):
        """
        Initializes recording vectors. Internal function
        """

        if self.n_celltypes > 1:
            #print "self.n_borders:",self.n_borders
            for n in range(self.n_celltypes):
                if self.n_borders[n] in self.gidlist[n]:
                    #print "np.shape(self.rec_v):",np.shape(self.rec_v)
                    #print "np.shape(self.cells):",np.shape(self.cells)
                    self.rec_v[n].record(self.cells[n][0].soma(0.5)._ref_v)                       

        
        if self.id == 0:  # only for first node and first cell
            
            # Voltage
            self.rec_v[0].record(self.cells[self.a_celltype[0]][0].soma(0.5)._ref_v) 
            
            # Stimuli
            self.rec_i = h.Vector()

            if (self.plays != []): 
                if (isinstance(self.plays[0], list) is False): 
                    self.rec_i.record(self.plays[0]._ref_i)
                else:
                    self.rec_i.record(self.plays[0][0]._ref_i)    
            
            self.rec_ich = h.Vector()
            if self.ic_holds != [] and (isinstance(self.ic_holds[0], list) is False):            
                self.rec_ich.record(self.ic_holds[0]._ref_i)
            
            self.rec_ics = h.Vector()
            if self.ic_starts != []:       
                self.rec_ics.record(self.ic_starts[0]._ref_i)
            
            self.rec_n = h.Vector() 
            
            if self.fluct_s[0] > 0:
                # Fluctuating input        
                self.rec_n.record(self.flucts[0]._ref_i)
                print "recording noise"
            elif (len(self.flucts) > 0) and (len(self.fluct_g_i0)>0):
                self.rec_n.record(self.flucts[0]._ref_g_i)
                print "recording g noise"
            else:
                print "nonoise"
            
            if hasattr(self.cells[self.a_celltype[0]][0], 'lkg2_noise'):
                if self.cells[self.a_celltype[0]][0].lkg2_noise > 0:
                    self.rec_n.record(self.cells[self.a_celltype[0]][0].fluct._ref_il)
                    print "recording tonic gaba noise"                
                
            self.rec_step = h.Vector()
            if self.ic_steps != []:     
                self.rec_step.record(self.ic_steps[0]._ref_i)                
            
            # Time
            self.rec_t = h.Vector()
            self.rec_t.record(h._ref_t)
            
    
    def run(self, tstop = 10*s, do_loadstate = True):
        """
        Starts the stimulation.
        """
        self.record()
        
        if self.first_run:

            if self.use_mpi: self.pc.set_maxstep(100)
            #self.pc.spike_compress(1) #test
    
            if self.use_multisplit:
                import multiprocessing
                
                Hines = h.CVode()
                Hines.active(0)
                
                h.load_file("parcom.hoc")
                p = h.ParallelComputeTool()
                
                if self.use_mpi:
                    cpus = multiprocessing.cpu_count() #32 #self.pc.nhost()
                else:
                    cpus = multiprocessing.cpu_count() #32   
                    
                p.change_nthread(cpus,1)    
                p.multisplit(1)
                print "Using multisplit, cpus:", cpus
                
            else:
                
                h.load_file("stdrun.hoc")
            
            if self.use_local_dt:
                h.cvode.active(1) 
                h.cvode.use_local_dt(1)   
            
            h.celsius = self.temperature  
            h.dt = self.dt/ms  # Fixed dt
            h.steps_per_ms = 1 / (self.dt/ms)
            
            if self.cells[self.a_celltype[0]] != []: 
                if hasattr(self.cells[self.a_celltype[0]][0], 'v_init'):
                    h.v_init = self.cells[self.a_celltype[0]][0].v_init  # v_init is supplied by cell itself!
                else:    
                    h.v_init = -60    
                    
            h.stdinit()  

            h.finitialize()
    
            if hasattr(self.cells[self.a_celltype[0]][0], 'load_states') and do_loadstate:
                m = md5.new()
                cell_exe_new = self.cell_exe[0]
                m.update(cell_exe_new)
                filename = './states_' + self.celltype[0] + '_' + m.hexdigest() + '_Population.b'
                self.cells[self.a_celltype[0]][0].load_states(filename)
        
        else:

            pass            
            
            
        if self.id == 0:
            import time
            t0 = time.time()

        if self.simstep == 0:
            if self.id == 0: print "Running without steps",
            
            if self.use_mpi:
                self.pc.psolve(tstop/ms)
            else:
                h.init()
                h.tstop = tstop/ms
                h.run()

        else:
            
            h.finitialize()
            cnt = 1
        
            #if self.id == 50: 
            #    print len(self.cells[1][0].nc), self.cells[1][0].nc[0].weight[0]
            #    print len(self.cells[0][0].nc_inh), self.cells[0][0].nc_inh[0].weight[0]
        
            h.t = 0
            while h.t < tstop/ms:
            
                if self.id == 0:
                    print "Running...",
                    if self.use_mpi:
                        past_time = self.pc.time()
                    
                h.continuerun(cnt*self.simstep/ms)
                if self.use_mpi: self.pc.barrier()
                
                if self.id == 0:
                    if self.use_mpi:
                        print "Simulated time =",h.t*ms, "s, Real time = ", (self.pc.time()-past_time), 's'
                    else:
                        print "Simulated time =",h.t*ms, "s"
                
                #if self.id == 0:
                #    print hpy.heap().byrcs
                cnt += 1

        if self.id == 0: print "psolve took ", time.time() - t0, "seconds"
       
        self.first_run = False
        
        self.barrier()   # wait for other nodes

        self.tstop = tstop      
               
        
    def get(self, t_startstop=[], i_startstop=[], N = []):
        """
        Gets the recordings.
        """
        
        if N == []:
            N = self.N
        
        if t_startstop == []:
            t_startstop = np.array([2, self.tstop])
        
        t_all_vec = []
        id_all_vec = []
                
        fmean = []
        fbase = []
        fmax = []
        fmstd = []
        fcvm = []
        fstdm = []
        gid_del = []
        f_cells_mean_all = []
        f_cells_base_all = []
        f_cells_cv_all = [] 
        f_cells_std_all = []
        
        fmeanA = []
        fmstdA = []
        fmaxA = []
        fcvmA = []
        fstdmA = []
        fbaseA = []
        fbstdA = []
        
        if self.id == 0: print "start gathering spikes"
        
        for n in range(self.n_celltypes):

            if self.use_mpi: 
                
                self.barrier()   # wait for other node
                t_vec = np.array(self.t_vec[n]).flatten()*ms - 1*ms # shift time because of output delay
                id_vec = np.array(self.id_vec[n]).flatten()
                
            else:
                
                t_vec = np.array([])
                id_vec = np.array([])
                print np.shape(self.t_vec)
                for i in self.gidlist[n]:
                    t_vec0 = np.array(self.t_vec[n][i]).flatten()*ms 
                    t_vec = np.append(t_vec, t_vec0).flatten()
                    id_vec = np.append(id_vec, np.ones(len(t_vec0))*i).flatten()                    

            fmean0, fmax0, fmstd0, fcvm0, fstdm0, gid_del0, f_cells_mean_all0, f_cells_cv_all0, f_cells_std_all0, fbase0, f_cells_base_all0 = self.get_fmean(t_vec, id_vec, t_startstop = t_startstop, gidlist = self.gidlist[n]) 
            fmean.append(fmean0); fmax.append(fmax0), fmstd.append(fmstd0), fcvm.append(fcvm0), fstdm.append(fstdm0), gid_del.append(gid_del0), f_cells_mean_all.append(f_cells_mean_all0), f_cells_cv_all.append(f_cells_cv_all0), f_cells_std_all.append(f_cells_std_all0)
            fbase.append(fbase0); f_cells_base_all.append(f_cells_base_all0)
            
            t_all_vec.append(self.do_gather(t_vec))
            id_all_vec.append(self.do_gather(id_vec))
            
        if (self.id == 0) and (self.no_fmean == False): 
            f_cells_mean_all = np.array(f_cells_mean_all).flatten()
            fmeanA = mean(f_cells_mean_all)  # compute mean of mean rate for all cells
            fmstdA = std(f_cells_mean_all) 
            fmaxA = max(f_cells_mean_all)
            
            f_cells_base_all = np.array(f_cells_base_all).flatten()
            fbaseA = mean(f_cells_base_all)  # compute mean of mean rate for all cells
            fbstdA = std(f_cells_base_all)
            
            f_cells_cv_all = np.concatenate((np.array(f_cells_cv_all)))
            f_cells_std_all = np.concatenate((np.array(f_cells_std_all)))
            
            fcvmA = mean(f_cells_cv_all)
            fstdmA = mean(f_cells_std_all)
            
            print "- ALL, fmean: ",fmeanA, "fmax: ",fmaxA, "Hz", "fmstd: ",fmstdA, "Hz", "fcvm: ",fcvmA, "fstdm: ",fstdmA, "Hz", "fbase: ",fbaseA, "Hz", "fbstd: ", fbstdA, "Hz"
        
        if self.id == 0: print "all spikes have been gathered"

        self.barrier()

        # do this here to have something to return
        voltage = []
        current = []
        time = []
                
        freq_times = []
        spike_freq = []
        gsyn = []
                 
        if self.id == 0:  # only for first node
        
            time = np.array(self.rec_t)*ms

            # use self.bin_width as bin width!
            freq_times = arange(0, time[-1], self.bin_width)

            voltage.append(np.array(self.rec_v[0])*mV)
            current = np.zeros(len(time))

            if len(np.array(self.rec_ics)) > 0:
                current = current + np.array(self.rec_ics) 
                
            if len(np.array(self.rec_ich)) > 0:
                current = current + np.array(self.rec_ich)
                
            if len(np.array(self.rec_i)) > 0:
                current = current + np.array(self.rec_i)   
                
            if len(np.array(self.rec_n)) > 0:
                current = current + np.array(self.rec_n)  
                print np.array(self.rec_n) 
                
            if len(np.array(self.rec_step)) > 0:
                current = current + np.array(self.rec_step)  

        else:
            time = [0]
        
        self.barrier()
        time = self.broadcast(time, fast = True)

        gsyn_in = []
        gsyn_in0 = []
        
        if 'gsyn_in' in self.method_interpol:
            
            gsyn_in = None
            if self.id == 0: print "- collecting gsyn_in"
            gsyn_in0 = np.zeros(len(time), dtype='d')
            if self.record_syn is not []:
                for i, j in enumerate(self.record_syn):
                    gsyn_in0 = gsyn_in0 + self.gsyn_in_fac[i] * np.array(j, dtype='d') 
            
            if self.use_mpi:
                count = len(time)
                
                #if self.id == 0: gsyn_in = np.empty(count*self.nhost, dtype='d')
                #self.comm.Gatherv(sendbuf=[gsyn_in0, MPI.DOUBLE], recvbuf=[gsyn_in, MPI.DOUBLE], root=0)
                
                gsyn_in = self.do_gather(gsyn_in0)
                
                if self.id == 0:
                    gsyn_in = np.reshape(gsyn_in, (self.nhost,count))
                    gsyn_in = sum(gsyn_in,0)
                
            else:
                gsyn_in = gsyn_in0
                        
        self.barrier()   # wait for other nodes
        
        if self.n_celltypes > 1:
            if self.id == 0: print "more than one celltype send voltage of first other cell to root"
            
            for n in range(1, self.n_celltypes):
                
                if self.use_pc == True:
                    
                    srclist = [None]*self.nhost
                    
                    if (self.n_borders[n] in self.gidlist[n]):
                        srclist[0] = np.array(self.rec_v[n])*mV
                        
                    destlist = self.pc.py_alltoall(srclist) 
                    
                    if self.id == 0:
                        idx = [i for i, x in enumerate(destlist) if x is not None]
                        if len(idx) > 1: raise ValueError('Error, too many vectors sent, should be one at a time!')
                        voltage.append(np.array(destlist[idx[0]]))
                    
                else:
                    
                    if self.id == 0:
                        if (self.n_borders[n] in self.gidlist[n]):  # first node has it, do not wait to receive it!
                            v_temp = np.array(self.rec_v[n])*mV
                        else:
                            v_temp = np.zeros(len(voltage[0]))
                            self.comm.Recv([v_temp, MPI.DOUBLE], source = MPI.ANY_SOURCE, tag=int(sum(N)+33))
                        
                        voltage.append(v_temp)
                    else:
                        if self.n_borders[n] in self.gidlist[n]:
                            voltage = np.array(self.rec_v[n])*mV 
                            self.comm.Ssend([voltage, MPI.DOUBLE], dest=0, tag=int(sum(N)+33))

        self.barrier()   # wait for other nodes     

        times = arange(0, time[-1], 1*ms)        
        gsyns = []
        if self.called_syn_out_all == True:
            
            for n in range(self.n_celltypes):
                gsyns.append([])
            
                if self.use_pc == True:

                    for i, gid in enumerate(self.global_gidlist[n]): 
                    
                        srclist = [None]*self.nhost
                        
                        if gid in self.gidlist[n]: #only one node does this
                            a = np.array(self.cells[n][self.gidlist[n].index(gid)].record['gsyn'])
                            c = np.zeros(int((1*ms)/self.dt))
                            temp = np.append(a, c).flatten()
                            temp = temp[int((1*ms)/self.dt):len(temp)+1]
                            gtemp = interp(times,time,temp)
                            
                            srclist[0] = gtemp # send to root only
    
                        destlist = self.pc.py_alltoall(srclist)    
                        
                        if self.id == 0:
                            idx = [i for i, x in enumerate(destlist) if x is not None]
                            if len(idx) > 1: raise ValueError('Error, too many vectors sent, should be one at a time!')
                            gsyns[n].append(np.array(destlist[idx[0]]))
                
                else:    
                                
                    for i, gid in enumerate(self.global_gidlist[n]): 
                    
                        if self.id == 0:
                            if gid in self.gidlist[n]:
                                a = np.array(self.cells[n][self.gidlist[n].index(gid)].record['gsyn'])
                                c = np.zeros(int((1*ms)/self.dt))
                                temp = np.append(a, c).flatten()
                                temp = temp[int((1*ms)/self.dt):len(temp)+1]
                                gtemp = interp(times,time,temp)
    
                            else:
                                gtemp = np.zeros(len(times))
                                self.comm.Recv([gtemp, MPI.DOUBLE], source = MPI.ANY_SOURCE, tag=int(gid))
                            
                            gsyns[n].append(np.array(gtemp))
                            
                        else:
                            if gid in self.gidlist[n]:
                                a = np.array(self.cells[n][self.gidlist[n].index(gid)].record['gsyn'])
                                c = np.zeros(int((1*ms)/self.dt))
                                temp = np.append(a, c).flatten()
                                temp = temp[int((1*ms)/self.dt):len(temp)+1]
                                gtemp = interp(times,time,temp)                             
                                #np.array(self.cells[n][self.gidlist[n].index(gid)].record['gsyn'])
                                self.comm.Ssend([gtemp, MPI.DOUBLE], dest=0, tag=int(gid))
                            
            if self.id == 0: print "root gathered synaptic output conductance" 
            
            
        self.barrier()   # wait for other nodes      
        
        times = arange(0, time[-1], 10*ms)
        
        w_mat = []
        winh_mat = []
                
        if self.stdp_used == True:
                    
            for n in range(self.n_celltypes):
                w_mat.append([])                
                
                for i, gid in enumerate(self.global_gidlist[n]): 
                
                    if self.id == 0:
                        
                        wall = []
                        
                        if gid in self.gidlist[n]:

                            walltemp = self.cells[n][self.gidlist[n].index(gid)].record['w']                            
                            if len(walltemp) > 0:
                                for l in range(len(walltemp)):
                                    wtemp = np.array(walltemp[l])
                                    wtemp = interp(times,time,wtemp)
                                    wall.append(wtemp)
                                    
                        else:
                            
                            while 1:
                                wtemp = np.zeros(len(times))
                                self.comm.Recv([wtemp, MPI.DOUBLE], source = MPI.ANY_SOURCE, tag=int(gid))
                                
                                if wtemp[0] == -1:
                                    break
                                else:
                                    wall.append(wtemp)
                            
                        w_mat[n].append(wall)
                        
                    else:
                        if gid in self.gidlist[n]:
                            walltemp = self.cells[n][self.gidlist[n].index(gid)].record['w']
                            
                            if len(walltemp) > 0:
                                for l in range(len(walltemp)):
                                    wtemp = np.array(walltemp[l])
                                    wtemp = interp(times,time,wtemp)
                                    self.comm.Ssend([wtemp, MPI.DOUBLE], dest=0, tag=int(gid))
                                    
                            wtemp = np.ones(len(times))*-1
                            self.comm.Ssend([wtemp, MPI.DOUBLE], dest=0, tag=int(gid))                               

            if self.id == 0: 
                print "root gathered synaptic input conductance" 


            self.barrier()   # wait for other nodes     

    
            for n in range(self.n_celltypes):
                    winh_mat.append([])
                    
                    for i, gid in enumerate(self.global_gidlist[n]): 
                    
                        if self.id == 0:
                            
                            wall = []
                            
                            if gid in self.gidlist[n]:
    
                                walltemp = self.cells[n][self.gidlist[n].index(gid)].record['w_inh']                            
                                if len(walltemp) > 0:
                                    for l in range(len(walltemp)):
                                        wtemp = np.array(walltemp[l])
                                        wtemp = interp(times,time,wtemp)
                                        wall.append(wtemp)
                                        
                            else:
                                
                                while 1:
                                    wtemp = np.zeros(len(times))
                                    self.comm.Recv([wtemp, MPI.DOUBLE], source = MPI.ANY_SOURCE, tag=int(gid))
                                    
                                    if wtemp[0] == -1:
                                        break
                                    else:
                                        wall.append(wtemp)
                                
                            winh_mat[n].append(wall)
                            
                        else:
                            if gid in self.gidlist[n]:
                                walltemp = self.cells[n][self.gidlist[n].index(gid)].record['w_inh']
                                
                                if len(walltemp) > 0:
                                    for l in range(len(walltemp)):
                                        wtemp = np.array(walltemp[l])
                                        wtemp = interp(times,time,wtemp)
                                        self.comm.Ssend([wtemp, MPI.DOUBLE], dest=0, tag=int(gid))
                                        
                                wtemp = np.ones(len(times))*-1
                                self.comm.Ssend([wtemp, MPI.DOUBLE], dest=0, tag=int(gid))
                                    
    
            if self.id == 0: 
                print "root gathered synaptic input conductance"  
  

            self.barrier()   # wait for other nodes     
            

        t_all_vec_vec = []
        id_all_vec_vec = []
        f_cells_mean = []
        
        if self.id == 0:  # only for first node
            
            for n in range(self.n_celltypes):
                
                ie = argsort(t_all_vec[n]) 
                t_all_vec_vec.append( t_all_vec[n][ie] )
                id_all_vec_vec.append( id_all_vec[n][ie].astype(int) ) # 
                
            print "all spikes have been sorted"

            if self.jitter > 0:  # add jitter!
                np.random.seed(40)
                x = np.random.normal(0, self.jitter, len(t_all_vec_vec[self.a_celltype[0]]))             
                t_all_vec_vec[self.a_celltype[0]] = t_all_vec_vec[self.a_celltype[0]] + x
           
            if self.delta_t > 0:
                t_all_vec_vec[self.a_celltype[0]] = t_all_vec_vec[self.a_celltype[0]] + self.delta_t
            
            gsyn = zeros(len(freq_times))
            
            if 'gsyn_in' in self.method_interpol:
                pass
            else:    
                bvec = ["syn" in st for st in self.method_interpol]
                if np.any(bvec):
                    
                    if (not hasattr(self, 'passive_target')) | (self.jitter > 0):  # if not already done in neuron via artificial cell
                        
                        [resp, _] = neuronpy.util.spiketrain.get_histogram(t_all_vec_vec[self.a_celltype[0]], bins = freq_times)
                        resp = np.concatenate((zeros(1),resp))
                        
                        Ksyn = syn_kernel(arange(0,10*self.syn_tau2,self.bin_width), self.syn_tau1, self.syn_tau2) 
                        Ksyn = np.concatenate((zeros(len(Ksyn)-1),Ksyn))
                        gsyn = np.convolve(Ksyn, resp, mode='same')
                        print "Generated gsyn by convolution with Ksyn"
                        self.nc_delay = 0   
                    
                    else:
                        gsyn = interp(freq_times,time,np.array(self.rec_g)) 
                
                spike_freq = np.zeros(len(freq_times))
                
                for j in self.a_celltype:
                    
                    #plt.figure('results_voltage') 
                    #ax99 = plt.subplot(2,1,1)
                    #ax99.plot(time,voltage[j])
                    
                    #plt.text(0.5, 1.1, r'CF=' + str(round(fmean,1)) + ',fmax=' + str(round(fmax,1)) + ',fmstd=' + str(round(fmstd,1)), transform=ax99.transAxes, fontsize=10, va='center', ha='center')
                    #plt.savefig("./figs/Pub/Voltage_" + str(self.pickle_prefix) + "_cell" + str(j) + "_N" + str(self.N[j]) + ".pdf", dpi = 300, transparent=True) # save it  
                    #plt.show()
                    #plt.clf()                    
                    
                    [num_spikes, _] = neuronpy.util.spiketrain.get_histogram(t_all_vec_vec[j], bins = freq_times)
                    
                    if isinstance(self.factor_celltype[j], ( int, long ) ):
                        f = self.factor_celltype[j] 
                    else:
                        f = self.factor_celltype[j][0]    
                    
                    spike_freq = spike_freq + f * np.concatenate((zeros(1),num_spikes)) / self.bin_width

        self.barrier()   # wait for other nodes
          
        #figure('1')
        #plot(time,np.array(self.rec_s1),'b', time,np.array(self.rec_s2),'r')
        #plt.show()
        
        return {'time':time, 'voltage':voltage, 'current':current, 'fmean':fmean, 'f_cells_mean':f_cells_mean,
        'gsyn':gsyn, 'freq_times':freq_times, 'spike_freq':spike_freq, 'gsyn_in':gsyn_in, 'fmeanA':fmeanA, 'fmaxA':fmaxA, 'fmstdA':fmstdA, 'fcvmA':fcvmA, 'fstdmA':fstdmA, 'fbstdA':fbstdA,
        't_all_vec_vec':t_all_vec_vec, 'id_all_vec_vec':id_all_vec_vec, 'gsyns':gsyns, 'w_mat':w_mat, 'winh_mat':winh_mat, 'fmax':fmax, 'fmstd':fmstd, 'fcvm':fcvm, 'fbaseA':fbaseA, 'fbase':fbase}
        
        
    def clean(self):
    
        self.pc.runworker()    
        self.pc.done()   
        
                 
    def compute_Transfer(self, stimulus, spike_freq, freq_times, t, noise_data_points, gsyn, gsyn_in, do_csd, t_qual, K_mat_old, t_startstop, inh_factor=[1]):

        stimulus0 = np.zeros(len(stimulus[0]))
        
        for a in self.a_celltype:
            # sum input to produce linear input that should be reconstructed!
        
            if (any(self.syn_inh_dist) > 0) and (any(self.syn_ex_dist) > 0):
                if  max(self.syn_inh_dist) == max(self.syn_ex_dist): # same signal through ex and inh
                    print "inh_factor = [0,1]"
                    inh_factor = [0,1]   
                    
            for ni in self.syn_ex_dist[a]:
                if ni != 0:
                    stimulus0 += inh_factor[ni-1] * stimulus[ni-1]
                    print "+ex:", ni-1

            for ni in self.syn_inh_dist[a]:
                if ni != 0:
                    stimulus0 -= inh_factor[ni-1] * stimulus[ni-1] #old: +nemax
                    print "-inh:", ni-1 #old: +nemax
                
        if (max(self.n_syn_ex) == 0) and (max(self.n_syn_inh) == 0): 
            stimulus0 += stimulus[0] 
            print "current"
            
        #if self.n_syn_ex[self.celltype_syn[0]] == 0:
        #    stimulus0 += stimulus[0] 
        
        # amplitude should not matter since filter amplitude is simply adjusted 
        #stimulus = stimulus0 #/len(self.syn_ex_dist)

        stimulus0 = stimulus0 / std(stimulus0) / 2
        
        # linear interpolation inside compute_Transfer !!!
        print "max(stimulus0):",max(stimulus0)
        results = compute_Transfer(spike_freq = spike_freq, freq_times = freq_times, 
            stimulus = stimulus0, t = t, noise_data_points = noise_data_points, gsyn = gsyn, gsyn_in = gsyn_in, do_csd = do_csd, t_kernel = 1*s,
            method_interpol = self.method_interpol, nc_delay = self.nc_delay, w_length = 3, t_qual = t_qual, K_mat_old = K_mat_old, t_startstop = t_startstop, give_psd = self.give_psd) # freq_wp not defined, use all frequencies
        
        # TEST:
        #VAF = results.get('VAFf_mat')
        #freq_used = results.get('freq_used')
        
        #iend = mlab.find(freq_used >= self.xmax)[0] 
        #err = 1-mean(VAF[1][0,1:iend-1])
        #print "err: ", err    
        
        return results
        
    
    def residuals_compute_Transfer(self, p, stimulus, spike_freq, freq_times, t, noise_data_points, gsyn, gsyn_in, do_csd, t_qual, K_mat_old, t_startstop, inh_factor):
        
        inh_factor_in = inh_factor[:]
        ip = 0
        for i, inhf in enumerate(inh_factor_in):
            if inhf < 0:
                inh_factor_in[i] = p[ip]
                ip += 1
        
        results = self.compute_Transfer(stimulus = stimulus, spike_freq = spike_freq, freq_times = freq_times, 
            t = t, noise_data_points = noise_data_points, gsyn = gsyn, gsyn_in = gsyn_in, 
            do_csd = do_csd, t_qual = t_qual, K_mat_old = K_mat_old, t_startstop = t_startstop, inh_factor = inh_factor_in)    
        
        VAF = results.get('VAFf_mat')
        freq_used = results.get('freq_used')
        
        iend = mlab.find(freq_used >= self.xmax)[0] 
        err = 1-mean(VAF[1][0,0:iend])
        print "inh_factor:", inh_factor_in, "err: ", err  
               
        return err    
                    
    #@profile            
    def fun_cnoise_Stim(self, t_stim = 10*s, sexp = 0, cutf = 0, do_csd = 1, t_qual = 0, freq_used = np.array([]), K_mat_old = np.array([]), inh_factor = [1], onf = None, equi = 0):
        """
        Stimulate cell with colored noise
        sexp = spectral exponent: Power ~ 1/freq^sexp
        cutf = frequency cutoff: Power flat (white) for freq <~ cutf  
        do_csd = 1: use cross spectral density function for computation
        """
        self.barrier()   # wait for other nodes
        
        filename = str(self.pickle_prefix) + "_results_pop_cnoise.p"
        filepath = self.data_dir + "/" + filename
        
        if self.id == 0: print "- filepath:", filepath 
        
        if self.do_run or (os.path.isfile(filepath) is False):

            tstart = 0;        
            fs = 1 / self.dt # sampling rate 
            fmax = fs / 2 # maximum frequency (nyquist)
            
            t_noise = arange(tstart, t_stim, self.dt) # create stimulus time vector, make sure stimulus is even!!!

            #print self.syn_ex_dist
            #print self.syn_inh_dist
            #exit()
            
            if (self.syn_ex_dist == []):
                for nt in range(self.n_celltypes): # loop over all cells
                    #print "nt", nt
                    if hasattr(self.cells[nt][0], 'input_vec'):
                        self.syn_ex_dist.append([1] * len(self.cells[nt][0].input_vec)) # default ex for all by default!!!
                    else:    
                        self.syn_ex_dist.append([1] * self.n_syn_ex[nt]) # default ex for all by default!!!
                
            #print self.syn_ex_dist
            
            if (self.syn_ex_dist[0] == []):
                nemax = 1
            else:
                nemax = max([item for sublist in self.syn_ex_dist for item in sublist])
            
            if (self.syn_inh_dist == []): # and (any(self.n_syn_inh) > 0)
                for nt in range(self.n_celltypes): # loop over all cells
                        self.syn_inh_dist.append([0] * self.n_syn_inh[nt]) # default no inh for all by default!!!
                    
            #print self.syn_inh_dist
            #exit()
            
            if (self.syn_inh_dist[0] == []):
                nimax = 0
            else:
                nimax = max([item for sublist in self.syn_inh_dist for item in sublist]) 
            
            #print "self.syn_inh_dist, self.syn_ex_dist", self.syn_inh_dist, self.syn_ex_dist
            
            n_noise = max([nemax,nimax]) # number of noise sources
            #print n_noise,nemax,nimax
            # create reproduceable input
            noise_data = []

            for nj in range(n_noise):
            
                if self.id == 0:  # make sure all have the same signal !!!
                    if len(freq_used) == 0: 
                        noise_data0 = create_colnoise(t_noise, sexp, cutf, self.seed+nj, onf = onf)
                    else:
                        noise_data0, _, _, _ = create_multisines(t_noise, freq_used)  # create multi sine signal
                else:
                    noise_data0 = np.empty(len(t_noise), dtype=np.float64)

                noise_data0 = self.broadcast(noise_data0, fast = True)    
                    
                noise_data.append(noise_data0)
                noise_data0 = []            
                   
            noise_data_points = len(noise_data[0]) 

            # Create signal weight vector inh_factor if it is not fully given
            if len(noise_data) > len(inh_factor):
                inh_factor = [inh_factor[0]] * len(noise_data)    
                print "inh_factor:", inh_factor

            #if equi:
                #pass
            #    tstop = t_stim
            
            if max(self.n_syn_ex) == 0: # this means current input
                
                self.set_IStim() # sets amp
                
                if self.fluct_s != []:
                    if self.fluct_s[self.a_celltype[0]] > 0:
                        if self.id == 0: print "- adding i fluct"
                        self.connect_fluct()
                    
                for i, m in enumerate(self.method_interpol):
                            if "syn" in m: self.method_interpol[i] = "syn " + str(self.syn_tau1/ms) + "/" + str(self.syn_tau2/ms) + "ms"
                            if "bin" in m: self.method_interpol[i] = "bin " + str(self.bin_width/ms) + "ms"
                
                stimulus = []
                for nj in range(len(noise_data)):
                    stimulus0, t, t_startstop = construct_Stimulus(noise_data[nj], fs, self.amp[self.a_celltype[0]], ihold = 0, delay_baseline = self.delay_baseline) # , tail_points = 0
                    stimulus.append(stimulus0)
                tstop = t[-1]
                
                self.set_IPlay2(stimulus, t)
                if self.id == 0: print "- starting colored noise transfer function estimation! with amp = " + str(np.round(self.amp[self.a_celltype[0]],4)) + ", ihold = " + str(np.round(self.ihold[self.a_celltype[0]],4)) + ", ihold_sigma = " + str(np.round(self.ihold_sigma,4)) + ", dt = " + str(self.dt) + " => maximum frequency = " + str(fmax) + "\r"    
                 
            else:

                self.give_freq = False
                ihold = self.set_i(self.ihold) # just sets amp, ihold should not change! 

                if 'gsyn_in' not in self.method_interpol:              
                    pass
                else:
                    self.g_syn_ex = [1]*len(self.N)
                    
                    
                if ((self.fluct_g_e0 != []) or (self.fluct_g_i0 != [])):
                    if ((self.fluct_g_e0[self.a_celltype[0]] > 0) or (self.fluct_g_i0[self.a_celltype[0]] > 0)):
                        if self.id == 0: print "- adding g fluct"
                        self.connect_gfluct(E_i=-65)
                
                stimulus = []
                for nj in range(len(noise_data)):
                    stimulus0, t, t_startstop = construct_Stimulus(noise_data[nj], fs, amp=1, ihold = 0, tail_points = 0, delay_baseline = self.delay_baseline) # self.amp
                    stimulus.append(stimulus0)
                
                noise_data = []    
                tstop = t[-1]
                                
                if self.N[self.a_celltype[0]] > 1:
                    self.set_IStim(ihold = [0]*self.n_celltypes, ihold_sigma = [0]*self.n_celltypes, random_start = True, tstart_offset = 1)
                    if self.id == 0: print "- add random start"
                
                #print "Enter Synplay()"
                self.set_SynPlay(stimulus, t, t_startstop = t_startstop)                
                #print "Exit Synplay()"

                if self.id == 0: print "- starting colored noise transfer function estimation with synaptic input! with amp = " + str(np.round(self.amp,4)) + ", ihold = " + str(np.round(self.ihold,4)) + ", ihold_sigma = " + str(np.round(self.ihold_sigma,4)) + ", dt = " + str(self.dt) + " => maximum frequency = " + str(fmax) + "\r"    
 
            amp_vec = []
            mag_vec = []  
            pha_vec = []
            freq_used = []
            ca = []
            SNR_mat = []
            VAFf_mat = []
            Qual_mat = []
            CF_mat = [] 
            VAF_mat = []
            stim = []
            stim_re_mat = []
            resp_mat = []
            current_re = []
            ihold1 = []
            tk = []
            K_mat = []
            gsyn_in = []
            fmean = []
            fmax = [] 
            fmstd = [] 
            fcvm = [] 
            fmeanA = []
            fmaxA = [] 
            fmstdA = [] 
            fcvmA = [] 
            t_all_vec_input_sorted = []
            id_all_vec_input_sorted = []
            
            if (self.id == 0) and (max(self.n_syn_ex) > 0):
                print range(self.n_celltypes), np.shape(self.t_all_vec_input)
                for l in range(self.n_celltypes): 
                    ie = argsort(self.t_all_vec_input[l]) 
                    t_all_vec_input_sorted.append( self.t_all_vec_input[l][ie] )
                    id_all_vec_input_sorted.append( self.id_all_vec_input[l][ie].astype(int) )
            
            #if (self.id == 0):            
            #    print self.g_syn_ex
            #    print np.array(self.g_syn_ex)>= 0
            
            #print "g_syn_ex:",self.g_syn_ex
            if np.array(np.array(self.g_syn_ex)>= 0).any():
                
                if hasattr(self.cells[self.a_celltype[0]][0], 'get_states') and equi:
                    print "- Equilibrate!"
                    self.run(tstop, do_loadstate = False)
                    m = md5.new()
                    cell_exe_new = self.cell_exe[0]
                    m.update(cell_exe_new)
                    filename = './states_' + self.celltype[0] + '_' + m.hexdigest() + '_Population.b'
                    self.cells[self.a_celltype[0]][0].get_states(filename)
                else:
                    self.run(tstop, do_loadstate = False)
                    
                i_startstop = []
                
                results = self.get(t_startstop, i_startstop) 
                time = results.get('time')
                current = results.get('current') 
                voltage = results.get('voltage') 
                fmean = results.get('fmean')  
                gsyn = results.get('gsyn') 
                freq_times = results.get('freq_times')
                spike_freq = results.get('spike_freq')
                t_all_vec_vec = results.get('t_all_vec_vec')
                id_all_vec_vec = results.get('id_all_vec_vec')
                gsyns = results.get('gsyns')
                gsyn_in = results.get('gsyn_in')
                
                fmax = results.get('fmax')
                fmstd = results.get('fmstd')
                fcvm = results.get('fcvm')
                
                fmeanA = results.get('fmeanA')  
                fmaxA = results.get('fmaxA')
                fmstdA = results.get('fmstdA')
                fcvmA = results.get('fcvmA')
                
                fbaseA = results.get('fbaseA') 
                fbase = results.get('fbase')
                fbstdA = results.get('fbstdA')
                
                
            else: # do not run, analyse input!!!
            
                time = t
                voltage = []
                for l in range(self.n_celltypes): 
                    voltage.append(np.zeros(len(t)))
                current = []
                
                freq_times = []
                spike_freq = []
                gsyn = []
                gsyn_in = []
                
                t_all_vec_vec = []
                id_all_vec_vec = []
                
                fmean = []
                fmax = []
                fmstd = []
                fcvm = []
                fstdm = []
                
                fmeanA = []
                fmaxA = []
                fmstdA = []
                fcvmA = []
                fbaseA = []
                fbase = []
                fbstdA = []
                
                if self.id == 0:
                    
                    current = self.n_train_ex
                
                    #t_all_vec = self.t_all_vec_input
                    #id_all_vec = self.id_all_vec_input

                    #ie = argsort(t_all_vec) 
                    #t_all_vec_vec.append( t_all_vec[ie] )
                    #id_all_vec_vec.append( id_all_vec[ie].astype(int) )
                    
                    t_all_vec_vec = t_all_vec_input_sorted
                    id_all_vec_vec = id_all_vec_input_sorted
                    
                    freq_times = arange(0, tstop, self.bin_width)
                    spike_freq = np.zeros(len(freq_times))
                    
                    for j in self.a_celltype:
                        
                        [num_spikes, _] = neuronpy.util.spiketrain.get_histogram(t_all_vec_vec[j], bins = freq_times)

                        if self.tau2_ex[0] > 0:
                            spike_freq = np.concatenate((zeros(1),num_spikes)) 
                            print "NOSYN TEST: start convolution with Ksyn"
                            Ksyn = syn_kernel(arange(0,10*self.tau2_ex[0],self.bin_width), self.tau1_ex[0], self.tau2_ex[0]) 
                            Ksyn = np.concatenate((zeros(len(Ksyn)-1),Ksyn))
                            spike_freq = np.convolve(Ksyn, spike_freq, mode='same')
                            print "NOSYN TEST: convolution finished"
                        else:

                            if isinstance(self.factor_celltype[j], ( int, long ) ):
                                f = self.factor_celltype[j] 
                            else:
                                f = self.factor_celltype[j][0]    
                            
                            spike_freq = spike_freq + f * np.concatenate((zeros(1),num_spikes)) / self.bin_width

                    fmean.append(self.fmean_input)
                    fmax.append(self.fmax_input) 
                    fmstd.append(self.fmstd_input) 
                    fcvm.append(self.fcvm_input) 
                    fstdm.append(self.fstdm_input)

                    if self.no_fmean == True:
                        fmean.append(ihold)
                        
                    #plt.figure('spike_freq') 
                    #plt.plot(freq_times, spike_freq)
                    #plt.savefig("./figs/Pub/Spike_freq_" + str(self.pickle_prefix) + ".pdf", dpi = 300, transparent=True) # save it  
                    #plt.clf()
                        
                    fmeanA = fmean[0]
                    fmaxA = fmax[0]
                    fmstdA = fmstd [0] 
                    fcvmA = fcvm[0]
                    fstdmA = fstdm[0]
                        
                        
            if self.id == 0: 
                
                if any([i<0 for i in inh_factor]):
                    
                    p0 = []
                    inhf_idx = []
                    for i, inhf in enumerate(inh_factor):
                        if inhf < 0: 
                            p0.append(0)  
                            inhf_idx.append(i)
                            
                    plsq = fmin(self.residuals_compute_Transfer, p0, args=(stimulus, spike_freq, freq_times, t, noise_data_points, gsyn, gsyn_in, do_csd, t_qual, K_mat_old, t_startstop, inh_factor))
                    p = plsq
                    
                    ip = 0
                    for i in inhf_idx:
                        inh_factor[i] = p[ip]
                        ip += 1
                            

                    print "Final inh_factor: ", inh_factor
                
                
                results = self.compute_Transfer(stimulus, spike_freq = spike_freq, freq_times = freq_times, 
                    t = t, noise_data_points = noise_data_points, gsyn = gsyn, gsyn_in = gsyn_in, 
                    do_csd = do_csd, t_qual = t_qual, K_mat_old = K_mat_old, t_startstop = t_startstop, inh_factor=inh_factor)
                        
                mag_vec, pha_vec, ca, freq, freq_used, fmean_all = results.get('mag_mat'), results.get('pha_mat'), results.get('ca_mat'), results.get('freq'), results.get('freq_used'), results.get('fmean')          
                SNR_mat, VAFf_mat, Qual_mat, CF_mat, VAF_mat = results.get('SNR_mat'), results.get('VAFf_mat'), results.get('Qual_mat'), results.get('CF_mat'), results.get('VAF_mat')  
                stim, resp_mat, stim_re_mat, tk, K_mat = results.get('stim'), results.get('resp_mat'), results.get('stim_re_mat'), results.get('tk'), results.get('K_mat')                     
        
        
            self.barrier()   # wait for other nodes
        
        
            if self.id == 0:
                
                if t_qual > 0:
                    #print t_startstop[0], t_startstop[0]/self.dt, (t_startstop[0]+t_qual)/self.dt
                    current_re = current[int(t_startstop[0]/self.dt):int((t_startstop[0]+t_qual)/self.dt)]
                    current_re = current_re[int(len(K_mat[self.a_celltype[0]])):int(len(current_re))-int(len(K_mat[self.a_celltype[0]]))]
                
                if len(self.i_holdrs) > 0:
                    ihold1 = self.i_holdrs[self.a_celltype[0]][0]
                else:
                    ihold1 = []
             
                for l in range(len(self.method_interpol)):  # unwrap 
                    pha_vec[l,:] = unwrap(pha_vec[l,:] * (pi / 180)) * (180 / pi)  # unwrap for smooth phase
                
                # only return fraction of actual signal, it is too long!!!     
                if time[-1] > self.tmax:  
                    imax = -1*int(self.tmax/self.dt)
                    time = time[imax:]; current = current[imax:]; gsyn = gsyn[imax:]; gsyn_in = gsyn_in[imax:]
                    for n in range(self.n_celltypes): 
                        voltage[n] = voltage[n][imax:]
                    
                if freq_times != []:   
                    if freq_times[-1] > self.tmax:
                        imax2 = where(freq_times > self.tmax)[0][0]  # for spike frequency   
                        freq_times = freq_times[0:imax2]; spike_freq = spike_freq[0:imax2] 
            
                bvec = ["_syn" in st for st in self.method_interpol]
                if np.any(bvec):
                    # normalize synaptic integration with others    
                    mag_vec[1,:]= mag_vec[0,0]*mag_vec[1,:]/mag_vec[1,0] 
                    
            if self.id == 0: print "start pickle"
             
            results = {'freq_used':freq_used, 'amp':amp_vec,'mag':mag_vec,'pha':pha_vec,'ca':ca,'voltage':voltage,'tk':tk,'K_mat':K_mat, 'ihold1': ihold1, 't_startstop':t_startstop, #'stimulus':stimulus,
                        'current':current,'t1':time,'freq_times':freq_times,'spike_freq':spike_freq, 'stim':stim, 'stim_re_mat':stim_re_mat, 'resp_mat':resp_mat, 'current_re':current_re, 'gsyn_in':gsyn_in, 'fmeanA':fmeanA, 'fmaxA':fmaxA, 'fmstdA':fmstdA, 'fcvmA':fcvmA, 'fbaseA':fbaseA, 'fbase':fbase, 'fbstdA':fbstdA,
                        'fmean':fmean,'method_interpol':self.method_interpol, 'SNR':SNR_mat, 'VAF':VAFf_mat, 'Qual':Qual_mat, 'CF':CF_mat, 'VAFs':VAF_mat, 'fmax':fmax, 'fmstd':fmstd, 'fcvm':fcvm, 'inh_factor':inh_factor, 't_all_vec_vec':t_all_vec_vec, 'id_all_vec_vec':id_all_vec_vec}       
        
            if self.id == 0:
                if self.dumpsave == 1:
                    pickle.dump( results, gzip.GzipFile( filepath, "wb" ) )
                    print "pickle done" 
            
                
                if self.plot_train:
                    
                    for a in self.a_celltype:

                        #i_start = mlab.find(t_all_vec_vec[a] >= 0)[0]
                        #i_stop = mlab.find(t_all_vec_vec[a] >= 5)[0]
                        
                        #t_all_cut = t_all_vec_vec[a][i_start:i_stop]
                        #id_all_cut = id_all_vec_vec[a][i_start:i_stop]
                        
                        t_all_cut = t_all_vec_vec[a]
                        id_all_cut = id_all_vec_vec[a]
                        
                        f_start_in = mlab.find(t_all_cut >= 0) 
                        f_stop_in = mlab.find(t_all_cut <= 10) 
                        
                        f_start = f_start_in[0]  
                        f_stop = f_stop_in[-1]+1      
                        use_spikes = t_all_cut[f_start:f_stop]
                        use_id = id_all_cut[f_start:f_stop]
                        
                        plt.figure('results_train') 
                        ax99 = plt.subplot(1,1,1)
                        ax99.plot(use_spikes,use_id,'|', ms=2)
                        plt.text(0.5, 1.1, r'CF=' + str(round(fmean,1)) + ',fmax=' + str(round(fmax,1)) + ',fmstd=' + str(round(fmstd,1)), transform=ax99.transAxes, fontsize=10, va='center', ha='center')
                        plt.savefig("./figs/Pub/Train_" + str(self.pickle_prefix) + "_cell" + str(a) + "_N" + str(self.N[a]) + ".pdf", dpi = 300, transparent=True) # save it  
                                     
                        plt.clf()
                        
                        if len(t_all_cut) > 0:
                        
                            tbin = 100*ms
                            tb = np.arange(0,t[-1],tbin)
                            [all_rate, _] = neuronpy.util.spiketrain.get_histogram(t_all_cut, bins = tb)
                            all_rate = np.concatenate((np.zeros(1),all_rate)) / self.N[a] / tbin
                
                            plt.figure('results_train2') 
                            plt.plot(tb,all_rate)
                            plt.savefig("./figs/Pub/PSTH_" + str(self.pickle_prefix) + "_cell" + str(a) + "_N" + str(self.N[a]) + ".pdf", dpi = 300, transparent=True) # save it  
                            plt.clf()
                        
                        plt.figure('results_noise') 
                        plt.plot(time,current)
                        plt.savefig("./figs/Pub/Noise_" + str(self.pickle_prefix) + "_cell" + str(a) + "_N" + str(self.N[a]) + ".pdf", dpi = 300, transparent=True) # save it  
                        plt.clf()
                    
                
                if self.plot_input:
                    
                    if len(t_all_vec_input_sorted[0]) > 0:
                        
                        i_start = mlab.find(t_all_vec_input_sorted[0] >= 0)[0]
                        i_stop = mlab.find(t_all_vec_input_sorted[0] >= 5)[0]
                    
                        t_all_cut = t_all_vec_input_sorted[0][i_start:i_stop]
                        id_all_cut = id_all_vec_input_sorted[0][i_start:i_stop]
                         
                        plt.figure('results_input') 
                        ax99 = plt.subplot(1,1,1)
                        ax99.plot(t_all_cut,id_all_cut,'|', ms=2)
                        plt.text(0.5, 1.1, r'fmean=' + str(round(self.fmean_input,1)) + ',fmax=' + str(round(self.fmax_input,1)) + ',fmstd=' + str(round(self.fmstd_input,1)) + ',fcvm=' + str(round(self.fcvm_input,1)) + ',fstdm=' + str(round(self.fstdm_input,1)), transform=ax99.transAxes, fontsize=10, va='center', ha='center')
                        plt.savefig("./figs/Pub/Input_" + str(self.pickle_prefix) + "_N" + str(self.N[self.a_celltype[0]]) + ".pdf", dpi = 300, transparent=True) # save it  
                        plt.clf()
                    

        else:
        
            if self.id == 0:
                results = pickle.load( gzip.GzipFile( filepath, "rb" ) )
                
                #print results
                #print {key:np.shape(value) for key,value in results.iteritems()}
                
                if self.minimal_dir: # save only info needed for plot
                
                    print {key:np.shape(value) for key,value in results.iteritems()}
                    
                    if "Fig6_pop_transfer_grc_syngr_nsyn4_cn_a1_noisesynlow_inhlow_adjfinh_varih_N100_CFo6.0_results_pop_cnoise.p" in filename:
                        results['ca'] = [] 
                        results['resp_mat'] = []
                        results['stim'] = []
                        results['current'] = []
                        results['tk'] = []
                        results['K_mat'] = []
                        results['freq_times'] = []
                        results['spike_freq'] = []
                        results['stim_re_mat'] = []
                        results['current_re'] = []
                        results['t_all_vec_vec'] = []
                        results['id_all_vec_vec'] = []  
                        results['gsyn_in'] = []
                        
                    elif ("Fig8_pop_transfer_none_synno_cn_cutf30_a1_noisesynlow_ih20_varih_N100_CFo-1_results_pop_cnoise.p" in filename) \
                         or ("Fig8_pop_transfer_none_synno_cn_cutf30_a10_noisesynlow_ih20_varih_N100_CFo-1_results_pop_cnoise.p" in filename) \
                         or ("Fig8_pop_transfer_grc_syngr_nsyn4_cn_cutf30_a1_noisesynlow_inhlow_adjfinh_varih_varinhn_N100_CFo9.0_results_pop_cnoise.p" in filename) \
                         or ("Fig8_pop_transfer_grc_syngr_nsyn4_cn_cutf30_a10_noisesynlow_inhlow_adjfinh_varih_varinhn_N100_is0.14_CFo9.0_results_pop_cnoise.p" in filename) \
                         :

                        results['ca'] = [] 
                        results['resp_mat'] = []
                        results['current'] = []
                        results['tk'] = []
                        results['K_mat'] = []
                        results['voltage'] = [] 
                        results['current_re'] = []
                        results['t_all_vec_vec'] = []
                        results['id_all_vec_vec'] = []
                        results['t1'] = []
                        results['gsyn_in'] = []
                        
                    elif ("Fig8_pop_transfer_none_synno_cn_cutf30_a1_noisesynlow_ih20_varih_N50_twopop_CFo-1_results_pop_cnoise.p" in filename) \
                         or ("Fig8_pop_transfer_none_synno_cn_cutf30_a10_noisesynlow_ih20_varih_N50_twopop_CFo-1_results_pop_cnoise.p" in filename) \
                         or ("Fig8_pop_transfer_grc_syngr_nsyn4_cn_cutf30_a1_noisesynlow_inhlow_adjfinh_varih_varinhn_N50_twopop_CFo9.0_results_pop_cnoise.p" in filename) \
                         or ("Fig8_pop_transfer_grc_syngr_nsyn4_cn_cutf30_a10_noisesynlow_inhlow_adjfinh_varih_varinhn_N50_is0.14_twopop_CFo9.0_results_pop_cnoise.p" in filename) \
                         or ("Fig8_pop_transfer_grc_syngr_nsyn4_cn_cutf5_a1_noisesynlow_inhlow_adjfinh_varih_varinhn_N100_CFo14.0_results_pop_cnoise.p" in filename) \
                         or ("Fig8_pop_transfer_grc_syngr_nsyn4_cn_cutf5_a1_noisesynlow_inhlow_adjfinh_varih_varinhn_N50_twopop_CFo14.0_results_pop_cnoise.p" in filename) \
                         :
 
                        results['ca'] = [] 
                        results['resp_mat'] = []
                        results['current'] = []
                        results['tk'] = []
                        results['K_mat'] = []
                        results['voltage'] = [] 
                        results['current_re'] = []
                        results['t_all_vec_vec'] = []
                        results['id_all_vec_vec'] = []
                        results['t1'] = []
                        results['gsyn_in'] = []
                        results['freq_times'] = []
                        results['spike_freq'] = []
                        
                    elif ("Fig4_pop_transfer_grc_cn_addn100_N[100]_CF[40]_amod[1]_results_pop_cnoise.p" in filename) \
                         or ("Fig4_pop_transfer_grc_cn_addn1_N[100]_CF[40]_amod[1]_results_pop_cnoise.p" in filename) \
                         or ("Fig4b_pop_transfer_grc_lowcf_cn_twopop_N[50, 50]_CF[0.0055, 0.0055]_amod[None, None]_results_pop_cnoise.p" in filename) \
                         or ("Fig4b_pop_transfer_grc_lowcf_cn_N[100]_CF[0.0055]_amod[None]_results_pop_cnoise.p" in filename) \
                         or ("Fig4b_pop_transfer_grc_lowcf_slownoise_cn_twopop_N[50, 50]_CF[0.0051, 0.0051]_amod[None, None]_results_pop_cnoise.p" in filename) \
                         or ("Fig4b_pop_transfer_grc_lowcf_slownoise_cn_N[100]_CF[0.0051]_amod[None]_results_pop_cnoise.p" in filename) \
                         :
 
                        results['ca'] = [] 
                        results['resp_mat'] = []
                        results['current'] = []
                        results['tk'] = []
                        results['K_mat'] = []
                        results['voltage'] = [] 
                        results['t_all_vec_vec'] = []
                        results['id_all_vec_vec'] = []
                        results['t1'] = []
                        results['gsyn_in'] = []
                        results['freq_times'] = []
                        results['spike_freq'] = []
                        
                    elif ("Fig2_pop_transfer_" in filename) \
                         :
                             
                        results['ca'] = [] 
                        results['resp_mat'] = []
                        results['current'] = []
                        results['t1'] = []
                        results['voltage'] = [] 
                        results['freq_times'] = []
                        results['spike_freq'] = []
                        results['current_re'] = []
                        results['t_all_vec_vec'] = []
                        results['id_all_vec_vec'] = []
                        results['gsyn_in'] = []
 
                    else:
                        results['ca'] = [] 
                        results['resp_mat'] = []
                        results['stim'] = []
                        results['current'] = []
                        results['tk'] = []
                        results['K_mat'] = []
                        results['t1'] = []
                        results['voltage'] = [] 
                        results['freq_times'] = []
                        results['spike_freq'] = []
                        results['stim_re_mat'] = []
                        results['current_re'] = []
                        results['t_all_vec_vec'] = []
                        results['id_all_vec_vec'] = []
                        results['gsyn_in'] = []

                    print {key:np.shape(value) for key,value in results.iteritems()}

                    pickle.dump( results, gzip.GzipFile( self.minimal_dir + "/" + filename, "wb" ) )                
                
            else:
                results = {'freq_used':[], 'amp':[],'mag':[],'pha':[],'ca':[],'voltage':[], 'tk':[],'K_mat':[], 'ihold1':[], 't_startstop':[], #'stimulus':[],
                    'current':[],'t1':[],'freq_times':[],'spike_freq':[], 'stim':[], 'stim_re_mat':[], 'current_re':[], 'gsyn_in':[], 'fmeanA':[], 'fmaxA':[], 'fmstdA':[], 'fcvmA':[], 'fbaseA':[], 'fbase':[], 'fbstdA':[],
                     'fmean':[],'method_interpol':self.method_interpol, 'SNR':[], 'VAF':[], 'Qual':[], 'CF':[], 'VAFs':[], 'fmax':[], 'fmstd':[], 'fcvm':[], 'inh_factor':[], 't_all_vec_vec':[], 'id_all_vec_vec':[]}  
        
        if self.id == 0:             

            if self.plot_train: 

                for a in self.a_celltype:
                    
                    t1 = results.get('t1') 
                    voltage = results.get('voltage') 
                    fmean = results.get('fmean')  
                    fmax = results.get('fmax')  
                    fmstd = results.get('fmstd') 
                    
                    
                    if results.has_key('t_all_vec_vec'):
                                
                        if len(results['t_all_vec_vec']) > 0:        
                            t_all_vec_vec = results.get('t_all_vec_vec') 
                            id_all_vec_vec = results.get('id_all_vec_vec') 
                            
                            t_all_cut = t_all_vec_vec[a]
                            id_all_cut = id_all_vec_vec[a]
                            
                            f_start_in = mlab.find(t_all_cut >= 0) 
                            f_stop_in = mlab.find(t_all_cut <= 10) 
                            
                            f_start = f_start_in[0]  
                            f_stop = f_stop_in[-1]+1      
                            use_spikes = t_all_cut[f_start:f_stop]
                            use_id = id_all_cut[f_start:f_stop]
                            
                            plt.figure('results_train') 
                            ax97 = plt.subplot(1,1,1)
                            ax97.plot(use_spikes,use_id,'|', ms=6)
                            plt.text(0.5, 1.1, r'CF=' + str(round(fmean,1)) + ',fmax=' + str(round(fmax,1)) + ',fmstd=' + str(round(fmstd,1)), transform=ax97.transAxes, fontsize=10, va='center', ha='center')
                            plt.savefig("./figs/Pub/Train_" + str(self.pickle_prefix) + "_cell" + str(a) + "_N" + str(self.N[a]) + ".pdf", dpi = 300, transparent=True) # save it  

                    
                    plt.figure('results_voltage') 
                    ax99 = plt.subplot(2,1,1)
                    ax99.plot(t1,voltage[a])
                    
                    t_noise = arange(0, t_stim, self.dt)
                    noise_data = create_colnoise(t_noise, sexp, cutf, 50, onf = onf)
                    stimulus, t, t_startstop = construct_Stimulus(noise_data, 1/self.dt, amp=1, ihold = 0, tail_points = 0, delay_baseline = self.delay_baseline) 
                    ax98 = plt.subplot(2,1,2)
                    ax98.plot(t[0:10/self.dt],stimulus[0:10/self.dt],color='k')
                    
                    plt.text(0.5, 1.1, r'CF=' + str(round(fmean,1)) + ',fmax=' + str(round(fmax,1)) + ',fmstd=' + str(round(fmstd,1)), transform=ax99.transAxes, fontsize=10, va='center', ha='center')
                    plt.savefig("./figs/Pub/Voltage_" + str(self.pickle_prefix) + "_cell" + str(a) + "_N" + str(self.N[a]) + ".pdf", dpi = 300, transparent=True) # save it  
                    plt.show()
                    plt.clf()
                        
        if (self.id == 0) and (do_csd == 1):
            Qual = results.get('Qual')            
            for i, ii in enumerate(self.method_interpol):
                        print "\n[QUAL:] Interpol:", ii, "SNR0:", Qual[i,0,0], "SNR_cutff:", Qual[i,0,1], "SNR_mean:", Qual[i,0,2], "\n VAF0:", Qual[i,1,0], "VAF_cutff:", Qual[i,1,1], "VAF_mean:", Qual[i,1,2],  "\n CF(subtracted):", Qual[i,2,0], "VAF(subtracted):", Qual[i,2,1]     
            
            VAF = results.get('VAF')
            freq_used = results.get('freq_used') 
            iend = mlab.find(freq_used >= self.xmax)[0] 
            print 'm(VAF)=' + str(np.mean(VAF[1][0,0:iend])) 
                
        self.barrier()   # wait for other nodes
            
        return results


#    def fun_ssine_Stim(self, freq_used = np.array([1, 10, 100, 1000])*Hz):
#        """
#        Compute impedance and/or transfer function using Single sine stimulation
#        Only compute transfer function if there is a steady state (resting) firing rate!
#        """
#        self.barrier()   # wait for other nodes
#        
#        filepath = "./data/" + str(self.pickle_prefix) + "_results_pop_ssine.p"
#        
#        if self.do_run or (os.path.isfile(filepath) is False):
#            
#            fs = 1 / self.dt # sampling rate 
#            fmax = fs / 2 # maximum frequency (nyquist)
#            
#            if self.id == 0: print "- starting single sine transfer function estimation! with amp = " + str(np.round(self.amp[a_celltype[0]],4)) + ", ihold = " + str(np.round(self.ihold[self.a_celltype[0]],4)) + ", dt = " + str(self.dt) + " => maximum frequency = " + str(fmax) + "\r"    
#    
#            if max(self.n_syn_ex) == 0:
#                self.set_IStim() 
#                
#                if self.fluct_s != []:
#                    if self.fluct_s[self.a_celltype[0]] > 0:
#                        if self.id == 0: print "- adding i fluct"
#                        self.connect_fluct()
#                    
#                for i, m in enumerate(self.method_interpol):
#                    if "syn" in m: self.method_interpol[i] = "syn " + str(self.syn_tau1/ms) + "/" + str(self.syn_tau2/ms) + "ms"
#                    if "bin" in m: self.method_interpol[i] = "bin " + str(self.bin_width/ms) + "ms"
#            
#            else:
#                self.give_freq = False
#                ihold = self.set_i(self.ihold) # just sets amp, ihold should not change! 
#                
#                if ((self.fluct_g_e0 != []) or (self.fluct_g_i0 != [])):
#                    if ((self.fluct_g_e0[self.a_celltype[0]] > 0) or (self.fluct_g_i0[self.a_celltype[0]] > 0)):
#                        if self.id == 0: print "- adding g fluct"
#                        self.connect_gfluct(E_i=-65)
#                        
#                #if ((self.fluct_std_e[self.a_celltype[0]] != []) or (self.fluct_std_i[self.a_celltype[0]] != [])):
#                #    if ((self.fluct_std_e[self.a_celltype[0]] > 0) or (self.fluct_std_i[self.a_celltype[0]] > 0)):
#                #        if self.id == 0: print "- adding g fluct"
#                #        self.connect_gfluct(E_i=-65)
#                
#                if 'gsyn_in' not in self.method_interpol:              
#                    pass
#                else:
#                    self.g_syn_ex = 1
#                
#            
#            for i, fu in enumerate(freq_used):
#                
#                if self.id == 0: print "- single sine processing frequency = " + str(fu)
#                
#                t, stimulus, i_startstop, t_startstop = create_singlesine(fu = fu, amp = self.amp[a_celltype[0]], ihold = 0, dt = self.dt, periods = 20, minlength = 2*s, t_prestim = 1*s)
#                tstop = t[-1]
#                
#                if i == 0: t_startstop_plot = t_startstop
#                    
#                if max(self.n_syn_ex) == 0:
#                    self.set_IPlay(stimulus, t)
#                else:
#                    self.set_SynPlay(stimulus, t)    
#                
#                if self.g_syn_ex >= 0: # should also be true for current input!!!
#              
#                    self.run(tstop)
#                                    
#                    if i == 0:   # do this here to have something to return
#                        
#                        # select first sinusoidal to plot, later
#                        voltage_plot = []
#                        current_plot = []
#                        time_plot = []
#                        freq_times_plot = []
#                        spike_freq_plot = []
#                        gsyn_plot = []
#                        
#                        # construct vectors
#                        amp_vec = zeros(len(freq_used)) # amplitude vector
#                        fmean_all = zeros(len(freq_used)) # mean firing frequency (all cells combined)
#                        fmean = zeros(len(freq_used)) # mean firing frequency (one cell)
#                        ca = zeros(len(freq_used), dtype=complex)
#                        
#                        # create matrix to hold all different interpolation methods:
#                        mag_vec = zeros((len(self.method_interpol),len(freq_used)))  # magnitude vector
#                        pha_vec = zeros((len(self.method_interpol),len(freq_used))) # phase vector 
#                        NI_vec = zeros((len(self.method_interpol),len(freq_used)))  # NI vector
#                        VAF_vec = zeros((len(self.method_interpol),len(freq_used)))  # VAF vector
#                    
#                    results = self.get(t_startstop, i_startstop) # t1 should be equal to t!!!
#                    time, voltage, current, fmean0, gsyn = results.get('time'), results.get('voltage'), results.get('current'), results.get('fmean'), results.get('gsyn')
#                    freq_times, spike_freq, t_all_vec_vec, id_all_vec_vec, gsyns = results.get('freq_times'), results.get('spike_freq'), results.get('t_all_vec_vec'), results.get('id_all_vec_vec'), results.get('gsyns')
# 
#                else:
#                    
#                    time = t
#                    voltage = []
#                    voltage.append(np.zeros(len(t)))
#                    current = stimulus
#                    
#                    freq_times = []
#                    spike_freq = []
#                    fmean0 = ihold
#                    gsyn = []
#                    gsyn_in = []
#                    
#                    t_all_vec_vec = []
#                    id_all_vec_vec = []
#                    
#                    
#                    if self.id == 0:
#                    
#                        t_all_vec = []
#                        t_all_vec.append([])
#                        t_all_vec[0] = np.concatenate(self.t_all_vec_input)
#                        
#                        id_all_vec = []
#                        id_all_vec.append([])
#                        id_all_vec[0] = np.concatenate(self.id_all_vec_input)
#                        
#                        ie = argsort(t_all_vec[0]) 
#                        t_all_vec_vec.append( t_all_vec[0][ie] )
#                        id_all_vec_vec.append( id_all_vec[0][ie].astype(int) ) # 
#                        
#                                            
#                        freq_times = arange(0, tstop, self.bin_width)
#                        [num_spikes, _] = neuronpy.util.spiketrain.get_histogram(t_all_vec_vec[0], bins = freq_times)
#                        spike_freq = np.concatenate((zeros(1),num_spikes)) / self.bin_width
#
#    
#                if self.id == 0:
#
#                    fmean[i] = fmean0[0]
#
#                    if i == 0:   
#                    
#                        # select first sinusoidal to plot
#                        voltage_plot = voltage
#                        current_plot = current
#                        time_plot = time
#                        freq_times_plot = freq_times
#                        spike_freq_plot = spike_freq
#                        gsyn_plot = gsyn
#                    
#                    
#                    for l in range(len(self.method_interpol)):
#                        
#                        if "bin" in self.method_interpol[l]:
#    
#                            # binning and linear interpolation
#                            stimulus_signal = stimulus[i_startstop[0]:i_startstop[1]] # cut out relevant signal
#                            t_input_signal = t[i_startstop[0]:i_startstop[1]] - t[i_startstop[0]]
#                            
#                            spike_freq_interp = interp(t, freq_times, spike_freq, left=0, right=0) # interpolate to be eqivalent with input, set zero at beginning and end!
#                            freq_out_signal_interp = spike_freq_interp[i_startstop[0]:i_startstop[1]] # cut out relevant signal
#                            vamp, mag_vec[l,i], pha_vec[l,i], fmean_all[i], _ = get_magphase(stimulus_signal, t_input_signal, freq_out_signal_interp, t_input_signal, method = "fft", f = fu)
#                            
#                            results = est_quality(t_input_signal, fu, freq_out_signal_interp, self.amp[a_celltype[0]]*mag_vec[l,i], pha_vec[l,i]/ (180 / pi), fmean_all[i])  
#                            NI_vec[l,i], VAF_vec[l,i] = results.get('NI'), results.get('VAF')
#                            print "-[bin] NI: " + str(NI_vec[l,i]) + ", VAF: " + str(VAF_vec[l,i])
#                            
#                        if "syn" in self.method_interpol[l]:
#                            
#                            # synaptic integration  
#                            dt_out = t_input_signal[2] - t_input_signal[1]
#                            shift = self.nc_delay/dt_out  # shift response by the nc delay to remove offset
#                            freq_out_signal_syn = gsyn[i_startstop[0]+shift:i_startstop[1]+shift] # cut out relevant signal
#                            
#                            vamp, mag_vec[l,i], pha_vec[l,i], fm, _ = get_magphase(stimulus_signal, t_input_signal, freq_out_signal_syn, t_input_signal, method = "fft", f = fu)
#                            
#                            results = est_quality(t_input_signal, fu, freq_out_signal_syn, self.amp[a_celltype[0]]*mag_vec[l,i], pha_vec[l,i]/ (180 / pi), fm)  
#                            NI_vec[l,i], VAF_vec[l,i] = results.get('NI'), results.get('VAF')
#                            print "-[syn] NI: " + str(NI_vec[l,i]) + ", VAF: " + str(VAF_vec[l,i])
#                            
#    
#                self.barrier()   # wait for other nodes
#                
#            #print "rest: " + str(vrest) + " freq_used:" + str(freq_used) + " amp_vec:" + str(amp_vec) + " mag_vec:" + str(mag_vec) + " pha_vec:" + str(pha_vec)
#            
#            if self.id == 0:
#                
#                for l in range(len(self.method_interpol)):  # unwrap 
#                    pha_vec[l,:] = unwrap(pha_vec[l,:] * (pi / 180)) * (180 / pi)  # unwrap for smooth phase
#                
#                # only return fraction of actual signal, it is too long!!!     
#                if time_plot[-1] > self.tmax:  
#                    imax = where(time_plot > self.tmax)[0][0]  # for voltage, current and time
#                    time_plot = time_plot[0:imax];  current_plot = current_plot[0:imax]; gsyn_plot = gsyn_plot[0:imax]
#                    for n in range(self.n_celltypes): 
#                        voltage_plot[n] = voltage_plot[n][0:imax]
#                        
#                if freq_times_plot != []:   
#                    if freq_times_plot[-1] > self.tmax:
#                        imax2 = where(freq_times_plot > self.tmax)[0][0]  # for spike frequency   
#                        freq_times_plot = freq_times_plot[0:imax2]; spike_freq_plot = spike_freq_plot[0:imax2] 
#            
#                # normalize synaptic integration with with first magnitude, may by syn itself! 
#                bvec = ["syn" in st for st in self.method_interpol]
#                if np.any(bvec):
#                    k = where(bvec)   
#                    mag_vec[k,:]= mag_vec[0,0]*mag_vec[k,:]/mag_vec[k,0]
#            
#            NI_vec = (freq_used, NI_vec)
#            VAF_vec = (freq_used, VAF_vec)
#            results = {'freq_used':freq_used, 'amp':amp_vec,'mag':mag_vec,'pha':pha_vec,'ca':ca,'voltage':voltage_plot, 't_startstop':t_startstop_plot,
#                'current':current_plot,'t1':time_plot,'freq_times':freq_times_plot,'spike_freq':spike_freq_plot,
#                'fmean':mean(fmean),'method_interpol':self.method_interpol, 'NI':NI_vec, 'VAF':VAF_vec}
#            
#            if self.id == 0:
#                pickle.dump( results, gzip.GzipFile( filepath, "wb" ) )
#            
#        else:
#        
#            if self.id == 0:
#                results = pickle.load( gzip.GzipFile( filepath, "rb" ) )
#            else:
#                results = {'freq_used':[], 'amp':[],'mag':[],'pha':[],'ca':[],'voltage':[], 't_startstop':[],
#                    'current':[],'t1':[],'freq_times':[],'spike_freq':[],
#                     'fmean':[],'method_interpol':self.method_interpol,'NI':[],'VAF':[]}                
#            
#        return results
        
    def get_RC(self, opt_plot):
        
        if self.id == 0:
            if "analytical" in opt_plot: # simplest case, only uses rm and tau, scaling necessary 
                exec self.cell_exe[self.a_celltype[0]]
                sim = Stimulation(cell, temperature = self.temperature)
                rm, cm, taum = sim.get_RCtau()
            else:
                rm = cm = taum = 0
            
            if "if" in opt_plot:
                Vrest = cell.soma(0.5).pas.e*mV
                Vth = cell.spkout.thresh*mV 
                Vreset = cell.spkout.vrefrac*mV
            else:
                Vreset = 0*mV; Vth = 1*mV; Vrest = 0*mV
                
            sim = None
            cell = None  
        else:
            rm = cm = taum = 0
            Vreset = 0*mV; Vth = 1*mV; Vrest = 0*mV
            
        return rm, cm, taum, Vreset, Vth, Vrest


    def fun_plot(self, currlabel="control", dowhat="cnoise", freq_used=np.array([]), cutf=10, sexp=0, t_stim=100*s, ymax=0, ax=None, SNR=None, VAF=None, t_qual=0, opt_plot=np.array([]), method_interpol_plot=[], do_csd = 1):

        SNR_switch = SNR
        VAF_switch = VAF
        
        rm, cm, taum, Vreset, Vth, Vrest = self.get_RC(opt_plot)
        
        if dowhat == "cnoise":
            
            if do_csd == 0:
                t_qual = 0; SNR_switch = 0; VAF_switch = 0

            results = self.fun_cnoise_Stim(t_stim = t_stim, cutf = cutf, sexp = sexp, t_qual = t_qual, freq_used = freq_used, do_csd = do_csd)
            
            freq_used, amp_vec, mag, pha, ca, voltage, current, t1 = results.get('freq_used'), results.get('amp'), results.get('mag'), results.get('pha'), results.get('ca'), results.get('voltage'), results.get('current'), results.get('t1') 
            freq_times, spike_freq, fmean, method_interpol, SNR, VAF, Qual = results.get('freq_times'), results.get('spike_freq'), results.get('fmean'), results.get('method_interpol'), results.get('SNR'), results.get('VAF'), results.get('Qual')                              
            stim, stim_re_mat, current_re, tk, K_mat_old = results.get('stim'), results.get('stim_re_mat'), results.get('current_re'), results.get('tk'), results.get('K_mat')
            
        elif dowhat == "ssine":
            
            results = self.fun_ssine_Stim(freq_used = freq_used0)
            
            freq_used, amp_vec, mag, pha, ca, voltage, current, t1 = results.get('freq_used'), results.get('amp'), results.get('mag'), results.get('pha'), results.get('ca'), results.get('voltage'), results.get('current'), results.get('t1') 
            freq_times, spike_freq, fmean, method_interpol, VAF = results.get('freq_times'), results.get('spike_freq'), results.get('fmean'), results.get('method_interpol'), results.get('VAF')                     
            tk = []
            K_mat_old = []

        # analyse
        if self.id == 0:
            
            print "Mean rate: " + str(fmean)
            
            # Turn it off if set to zero
            if SNR_switch == 0: SNR = None
            if VAF_switch == 0: VAF = None 

            
            if t_qual > 0:
                
                plt.figure("Reconstruct")
                
                ax1 = subplot(2,1,1)
            
                ax1.plot(np.arange(len(stim))*dt-1, current_re*1e3, 'b', linewidth=1)  
                ax1.plot(np.arange(len(stim))*dt-1, (stim)*1e3, 'k-', linewidth=1)
                ax1.plot(np.arange(len(stim))*dt-1, (stim_re_mat[0,:])*1e3, 'r', linewidth=1, alpha=1)
                                
                #adjust_spines(ax1, ['left','bottom'], d_out = 10) 
                #ax1.axis(xmin=0, xmax=1) 
                
                #ax1.axis(ymin=8.3, ymax=10.7)
                #ax1.yaxis.set_ticks(array([8.5,9,9.5,10,10.5]))
                #ax1.set_title("Reconstruction") 
                
                #ax1.set_xlabel("s") 
                #ax1.set_ylabel("pA")
                
                #ax1.text(0.15, 10.7, "Input current", color=color3, fontsize = 8)
                #ax1.text(0.8, 10.7, "Signal", color="#000000", fontsize = 8)
                #ax1.text(0.0, 8.2, "Reconstruction", color=color2, fontsize = 8)
                
                ax2 = subplot(2,1,2)
                ax2.plot(tk, K_mat_old[0], 'k', linewidth=1)  
 
            
                self.save_plot(directory = "./figs/dump/", prefix = "reconstruct")
            
            plt.figure("Transfer")
            
            currtitle = currlabel + " pop " + dowhat + ", " + self.celltype[self.a_celltype[0]]         
                   
            ax = plot_transfer(currtitle, freq_used, mag, pha, t1, current, voltage[self.a_celltype[0]], freq_times, spike_freq, taum, fmean, self.ihold, rm, Vreset, Vth, Vrest, method_interpol, method_interpol_plot, SNR = SNR, VAF = VAF, ymax = self.ymax, ax = self.ax, linewidth = self.linewidth, color_vec = self.color_vec, alpha = self.alpha, opt_plot = opt_plot)            
            
            suptitle("Population transfer function of " + str(self.N[self.a_celltype[0]]) + " " + self.celltype[self.a_celltype[0]] + ", amp: " + str(np.round(self.amp[self.a_celltype[0]],4)) + ", amod: " + str(self.amod) + ", ih: " + str(np.round(self.ihold,4)) + ", ih_s: " + str(np.round(self.ihold_sigma,4)) + ", fm: " + str(np.round(fmean,2)) + ", fl_s: " + str(self.fluct_s)) 
    
        return VAF, SNR, ax, tk, K_mat_old       
        

    def save_plot(self, directory = "./figs/dump/", prefix = " "):
        
        if pop.id == 0:
            
            from datetime import datetime
            idate = datetime.now().strftime('%Y%m%d_%H%M')  # %S
            savefig(directory + idate + "-pop_transfer_" + prefix + "_" + self.celltype[self.a_celltype[0]] + "_N" + str(self.N[self.a_celltype[0]]) + "_ihold" + str(np.round(self.ihold,4)) + "_amp" + str(np.round(self.amp[self.a_celltype[0]],4)) + ".pdf", dpi = 300)  # save it

                 
    def do_pca_ica(self, t_analysis_delay=0, t_analysis_stop=1, time=0, signals=0, output_dim=10, n_processes=32, n_chunks=32, do_ica=1, n_celltype = 0):
    
        if self.use_mpi:
            
            filepath = self.data_dir + "/" + str(self.pickle_prefix) + "_results_pop_pca_ica.p"
            
            if self.do_run or (os.path.isfile(filepath) is False):
                
                # PCA
                
                # remove beginning
                dt = time[2]-time[1]
                t = time[int(t_analysis_delay/dt):int(t_analysis_stop/dt)] 
                pca_mat = np.array(signals[n_celltype]).T[int(t_analysis_delay/dt):int(t_analysis_stop/dt),:]
                
                node = mdp.nodes.PCANode(output_dim=output_dim, svd=True)
                
                # pad with zeros to be able to split into chunks!
                n_add = n_chunks-np.remainder(np.shape(pca_mat)[0],n_chunks)
                mat_add = np.zeros((n_add, np.shape(pca_mat)[1]))
                pca_mat_add = np.concatenate((pca_mat, mat_add))
                pca_mat_iter = np.split(pca_mat_add, n_chunks) 
    
                flow = mdp.parallel.ParallelFlow([node])
                
                start_time = ttime.time()
                
                with mdp.parallel.ProcessScheduler(n_processes=n_processes, verbose=True) as scheduler:
                    flow.train([pca_mat_iter], scheduler=scheduler) # input has to be list, why??
                    
                process_time = ttime.time() - start_time
                
                s = np.array(flow.execute(pca_mat_iter))
                s = s[0:len(t),:]  # resize to length of t!
                
                #print "node.d: ",node.d
                var_vec = node.d/sum(node.d)
                print 'Explained variance (', 0, ') : ',  round(node.explained_variance,4)
                print 'Variance (' , 0, ') : ', var_vec
                print 'Time to run (' , 0, ') : ', process_time
                
                s2 = []
                if do_ica:
                # ICA
                    #s2 = mdp.fastica(s)
                    ica = mdp.nodes.FastICANode() #CuBICANode()
                    ica.train(s)
                    s2 = ica(s)
                
                results = {'t':t, 'pca':s,'pca_var':var_vec,'pca_var_expl':round(node.explained_variance,4), 'ica':s2}
                
                if self.id == 0:
                    if self.dumpsave == 1:
                        pickle.dump( results, gzip.GzipFile( filepath, "wb" ) )
            
            else:
            
                if self.id == 0:
                    results = pickle.load( gzip.GzipFile( filepath, "rb" ) ) 
                
        else:
            
            # remove beginning
            dt = time[2]-time[1]
            t = time[int(t_analysis_delay/dt):int(t_analysis_stop/dt)] 
            pca_mat = np.array(signals[n_celltype]).T[int(t_analysis_delay/dt):int(t_analysis_stop/dt),:]
            
            node = mdp.nodes.PCANode(output_dim=output_dim, svd=True)

            start_time = ttime.time()
            
            node.train(pca_mat)
            s = node(pca_mat)
            
            process_time = ttime.time() - start_time  
            #print "node.d: ",node.d
            var_vec = node.d/sum(node.d)
            print 'Explained variance (', 0, ') : ',  round(node.explained_variance,4)
            print 'Variance (' , 0, ') : ', var_vec
            print 'Time to run (' , 0, ') : ', process_time
            
            s2 = []
            if do_ica:
            # ICA
                #s2 = mdp.fastica(s)
                ica = mdp.nodes.FastICANode() #CuBICANode()
                ica.train(s)
                s2 = ica(s)
            
            results = {'t':t, 'pca':s,'pca_var':var_vec,'pca_var_expl':round(node.explained_variance,4), 'ica':s2}

        return results
        
        
    def net_run(self, tstop, simprop = "default", t_analysis_delay=0, t_analysis_stop=1, stim_start=0):

        freq_times = []
        t_all_vec_vec = []
        id_all_vec_vec = []
        gsyns = []
        w_mat = []
        winh_mat = []
        time = []
        voltage = []
        current = []
        
        filepath = self.data_dir + "/" + str(self.pickle_prefix) + "_results_pop_randomnet.hdf5"
        
        if self.do_run or (os.path.isfile(filepath) is False):
            
            self.run(tstop)
            
            self.no_fmean = True
            results = self.get()
            
            time, voltage, current, fmean, gsyn = results.get('time'), results.get('voltage'), results.get('current'), results.get('fmean'), results.get('gsyn')
            freq_times, spike_freq, t_all_vec_vec, id_all_vec_vec, gsyns, w_mat, winh_mat = results.get('freq_times'), results.get('spike_freq'), results.get('t_all_vec_vec'), results.get('id_all_vec_vec'), results.get('gsyns'), results.get('w_mat'), results.get('winh_mat')
            
            if self.id == 0:
                if self.dumpsave == 1:
                    #pickle.dump( results, open( filepath, "wb" ) ) # gzip.GzipFile
                    
                    print "- Saving", filepath
                    
                    f = h5py.File(filepath, 'w')
                    f.create_dataset('time', data=time, compression='gzip', shuffle=True)
                    f.create_dataset('voltage', data=np.array(voltage), compression='gzip', shuffle=True)
                    f.create_dataset('current', data=current, compression='gzip', shuffle=True)
                    f.create_dataset('freq_times', data=freq_times, compression='gzip', shuffle=True)
                                        
                    #f.create_dataset('t_all_vec_vec', data=np.array(t_all_vec_vec), compression='lzf', shuffle=True)
                    #f.create_dataset('id_all_vec_vec', data=np.array(id_all_vec_vec), compression='lzf', shuffle=True)
                    #f.create_dataset('gsyns', data=np.array(gsyns), compression='lzf', shuffle=True)

                    for i in range(len(self.N)):
                        subgroup = f.create_group("cell" + str(i))
                        subgroup.create_dataset('t_all_vec_vec', data=t_all_vec_vec[i], compression='gzip', shuffle=True)
                        subgroup.create_dataset('id_all_vec_vec', data=id_all_vec_vec[i], compression='gzip', shuffle=True)
                        subgroup.create_dataset('g', data=gsyns[i], compression='gzip', shuffle=True)

                        #for j in range(len(gsyns[i])):
                        #    subsubgroup = subgroup.create_group("gsyn" + str(j))
                        #    subsubgroup.create_dataset('g', data=gsyns[i][j], compression='lzf', shuffle=True)
                        
                    f.close()   
                    print "- Save finished"
                    
                    #filename = slugify(simprop)

                    #syn_grc = np.array(gsyns[0])
                    
                    #import scipy
                    #from scipy import io
                    
                    #print "Saving .mat"
                    #data = {}
                    #data['syn_grc'] = syn_grc[:,int(t_analysis_delay/self.bin_width):int(t_analysis_stop/self.bin_width)]
                    #data['time'] = freq_times[int(t_analysis_delay/self.bin_width):int(t_analysis_stop/self.bin_width)]-stim_start
                    #scipy.io.savemat('./figs/' + filename + '.mat',data)
                        
        else:
                
            if self.id == 0:
                #results = pickle.load( open( filepath, "rb" ) ) #gzip.GzipFile
                f = h5py.File(filepath, 'r')
                
                time = np.array(f['time'])
                voltage = np.array(f['voltage'])
                current = np.array(f['current'])
                freq_times = np.array(f['freq_times'])
                
                
                for i in range(len(self.N)):
                    t_all_vec_vec.append(np.array(f['/cell' + str(i) + '/t_all_vec_vec'])) 
                    id_all_vec_vec.append(np.array(f['/cell' + str(i) + '/id_all_vec_vec'])) 
                    gsyns.append(np.array(f['/cell' + str(i) + '/g'])) 
                    
                    #gsyns.append([])
                    #for j in range(self.N[i]):
                    #    gsyns[i].append(np.array(f['/cell' + str(i) + '/gsyn' + str(j) + '/g' ])) 

                f.close()
                
        return time, voltage, current, t_all_vec_vec, id_all_vec_vec, gsyns, freq_times, w_mat, winh_mat 

        
    def delall(self):        
        
        if self.use_mpi: 
            self.pc.gid_clear()
            print "- clearing gids"
        else:
            pass
            #h.topology()            
            #for sec in h.allsec():
            #    print "- deleting section:", sec.name()
            #    #h("%s{delete_section()}"%sec.name())
            #    sec.push()
            #    h.delete_section()
            #h.topology()
            
        for n in range(self.n_celltypes): 
            for m in self.cells[n]:
                m.destroy()
                del m 
        del self.cells
        del self.nc_vecstim
        del self.netcons
        del self.nclist
        print h.topology()   
        
        
    def delrerun(self):        
        
        del self.nc_vecstim
        del self.netcons
        del self.nclist
        del self.vecstim
        del self.spike_vec
        del self.ST_stims
        del self.PF_stims
        
        self.netcons = [] 
        self.nclist = []
        self.nc_vecstim = []
        self.vecstim = []
        self.spike_vec = []
        self.ST_stims = []
        self.PF_stims = []
        
        self.t_vec = []
        self.id_vec = []
        self.rec_v = []
                
        for n in range(self.n_celltypes):
            if self.use_mpi:
                self.t_vec.append(h.Vector()) # np.array([0])
                self.id_vec.append(h.Vector()) # np.array([-1], dtype=int)
            else:
                self.t_vec.append([])
                
            self.rec_v.append(h.Vector())

            for cell in self.cells[n]:
                self.t_vec[n].append(h.Vector())
                cell.nc_spike.record(self.t_vec[n][-1])                

        self.flucts = []            # Fluctuating inputs on this host
        self.noises = []            # Random number generators on this host
        self.plays = []             # Play inputs on this host
        self.rec_is = []
        self.trains = []     
        
        self.ic_holds = []
        self.i_holdrs = []
        self.i_holds = []
        self.ic_starts = [] 
        self.vc_starts = []
        self.ic_steps = []
        self.tvecs = []
        self.ivecs = []  
        self.noises = []
        self.record_syn = []
        self.id_all_vec_input = []
        self.t_all_vec_input = []
        self.syn_ex_dist = []
        self.syn_inh_dist = []

            
# test code
if __name__ == '__main__':
    
    # mpiexec -f ~/machinefile -enable-x -n 96 python Population.py --noplot
    
    from Stimulation import *
    from Plotter import *
    from Stimhelp import *

    from cells.IfCell import *
    import scipy
    from scipy import io
            
    dt = 0.1*ms
    dt = 0.025*ms
    
    do_run = 1
    if results.norun:  # do not run again use pickled files!
        print "- Not running, using saved files"
        do_run = 0
    
    
    do = np.array(["transfer"])
    opts = np.array(["if_cnoise", "grc_cnoise"]) #ssine  
    #opts = np.array(["if_cnoise"]) #ssine
    #opts = np.array(["if_recon"]) #ssine
    opts = np.array(["if_syn_CFvec"]) 
    #opts = np.array(["prk_cnoise"])
    opts = np.array(["if_cnoise", "if_ssine"]) #ssine 
    opts = np.array(["if_ssine"]) #ssine 
    opts = np.array(["grc_cnoise_addn_cn_", "grc_cnoise_cn_", "grc_cnoise_addn_cn_a01"]) 
    opts = np.array(["grc_cnoise_addn100_cn_", "grc_cnoise_addn_cn_", "grc_cnoise_cn_"]) 
    opts = np.array(["grc_cnoise_addn100_cn_"])
    opts = np.array(["grc_cnoise_addn100_"])
    opts = np.array(["grc_cnoise_addn_cn_"])
    #opts = np.array(["grc_cnoise"])
    #opts = np.array(["grc_cnoise_cn", "grc_cnoise_addn_cn"]) 
    #opts = np.array(["if_cnoise_addn", "if_cnoise"]) 
    
    do = np.array(["timeconst"])
    
    #do = np.array(["transfer"])
    #opts = np.array(["grc_cnoise_syn"])
    #opts = np.array(["grc_recon_syn"])
    
    #do = np.array(["prk_test"])
    
        
    if "prk_test" in do:
        
        import multiprocessing
        from Purkinje import Purkinje
        cell = Purkinje()  

        # set up recording
        # Time
        rec_t = h.Vector()
        rec_t.record(h._ref_t)
        
        # Voltage
        rec_v = h.Vector()
        rec_v.record(cell.soma(0.5)._ref_v)

        tstop = 500
        v_init = -60
        
        stim = h.IClamp(cell.soma(0.5))
        stim.amp = 0.0/nA
        stim.delay = 1
        stim.dur = 1000
        
        cpu = multiprocessing.cpu_count()
        h.load_file("parcom.hoc")
        p = h.ParallelComputeTool()
        p.change_nthread(cpu,1)
        p.multisplit(1)
        print 'cpus:', cpu
        
        h.load_file("stdrun.hoc")
        h.celsius = 37              
        h.init()
        h.tstop = tstop
        dt = 0.025 # ms
        h.dt = dt
        h.steps_per_ms = 1 / dt     
        h.v_init = v_init  
        
        h.finitialize()
        h.run()
    
        t1 = np.array(rec_t)
        voltage = np.array(rec_v)
        s, spike_times = get_spikes(voltage, -20, t1)

        print 1000/diff( spike_times)

        plt.figure()
        plt.subplot(2,1,1)
        plt.plot(t1, voltage)
                
        plt.show()


    if "transfer" in do:
        
        # SET DEFAULT VALUES FOR THIS PLOT
        fig_size =  [11.7, 8.3]
        params = {'backend': 'ps', 'axes.labelsize': 9, 'axes.linewidth' : 0.5, 'title.fontsize': 8, 'text.fontsize': 9,
            'legend.borderpad': 0.2, 'legend.fontsize': 8, 'legend.linewidth': 0.1, 'legend.loc': 'best', # 'lower right'    
            'legend.ncol': 4, 'xtick.labelsize': 8, 'ytick.labelsize': 8, 'text.usetex': False, 'figure.figsize': fig_size}
        rcParams.update(params)        
        
        
        freq_used0 = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 35, 40, 45, 50, 55, 60, 65, 70, 80, 100, 1000])*Hz
        #freq_used0 = np.concatenate((arange(0.1, 1, 0.1), arange(1, 501, 1) ))
        freq_used0 = np.array([1, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 64, 66, 68, 70, 72, 74, 76, 78, 80, 82, 84, 86, 88, 90, 92, 94, 96, 98, 100, 200, 400, 600, 800, 1000])
                            
        SNR = None 
        NI = None
        VAF = None
        
        t_stim = 1000*s  # only for cnoise     
        
        opt_plot = np.array(["only_mag","normalize", "dB"]) #  
        #opt_plot = np.array(["normalize", "dB"]) # 
        
        color_vec = (np.array(["Red", "Blue", "HotPink", "Indigo"]), np.array(["Blue", "Orange", "HotPink", "Indigo"]))
        #color=cm.jet(1.*i/x)
        
        method_interpol = np.array(['bin','syn'])   
        method_interpol = np.array(['bin'])   
        
        for i, o in enumerate(opts):
            
            dt = 0.025*ms
            bin_width = 5*ms
            bin_width = dt
            jitter = 0*ms
            
            n_syn_ex = [0] 
            g_syn_ex = [1]
            noise_syn = 0 
            inh_hold = 0 
            n_syn_inh = [0] 
            g_syn_inh = [1]
            tau1_ex = 0
            tau2_ex = 10*ms
            tau1_inh = 0
            tau2_inh = 100*ms
            
            cutf = 20
            sexp = -1

            cutf = 0
            sexp = 0
            
            ihold = [10]
            amod = 0.1 # relative value
            give_freq = True
            
            anoise = [0]
            fluct_tau = 0*ms 
            
            N = [100]
            
            amp = 0 # absolute value
            fluct_s = [0] # absolute value 0.0008
            ihold_sigma = [0] # 0.01 absolute value
            
            CF_var = [[5,10,20]]
            CF_var = False
            
            syn_tau1 = 5*ms
            syn_tau2 = 5*ms
            
            do_csd = 1
            
            if "if" in o:
                
                do_csd = 1
                
                color_vec = (np.array(["Blue"]), np.array(["Blue"]))
                #color_vec = (np.array(["Red"]), np.array(["Red"]))
                
                cellimport = []
                celltype = ["IfCell"]
                #cell_exe = ["cell = IfCell()"]
                #cell_exe = ["cell = IfCell(e = -70*mV, thresh = -69*mV, vrefrac = -70*mV)"]  
                #cell_exe = ["cell = IfCell(e = 0*mV, thresh = 1*mV, vrefrac = 0*mV)"]
                
                # Brunel
                #cell_exe = ["cell = IfCell(C = 0.0005 *uF, R = 40*MOhm, e = -70*mV, thresh = -50*mV, vrefrac = -56*mV); cell.add_resonance(tau_r = 100*ms, gr = 0.025*uS)"]  
                
                #cell_exe = ["cell = IfCell(C = 0.0001*uF, R = 40*MOhm, sigma_C = 0.2, sigma_R = 0.2)"] 
                #cell_exe = ["cell = IfCell(C = 0.0001*uF, R = 40*MOhm)"] # tau = 4 ms
                #cell_exe = ["cell = IfCell(C = 0.0001*uF, R = 40*MOhm, s_reset_noise = 0*mV)"] # tau = 4 ms
                
                #GrC resting: 737 MOhm, 2.985e-06 uF   tau: 0.0022 s
                #GrC transfer fit: tau: 0.027 s => with 2.985e-06 uF, R = 0.027/2.985e-12 = 9045 MOhm
                
                #cell_exe = ["cell = IfCell(C = 2.985e-06*uF, R = 9045*MOhm)"] 
                
                thresh = -41.8    
                R = 5227*MOhm
                #tau_passive = 3e-06*5227 = 15.7ms
            
                cell_exe = ["cell = IfCell(C = 3.0e-06*uF, R = " + str(R) + ", e = -71.5*mV, thresh =" + str(thresh) + ", vrefrac = -71.5*mV)"]
                
                prefix = "if_tf"
 
                istart = 0 
                istop = 0.01
                di = 0.00001
                
                
                syn_tau1 = 10*ms
                syn_tau2 = 10*ms
                
                # Indirect
                give_freq = True
                ihold = [40]
                amod = 1 # relative value
                
                anoise = [0] 
                fluct_tau = 0*ms 
                
                #anoise = 0.1
                #fluct_tau = 100*ms
                
#                # Direct
#                give_freq = False
#                ihold = [0.00569223341176]
#                amod = None
#                amp = 7.31353725e-06
#                
#                anoise = None
#                fluct_s = [3.65676863e-06] 
#                fluct_tau = 0*ms
#                
#                # Low CF, No low noise
#                N = [10000]
#                give_freq = False
#                ihold = [0.004]
#                ihold_sigma = [0.1/2] # 0.1/2 0.01 realtive value
#                amod = None
#                amp = 0.0021
#                
#                anoise = None
#                fluct_s = [0.00]  # .005
#                fluct_tau = 0*ms
                
                
#                # Low CF, With low noise
#                N = [10000]
#                give_freq = False
#                ihold = [0.002]
#                ihold_sigma = [0.1/2] # 0.1/2 0.01 realtive value
#                amod = None
#                amp = 0.001
#                
#                anoise = None
#                fluct_s = [0.002]  # .005
#                fluct_tau = 100*ms
                
            if "resif" in o:
                
                do_csd = 1
                
                color_vec = (np.array(["Blue"]), np.array(["Blue"]))
                #color_vec = (np.array(["Red"]), np.array(["Red"]))
                
                cellimport = []
                celltype = ["IfCell"]
                
                gr = 5.56e-05*uS 
                tau_r = 19.6*ms
                R = 5227*MOhm
                delta_t = 4.85*ms
                thresh = (0.00568*nA * R) - 71.5*mV # 
                thresh = -41.8          
    
                cellimport = []
                celltype = "IfCell"
                cell_exe = "cell = IfCell(C = 3e-06*uF, R = " + str(R) + ", e = -71.5*mV, thresh =" + str(thresh) + ", vrefrac = -71.5*mV, dgk =" + str(gr) + ", egk = -71.5*mV, ctau =" + str(tau_r) + ")"

                prefix = "resif_tf"
 
                istart = 0 
                istop = 0.01
                di = 0.00001
                
                
                syn_tau1 = 10*ms
                syn_tau2 = 10*ms
                
                # Indirect
                give_freq = True
                ihold = [40]
                amod = 1 # relative value
                
                anoise = [0] 
                fluct_tau = 0*ms 
                dt = 0.1*ms
                
            
                
            if "if_syn" in o:
                
                N = [1]    
                ihold = [40] 
                amod = 1 # relative value
                
                prefix = "if_syntf"   
                
                n_syn_ex = 1 

                g_syn_ex = 0 
                
                noise_syn = 0

                fluct_tau = 0*ms 
                
                freq_used = np.array([])
                
                tau1_ex=0*ms
                tau2_ex=10*ms
                
                anoise = [0]

                
            if "grc" in o:
                
                color_vec = (np.array(["Blue"]), np.array(["Blue"]))

                cellimport = ["from GRANULE_Cell import Grc"]
                celltype = ["Grc"]
                cell_exe = ["cell = Grc(np.array([0.,0.,0.]))"]    
        
                prefix = "grc_tf"    

                istart = 0 
                istop = 0.1
                di = 0.01
                
                syn_tau1 = 10*ms
                syn_tau2 = 10*ms
                
                # Indirect
                give_freq = True
                ihold = [40]
                amod = 1 # relative value
                
                anoise = [0]
                fluct_tau = 0*ms 
                
                #anoise = 0.1
                #fluct_tau = 100*ms
                
#                # Direct
#                give_freq = False
#                ihold = [0.0058021085712642992]
#                amod = None
#                amp = 7.31353725e-06
#                
#                anoise = None
#                fluct_s = [3.65676863e-06] 
#                fluct_tau = 0*ms
#                
#                # Low CF, No low noise
#                N = [50]
#                give_freq = False
#                ihold = [0.0049]
#                ihold_sigma = [0.1/2] # 0.1/2 0.01 realtive value
#                amod = None
#                amp = 0.0021
#                
#                anoise = None
#                fluct_s = [0.00]  # .005
#                fluct_tau = 0*ms
#                
#                
#                # Low CF, With low noise
#                N = [10000]
#                give_freq = False
#                ihold = [0.003]
#                ihold_sigma = [0.1/2] # 0.1/2 0.01 realtive value
#                amod = None
#                amp = 0.001
#                
#                anoise = None
#                fluct_s = [0.002]  # .005
#                fluct_tau = 100*ms
                
            
            use_multisplit = False
            use_mpi = True
            simstep = 1*s
                
            if "prk" in o:
                
                N = [1]    
                ihold = [60] 
                
                color_vec = (np.array(["Blue"]), np.array(["Blue"]))

                cellimport = ["from Purkinje import Purkinje"]
                celltype = ["Prk"]
                cell_exe = ["cell = Purkinje()"]    
        
                prefix = "prk_tf"    

                temperature = 37

                istart = 0 
                istop = 0.1
                di = 0.005
            
                use_multisplit = True
                use_mpi = False
                
                t_stim = 5*s  # only for cnoise 
                simstep = 1*s


            if "grc_syn" in o:
                
                N = [1]    
                ihold = [125] 
                amod = 1 # relative value
                
                prefix = "grc_syntf"    
                
                cutf = 20
                sexp = -1
                
                cutf = 0
                sexp = 0
                                              
                n_syn_ex = 1 
                g_syn_ex = -1
                noise_syn = 1

                n_syn_inh = -1
                inh_hold = 0
                g_syn_inh = 0
                
                fluct_tau = 0*ms 
                
                freq_used = np.array([])
                
                anoise = 0
                
            
            if "_addn" in o:
                
                anoise = [6] # RESPONSIBLE FOR FILTERING EFFECT!!!
                fluct_tau = 1*ms 
                prefix = prefix + "_addn"
                color_vec = (np.array(["Red"]), np.array(["Red"]))
                
            if "_addn100" in o:
                
                anoise = [2] # RESPONSIBLE FOR FILTERING EFFECT!!!
                fluct_tau = 100*ms 
                prefix = prefix + "100"
                color_vec = (np.array(["Green"]), np.array(["Green"]))
                
            if "_cn_" in o:
                
                cutf = 20
                sexp = -1
                prefix = prefix + "_cn"
                
            if "_a01" in o:
                
                amod=0.1
                prefix = prefix + "_a01"
                

            
            plt.figure(i)
            pickle_prefix = "Population.py_" + prefix
            
            #comm = MPI.COMM_WORLD
            #comm.Barrier()   # wait for other nodes
            
            pop = Population(cellimport = cellimport,  celltype = celltype, cell_exe = cell_exe, N = N, temperature = 37, ihold = ihold, ihold_sigma = ihold_sigma, amp = amp, amod = amod, give_freq = give_freq, do_run = do_run, pickle_prefix = pickle_prefix, istart = istart, istop = istop, di = di, dt = dt)      
            pop.bin_width = bin_width
            pop.jitter = jitter
            pop.anoise = anoise
            pop.fluct_s = fluct_s 
            pop.fluct_tau = fluct_tau 
            pop.method_interpol = method_interpol  
            pop.no_fmean = False
            pop.CF_var = CF_var
            
            pop.tau1_ex=tau1_ex
            pop.tau2_ex=tau2_ex
            pop.tau1_inh=tau1_inh
            pop.tau2_inh=tau2_inh
            
            pop.n_syn_ex = n_syn_ex 
            pop.g_syn_ex = g_syn_ex 
            
            pop.noise_syn = noise_syn 
            pop.inh_hold = inh_hold 
            pop.n_syn_inh = n_syn_inh 
            pop.g_syn_inh = g_syn_inh
            
            pop.force_run = False
            pop.use_multisplit = use_multisplit
            pop.use_mpi = use_mpi
            pop.simstep = simstep
            pop.use_local_dt = False
            pop.syn_tau1 = syn_tau1
            pop.syn_tau2 = syn_tau2
            pop.plot_input = False
            
            
            if n_syn_inh == -1:
                pop.connect_gfluct(g_i0=g_syn_inh)
                
            #pop.test_mod(n_syn_ex = n_syn_ex, g_syn_ex = g_syn_ex, noise_syn = noise_syn, inh_hold = inh_hold, n_syn_inh = n_syn_inh, g_syn_inh = g_syn_inh, do_plot = True)
             
            if "ssine" in o:
                pop.color_vec = color_vec
                #pop.color_vec = (np.array(["Red", "Orange", "HotPink", "Indigo"]), np.array(["Red", "Orange", "HotPink", "Indigo"]))     
                pop.fun_plot(currlabel = "control", dowhat = "ssine", freq_used = freq_used0, opt_plot = opt_plot)

                pop.save_plot(directory = "./figs/dump/")          
           
            if "cnoise" in o:
                
                freq_used = np.array([])
                pop.color_vec = color_vec
                #pop.color_vec = (np.array(["Blue", "Green", "DimGray", "DarkGoldenRod"]), np.array(["Blue", "Green", "DimGray", "DarkGoldenRod"]))    
                pop.fun_plot(currlabel = "control", dowhat = "cnoise", t_stim = t_stim, cutf = cutf, sexp = sexp, t_qual = 0, opt_plot = opt_plot, freq_used = freq_used, do_csd = do_csd)
            
                pop.save_plot(directory = "./figs/dump/")  
                
                
            if "recon" in o:
                
                pop.color_vec = color_vec 
                #VAF, SNR, ax, tk, K_mat_old = pop.fun_plot(currlabel = "control", dowhat = "cnoise", t_stim = t_stim, cutf = cutf, sexp = sexp, t_qual = 0, opt_plot = opt_plot, n_syn_ex = n_syn_ex, g_syn_ex = g_syn_ex, noise_syn = noise_syn, inh_hold = inh_hold, n_syn_inh = n_syn_inh, g_syn_inh = g_syn_inh, SNR=0, freq_used = freq_used)
                
                # RECONSTRUCT!
                freq_used = np.array([9, 47, 111, 1000])*Hz
                t_stim = 10*s

                tk = arange(0,0.8192*2,pop.dt)
                K_mat_old = zeros((len(method_interpol),len(tk)), dtype=complex)
                
                if pop.id == 0:

                    sigma = 0.1e-3
                    a=0.1
                    t0 = tk[floor(len(tk)/2)]
                    K_mat_old[0] = gauss_func(tk, a, t0, sigma)
                    
                K_mat_old = np.array([])

                results = pop.fun_cnoise_Stim(t_stim = t_stim, cutf = cutf, sexp = sexp, t_qual = 5, n_syn_ex = n_syn_ex, g_syn_ex = g_syn_ex, noise_syn = noise_syn, inh_hold = inh_hold, n_syn_inh = n_syn_inh, g_syn_inh = g_syn_inh, freq_used = freq_used, K_mat_old = K_mat_old, seed = 311)
                
                freq_used, amp_vec, mag, pha, ca, voltage, current, t1 = results.get('freq_used'), results.get('amp'), results.get('mag'), results.get('pha'), results.get('ca'), results.get('voltage'), results.get('current'), results.get('t1') 
                freq_times, spike_freq, fmean, method_interpol, SNR, VAF, Qual = results.get('freq_times'), results.get('spike_freq'), results.get('fmean'), results.get('method_interpol'), results.get('SNR'), results.get('VAF'), results.get('Qual')                              
                stim, resp_mat, stim_re_mat = results.get('stim'), results.get('resp_mat'), results.get('stim_re_mat')
                
                if pop.id == 0:
                
                    plt.figure('Reconstruct')
                    axR0 = plt.subplot(4,1,1)
                    axR1 = plt.subplot(4,1,2)
                    axR2 = plt.subplot(4,1,3)
                    axR3 = plt.subplot(4,1,4)
                    
                    axR0.plot(np.arange(len(stim))*pop.dt, resp_mat[0,:])
                    axR0.axis(xmin=0.9, xmax=1)
                    #axR0.plot(t1, voltage[0])
                    axR1.plot(np.arange(len(stim))*pop.dt, stim, 'b')
                    axR1.axis(xmin=0.9, xmax=1)
                    axR2.plot(np.arange(len(stim))*pop.dt, stim_re_mat[0,:], 'r')
                    axR2.axis(xmin=0.9, xmax=1)
                    axR3.plot(tk, K_mat_old[0])
                    plt.savefig("./figs/dump/Reconstruct.pdf", dpi = 300, transparent=True) # save it
            
            pop = None
            
            
    plt.show()
    
   
    if "timeconst" in do:
        
        from lmfit import minimize, Parameters
        
        # SET DEFAULT VALUES FOR THIS PLOT
        fig_size =  [11.7, 8.3]
        params = {'backend': 'ps', 'axes.labelsize': 9, 'axes.linewidth' : 0.5, 'title.fontsize': 8, 'text.fontsize': 9,
            'legend.borderpad': 0.2, 'legend.fontsize': 8, 'legend.linewidth': 0.1, 'legend.loc': 'best', # 'lower right'    
            'legend.ncol': 4, 'xtick.labelsize': 8, 'ytick.labelsize': 8, 'text.usetex': False, 'figure.figsize': fig_size}
        rcParams.update(params)     
        
        dt = 0.025*ms
                    
        prefix = "timeconst"
        pickle_prefix = "Population.py_" + prefix
        
        stimtype = "inh_50ms_20ms"
                    
        if stimtype == "ex_20ms":
            
            trun = 2.9
            tstart = 1.8
            tstop = 2.7

            celltype = ["IfCell"]
            cell_exe = ["cell = IfCell(C = 0.0001*uF, R = 200*MOhm)"]
            N = [5000]
                    
            pop = Population(celltype = celltype, cell_exe = cell_exe, N = N, temperature = 0, do_run = do_run, pickle_prefix = pickle_prefix, dt = dt)      
            
            pop.method_interpol = np.array(["bin", "syn"])
            pop.method_interpol = np.array(["bin"])
            
            modulation_vec = pop.set_PulseStim(start_time=[100*ms], dur=[3000*ms], steadyf=[100*Hz], pulsef=[150*Hz], pulse_start=[2000*ms], pulse_len=[500*ms], weight0=[1*nS], tau01=[0*ms], tau02=[20*ms], weight1=[0*nS], tau11=[0*ms], tau12=[1*ms])
            
            params = Parameters()
            params.add('amp', value=0.1)
            params.add('shift', value=10)
            params.add('tau1', value=1, vary=False) # alpha! 
            params.add('tau2', value=20*ms) 
            
            
        if stimtype == "ex_gr":
            
            trun = 6.9
            tstart = 4.8
            tstop = 6.5

            cellimport = ["from GRANULE_Cell import Grc"]
            celltype = ["Grc"]
            cell_exe = ["cell = Grc(np.array([0.,0.,0.]))"]
            N = [4096*10]
                    
            pop = Population(cellimport = cellimport, celltype = celltype, cell_exe = cell_exe, N = N, temperature = 37, do_run = do_run, pickle_prefix = pickle_prefix, dt = dt)      
            
            pop.method_interpol = np.array(["bin", "syn"])
            pop.method_interpol = np.array(["bin"])
            
            modulation_vec = pop.set_PulseStim(start_time=[100*ms], dur=[7000*ms], steadyf=[20*Hz], pulsef=[30*Hz], pulse_start=[5000*ms], pulse_len=[500*ms])
            
            params = Parameters()
            params.add('amp', value=0.1)
            params.add('shift', value=10)
            params.add('tau1', value=1, vary=False) # alpha! 
            params.add('tau2', value=20*ms) 
            
            
        if stimtype == "inh_50ms_20ms":
            
            trun = 2.9
            tstart = 1.8
            tstop = 2.7
            
            celltype = ["IfCell", "IfCell"]
            cell_exe = ["cell = IfCell()", "cell = IfCell()"]
        
            N = [10000,10000]
            
            pop = Population(celltype = celltype, cell_exe = cell_exe, N = N, temperature = 0, do_run = do_run, pickle_prefix = pickle_prefix, dt = dt)  
            
            pop.method_interpol = np.array(["bin", "syn"])
            pop.method_interpol = np.array(["bin"])
            
            modulation_vec = pop.set_PulseStim(start_time=[100*ms,100*ms], dur=[3000*ms,3000*ms], steadyf=[100*Hz,50*Hz], pulsef=[100*Hz,80*Hz], pulse_start=[2000*ms,2000*ms], pulse_len=[500*ms,500*ms], weight0=[1*nS,1*nS], tau01=[1*ms,1*ms], tau02=[20*ms,20*ms], weight1=[0,0], tau11=[0*ms,0*ms], tau12=[1*ms,1*ms])

            pop.connect_cells(conntype='inh', weight=0.001, tau=50)
            
            params = Parameters()
            params.add('amp', value=-0.1)
            params.add('shift', value=10)
            params.add('tau1', value=1, vary=False) # alpha! 
            params.add('tau2', value=20*ms)
            
            
        if stimtype == "inh_gr":

            trun = 9.9                
            tstart = 4.8
            tstop = 8
            
            cellimport = ["from GRANULE_Cell import Grc", "from templates.golgi.Golgi_template import Goc"]
            celltype = ["Grc","Goc_noloop"]
            cell_exe = ["cell = Grc(np.array([0.,0.,0.]))","cell = Goc(np.array([0.,0.,0.]))"]
            N = [100,4]
            #N = [4096, 27]
            #N = [4096*5, 27*5]

            pop = Population(cellimport = cellimport, celltype = celltype, cell_exe = cell_exe, N = N, temperature = 37, do_run = do_run, pickle_prefix = pickle_prefix, dt = dt)  
            
            pop.method_interpol = np.array(["bin", "syn"])
            pop.method_interpol = np.array(["bin"])
            
            modulation_vec = pop.set_PulseStim(start_time=[100*ms,100*ms], dur=[9800*ms,9800*ms], steadyf=[60*Hz,10*Hz], pulsef=[60*Hz,22*Hz], pulse_start=[5000*ms,5000*ms], pulse_len=[1500*ms,1500*ms])

            pop.connect_cells(conntype='inh_gr', weight = 0.3)
            
            params = Parameters()
            params.add('amp', value=-0.1)
            params.add('shift', value=10)
            params.add('tau1', value=1, vary=False) # alpha! 
            params.add('tau2', value=20*ms)
            
            
        if stimtype == "inh_50ms_curr":
        
            trun = 2.9
            tstart = 1.8
            tstop = 2.8
            
            celltype = ["IfCell", "IfCell"]
            cell_exe = ["cell = IfCell()", "cell = IfCell()"]
            
            N = [1000,1000]
            
            give_freq = True
                            
            istart = 0 
            istop = 0.2
            di = 0.01
            
            ihold = [100, 50] 
            ihold_sigma = [0.01, 0.01] # relative sigma
             
            pop = Population(celltype = celltype, cell_exe = cell_exe, N = N, temperature = 0, ihold = ihold, ihold_sigma = ihold_sigma, give_freq = give_freq, do_run = do_run, pickle_prefix = pickle_prefix, istart = istart, istop = istop, di = di, dt = dt)  
            
            pop.method_interpol = np.array(["bin", "syn"])
            pop.method_interpol = np.array(["bin"])
            
            tstep = 2 
            tdur = 0.5
            
            istep = [100,100]
            current1 = np.concatenate(([ihold[1]*np.ones(round((tstep)/pop.dt)), istep[1]*np.ones(round(tdur/pop.dt)),ihold[1]*np.ones(round((trun-tstep-tdur)/pop.dt)) ])) 
            
            pop.set_IStim()
            pop.set_IStep(istep = istep, istep_sigma = [0.01,0.01], tstep = tstep, tdur = tdur)
            
            pop.connect_cells(conntype='inh', weight=0.0003, tau=50)
                            
            pop.fluct_s = [0.02,0.05]
            pop.connect_fluct() 
            
            params = Parameters()
            params.add('amp', value=-0.1)
            params.add('shift', value=10)
            params.add('tau1', value=1, vary=False) # alpha! 
            params.add('tau2', value=20*ms)
            
            
        if stimtype == "inh_gr_curr":
            
            trun = 9.9                
            tstart = 4.8
            tstop = 8
            
            cellimport = ["from GRANULE_Cell import Grc", "from templates.golgi.Golgi_template import Goc"]
            celltype = ["Grc","Goc_noloop"]
            cell_exe = ["cell = Grc(np.array([0.,0.,0.]))","cell = Goc(np.array([0.,0.,0.]))"]
            N = [100,4]
            N = [4096, 27]
            N = [4096*10, 27*10]            

            give_freq = True
            
            # GRC                
            #istart = 0 
            #istop = 0.1
            #di = 0.01
            
            #GOC
            istart = 0 
            istop = 0.5
            di = 0.02
                           
            ihold = [100, 10] 
            ihold_sigma = [0, 0] # relative sigma
             
            pop = Population(cellimport = cellimport, celltype = celltype, cell_exe = cell_exe, N = N, temperature = 37, ihold = ihold, ihold_sigma = ihold_sigma, give_freq = give_freq, do_run = do_run, pickle_prefix = pickle_prefix, istart = istart, istop = istop, di = di, dt = dt)  
            
            pop.method_interpol = np.array(["bin", "syn"])
            pop.method_interpol = np.array(["bin"])
            
            tstep = 5 
            tdur = 2
            
            istep = [100,50]
            current1 = np.concatenate(([ihold[1]*np.ones(round((tstep)/pop.dt)), istep[1]*np.ones(round(tdur/pop.dt)),ihold[1]*np.ones(round((trun-tstep-tdur)/pop.dt)) ])) 
            
            pop.set_IStim()
            pop.set_IStep(istep = istep, istep_sigma = [0,0], tstep = tstep, tdur = tdur)
            
            pop.connect_cells(conntype='inh_gr', weight = 0.4)
                            
            pop.fluct_s = [0.05,2]
            pop.connect_fluct() 
            
            params = Parameters()
            params.add('amp', value=-0.1)
            params.add('shift', value=10)
            params.add('tau1', value=1, vary=False) # alpha! 
            params.add('tau2', value=20*ms)               
           
        
        pop.run_steps(trun)
        
        self.no_fmean = True
        results = pop.get()
        time, voltage, current, fmean, gsyn = results.get('time'), results.get('voltage'), results.get('current'), results.get('fmean'), results.get('gsyn')
        freq_times, spike_freq, t_all_vec_vec, id_all_vec_vec, gsyns = results.get('freq_times'), results.get('spike_freq'), results.get('t_all_vec_vec'), results.get('id_all_vec_vec'), results.get('gsyns')
        
        if pop.id == 0:
            
            bin_width = 1*ms
            freq_times = arange(0, time[-1], bin_width)
            [num_spikes, _] = neuronpy.util.spiketrain.get_histogram(t_all_vec_vec[0], bins = freq_times)
            spike_freq = np.concatenate((zeros(1),num_spikes)) / bin_width / N[0]
            
            
            if "inh" in stimtype: # generate input current, to complicated to get it out
            
                if "curr" in stimtype:
                    time1 = np.arange(0, trun, pop.dt)
                                                                   
                    r_mod = interp(freq_times, time1, current1, left=0, right=0)
                    
                    [num_spikes, _] = neuronpy.util.spiketrain.get_histogram(t_all_vec_vec[1], bins = freq_times)
                    spike_freq1 = np.concatenate((zeros(1),num_spikes)) / bin_width / N[1]
                else:
                    r_mod = interp(freq_times, modulation_vec[1][0], modulation_vec[1][1], left=0, right=0)
                
                    [num_spikes, _] = neuronpy.util.spiketrain.get_histogram(t_all_vec_vec[1], bins = freq_times)
                    spike_freq1 = np.concatenate((zeros(1),num_spikes)) / bin_width / N[1]
                    
            elif "ex" in stimtype:
                r_mod = interp(freq_times, modulation_vec[0][0], modulation_vec[0][1], left=0, right=0) 


            def modelfun(amp, shift, tau1, tau2, bin_width, r_mod):
                
                tau1 = tau1
                tau2 = tau2
            
                t1 = np.arange(0,10*tau2,bin_width)
                K = amp*syn_kernel(t1, tau1, tau2) 
                K = np.concatenate((np.zeros(len(K)-1),K))
                t2 = np.arange(0,len(K)*bin_width,bin_width)
                
                model = np.convolve(K, r_mod, mode='same') + shift
                
                return model

                            
            def residual(params, r_mod, data=None, bin_width=1*ms, tstart=0, tstop=3):
                
                amp = params['amp'].value
                shift = params['shift'].value
                tau1 = params['tau1'].value
                tau2 = params['tau2'].value
                
                model = modelfun(amp, shift, tau1, tau2, bin_width, r_mod)
                
                return (data[int(tstart/bin_width):int(tstop/bin_width)]-model[int(tstart/bin_width):int(tstop/bin_width)])
                
            
            result = minimize(residual, params, args=(r_mod, spike_freq, bin_width, tstart, tstop))
            
            print "chisqr: ", result.chisqr
            print 'Best-Fit Values:'
            for name, par in params.items():
                print '  %s = %.4f +/- %.4f ' % (name, par.value, par.stderr)
            
            amp = params['amp'].value
            shift = params['shift'].value
            tau1 = params['tau1'].value
            tau2 = params['tau2'].value
                
            model = modelfun(amp, shift, tau1, tau2, bin_width = bin_width, r_mod = r_mod)               
            
            
            if "ex" in stimtype:
                plt.figure(0)
                plt.plot(freq_times[int(0.5/bin_width):int(trun/bin_width)], spike_freq[int(0.5/bin_width):int(trun/bin_width)], freq_times[int(0.5/bin_width):int(trun/bin_width)], model[int(0.5/bin_width):int(trun/bin_width)])
                plt.figure(1)
                plt.plot(time, voltage[0]), freq_times, r_mod, time, current
                #plt.figure(100)                
                #plt.plot(t_all_vec_vec[0],id_all_vec_vec[0],'k|')
                #plt.savefig("./figs/dump/taufit_" + str(stimtype) + "_spikes.pdf", dpi = 300)  # save it 
            
            else:
                plt.figure(0)
                plt.plot(freq_times[int(0.5/bin_width):int(trun/bin_width)], spike_freq1[int(0.5/bin_width):int(trun/bin_width)], freq_times[int(0.5/bin_width):int(trun/bin_width)], spike_freq[int(0.5/bin_width):int(trun/bin_width)], freq_times[int(0.5/bin_width):int(trun/bin_width)], model[int(0.5/bin_width):int(trun/bin_width)])
                plt.figure(1)
                plt.plot(time, voltage[0], time, voltage[1], freq_times, r_mod, time, current)
                plt.figure(100)                
                #plt.plot(t_all_vec_vec[0],id_all_vec_vec[0],'k|')
                #plt.plot(t_all_vec_vec[1],id_all_vec_vec[1],'b|')
                #plt.savefig("./figs/dump/taufit_" + str(stimtype) + "_spikes.pdf", dpi = 300)  # save it 
                
            
            plt.figure(0)
            plt.title('Fit: ' + str(stimtype) + ', tau1=' + str(tau1) +  ' tau2=' + str(tau2))
            plt.savefig("./figs/dump/taufit_" + str(stimtype) + "_rate.png", dpi = 300)  # save it  
            
            plt.figure(1)
            plt.savefig("./figs/dump/taufit_" + str(stimtype) + "_voltage.png", dpi = 300)  # save it 
            
            
    plt.show()