Cycle skipping in ING Type 1 / Type 2 networks (Tikidji-Hamburyan & Canavier 2020)

 Download zip file 
Help downloading and running models
Accession:259366
"All-to-all homogeneous networks of inhibitory neurons synchronize completely under the right conditions; however, many modeling studies have shown that biological levels of heterogeneity disrupt synchrony. Our fundamental scientific question is “how can neurons maintain partial synchrony in the presence of heterogeneity and noise?” A particular subset of strongly interconnected interneurons, the PV+ fast spiking basket neurons, are strongly implicated in gamma oscillations and in phase locking of nested gamma oscillations to theta. Their excitability type apparently varies between brain regions: in CA1 and the dentate gyrus they have type 1 excitability, meaning that they can fire arbitrarily slowly, whereas in the striatum and cortex they have type 2 excitability, meaning that there is a frequency threshold below which they cannot sustain repetitive firing. We constrained the models to study the effect of excitability type (more precisely bifurcation type) in isolation from all other factors. We use sparsely connected, heterogeneous, noisy networks with synaptic delays to show that synchronization properties, namely the resistance to suppression and the strength of theta phase to gamma amplitude coupling, are strongly dependent on the pairing of excitability type with the type of inhibition. ..."
Reference:
1 . Tikidji-Hamburyan RA, Canavier CC (2020) Shunting Inhibition Improves Synchronization in Heterogeneous Inhibitory Interneuronal Networks with Type 1 Excitability Whereas Hyperpolarizing Inhibition is Better for Type 2 Excitability. eNeuro [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type:
Brain Region(s)/Organism:
Cell Type(s): Abstract single compartment conductance based cell;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON;
Model Concept(s):
Implementer(s): Tikidji-Hamburyan, Ruben [ruben.tikidji.hamburyan at gmail.com] ;
"""
/***********************************************************************************************************\

 This a NEURON + Python 2 script associated with the paper:

 Shunting Inhibition Improves Synchronization in Heterogeneous Inhibitory Interneuronal Networks with Type 1 Excitability Whereas Hyperpolarizing Inhibition is Better for Type 2 Excitability

 eNeuro (in press)

 Copyright: Ruben Tikidji-Hamburyan <rath@gwu.edu> <rtikid@lsuhsc.edu> Aug.2018 - ....

\************************************************************************************************************/
"""
import numpy as np
from scipy import optimize
from scipy import integrate
import scipy.stats as sps
import sys,os,csv,threading
import random as rnd
try:
	import cPickle as pickle
except:
	import pickle
from datetime import datetime
import time
import fcntl

from neuron import h
h.load_file("stdgui.hoc")




###### Abbreviations:
Abbreviations=(
	( 'I',       'I',   "Current"),
	( 'TI',      'TI',  "Total Current"),
	( 'TSI',     'TSI', "Total Synaptic Current"),
	( 'MTI',     'MTI', "Mean Total Curent"),
	( 'MTSI',    'MTSI',"Mean Total Synaptic Current"),
	( 'G',       'G',   "Condunctance"),
	( 'TG',      'TG',  "Total Conductance"),
	( 'MTG',     'MTG', "Mean Total Conducntance"),
	( 'FR',      'FR',  "Firing Rate"),
	( 'NORM',   'NORM',"Normal distribution"),
	( 'LOGN',   'LOGN',"Log Normal Distribution"),
	( 'Mstate', 'm',   "Sodium Activation variable"), 
	( 'Hstate', 'h',   "Sodium Inctivation variable"),
	( 'Nstate', 'n',   "Potassium Activation variable"),
	( 'ucon',   'ucon',"Uniform distribution of number connection per cell"),
	( 'ncon',   'ncon',"Normal distribution of number connection per cell"),
	( 'bcon',   'bcon',"Binomial distribution of number connection per cell"),
	( 'ON',     'True',"Turn ON some parameter"),
	( 'OFF',     'False',"Turn OFF some parameter"),
)

for ab,val,meaning in Abbreviations:	
	print "Applay abbreviation % 8s for % 8s for: "%(ab,val)+meaning,
	try:
		exec ab+'=\''+val+'\''
	except:
		print "Fail!"
		exit(1)
	print "DONE"


###### Paramters:
methods		= {
	'ncell'		: 300,			# number of neurons in population
	'ncon'		: ('b',0.133),	# number of input connections per neuron
								# constant or uniform distribution(from, to) or normalized uniform distribution (mena, stder, ncon-norm)
								# OR uniform distribution  ('u', from, {to},    {{ncon-norm}})
								# OR normal distribution   ('n', mean, {stdev}, {{ncon-norm}})
								# OR binomial distribution ('b', prob,           {ncon-norm} )
	'neuron'	: {
		'Vinit'		: (-50.,20),#(-51.86007190636312,20),	# Constant or (mean,stdev) or [value for each neuron] or string or file name there values for each neuron are contained.
		'Type'      : 1,
		'Iapp'      : None,
		'Istdev'	: None,						# ---/---/---

	},
	'synapse'	: {
		#'weight'	: 0.75e-2,					# Synaptic conductance
		'weight'	: 0.03e-2,					# Synaptic conductance
		'delay'		: 0.8,						# Axonal Delays
		'gsynscale'	: 1.0,						# Conductance caramelization
		'tau1'		: 1.0,
		'tau2'		: 3.0,
		'Esyn'		: -75.0,					# Synaptic reversal potential
												# Constant or (mean,stdev) or [value for each neuron] or string or file name there values for each neuron are contained.
		'synscaler'	: None,

	},		
	'R2'		: True,
	'maxFreq'	: 200.0,		# max frequency
	'peakDetec' : True,			# Turn on/off peak detector
	'gkernel'	: (3.,25.),#(5.0,25.0),#(10.0,50.0),	# Kernel and size (5,25),#
	"netFFT"	: False,#True,#False,#True,#False,		# Turn on/off network FFT
	"nrnFFT"	: False,#True,#False,		# Turn on/off neuron FFT
	'netISI'	: 30001,			# max net ISI
	'nrnISI'	: 30001,			# max neuron ISI
	'cliptrn'	: False,#1000,#False alse,#500,#False,	# Clip transience for first n ms or False
	'traceView'	: False,#'n',
	'tV-synmax' : False,
	'traceWidth': 55.0,			#
	'tracetail'	: 'mean total conductance',#'conductance',#'mean total conductance',#'conductance',#'current',#'conductance',#'firing rate',#'total conductance', #'total current' 'total current',
	'patview'	: True,			# Turn on/off Pattern view
	'gui'		: True,
	'git'		: False,		# Turn on/off git core (Never turn-on at the head node!!!!!!!)
	'gif'		: False,		# Generate gif instead pop up on a screan.
	'corefunc'	: (4,8,64),
	'coreindex'	: False,			# Turn on/off Core indexing
	'corelog'	: 'network',
	'noexit'	: False,
	'GPcurve'	: False,
	'IGcurve'	: False,
	'Conn-rec'	: False,
	'Conn-stat'	: False,
	'G-Spikerate'	: False,
	
	'Gtot-dist'	: False,
	'Gtot-rec'	: False,#True,	#record all gtotal in neurons
	'Gtot-stat'	: True, #False, #record gtot statistic
	'sycleprop'	: False,#True,
	'external'	: False,
	'extprop'	: 0.5,				# Calculate probability to fire after external input
	'timestep'	: 0.025,#0.005,
	'sortbysk'	: False, #'ST',#'F',#'GT',#'NC',#'GT','G',#'I', #'F',#'I',#'F',#False,			#Do not use
	'taunorm'	: False,#True,#False,#True,
	'nstart'	: False, #(900.,0.1e-5,1000),<Noise for paper  #False,#(900.,0.2e-5,1000),#False,#(200,0.000002,900),#False, (delay, ampl, dur)
	'cliprst'	: False,#10,#False,#20,
	'T&S'		: False,#True,
	'lastspktrg': True,
	'fullrast'	: True,
	'gtot-dist'	: 'LOGN', #LOGN - lognormal, 'NORM' - normal
	'gsyn-dist' : 'LOGN',#'NORM',#'LOGN', #same
	'cycling'	: False, #4,False
	'popfr'		: False,#True,	#calculate population firing rate
	'cmd-file'	: 'network.start',
	'preview'	: True,
	'2cintercon': False,#True,
	'2clrs-stat': False,#True,
	'tv'		: (0., 500.),
	'tstop'		: 10001,
	'jitter-rec': False,#True,
	'pop-pp-view':False,#True,
	'N2NHI'     : False, #True,		#neuron to network harmonic index
	'N2NHI-netISI' :False,	#the same index but using netowkr ISI to get network harmonics (it is slow)
	'vpop-view' : False,
	'CtrlISI'   : False,#{'bin'   : 5.,'max'   : 120.,},# ISI histogram with controled bin and max
	'nrnFRhist' : False, #Neuron Firing rate histogram
}



class neuron:
	def __init__(self):
		self.soma = h.Section()
		if checkinmethods('/neuron/L'):
			self.soma.L		= methods["neuron"]["L"]
		else:
			self.soma.L		= 100.
		if checkinmethods('/neuron/diam'):
			self.soma.diam	= methods["neuron"]["diam"]
		else:
			self.soma.diam	= 10./np.pi
		if checkinmethods('/neuron/nseg'):
			self.soma.nseg	= int(methods["neuron"]["nseg"])
		else:
			self.soma.nseg	= 1
		if checkinmethods('/neuron/cm'):
			self.soma.cm	= float(methods["neuron"]["cm"])
		self.soma.insert('type21')
		self.soma(0.5).type21.type21 = 1 #default type
		self.type21 = 1
		if checkinmethods('/neuron/set'):
#			print "=================================="
#			print "===      SETTING PARAMETERS    ==="
#			print "  >  Cells in the Cluster A      :",countA
			for p in methods["neuron"]["set"]:
				exec "self.soma(0.5)."+p+" = {}".format(methods["neuron"]["set"][p])
			
		self.soma(0.5).v = -67.
		self.isyn	= h.Exp2Syn(0.5, sec=self.soma)
		self.isyn.e		= -75.0
		self.isyn.tau1	= 2.0
		self.isyn.tau2	= 10.0
		######## Recorders ##########
		self.spks	= h.Vector()
		self.sptr	= h.APCount(.5, sec=self.soma)
		self.sptr.thresh = 0.#25
		self.sptr.record(self.spks)
		#self.sptr = h.NetCon(self.soma(0.5)._ref_v,None,sec=self.soma)
		#self.sptr.threshold = 25.
		#self.sptr.record(self.spks)
		if checkinmethods('gui'):
			self.volt	= h.Vector()
			self.volt.record(self.soma(0.5)._ref_v)
			if checkinmethods("neuron/record/current"):
				self.isyni	= h.Vector()
				self.isyni.record(self.isyn._ref_i)
			if checkinmethods('neuron/record/conductance') or checkinmethods('pop-pp-view'):
				self.isyng	= h.Vector()
				self.isyng.record(self.isyn._ref_g)
			if checkinmethods('traceView') or checkinmethods('pop-pp-view'):
				self.svar   = h.Vector()
				self.svar.record(self.soma(0.5)._ref_n_type21)
		elif checkinmethods('get-steadystate'):
			self.volt	= h.Vector()
			self.volt.record(self.soma(0.5)._ref_v)
		if checkinmethods('sinmod'):
			self.sin = h.sinIstim(0.5, sec=self.soma)
			if type(methods['sinmod']) is dict:
				for name in methods['sinmod']:
					exec "self.sin."+name+"= {}".format(methods['sinmod'][name])
		if checkinmethods('singmod'):
			self.sing = h.sinGstim(0.5, sec=self.soma)
			if type(methods['singmod']) is dict:
				for name in methods['singmod']:
					exec "self.sing."+name+"= {}".format(methods['singmod'][name])

		######## Registrations ###### 
		self.gsynscale	= 0.0
		self.concnt		= 0.0
		self.gtotal		= 0.0
		self.tsynscale	= 1.0

	def setparams(self, 
			V=None, N=-1, Type=1,
			Iapp = None, Insd = None, delay = None, duration = None, period = None,
			SynE = None, SynT1 = None, SynT2 = None,
			):
		if not    V is None : self.soma(0.5).v             = V
		if not    N is None : self.soma(0.5).type21.ninit  = N
		if not Type is None : self.soma(0.5).type21.type21 = self.type21 = Type
		########
		if not Iapp is None or not Insd is None:
			self.innp	= h.InNp(0.5, sec=self.soma)
			self.rnd	= h.Random(np.random.randint(0,32562))
			if checkinmethods("/neuron/rnd") and ( methods["neuron"]["rnd"] == "u"  or methods["neuron"]["rnd"] == "U" ):
				self.rnd.uniform(0.,1.)
			self.innp.noiseFromRandom(self.rnd)
			self.innp.dur	= 1e9 if duration == None else duration
			self.innp.delay	= 0   if delay == None else delay
			self.innp.per	= 0.1 if period == None else period
			self.innp.mean	= 0.0 if Iapp == None else Iapp
			self.innp.stdev	= 0.0 if Insd == None else Insd
			if checkinmethods("/neuron/rnd") and ( methods["neuron"]["rnd"] == "u"  or methods["neuron"]["rnd"] == "U" ):
				self.innp.stdev = -self.innp.stdev - self.innp.mean
			if methods['gui']:
				self.inoise	= h.Vector()
				self.inoise.record(self.innp._ref_i)
			elif checkinmethods("rawdata") and type(methods["rawdata"]) is str:
				self.inoise	= h.Vector()
				self.inoise.record(self.innp._ref_i)
		########
		if not SynE  is None: self.isyn.e		= SynE
		if not SynT1 is None: self.isyn.tau1	= SynT1
		if not SynT2 is None: self.isyn.tau2	= SynT2
		########
			
			

	def addnoise(self,Iapp=0.,Insd=0.,delay=0.,dur=0.,per=0.1):
		self.andnoise = h.InNp(0.5, sec=self.soma)
		self.andrnd	= h.Random(np.random.randint(0,32562))
		self.andnoise.noiseFromRandom(self.andrnd)
		self.andnoise.mean  = Iapp
		self.andnoise.stdev = Insd
		self.andnoise.delay	= delay
		self.andnoise.per	= per
		self.andnoise.dur	= dur
		self.iandnoise	    = h.Vector()
		self.iandnoise.record(self.andnoise._ref_i)

#class symulation:
	#def __init___(self,params):
		#if params.get("a",False):

def getType21(nid):
	for i in 'type21,gna,ena,gk,ek,gl,el,n0,sn,t0,st,v12,sl,a,b'.split(","):
		exec "{0} = neurons[{1}].soma(0.5).type21.{0}".format(i,nid)
	s  = neurons[nid].soma.L * neurons[0].soma.diam * np.pi * 1e-8 # cm2
	es = neurons[nid].isyn.e
	return type21,gna,ena,gk,ek,gl,el,n0,sn,t0,st,v12,sl,a,b,s,es


def getnulls(nid,vmin,vmax,gsyn,inoise,ibias,gstim = 0.,estim = 0.):
	type21,gna,ena,gk,ek,gl,el,n0,sn,t0,st,v12,sl,a,b,s,es = getType21(nid)
	#DB>>
	if checkinmethods("local-parameters"):
		print "======== THE SET ======="
		for i in 'type21,gna,ena,gk,ek,gl,el,n0,sn,t0,st,v12,sl,a,b,s,gsyn,es,inoise,ibias'.split(","):
			print " > % 6s          = "%i,eval(i)
	#<<DB
	
	## N-nullcline
	vx=np.linspace(vmin,vmax,200)
	ninf = lambda v:n0 + sn/(1.+np.exp(-(v-v12)/sl ))
	ntau = lambda v:t0 + st/(1.+np.exp( (v+40.)/20.))
	n0c  = np.dstack( (vx, ninf(vx)) )[0]

	## V - null	
	dvdt = lambda n,v,I:I+gl*(el-v)+gna*(1./(1.+np.exp(-(v+40.)/9.5)))**3*(b*n+a)*(ena-v)+gk*n**4*(ek-v)	
	def vfun(vx,I):
		"""
		Solves 4th order algebraic equation and returns null-cline
		"""
		if type(I) is float or type(I) is int:
			I = np.ones(vx.shape[0])*I 
		res = []
		for vp,ip in zip(vx,I):
			try:
				n=optimize.fsolve(dvdt,0.5,args=(vp,ip),xtol=0.01)[0]	
			except: return np.array([[],[]])
			res.append((vp,n))
		return np.array(res)	

	v0c  = vfun(vx,( -ibias  - gstim*(vx-estim)                )*1e-3 /s )
	v0n  = vfun(vx,( -inoise - gstim*(vx-estim) - gsyn*(vx- es))*1e-3 /s )
	if type21 == 2:
		vXn = v0c[ np.where(v0c[:,0] > -40.) ]
		try:
			d2vdt2 = np.polyder(np.polyfit(vXn[:,0], vXn[:,1], 3),m=2)
			vnX = vXn[ np.argmax(d2vdt2) ]
		except:
			vnX = vXn[ np.argmax(vXn[:,1]) ]
		#DB>>
		#print " > % 9s       = "%("Vinit"),vnX
		#<<DB
		dvdt = lambda n,v:-ibias*1e-3/s+gl*(el-v)+gna*(1./(1.+np.exp(-(v+40.)/9.5)))**3*(b*n+a)*(ena-v)+gk*n**4*(ek-v)
		def rhs(t,Y): return [-dvdt(Y[1],Y[0]),-(ninf(Y[0])-Y[1])/ntau(Y[0])]
		slv = integrate.ode(rhs).set_initial_value(vnX, 0)#.set_integrator('zvode', method='bdf')
		thc = [ slv.integrate(slv.t+0.01) ]
		vbl,vbr = vx[0],vnX[0]+0.1
		while slv.successful() and slv.t < 100. and vbl < thc[-1][0] < vbr and 0. < thc[-1][1] < 1:
			thc.append( slv.integrate(slv.t+0.01) )
		thc = np.array(thc)
		#========#
		vXn = v0n[ np.where(v0n[:,0] > -40.) ]
		try:
			d2vdt2 = np.polyder(np.polyfit(vXn[:,0], vXn[:,1], 3),m=2)
			vnX = vXn[ np.argmax(d2vdt2) ]
		except:
			vnX = vXn[ np.argmax(vXn[:,1]) ]
		#DB>>
		#print " > % 9s       = "%("Vinit"),vnX
		#<<DB
		dvdt = lambda n,v:(-inoise-gsyn*(v-es))*1e-3/s+gl*(el-v)+gna*(1./(1.+np.exp(-(v+40.)/9.5)))**3*(b*n+a)*(ena-v)+gk*n**4*(ek-v)
		def rhs(t,Y):return [-dvdt(Y[1],Y[0]),-(ninf(Y[0])-Y[1])/ntau(Y[0])]
		slv = integrate.ode(rhs).set_initial_value(vnX, 0)#.set_integrator('zvode', method='bdf')
		thn = [ slv.integrate(slv.t+0.01) ]
		vbl,vbr = vx[0],vnX[0]+0.1
		while slv.successful() and slv.t < 100. and vbl < thn[-1][0] < vbr and 0. < thn[-1][1] < 1:
			thn.append( slv.integrate(slv.t+0.01) )
		thn = np.array(thn)
		
	else:
		thc = None
		thn = None

	return n0c,v0c,v0n,thc,thn,type21

def onclick1(event):
	if not hasattr(onclick1,"aix"):
		onclick1.aix=zooly.add_subplot(111)
	onclick1.et = event.xdata
	
	### BUG
	onclick1.tl, onclick1.tr = onclick1.et-methods['traceWidth'], onclick1.et+methods['traceWidth']
	onclick1.idx, = np.where( (t > onclick1.tl) * (t < onclick1.tr))
	
	if not hasattr(onclick1,"marks"):
		onclick1.marks = []
		onclick1.marks.append( p.plot([onclick1.tl,onclick1.tl],[-80,30],"r--",lw=2)[0] )
		onclick1.marks.append( p.plot([onclick1.tr,onclick1.tr],[-80,30],"r--",lw=2)[0] )
	else:
		onclick1.marks[0].set_xdata([onclick1.tl,onclick1.tl])
		onclick1.marks[1].set_xdata([onclick1.tr,onclick1.tr])
	
	if not hasattr(onclick1,"lines"):
		onclick1.lines = []
		for n in neurons:
			volt = np.array(n.volt)
			onclick1.lines.append(onclick1.aix.plot(t[onclick1.idx],volt[onclick1.idx])[0])
	else:
		vmin,vmax = 1000,-1000
#		for ind,n in map(None,xrange(methods["ncell"]),neurons):
		for ind,n in enumerate(neurons):
			volt = np.array(n.volt)
			if vmin > volt[onclick1.idx].min():vmin = volt[onclick1.idx].min()
			if vmax < volt[onclick1.idx].max():vmax = volt[onclick1.idx].max()
			onclick1.lines[ind].set_xdata(t[onclick1.idx])
			onclick1.lines[ind].set_ydata(volt[onclick1.idx])
			onclick1.lines[ind].set_linewidth(1)
			onclick1.lines[ind].set_ls("-")

		onclick1.aix.set_xlim(onclick1.tl,onclick1.tr)
		#print vmin,"---",vmax
		onclick1.aix.set_ylim(vmin,vmax)
	if hasattr(moddyupdate,"lines"):
		del moddyupdate.lines
	mainfig.canvas.draw()
	zoolyupdate(vindex)
	moddyupdate(moddyupdate.idx)



def neuronsoverview(event):
	global vindex
	if event.key == "up":
		vindex += 1
		if vindex >= methods["ncell"] : vindex = methods["ncell"] -1
	elif event.key == "down":
		vindex -= 1
		if vindex < 0 : vindex = 0
	elif event.key == "home": vindex = 0
	elif event.key == "end" : vindex = methods["ncell"] -1
	elif event.key == "/"   : vindex = methods["ncell"]/2
	elif event.key == "pageup":
		vindex += 10
		if vindex >= methods["ncell"] : vindex = methods["ncell"] -1
	elif event.key == "pagedown":
		vindex -= 10
		if vindex < 0 : vindex = 0
	nsorted = methods['sortbysk'] == 'I'  or methods['sortbysk'] == 'G' or methods['sortbysk'] == 'NC' or\
	          methods['sortbysk'] == 'GT' or methods['sortbysk'] == 'ST'or methods['sortbysk'] == 'N'  or\
	          methods['sortbysk'] == 'T'  or methods['sortbysk'] == 'FR'
	if nsorted:
		  print vindex,"->",nindex[vindex][1],"("+methods['sortbysk']+")=",nindex[vindex][0]
		  vtrace.set_ydata( np.array(neurons[nindex[vindex][1]].volt)[tproc:tproc+tprin.size])
	else:
		vtrace.set_ydata( np.array(neurons[vindex].volt)[tproc:tproc+tprin.size])
	
	if methods['tracetail'] == 'conductance':
		if nsorted :
			xvcrv.set_ydata( np.array(neurons[nindex[vindex][1]].isyng)[tproc:tproc+tprin.size]*1e5 )
		else:
			xvcrv.set_ydata( np.array(neurons[vindex].isyng)[tproc:tproc+tprin.size]*1e5 )
	elif methods['tracetail'] == 'current':
		if nsorted :
			xvcrv.set_ydata( np.array(neurons[nindex[vindex][1]].isyni)[tproc:+tprin.size]*1e5 )
		else:
			xvcrv.set_ydata( np.array(neurons[vindex].isyni)[tproc:+tprin.size]*1e5 )
	if event.key == "H":
		if nsorted:
			p.plot(tprin,np.array(neurons[nindex[vindex][1]].volt)[tproc:tprin.size+tproc],"-")
		else:
			p.plot(tprin,np.array(neurons[vindex].volt)[tproc:tprin.size+tproc],"-")
	mainfig.canvas.draw()
	if checkinmethods('traceView'):
		if nsorted:
			zoolyupdate(nindex[vindex][1])
		else: 
			zoolyupdate(vindex)
		moddyupdate(moddyupdate.idx)

def positiveGauss(mean,stdev):
	result = -1
	while result < 0:
		result = mean + np.random.randn()*stdev
	return result

def checkinmethods(item, dirtree = methods):
	def getsubitems(item):
		items = item.split("/")
		if items[ 0] == "" and len(items) !=1 : items = items[1:]
		if items[-1] == "" and len(items) !=1 : items = items[:-1]
		return items[0],"/".join(items[1:])
	item,subitems = getsubitems(item)
	if subitems != "":
		if not item in dirtree : return False
		if not type(dirtree[item]) is dict: return False
		return checkinmethods(subitems,dirtree[item])
	else:
		if not item in dirtree : return False
		if not ( (type(dirtree[item]) is bool or type(dirtree[item]) is int) ):
			if dirtree[item] is None: return False
			else: return True
		return bool(dirtree[item])

def ggap_var(t,t0,t1,r0,r1):
	if t < t0:
		for gj in gapjuctions: gj[0].r, gj[1].r = r0, r0
	elif t > t1:
		for gj in gapjuctions: gj[0].r, gj[1].r = r1, r1
	else :
		r = (r1-r0)*(t-t0)/(t1-t0)+r0
		for gj in gapjuctions: gj[0].r, gj[1].r = r, r
	#DB>>
#	print "ggap_var was called with parameters", t,t0,t1,r0,r1
#	exit(0)
	#<<DB
	
def getNdist(prm, nrntype=None):
	if type(prm) is float or type(prm) is int:
		return prm*np.ones(methods["ncell"])
	elif type(prm) is tuple:
		if len(prm) > 1:
			if type(prm[0]) is not str:
				return prm[0]+np.random.randn(methods["ncell"])*prm[1]
			else:
				if   prm[0][0] == "n" or prm[0][0] == "N":
					if   len(prm) == 3:return prm[1]+np.random.randn(methods["ncell"])*prm[2]
					elif len(prm) == 2:return prm[1]
					else:
						print "ERROR: normal distribution should have mean and std parameters ('n',mean,std)!"
						exit(1)
				elif prm[0][0] == "u" or prm[0][0] == "U":
					if   len(prm) == 3:return prm[1]+np.random.rand(methods["ncell"])*(prm[2]-prm[1])
					else:
						print "ERROR: normal distribution should have two parameters left and right borders ('u',min,max)!"
						exit(1)
				#elif prm[0][0] == "l" or prm[0][0] == "L":
				#elif prm[0][0] == "s" or prm[0][0] == "S": #shifted lognormal
				#elif prm[0][0] == "t" or prm[0][0] == "T": #Truncated normal
					#if   len(prm) == 4:
						#return prm[1]+np.random.rand(methods["ncell"])*(prm[2]-prm[1])
				elif prm[0][0] == "q" or prm[0][0] == "Q": #prameter depends upon the neuron type!
					if   len(prm) != 3:
						print "ERROR: type dependnet paramter should be  ('q',param_for_type1,param_for_type2)!"
						exit(1)
					return [ prm[2] if t == 2 else prm[1] for t in nrntype ]
				else:
					print "ERROR: unknown distribution type for parameter {}!".format(prm)
					exit(1)
		else:
			return prm[0]*np.ones(methods["ncell"])
	elif type(prm) is list:
		if len(prm) == methods['ncell'] :
			return np.array(prm)
		else: 
			return [ None for i in xrange(methods["ncell"]) ]
	elif type(prm) is str:
		return np.genfromtxt(prm)
		print "  > Read Vinit from file         :",prm
	elif prm is None:
		return[ None for i in xrange(methods["ncell"]) ]

#elif methods["delay-dist"] == "NORM":
						#### Trancated normal
						#dx = positiveGauss(methods['synapse']['delay'][0],methods['synapse']['delay'][1])
						#if len(methods['synapse']['delay']) < 3:
							#while(dx < methods['timestep']*2):
								#dx = positiveGauss(methods['synapse']['delay'][0],methods['synapse']['delay'][1])
						#else:
							#while(dx < methods['synapse']['delay'][2]):
								#dx = positiveGauss(methods['synapse']['delay'][0],methods['synapse']['delay'][1])
							
					#elif methods["delay-dist"] == "LOGN":
						#### Lognormal
						#dlym,dlys=methods['synapse']['delay'][0]-methods['timestep']*2.,methods['synapse']['delay'][1]
						#if len(methods['synapse']['delay']) < 3:
							#dx = np.random.lognormal(mean=np.log(dlym/np.sqrt(1.+dlys**2/dlym**2)),sigma=np.sqrt(np.log(1.+dlys**2/dlym**2)))+methods['timestep']*2
						#else:
							#dx = np.random.lognormal(mean=np.log(dlym/np.sqrt(1.+dlys**2/dlym**2)),sigma=np.sqrt(np.log(1.+dlys**2/dlym**2)))
							#while dx < methods['synapse']['delay'][2]:
								#dx = np.random.lognormal(mean=np.log(dlym/np.sqrt(1.+dlys**2/dlym**2)),sigma=np.sqrt(np.log(1.+dlys**2/dlym**2)))
					#elif methods["delay-dist"] == "LOGN-SHIFT":
						#### Lognormal
						#dlym,dlys=methods['synapse']['delay'][0]-methods['timestep']*2.,methods['synapse']['delay'][1]
						#if len(methods['synapse']['delay']) < 3:
							#dx = np.random.lognormal(mean=np.log(dlym/np.sqrt(1.+dlys**2/dlym**2)),sigma=np.sqrt(np.log(1.+dlys**2/dlym**2)))+methods['timestep']*2
						#else:
							#dx = np.random.lognormal(mean=np.log(dlym/np.sqrt(1.+dlys**2/dlym**2)),sigma=np.sqrt(np.log(1.+dlys**2/dlym**2)))+methods['synapse']['delay'][2]
					#elif methods["delay-dist"] == "DIST":
						#dmin = methods['synapse']['delay'][0]
						#dinciment = methods['synapse']['delay'][1] if len(methods['synapse']['delay']) >= 2 else 0.
						#dx = dmin + dinciment*float(abs(pre-x[0]))		
			
	
#===============================================#
#               MAIN PROGRAMM                   #
#===============================================#
if __name__ == "__main__":
	if len(sys.argv) > 1:
		def setmethod(arg):
			global methods
			if not "=" in arg: return
			try:
				name,value = arg.split("=")
			except: 
				print "ERROR! Parameter %s has not = symbol"%arg
				exit(1)
			if not "/" in name: return
			if name[0] != '/' : return
			names = name.split("/")
			if len(names) > 2:
				name = "methods"
				for nm in names[1:-1]:
					exec "inmethods=\'"+nm+"\' in "+name
					if inmethods :
						exec "inmethods=type("+name+"[\'"+nm+"\']) is dict"
						if inmethods :
							name += "[\'"+nm+"\']"
							continue
						else:
							exec "inmethods=type("+name+"[\'"+nm+"\']) is bool or type("+name+"[\'"+nm+"\']) is None"
							if inmethods :
								name += "[\'"+nm+"\']"
								exec name+"={}"
							else:
								sys.stdout.write("method item %s of %s isn't dict\n"%(name,nm))
					else:
						name += "[\'"+nm+"\']"
						exec name+"={}"
					
			cmd = "methods" + reduce(lambda x,y: x+"[\'"+y+"\']", names[1:],"") + "=" + value
			try:
				exec cmd
			except: 
				#cmd = "methods" + reduce(lambda x,y: x+"[\'"+y+"\']", names[1:],"") + "=" + "\'\\\'"+value+"\\\'\'"
				cmd = "methods" + reduce(lambda x,y: x+"[\'"+y+"\']", names[1:],"") + "=" + "\'"+value+"\'"
				exec cmd
		def readfromsimdb(simdb,ln):
			rec = None
			with open(simdb) as fd:
				for il,l in enumerate(fd.readlines()):
					if il == ln: rec = l
			if rec == None:
				sys.stderr.write("Cannot find line %d in the file %s\n"%(ln,simdb))
				exit(1)
			for itm in rec.split(":"):
				#DB>>
				print itm
				#<<DB
				setmethod(itm)
		simdbrec=None
		for arg in sys.argv:
			if arg[:len('--readsimdb=')] == '--readsimdb=':
				simdbrec=arg[len('--readsimdb='):]
			#if arg == '-h' or arg == '-help' or arg == '--h' or arg == '--help':
				#print __doc__
				#print 
				#print "USAGE: nrngui -nogui -python network.py [parameters]"
				#print "\nPARAMETERS:"
				#print "-n=          number of neurons in population"
				#print "-c=          number of connections per neuron"
				#print "-Iapp=       apply current. Use scaling factor 1e-5 to get nA"
				#print "             Iapp may be a constant or mean,standard deviation across population."
				#print "-Istd=       amplitude of noise. Should be scaled by 1e-5 to get nA"
				#print "-gui=ON|OFF  Turn on/off gui and graphs"
				#print "-F=          Set up neuron dynamics scale factor"
				#print "             F may be a constant or mean,standard deviation across population."
				#print "-gsyn=       conductance of single synapse. Use scaling factor 1e-5 to get nS"
				#print "             gsyn may be a constant or mean,standard deviation for all synapses in model"
				#print "-gsynscale=  total synaptic conductance.  Use scaling factor 1e-5 to get nS"
				#print "             gsynscale may be a constant or mean,standard deviation for all neurons within population"
				#print "-tau1=       rising time constant in ms"
				#print "-tau2=       falling time constant in ms"
				#print "-Esyn=       synaptic reversal potential in mV"
				#print "-taunorm=0|1 On or Off normalization by space under the curve"
				#print "-tsynscaler= scaling coefficient for synaptic time constants"
				#print "-delay=      axonal delay in ms"
				#print "             delay may be a constant or mean,standard deviation for all synapses in model"
				#print "-view        limits simulation and save memory"
				
				#exit(0)
			
		if not simdbrec is None:
			simdbrec = simdbrec.split(":")
			if len(simdbrec) < 2:
				sys.stderr.write("Error format --readsimdb=file:record\n")
			readfromsimdb( simdbrec[0], int(simdbrec[1]) ) 
		for arg in sys.argv: setmethod(arg)
			

	if 'cmd-file' in methods:
		if not type(methods['cmd-file']) is str: methods['cmd-file'] = 'network.start'
	else:
		methods['cmd-file'] = 'network.start'
	with open(methods['cmd-file'],"w") as fd:
		for arg in sys.argv: fd.write("%s "%arg)
	if methods['taunorm']:
		from norm_translation import getscale
		nFactor = getscale(1.0,3.0,methods['synapse']['tau1'],methods['synapse']['tau2'])
		if type(methods["synapse"]["weight"]) == tuple or type(methods["synapse"]["weight"]) == list:
			methods["synapse"]["weight"]		= (methods["synapse"]["weight"][0]*nFactor, methods["synapse"]["weight"][1]*nFactor)
		else:
			methods["synapse"]["weight"] *= nFactor
	if checkinmethods('preview'):
		methods['tstop'] = methods['tv'][1]

###DB>
	print "=================================="
	print "==       ::  METHODS  ::        =="
	def dicprn(dic, space):
		for nidx,name in enumerate(sorted([ x for x in dic ])):
			if type(dic[name]) is dict:
				rep = "%s%s\\ %s "%(space,"v-" if nidx==0 else "|-", name)
				print rep
				dicprn(dic[name], space+"  ")
			else:
				rep = "%s%s> %s "%(space,"`-" if nidx==len(dic)-1 else "|-", name)
				if len(rep) < 31:
					for x in xrange(31-len(rep)):rep += " "
				if type(dic[name]) is str:
					print rep," : ","\'%s\'"%dic[name]
				else:
					print rep," : ",dic[name]
	dicprn(methods,' ')
	print "==================================\n"
###<DB

	
	if methods["gui"]:
		import matplotlib
		if checkinmethods("/BackEnd"):
			matplotlib.use(methods["BackEnd"],warn=False, force=True)
		import matplotlib.pyplot as plt
		matplotlib.rcParams["savefig.directory"] = ""
		#cmap = matplotlib.cm.get_cmap('jet')
		#cmap = matplotlib.cm.get_cmap('plasma')
		#cmap = matplotlib.cm.get_cmap('autumn')
		#cmap = matplotlib.cm.get_cmap('gist_rainbow')
		cmap = matplotlib.cm.get_cmap('rainbow')
		print "=================================="
		print "===        GUI turned ON       ==="
		print "==================================\n"
	
	h.tstop 	= float(methods['tstop'])
	#h.v_init 	= V
	h.dt		= float(methods["timestep"])
	if checkinmethods('temperature'):
		temp = float(h.celsius)
		h.celsius = methods['temperature']
		print "=================================="
		print "===       SET TEMPERATURE      ==="
		print " > set temperature"
		print " \> from {} to {} celsius degree  ".format(temp,methods['temperature'])
		print "==================================\n"
		

	if checkinmethods('simvar'):
		if not type(methods['simvar']) is dict: methods['simvar'] = False
		elif not "type" in methods['simvar']: methods['simvar'] = False
		elif not type(methods['simvar']["type"]) is str: methods['simvar'] = False
		elif not (methods['simvar']["type"] == 'n' or methods['simvar']["type"] == 'c' or methods['simvar']["type"] == 'g'  ): methods['simvar'] = False
		elif not "var" in methods['simvar']: methods['simvar'] = False
		elif not type(methods['simvar']["var"]) is str: methods['simvar'] = False
		elif not "a0" in methods['simvar']: methods['simvar'] = False
		elif not (type(methods['simvar']["a0"]) is float or type(methods['simvar']["a0"]) is int): methods['simvar'] = False
		elif not "a1" in methods['simvar']: methods['simvar'] = False
		elif not (type(methods['simvar']["a1"]) is float or type(methods['simvar']["a1"]) is int): methods['simvar'] = False
		elif not "t0" in methods['simvar']: methods['simvar']['t0'] = methods['tv'][0]
		elif not (type(methods['simvar']["t0"]) is float or type(methods['simvar']["t0"]) is int): methods['simvar'] = False
		elif not "t1" in methods['simvar']: methods['simvar']['t1'] = methods['tv'][1]
		elif not (type(methods['simvar']["t1"]) is float or type(methods['simvar']["t1"]) is int): methods['simvar'] = False
	if methods['tracetail'] == 'R2':
		if not checkinmethods("cont-R2"):
			print "Need /cont-R2= parameter with period of R2 sliding window"
			exit(1)	

	#### Save mamory, record only what is needed
	if checkinmethods('gui'):
		print "=================================="
		print "===          RECORDER          ==="
		methods['neuron']['record'] = {}
		if methods['tracetail'] == 'total conductance' or methods['tracetail'] == 'mean total conductance' or\
		   methods['tracetail'] == 'TG'                or methods['tracetail'] == 'MTG'                    or\
		   methods['tracetail'] == 'conductance'       or\
		   checkinmethods('traceView'):
			   methods['neuron']['record']['conductance'] = True
			   print " > RECORD                       : cunductance"
		
		if methods['tracetail'] == 'total current'               or methods['tracetail'] == 'TI'   or\
		   methods['tracetail'] == 'mean total current'          or methods['tracetail'] == 'MTI'  or\
		   methods['tracetail'] == 'total synaptic current'      or methods['tracetail'] == 'TSI'  or\
		   methods['tracetail'] == 'mean total synaptic current' or methods['tracetail'] == 'MTSI' or\
		   methods['tracetail'] == 'current'                     or\
		   checkinmethods('spectrogram'):
			   methods['neuron']['record']['current'] = True
			   print " > RECORD                       : current"
		if methods['tracetail'] == 'LFP'                         or methods['tracetail'] == 'p2eLFP' or\
		   checkinmethods('PAC-Canolty06') or checkinmethods('PAC-Tort10') or checkinmethods('PAC-VS'):
			methods[ "peakDetec" ] = True
		print "==================================\n"

	#### Check that PAC is used with sinmod
	if ( checkinmethods('PAC-Canolty06') or checkinmethods('PAC-Tort10') or checkinmethods('PAC-VS') ) and not ( checkinmethods('sinmod') or checkinmethods('singmod') ):
		print "Need /sinmod or /singmod to calulate PAC index(s)"
		exit(1)
	
	#### Create Neurons and setup noise and Iapp
	print "=================================="
	print "===        Create Neurons      ==="
	neurons = [ neuron() for x in xrange(methods["ncell"]) ]
	
	if   type(methods["neuron"]["Type"]) is int:
		xT = [ methods["neuron"]["Type"] for i in xrange(methods["ncell"]) ]
	elif type(methods["neuron"]["Type"]) is float:
		xT = [ 2 if i < methods["neuron"]["Type"] else 1 for i in np.random.random(methods["ncell"]) ]
	else:
		print "ERROR: unknown type of /neuron/Type parameter!"
		exit(1)
	

	xV      = getNdist( methods["neuron"]["Vinit"] )
	xEsyn   = getNdist( methods['synapse']['Esyn'] , nrntype=xT)
	xIapp   = getNdist( methods["neuron"]["Iapp"]  )
	xIstdev = getNdist( methods["neuron"]["Istdev"])
	#DB>>
	#for dbt,dbe in zip(xT,xEsyn): print dbt,dbe
	#exit(0)
	#<<DB

	
	for n,i in zip(neurons,xrange(methods["ncell"])):
		if not methods['synapse']['synscaler'] is None:
			if type(methods['synapse']['synscaler']) is float or type(methods['synapse']['synscaler']) is int:
				n.tsynscale = float(methods['synapse']['synscaler'])
				xTau1,xTau2 = methods['synapse']['tau1'] * n.tsynscale, methods['synapse']['tau2'] * n.tsynscale
			elif type(methods['synapse']['synscaler']) is list or type(methods['synapse']['synscaler']) is tuple:
				if len(methods['synapse']['synscaler']) >= 2:
					n.tsynscale = -1.0
					while( n.tsynscale < 0.0 ):
						n.tsynscale = methods['synapse']['synscaler'][0] + np.random.randn()*methods['synapse']['synscaler'][1]
					xTau1,xTau2 = methods['synapse']['tau1'] * n.tsynscale, methods['synapse']['tau2'] * n.tsynscale
				else :
					n.tsynscale = float(methods['synapse']['synscaler'][0])
					xTau1,xTau2 = methods['synapse']['tau1'] * n.tsynscale, methods['synapse']['tau2'] * n.tsynscale
			else:
				xTau1,xTau2 = methods['synapse']['tau1'],methods['synapse']['tau2']
		else:
			xTau1,xTau2 = methods['synapse']['tau1'],methods['synapse']['tau2']

		n.setparams(
			V=xV[i], Type = xT[i],
			SynT1=xTau1, SynT2=xTau2, SynE=xEsyn[i], 
			Iapp = xIapp[i] if xIapp[i] is None else -1.*xIapp[i], Insd=xIstdev[i]
			)
	print "==================================\n"
	
	if checkinmethods("neuron/distribution"):
		if type(methods["neuron"]["distribution"]) is not dict:
			print "=================================="
			print "===           ERROR            ==="
			print "=== neuron/distribution should ==="
			print "=== be a set of parameters and ==="
			print "=== parameters def.            ==="
			print "==================================\n"
			exit(1)
		#>> Init neurons
		h.finitialize()
		#>> do not init as a specific type
		for n in neurons: n.soma(0.5).type21.type21 = 0
		#>> set new distribution parameters
		for k,v in methods["neuron"]["distribution"].items():
			pl = getNdist( v , nrntype=xT)
			if k == "F":
				print "=================================="
				print "===   Setting Frequency scale  ==="
				print " > F                             :",v
				methods["neuron"]["F"] = pl
				for n,p in zip(neurons,pl):
					n.soma(0.5).type21.gna	*= p
					n.soma(0.5).type21.gk	*= p
					n.soma(0.5).type21.gl	*= p
					n.soma(0.5).type21.st	/= p
					n.soma(0.5).type21.t0	/= p
					n.innp.mean				*= p
					n.innp.stdev			*= np.sqrt(p)
				print " > check in methods F            :",checkinmethods("neuron/F")
				print "==================================\n"
			else:
				for n,p in zip(neurons,pl):
					try:
						exec "n.soma(0.5).type21."+k+" = p"
					except: 
						print "Cannot set parameter",k

	if checkinmethods('nstart'):
		if type(methods['nstart']) is list or type(methods['nstart']) is tuple:
			methods['nstart'] = {
				'delay'    : methods['nstart'][0],
				'Istdev'   : methods['nstart'][1],
				'duration' : methods['nstart'][2],
			}
		if not checkinmethods('nstart/Iapp'     ): methods['nstart']['Iapp'    ] = 0.
		if not checkinmethods('nstart/period'   ): methods['nstart']['period'  ] = 0.1
		if not checkinmethods('nstart/delay'    ): methods['nstart']['delay'   ] = 0.
		if not checkinmethods('nstart/duration' ): methods['nstart']['duration'] = 1e9
			
		if not checkinmethods('nstart/Istdev'):
			raise RuntimeError("/nstart/Istdev isn't set up")
		for n in neurons:
			n.addnoise(\
				Iapp  = methods['nstart']['Iapp'],\
				Insd  = methods['nstart']['Istdev'],\
				delay = methods['nstart']['delay'],\
				dur   = methods['nstart']['duration'],\
				per   = methods['nstart']['period'] )


	#DB>>
	if checkinmethods("DB-nrn"):
		h.finitialize()
		h.fcurrent()
		h.frecord_init()
		for n in neurons:
			print n.soma(0.5).type21.type21
			for x in "gl,el,n0,v12,sl,t0,st,v0,sg".split(","):
				print "    ","% 3s"%x,"=",eval("n.soma(0.5).type21."+x)
			print "-----"
		exit(0)
	#<<DB

	t = h.Vector()
	t.record(h._ref_t)


	#### Create Connection List:
	if checkinmethods("ncon"):
		def CreateConnectionList():
			def CreateFixNumberOrRange(n0,n1=None):
				OUTList = [ [] for x in xrange(methods["ncell"])]
				for toid in xrange(methods["ncell"]):
					from_excaption = [ 0 for x in xrange(methods["ncell"]) ]
					from_excaption[toid] = 1
					upcnt = 0
					total = 0
					if not n1 is None:
						neurons[toid].concnt = int(np.random.random()*(n1-n0) + n0)
					else:
						neurons[toid].concnt = n0
					while  upcnt < 10000*methods["ncell"]:
						upcnt += 1
						fromid = rnd.randint(0, methods["ncell"]-1)
						if from_excaption[fromid] == 1 : continue
						upcnt  = 0
						total += 1
						from_excaption[fromid] = 1
						OUTList[toid].append(fromid)
						if total >= neurons[toid].concnt :break
					else:
						sys.stderr.write("Couldn't obey connections conditions\nNeuron:%d TOTLA:%d CURRENT:%d\n"%(toid,n0,total))
						for x in OUTList:
							sys.stderr.write("ID:%d #%d\n"%(x[0],x[1]))
						sys.exit(1)
				return OUTList
			def CreateNormalDistribution(mean,stdev=0.):
				OUTList = [ [] for x in xrange(methods["ncell"])]
				for toid in xrange(methods["ncell"]):
					from_excaption = [ 0 for x in xrange(methods["ncell"]) ]
					from_excaption[toid] = 1
					upcnt = 0
					total = 0
					neurons[toid].concnt = int( positiveGauss(mean,stdev) )
					while  upcnt < 10000*methods["ncell"]:
						upcnt += 1
						fromid = rnd.randint(0, methods["ncell"]-1)
						if from_excaption[fromid] == 1 : continue
						upcnt  = 0
						total += 1
						from_excaption[fromid] = 1
						OUTList[toid].append(fromid)
						if total >= neurons[toid].concnt :break
					else:
						sys.stderr.write("Couldn't obey connections conditions\nNeuron:%d TOTLA:%d CURRENT:%d\n"%(toid,n0,total))
						for x in OUTList:
							sys.stderr.write("ID:%d #%d\n"%(x[0],x[1]))
						sys.exit(1)
				return OUTList
			def CreateBinomialDistribution(prob):
				OUTList = [ [] for x in xrange(methods["ncell"])]
				for toid in xrange(methods["ncell"]):
					for fromid in xrange(methods["ncell"]):
						if fromid == toid: continue
						if np.random.random() > prob : continue
						OUTList[toid].append(fromid)
						neurons[toid].concnt += 1
					#DB>>
					#print OUTList[toid]
					#<<DB
				return OUTList
#>>			
			#DB>>
			#print type(methods["ncon"]),methods["ncon"]
			#<<DB
			if type(methods["ncon"]) is int:
				#DB>>
				#print "ncon - int"
				#<<DB
				return CreateFixNumberOrRange(methods["ncon"])
			elif type(methods["ncon"]) is tuple or type(methods["ncon"]) is list:
				#DB>>
				#print "Ncon tuple or list"
				#<<DB
				if type(methods["ncon"][0]) is int:
					if len(methods["ncon"]) > 2:
						methods["normalize-weight-by-ncon"] = methods["ncon"][2]
						return CreateFixNumberOrRange(methods["ncon"][0],methods["ncon"][1])
					elif len(methods["ncon"]) > 1:
						return CreateFixNumberOrRange(methods["ncon"][0],methods["ncon"][1])
					else:
						return CreateFixNumberOrRange(methods["ncon"][0])
				elif type(methods["ncon"][0]) is str:
					if methods["ncon"][0] == "u":
						#print "  > Uniform Distribution "
						if len(methods["ncon"]) > 3:
							methods["normalize-weight-by-ncon"] = methods["ncon"][3]
							return CreateFixNumberOrRange(methods["ncon"][1],methods["ncon"][2])
						elif len(methods["ncon"]) > 2:
							return CreateFixNumberOrRange(methods["ncon"][1],methods["ncon"][2])
						elif len(methods["ncon"]) > 1:
							return CreateFixNumberOrRange(methods["ncon"][1])
						else:
							print "ERROR in ncon parameter:\nUSAGE of uniform distribution: /ncom=('u',n-from, {n-to}, {{norm-by}})"
							sys.exit(1)
				
					if methods["ncon"][0] == "n":
						if len(methods["ncon"]) > 3:
							methods["normalize-weight-by-ncon"] = methods["ncon"][3]
							return CreateNormalDistribution(methods["ncon"][1],methods["ncon"][2])
						elif len(methods["ncon"]) > 2:
							return CreateNormalDistribution(methods["ncon"][1],methods["ncon"][2])
						elif len(methods["ncon"]) > 1:
							return CreateFixNumberOrRange(methods["ncon"][1])
						else:
							print "ERROR in ncon parameter:\nUSAGE for noormal distribution: /ncom=('n',mean, {stdev}, {{norm-by}})"
							sys.exit(1)

					if methods["ncon"][0] == "b":
						if len(methods["ncon"]) > 2:
							methods["normalize-weight-by-ncon"] = methods["ncon"][1]
							return CreateBinomialDistribution(methods["ncon"][1])
						elif len(methods["ncon"]) > 1:
							return CreateBinomialDistribution(methods["ncon"][1])
						else:
							print "ERROR in ncon parameter:\nUSAGE for binomial distribution: /ncom=('b',prob,{norm-by})"
							sys.exit(1)

		print "=================================="
		print "===        Map Connections     ==="
		if checkinmethods('connectom'):
			print "  > Try to Read Connectom file :", methods['connectom']
			if os.access(methods['connectom'],os.R_OK):
				print "  > File is accessible         : "
				with open(methods['connectom'],"r") as fd:
					xncell = pickle.load(fd)
					xnconn = pickle.load(fd)
					if xncell != methods['ncell'] or xnconn != methods['ncon']:
						print "  > File has different numbers : "
						print "  > n cell                     : ", xncell,"|",methods['ncell']
						print "  > n connection               : ", xnconn,"|",methods['ncon']
						OUTList = CreateConnectionList()
					else:
						print "  > Read Connection Map        : ",
						OUTList = pickle.load(fd)
						for n,cpn in zip(neurons,pickle.load(fd)):
							n.concnt = cpn
						if not checkinmethods("normalize-weight-by-ncon"):
							methods["normalize-weight-by-ncon"] = pickle.load(fd)
						print "Successfully"
					
			elif not os.access(methods['connectom'],os.F_OK):
				print "  > File dos not exist         : try to create"
				OUTList = CreateConnectionList()
				with open(methods['connectom'],"w") as fd:
					pickle.dump(methods['ncell'],fd)
					pickle.dump(methods['ncon'],fd)
					pickle.dump(OUTList,fd) 
					pickle.dump([ n.concnt for n in neurons ],fd)
					if checkinmethods("normalize-weight-by-ncon"):
						pickle.dump(methods["normalize-weight-by-ncon"],fd)
					else:
						pickle.dump(False,fd)
			else:
				print
				print "============= ERROR ============="
				print " > Cannot create file \'{}\' ".format(methods['connectom'])
				print
				exit(0)
		else:
			print " > Generate connections         :",
			OUTList = CreateConnectionList()
			print "Successfully"
		print "==================================\n"

	#DB>
		#for i in OUTList:
			#print len(i)
			#for j in i:	print "%03d"%(j),
			#print
		#sys.exit(0)
	#<DB
	if checkinmethods('cycling'):
		print "=================================="
		print "===      Cycles counting       ==="
		mat=np.matrix( np.zeros((methods["ncell"],methods["ncell"])) )
#		for i,vec in map(None,xrange(methods["ncell"]),OUTList):
		for i,vec in enumerate(OUTList):
			mat[i,vec]=1
		kx = []
		for cnt in xrange(methods['cycling']):
			kx.append(np.trace(mat)/methods["ncell"])
			mat = mat.dot(mat)
		print " > Cyclopedic numbers           : ",kx
		print "==================================\n"
		methods['cycling-result'] = kx
		del mat
		
		
	

	#### Create Conneactions:
	if checkinmethods("ncon"):
		print "=================================="
		print "===    Make the Connections    ==="
		print "==================================\n"
		connections = []
		if not checkinmethods("gtot-dist"):  methods["gtot-dist"]  = "NORM"
		if not checkinmethods("gsyn-dist"):  methods["gsyn-dist"]  = "NORM"
		if not checkinmethods("delay-dist"): methods["delay-dist"] = "NORM"
#		for x in map(None,xrange(methods["ncell"]),OUTList):
		for x in enumerate(OUTList):
			if checkinmethods("neuron/F"):
				gx = methods["neuron"]["F"][x[0]]
			elif type(methods['synapse']['gsynscale']) is int or type(methods['synapse']['gsynscale']) is float:
				gx = float(methods['synapse']['gsynscale'])
			elif type(methods['synapse']['gsynscale']) is tuple :
				#DB>>
				#print "DB: TUPLE size=",len(methods['synapse']['gsynscale'])
				#print "DB: DIST",methods["gtot-dist"]
				#exit(0)
				#<<DB
				if methods['synapse']['gsynscale'][1] <= 0.: 
					gx  = methods['synapse']['gsynscale'][0] 
					gx *= 1. if len(methods['synapse']['gsynscale']) < 3 else methods['synapse']['gsynscale'][2]
				elif methods["gtot-dist"] == "NORM":
					### Trancated normal
					gx = positiveGauss(methods['synapse']['gsynscale'][0],methods['synapse']['gsynscale'][1])
				elif methods["gtot-dist"] == "LOGN":
					### Lognormal
					if len(methods['synapse']['gsynscale']) != 2 and len(methods['synapse']['gsynscale']) != 3:
						print "ERROR: wrong scaler size!\n/synapse/gsynscale should have 2 or more parameters"
						exit(1)
					elif len(methods['synapse']['gsynscale']) == 2:
						gsymtotm,gsymtots,gsyntotsc=methods['synapse']['gsynscale'],1.
						##DB>>
						#print gsymtotm,gsymtots,gsyntotsc
						#exit(0)
						##<<DB
					elif len(methods['synapse']['gsynscale']) == 3:
						gsymtotm,gsymtots,gsyntotsc=methods['synapse']['gsynscale']
						##DB>>
						#print gsymtotm,gsymtots,gsyntotsc
						#exit(0)
						##<<DB
					gx = gsyntotsc*np.random.lognormal(mean=np.log(gsymtotm/np.sqrt(1.+gsymtots**2/gsymtotm**2)),sigma=np.sqrt(np.log(1.+gsymtots**2/gsymtotm**2)))
					#DB>>
					#print "DB: gx=",gx
					#<<DB
			else:
				print "ERROR: wrong type of /synapse/gsynscale"
				exit(1)
				
			neurons[x[0]].gsynscale = gx
			for pre in x[1]:
				if type(methods['synapse']['delay']) is tuple :
					if methods['synapse']['delay'][1] <= 0: dx = float(methods['synapse']['delay'][0])
					elif methods["delay-dist"] == "NORM":
						### Trancated normal
						dx = positiveGauss(methods['synapse']['delay'][0],methods['synapse']['delay'][1])
						if len(methods['synapse']['delay']) < 3:
							while(dx < methods['timestep']*2):
								dx = positiveGauss(methods['synapse']['delay'][0],methods['synapse']['delay'][1])
						else:
							while(dx < methods['synapse']['delay'][2]):
								dx = positiveGauss(methods['synapse']['delay'][0],methods['synapse']['delay'][1])
							
					elif methods["delay-dist"] == "LOGN":
						### Lognormal
						dlym,dlys=methods['synapse']['delay'][0]-methods['timestep']*2.,methods['synapse']['delay'][1]
						if len(methods['synapse']['delay']) < 3:
							dx = np.random.lognormal(mean=np.log(dlym/np.sqrt(1.+dlys**2/dlym**2)),sigma=np.sqrt(np.log(1.+dlys**2/dlym**2)))+methods['timestep']*2
						else:
							dx = np.random.lognormal(mean=np.log(dlym/np.sqrt(1.+dlys**2/dlym**2)),sigma=np.sqrt(np.log(1.+dlys**2/dlym**2)))
							while dx < methods['synapse']['delay'][2]:
								dx = np.random.lognormal(mean=np.log(dlym/np.sqrt(1.+dlys**2/dlym**2)),sigma=np.sqrt(np.log(1.+dlys**2/dlym**2)))
					elif methods["delay-dist"] == "LOGN-SHIFT":
						### Lognormal
						dlym,dlys=methods['synapse']['delay'][0]-methods['timestep']*2.,methods['synapse']['delay'][1]
						if len(methods['synapse']['delay']) < 3:
							dx = np.random.lognormal(mean=np.log(dlym/np.sqrt(1.+dlys**2/dlym**2)),sigma=np.sqrt(np.log(1.+dlys**2/dlym**2)))+methods['timestep']*2
						else:
							dx = np.random.lognormal(mean=np.log(dlym/np.sqrt(1.+dlys**2/dlym**2)),sigma=np.sqrt(np.log(1.+dlys**2/dlym**2)))+methods['synapse']['delay'][2]
					elif methods["delay-dist"] == "DIST":
						dmin = methods['synapse']['delay'][0]
						dinciment = methods['synapse']['delay'][1] if len(methods['synapse']['delay']) >= 2 else 0.
						dx = dmin + dinciment*float(abs(pre-x[0]))
					elif methods["delay-dist"] == "RING":
						dmin = methods['synapse']['delay'][0]
						dinciment = methods['synapse']['delay'][1] if len(methods['synapse']['delay']) >= 2 else 0.
						dist = min([float(abs(pre-x[0])),float(abs(pre-x[0]-methods['ncell'])),float(abs(pre-x[0]+methods['ncell']))])
						dx = dmin + dinciment*dist
					elif methods["delay-dist"] == "UNIFORM":
						dmin = methods['synapse']['delay'][0]
						dmax = methods['synapse']['delay'][1]
						dx = dmin + np.random.rand()*(dmax-dmin)
						
				else:
					dx = float(methods['synapse']['delay'])
				if checkinmethods("gtot-set"):
					wx=1.
				else:
					if type(methods['synapse']['weight']) is tuple :
						if methods['synapse']['weight'][1] <= 0: wx = methods['synapse']['weight'][0]
						elif methods["gsyn-dist"] == "NORM":
							#### Trancated normal
							wx = methods['synapse']['weight'][1]*np.random.randn()+methods['synapse']['weight'][0]
							while wx < 0.0 : wx = methods['synapse']['weight'][1]*np.random.randn()+methods['synapse']['weight'][0]
						elif methods["gsyn-dist"] == "LOGN":
							### Lognormal
							wm,ws=methods['synapse']['weight']
							wx = np.random.lognormal(mean=np.log(wm/np.sqrt(1.+ws**2/wm**2)),sigma=np.sqrt(np.log(1.+ws**2/wm**2)))
							#wx = np.random.lognormal(mean=np.log(wm**2/np.sqrt(wm**2+ws**2)),sigma=np.sqrt(np.log(1.+ws**2/wm**2)))
						elif methods["gsyn-dist"] == "DIST":
							wmin = methods['synapse']['weight'][0]
							winciment = methods['synapse']['weight'][1] if len(methods['synapse']['weight']) >= 2 else 0.
							wx = wmin + winciment*float(abs(pre-x[0]))
						elif methods["gsyn-dist"] == "RING":
							wmin = methods['synapse']['weight'][0]
							winciment = methods['synapse']['weight'][1] if len(methods['synapse']['weight']) >= 2 else 0.
							wist = min([float(abs(pre-x[0])),float(abs(pre-x[0]-methods['ncell'])),float(abs(pre-x[0]+methods['ncell']))])
							wx = wmin + winciment*wist
					else:
						wx = float(methods['synapse']['weight'])
					#if type(methods["ncon"]) is tuple or type(methods["ncon"]) is list:
						#if len(methods["ncon"]) > 2:
							#wx *= float(methods["ncon"][2])/float(neurons[x[0]].concnt)
					if checkinmethods("normalize-weight-by-ncon"):
						#DB>>
						#print "Norm by Factor",float(methods["normalize-weight-by-ncon"]),float(neurons[x[0]].concnt),float(methods["normalize-weight-by-ncon"])/float(neurons[x[0]].concnt)
						#<<DB

						wx *= float(methods["normalize-weight-by-ncon"])/float(neurons[x[0]].concnt)
					if methods['taunorm'] and not methods['synapse']['synscaler'] is None:
						#DB print "norm by factor",1./neurons[x[0]].tsynscale
						wx /= neurons[x[0]].tsynscale
				#####DB>>
				#print "DB:gx=",gx,"dx=",dx,"wx=",wx
				#####<<DB
				#connections.append( (h.NetCon(neurons[pre].soma(0.5)._ref_v,neurons[x[0]].isyn,
						#0., dx, gx*wx,
						#sec=neurons[pre].soma),pre,x[0]) )
				#connections.append( (h.NetCon(neurons[pre].soma(0.5)._ref_v,neurons[x[0]].isyn,
						#25., dx, gx*wx,
						#sec=neurons[pre].soma),pre,x[0]) )
				connections.append( (h.NetCon(neurons[pre].soma(0.5)._ref_v,neurons[x[0]].isyn,
						0., dx, gx*wx,
						sec=neurons[pre].soma),pre,x[0]) )
				neurons[x[0]].gtotal += gx*wx
		if checkinmethods('Conn-alter'):
			print "================================================"
			print "===        ALTER CONNECTIONS SETTINGS        ==="
			n_alter = methods['Conn-alter']['n']      if checkinmethods('Conn-alter/n')      else 1
			d_alter = methods['Conn-alter']['delay']  if checkinmethods('Conn-alter/delay')  else 0.8
			w_alter = methods['Conn-alter']['weight'] if checkinmethods('Conn-alter/weight') else 0.1e-2
			alter = range(len(connections))
			for i in xrange(n_alter):
				Xalter = alter[np.random.randint(len(alter))]
				connections[Xalter][0].weight[0] = w_alter
				connections[Xalter][0].delay     = d_alter
				print " > %03d -> %03d                                 : were altered"%(connections[Xalter][1],connections[Xalter][2])
				alter.remove(Xalter)
			print "Total number of connections                   :",len(connections)
			print "Number of altered connections                 :",n_alter
			print "Procentage                                    :",n_alter*100/len(connections)
			print "================================================\n"
				
		if checkinmethods('Conn-rec'):
			methods['Conn-rec-results'] = [ (n[1],n[2],n[0].weight[0],n[0].delay) for n in connections ]
		if checkinmethods('Conn-stat'):
			print "================================================"
			print "===           Connections Statistics         ==="
			statn = np.array( [ float(len(o)) for o in OUTList ] )
			meann,stdrn = np.mean(statn),np.std(statn)
			minin,maxin = np.min(statn), np.max(statn)
			statw = np.array( [ n[0].weight[0] for n in connections ] )
			meanw,stdrw = np.mean(statw),np.std(statw)
			miniw,maxiw = np.min(statw), np.max(statw)
			statd = np.array( [ n[0].delay for n in connections ] )
			meand,stdrd = np.mean(statd),np.std(statd)
			minid,maxid = np.min(statd), np.max(statd)
			methods['Conn-stat-results'] = {
				'ncon': {
					'mean':meann, 'stdr':stdrn, 'min':minin, 'max':maxin
				},
				'weight':{
					'mean':meanw, 'stdr':stdrw, 'min':miniw, 'max':maxiw
				},
				'delay':{
					'mean':meand, 'stdr':stdrd, 'min':minid, 'max':maxid
				}
			}
			print " > Number min / max / mean / stdev / CV       :",minin,"/",maxin,"/",meann,"/",stdrn,"/",stdrn/meann
			print " > Weight min / max / mean / stdev / CV       :",miniw,"/",maxiw,"/",meanw,"/",stdrw,"/",stdrw/meanw
			print " > Delay  min / max / mean / stdev / CV       :",minid,"/",maxid,"/",meand,"/",stdrd,"/",stdrd/meand
			print "================================================\n"
		#DB>>
		#plt.figure(0)
		#w=np.array([c[0].weight[0] for c in connections])
		#print np.mean(w), np.std(w)
		#plt.hist(w,bins=50,range=(0,1e-6))
		#plt.show()
		#exit(0)
		#<<DB
			
	
	#### Create gapjunctions:
	if checkinmethods('gapjunction'):
		if   not 'ncon' in methods['gapjunction']            : GJList = OUTList
		elif methods['gapjunction']['ncon'] is None          : GJList = OUTList
		elif not type(methods['gapjunction']['ncon']) is int : GJList = OUTList
		elif not methods['gapjunction']['ncon'] > 0          : GJList = OUTList
		else:
			GJList = [ [] for x in xrange(methods['ncell'])]
			gjncon = methods['gapjunction']['ncon']
			print "=================================="
			print "===       Map Gap-junctions     ==="
			print "==================================\n"
			for toid in xrange(methods['ncell']):
				from_excaption = [ 0 for x in xrange(methods['ncell']) ]
				from_excaption[toid] = 1
				upcnt = 0
				total = 0
				if type(gjncon) is tuple or type(gjncon) is list:
					if len(gjncon) > 1:
						neurons[toid].g_jcnt = int(np.random.random()*(gjncon[1]-gjncon[0]) + gjncon[0])
					else:
						neurons[toid].g_jcnt = gjncon[0]
				else:
					neurons[toid].g_jcnt = gjncon
				while  upcnt < 10000*methods['ncell']:
					upcnt += 1
					fromid = rnd.randint(0, methods['ncell']-1)
					if from_excaption[fromid] == 1 : continue
					upcnt  = 0
					total += 1
					from_excaption[fromid] = 1
					GJList[toid].append(fromid)
					if total >= neurons[toid].g_jcnt :break
				else:
					sys.stderr.write("Couldn't obey connections conditions\nNeuron:%d TOTLA:%d CURRENT:%d\n"%(toid,gjncon,total))
					for x in GJList:
						sys.stderr.write("ID:%d #%d\n"%(x[0],x[1]))
					sys.exit(1)
		
		print "=================================="
		print "===    Make the Gap-Junctions   ==="
		print "==================================\n"
		gapjuctions = []
		for cellid,gjlst in enumerate(GJList):
			for preidx in gjlst:
				gj0,gj1 = h.gap(0.5, sec=neurons[cellid].soma), h.gap(0.5, sec=neurons[preidx].soma)
				h.setpointer(neurons[preidx].soma(.5)._ref_v, 'vgap', gj0)
				h.setpointer(neurons[cellid].soma(.5)._ref_v, 'vgap', gj1)
				gj0.r, gj1.r = methods['gapjunction']['r'],methods['gapjunction']['r']
				gapjuctions.append( (gj0,gj1,neurons[cellid].soma,neurons[preidx].soma) )

	if checkinmethods('simvar'):
		print "=================================="
		print "===      SIMVAR was found!     ==="
		print "==================================\n"
		simvars = []
		if methods['simvar']['type'] == 'n':
			for n in neurons:
				sv = h.variator(0.5, sec=n.soma)
				exec "h.setpointer(n.soma(0.5)."+methods['simvar']["var"]+", \'var\',sv)"
				simvars.append(sv)
			simvarrec = h.Vector()
			exec "simvarrec.record(neurons[0].soma(0.5)."+methods['simvar']["var"]+")"
		#elif methods['simvar']['type'] == 'c':
			#for n in neurons:
				#sv = h.variator()
				#exec "h.setpointer(n."+methods['simvar']["var"]+", \'var\',sv)"
				#simvars.append(sv)
			#simvarrec = h.Vector()
			#exec "simvarrec.record(neurons[0]."+methods['simvar']["var"]+")"
		elif methods['simvar']['type'] == 'g':
			for g0,g1,s0,s1 in gapjuctions:
				sv = h.variator(0.5, sec=s0)
				#DB>>
				#print "h.setpointer(g0."+methods['simvar']["var"]+", \'var\',sv)"
				#<<DB
				exec "h.setpointer(g0."+methods['simvar']["var"]+", \'var\',sv)"
				simvars.append(sv)
				sv = h.variator(0.5, sec=s1)
				exec "h.setpointer(g1."+methods['simvar']["var"]+", \'var\',sv)"
				simvars.append(sv)
			simvarrec = h.Vector()
			exec "simvarrec.record(gapjuctions[0][0]."+methods['simvar']["var"]+")"
		for sv in simvars:
			sv.a0 = methods['simvar']['a0']
			sv.a1 = methods['simvar']['a1']
			sv.t0 = methods['simvar']['t0']
			sv.t1 = methods['simvar']['t1']
	
	if checkinmethods('external'):
		ex_netstim	= h.NetStim(.5, sec = neurons[0].soma)
		if type(methods['external']) is list:
               #              0      1          2         3     4    5    6    7                8
               #/external=\(Start,interval,spike-count,weight,Esyn,Tau1,Tau2,delay,probability of connections\)
			if len(methods['external']) < 8:
			   methods['external'] = {
					'start'       :methods['external'][0],
					'interval'    :methods['external'][1],
					'count'       :methods['external'][2],
					'weight'      :methods['external'][3],
					'E'           :methods['external'][4],
					'tau1'        :methods['external'][5],
					'tau2'        :methods['external'][6]
			   }
			elif len(methods['external']) < 9:
			   methods['external'] = {
					'start'       :methods['external'][0],
					'interval'    :methods['external'][1],
					'count'       :methods['external'][2],
					'weight'      :methods['external'][3],
					'E'           :methods['external'][4],
					'tau1'        :methods['external'][5],
					'tau2'        :methods['external'][6],
					'delay'       :methods['external'][8]
			   }
			elif len(methods['external']) >= 9:
			   methods['external'] = {
					'start'       :methods['external'][0],
					'interval'    :methods['external'][1],
					'count'       :methods['external'][2],
					'weight'      :methods['external'][3],
					'E'           :methods['external'][4],
					'tau1'        :methods['external'][5],
					'tau2'        :methods['external'][6],
					'delay'       :methods['external'][8],
					'p'           :methods['external'][9]
			   }
			if type(methods['external']['delay']) is list or type(methods['external']['delay']) is tuple:
				if len(methods['external']['delay']) < 2:
					methods['external']['delay'] = methods['external']['delay'][0]
				else:
					methods['external']['delay'] = {
						'mean'  : methods['external']['delay'][0],
						'stdev' : methods['external']['delay'][1]
					}
		if not checkinmethods('external/start'   ): methods['external']['start'   ] = methods["tstop"]/3.
		if not checkinmethods('external/interval'): methods['external']['interval'] = methods["tstop"]/6.
		if not checkinmethods('external/count'   ): methods['external']['count'   ] = 1
		if not checkinmethods('external/E'       ): methods['external']['E'       ] = 0
		if not checkinmethods('external/tau1'    ): methods['external']['tau1'    ] = 0.8
		if not checkinmethods('external/tau2'    ): methods['external']['tau2'    ] = 1.2
		if not checkinmethods('external/weight'  ): methods['external']['weight'  ] = 0.
		if not checkinmethods('external/delay'   ): methods['external']['delay'   ] = 1.
		print "================================================"
		print "===              External Input              ==="
		print " > Start                                      :", methods['external']['start']
		print " > Interval                                   :", methods['external']['interval']
		print " > Count                                      :", methods['external']['count']
		if  checkinmethods('external/p'):
			print " > P                                          :", methods['external']['p']
		print " > Reversal potential                         :", methods['external']['E']
		print " > Tau 1                                      :", methods['external']['tau1']
		print " > Tau 2                                      :", methods['external']['tau2']
		print " > Weight                                     :", methods['external']['weight']
		if checkinmethods('external/delay/mean') or checkinmethods('external/delay/stdev'):
			print " > Delay                                    "
			if checkinmethods('external/delay/mean'):
				print "   > mean                                     :", methods['external']['delay']['mean']
			if checkinmethods('external/delay/stdev'):
				print "   > stdev                                    :", methods['external']['delay']['stdev']
		else:
			print " > Delay                                      :", methods['external']['delay']
		print "================================================\n"
		ex_netstim.start	= methods['external']['start'   ]
		ex_netstim.noise	= 0
		ex_netstim.interval	= methods['external']['interval'] 
		ex_netstim.number	= methods['external']['count']
		ex_syn,ex_netcon = [],[]
		for n in neurons:
			if  checkinmethods('external/p'):
				if rnd.random() > methods['external']['p']: continue
			ex_syn_new = h.Exp2Syn(0.5, sec=n.soma)
			ex_syn_new.e	= methods['external']['E'   ]
			ex_syn_new.tau1	= methods['external']['tau1'] if checkinmethods('external/tau1') else 0.8
			ex_syn_new.tau2	= methods['external']['tau2'] if checkinmethods('external/tau2') else 1.2
			ex_syn.append(ex_syn_new)
			if checkinmethods('external/delay/mean') and checkinmethods('external/delay/stdev'):
				exdelay = -1.0
				while exdelay <= 0.0 : exdelay = methods['external']['delay']['mean']+np.random.randn()*methods['external']['delay']['stdev']
			elif checkinmethods('external/delay/mean'):
				exdelay = methods['external']['delay']['mean']
			else :
				exdelay = methods['external']['delay']
			ex_netcon_new	= h.NetCon(ex_netstim, ex_syn_new, 1,exdelay , methods['external']['weight'], sec = n.soma)
			ex_netcon.append(ex_netcon_new)

	if checkinmethods("wmod"):
		print "=================================="
		print "===      Weight Modulator      ==="
		if not checkinmethods("wmod/scale"          ):methods["wmod"]["scale"        ] = 2.
		if not checkinmethods("wmod/time-points"    ):methods["wmod"]["time-points"  ] = [0., methods['tstop']]
		if not checkinmethods("wmod/weight-points"  ):methods["wmod"]["weight-points"] = [ methods["synapse"]["weight"], methods["wmod"]["scale"]*methods['synapse']["weight"] ]
		print " > Scale                        : ",methods["wmod"]["scale"        ]
		print " > Time points                  : ",methods["wmod"]["time-points"  ]
		print " > Weight points                : ",methods["wmod"]["weight-points"]
		wmodT, wmodW = h.Vector(), h.Vector()
		wmodT.from_python(methods["wmod"]["time-points"  ])
		wmodW.from_python(methods["wmod"]["weight-points"  ])
		for c,pre,post in connections:
			wmodW.play(c._ref_weight[0],wmodT,1)
		print "==================================\n"
			

	if checkinmethods("imod") and checkinmethods("neuron/Iapp"):
		print "=================================="
		print "===     Current Modulator      ==="
		if not checkinmethods("imod/scale"          ):methods["imod"]["scale"         ] = 2.
		if not checkinmethods("imod/time-points"    ):methods["imod"]["time-points"   ] = [0.                           , methods['tstop']]
		if not checkinmethods("imod/current-points" ):methods["imod"]["current-points"] = [-1.*methods["neuron"]["Iapp"], -1.*methods["imod"]["scale"]*methods["neuron"]["Iapp"] ]
		print " > Time Points                  : ",methods["imod"]["time-points"   ]
		print " > Current Points               : ",methods["imod"]["current-points"]
		imodAll = []
		
		for idx, n in enumerate(neurons):
			imodt, imodw = h.Vector(), h.Vector()
			imodt.from_python(methods["imod"]["time-points"     ])
			imodw.from_python(methods["imod"]["current-points"  ])
			imodw.play(n.innp._ref_mean,imodt,1)
			imodAll.append( (imodt,imodw) )
			#++++
		print "==================================\n"


	##DB>>
	#for n in neurons:
		#print n.isyn.e
	#exit(0)
	##<<DB

	print "=================================="
	print "===           RUN              ==="
	npc = h.ParallelContext()
	if checkinmethods("ncon"):
		mindel = np.array([ x[0].delay for x in connections ] )
		mindel = np.min(mindel)
		if mindel > h.dt*2:
			if type(methods['corefunc']) is int:
				npc.nthread(methods['corefunc'])
				sys.stderr.write( " > Setup                            : %g core\n"%(methods['corefunc']) )
				print " > Setup                        : %g core"%(methods['corefunc']) 
			else:
				#### Setup parallel context if there are delays
				if not os.path.exists("/etc/beowulf") and os.path.exists("/sysini/bin/busybox"):
					#I'm not on head node. I can use all cores (^-^)
					methods['corefunc'] = methods['corefunc'][2]
					npc.nthread(methods['corefunc'])
					sys.stderr.write( " > Setup                        : %g core\n"%(methods['corefunc']) )
					print " > Setup                        : %g core"%(methods['corefunc']) 
				elif os.path.exists("/etc/beowulf"):
					#I'm on head node. I grub only half (*_*)
					methods['corefunc'] = methods['corefunc'][1]
					npc.nthread(methods['corefunc'])
					sys.stderr.write( " > Setup                        : %g core\n"%(methods['corefunc']) )
					print " > Setup                        : %g core"%(methods['corefunc']) 
				else:
					#I'm on Desktop (-.-)
					methods['corefunc'] = methods['corefunc'][0]
					npc.nthread(methods['corefunc'])
					sys.stderr.write( " > Setup                        : %g cores\n"%(methods['corefunc']) )
					print " > Setup                        : %g core"%(methods['corefunc']) 
		else:
			#I'm in (_!_)
			methods['corefunc'] = 0
			#npc.nthread(methods['corefunc'])
			sys.stderr.write( " > Setup                        : %g cores\n"%(methods['corefunc']) )
			print " > Setup                        : %g core"%(methods['corefunc']) 

	if checkinmethods("cvode"):
		cvode = h.CVode()
		cvode.active(1)
		print " > CVODE                    : ON"
		
	h.finitialize()
	h.fcurrent()
	h.frecord_init()
		
	while h.t < methods['tstop']:h.fadvance()

	print "==================================\n"
		
	

	print "=================================="
	print "===          Analysis          ==="
	print "==================================\n"
	
	t = np.array(t)
	if checkinmethods('gui'):
		plt.rc('text', usetex = True )
		plt.rc('font', family = 'serif')
		plt.rc('svg', fonttype = 'none')
		mainfig = plt.figure(1)
		if checkinmethods("MainFigTitle"):
			mainfig.suptitle(methods["MainFigTitle"])
		if checkinmethods('traceView'):
			cid = mainfig.canvas.mpl_connect('button_press_event', onclick1)
		nplot = 411
		if checkinmethods('simvar')      : nplot += 100
		if checkinmethods('spectrogram') : nplot += 100
		p = plt.subplot(nplot)
			
		tprin=np.array(t)
		tprin = tprin[ np.where( tprin < methods['tv'][1] ) ]
		if methods['tv'][0] <= 0.:
			tproc = 0
		else:
			tproc = tprin[ np.where( tprin < methods['tv'][0] ) ]
			tproc = tproc.shape[0]
		tprin = tprin[tproc:]
		vindex = (methods["ncell"]-1)/2
		##PLOT VOLTAGE>>
		vtrace, = plt.plot(tprin,np.array(neurons[vindex].volt)[tproc:tprin.size+tproc],"k")
		plt.ylim(ymax=40.)
		mainfig.canvas.mpl_connect('key_press_event',neuronsoverview)
		plt.ylabel("Voltage (mV)", fontsize=16)
		if methods["external"]:
			ex0 = methods["external"]['start']
			ex1 = methods["external"]['interval']
			for ex2 in xrange(methods["external"]['count']):
				plt.plot([ex0+ex1*ex2,ex0+ex1*ex2],[0,30],"r--")
		plt.subplot(nplot+1,sharex=p)
		nurch = np.arange(1,methods["ncell"]+1,1)
		if checkinmethods('sortbysk'):
			if methods['sortbysk'] == 'I':
				nindex = [ (-neurons[i].innp.mean,i) for i in xrange(methods["ncell"])]
				nindex.sort()
				#DB>>
				#print nindex
				#<<DB
				for i in xrange(methods["ncell"]):
					nurch[nindex[i][1]]=i
				#DB>>
				#for i in xrange(methods["ncell"]):
				#	print i,'->',nurch[i],'=',-neurons[nurch[i]].innp.mean,"|",-neurons[i].innp.mean
				#exit(0)
				#<<DB
			if methods['sortbysk'] == 'N':
				nindex = [ (-neurons[i].innp.stdev,i) for i in xrange(methods["ncell"])]
				nindex.sort()
				#print nindex
				for i in xrange(methods["ncell"]):
					nurch[nindex[i][1]]=i
			if methods['sortbysk'] == 'T':
				nindex = [ (neurons[i].type21,i) for i in xrange(methods["ncell"])]
				nindex.sort()
				#print nindex
				for i in xrange(methods["ncell"]):
					nurch[nindex[i][1]]=i
			if methods['sortbysk'] == 'G':
				nindex = [ (neurons[i].gsynscale,i) for i in xrange(methods["ncell"])]
				nindex.sort()
				#print nindex
				for i in xrange(methods["ncell"]):
					nurch[nindex[i][1]]=i
			if methods['sortbysk'] == 'NC':
				nindex = [ (-neurons[i].concnt,i) for i in xrange(methods["ncell"])]
				nindex.sort()
				#print nindex
				for i in xrange(methods["ncell"]):
					nurch[nindex[i][1]]=i
			if methods['sortbysk'] == 'GT':
				nindex = [ (neurons[i].gtotal,i) for i in xrange(methods["ncell"])]
				nindex.sort()
				#print nindex
				for i in xrange(methods["ncell"]):
					nurch[nindex[i][1]]=i
			if methods['sortbysk'] == 'ST':
				nindex = [ (neurons[i].tsynscale,i) for i in xrange(methods["ncell"])]
				nindex.sort()
				#print nindex
				for i in xrange(methods["ncell"]):
					nurch[nindex[i][1]]=i
			if methods['sortbysk'] == 'FR':
				nindex = [ (neurons[i].spks.size(),i) for i in xrange(methods["ncell"])]
				nindex.sort()
				#print nindex
				for i in xrange(methods["ncell"]):
					nurch[nindex[i][1]]=i
				
			#print nurch
	pmean, fmean = 0., 0
	pcnt,  fcnt  = 0 , 0
	nrnfr        = []  


	meancur = np.zeros(t.size)
	if checkinmethods('spectrogram'):
		populationcurrent = np.zeros(t.size)
	spbins  = np.zeros( int(np.floor(methods['tstop']))+1 )
	specX	= np.arange(spbins.size, dtype=float)
	specX	*= 1000.0/methods['tstop']
	#pnum	= specX.size/2
	pnum 	= int(200.*methods['tstop']/1000.0)
	specX	= specX[:pnum]
	if checkinmethods("nrnFFT"):
		specN	= np.zeros(pnum)
#	specV	= np.zeros(t.size())

	#NEED ISI for indices 
	if checkinmethods('N2NHI-netISI') or checkinmethods('N2NHI'):
		cg_nrnisi = []

	if 10 < methods["nrnISI"] <= 3000:
		isi		= np.zeros(methods["nrnISI"])
	if checkinmethods('coreindex'):
		coreindex = [0.0, 0.0]

	if checkinmethods('gui'):
		rast = []
		#xrast = np.array([])
		#yrast = np.array([])
	if checkinmethods('jitter-rec'):
		jallspikes = np.zeros(int((methods['tstop']-methods['cliptrn'])/methods['timestep'])\
		                       if checkinmethods('cliptrn') else \
		                      int(methods['tstop']/methods['timestep']) ) 
	
	analdur = (methods["tstop"]-methods['cliptrn']) if checkinmethods('cliptrn') else methods["tstop"]

	if checkinmethods("FI-curve"):
		FIcurve = []
	#for (idx,n) in map(None,xrange(methods["ncell"]),neurons):
	for idx,n in enumerate(neurons):
		n.spks = np.array(n.spks)
		if checkinmethods('gui'):
			if not methods['cliprst']:
				rast += [ (st,nurch[idx]) for st in n.spks if methods['tv'][0] < st < methods['tv'][1] ]
			elif nurch[idx]%methods['cliprst'] == 0:
				rast += [ (st,nurch[idx]) for st in n.spks if methods['tv'][0] < st < methods['tv'][1] ]
			#spk = n.spks[ np.where (n.spks < methods['tv'][1]) ]
			#if methods['tv'][0] > 0:
				#spk = spk[ np.where (spk > methods['tv'][0]) ]
			
				#xrast = np.append(xrast,spk)
				#yrast = np.append(yrast,np.repeat(,spk.size))
			##elif idx%methods['cliprst'] == 0:
			#elif nurch[idx]%methods['cliprst'] == 0:
				#xrast = np.append(xrast,spk)
				#yrast = np.append(yrast,np.repeat(nurch[idx],spk.size))
		
		if checkinmethods('cliptrn'):
			fstidx = np.where(n.spks > methods['cliptrn'] )[0]
			if len(fstidx) < 2:
				aisi = None
			else:
				fstidx = fstidx[0]
				aisi = n.spks[fstidx+1:] - n.spks[fstidx:-1]
		else:
			aisi = n.spks[1:] - n.spks[:-1]
		if checkinmethods('coreindex'):
			coreindex[0] += np.sum((aisi[1:] - aisi[:-1])/aisi[:-1])
			coreindex[1] += aisi.size - 1
		if 10 < methods["nrnISI"] <= 3000:
			for i in aisi[ np.where(aisi < methods["nrnISI"]) ]:
				isi[ int(np.floor(i)) ] += 1.0
		if not aisi is None  and( checkinmethods('N2NHI-netISI') or checkinmethods('N2NHI') ):
			cg_nrnisi += aisi.tolist()
		if not aisi is None:
			pmean += np.sum(aisi)
			pcnt  += aisi.shape[0]
			fr    = ( float(n.spks.shape[0] - fstidx)*1000./analdur ) if checkinmethods('cliptrn') else float(n.spks.shape[0])*1000./analdur
			if checkinmethods("FI-curve"):
				FIcurve.append( [ -n.innp.mean,fr ] )
			fmean += fr
			fcnt  += 1
			nrnfr.append(fr)
		if checkinmethods('gui'):
			if methods['tracetail'] == 'total current' or methods['tracetail'] == 'mean total current' or methods['tracetail'] == 'TI' or methods['tracetail'] == 'MTI':
				meancur += np.array(n.isyni.x) + np.array(n.inoise.x)
			elif methods['tracetail'] == 'total synaptic current' or methods['tracetail'] == 'mean total synaptic current' or methods['tracetail'] == 'TSI' or methods['tracetail'] == 'MTSI':
				meancur += np.array(n.isyni.x)
			elif methods['tracetail'] == 'total conductance' or methods['tracetail'] == 'mean total conductance' or methods['tracetail'] == 'TG' or methods['tracetail'] == 'MTG':
				meancur += np.array(n.isyng.x)
			
		if checkinmethods('spectrogram'):
			populationcurrent += np.array(n.isyni.x) + np.array(n.inoise.x)
			
		spn	= np.zeros(spbins.size)
		for sp in n.spks:
			spbins[ int( np.floor(sp) ) ] +=1
			spn[ int( np.floor(sp) ) ] +=1
		
			if checkinmethods('jitter-rec'):
				jps = int( (sp - methods['cliptrn'])/methods["timestep"] ) if checkinmethods('cliptrn') else int( sp/methods["timestep"] )
				jallspikes[jps] += 1
				
				
		if checkinmethods('cliptrn'):
			spn = spn[methods['cliptrn']:]
		if checkinmethods("nrnFFT"):
			fft = np.abs( np.fft.fft(spn)*1.0/methods['tstop'] )**2
			specN += fft[:pnum]
	
	methods["nrnPmean"] = None if pcnt < 1 else float(pmean)/float(pcnt)
	methods["nrnFmean"] = None if fcnt < 1 else float(fmean)/float(fcnt)
	methods["nrnfr"]    = nrnfr[:]
	print "=================================="
	print "===           Neurons          ==="
	print " > mean Period      (ms)        :",methods["nrnPmean"]
	print " > mean Firing Rate (Hz)        :",methods["nrnFmean"]
	print "==================================\n"

	if checkinmethods("FI-curve"):
		np.savetxt(methods["FI-curve"],np.array(sorted(FIcurve)) )
	
	if checkinmethods('jitter-rec') and checkinmethods('cliptrn'):
		jallspikes = jallspikes[ np.where( jallspikes > methods['cliptrn'] ) ]
		
	if checkinmethods('gui'):
		if not checkinmethods('rstmark'    ):methods['rstmark'    ]="."
		if not checkinmethods('rstmarksize'):methods['rstmarksize']=5
		#PLOT RASTER>>
		#plt.plot(xrast,yrast,"k"+methods['rstmark'],mew=0.75,ms=methods['rstmarksize'])#,ms=10)
		rast = np.array(rast)
		if rast.shape[0] > 2:
			plt.plot(rast[:,0],rast[:,1],"k"+methods['rstmark'],mew=0.75,ms=methods['rstmarksize'])#,ms=10)
		
		if methods['fullrast']	: plt.ylim(ymin=0,ymax=methods["ncell"])
		else			: plt.ylim(ymin=0)
	if checkinmethods("nrnFFT"):
		specN /= float(methods["ncell"])#	specV /= float(methods["ncell"])
		methods["nrnFFT-results"] = { 'sectrum':specN[:pnum], 'freq':specX }
	
	if checkinmethods('cliptrn'):
		spbins = spbins[methods['cliptrn']:]
	
	if checkinmethods('popfr'):
		popfr = np.mean(spbins)
		print "=================================="
		print "===       MEAN FIRING RATE     ==="
		print "  > MFR =           ",popfr
		print "==================================\n"
		methods['popfr-results'] = popfr

	if checkinmethods("netFFT") or checkinmethods("nrnFFT"):
		print "=================================="
		print "===            FFT             ==="
		print "==================================\n"
		fft = np.abs( np.fft.fft(spbins)*1.0/methods['tstop'] )**2
		methods["netFFT-results"] = { 'sectrum':fft[:pnum], 'freq':specX }

	##EN
	#probscale = np.zeros(methods["ncell"] + 1)
	#probscale[0] = 1./float(methods["ncell"] + 1)
	#for x in range(1,methods["ncell"] + 1):
		#probscale[x] = probscale[0]*probscale[x-1]
	#pspbin = np.array([ probscale[int(x)] for x in  spbins] )
	#en = np.sum( (-1)*pspbin*np.log(pspbin) )
	
	if checkinmethods('coreindex'):
		coreindex = coreindex[0]/coreindex[1]
		print coreindex, 1./(1.+ abs(coreindex))
		sys.exit(0)
		#with open("coreindex.csv","w") as fd:
			#for i in coreindex: fd.write("%g\n"%i)
		#coreindex = np.corrcoef(coreindex[:-1],y=coreindex[1:])[0][1]
	
	#external stimulation index
	if checkinmethods('external') and checkinmethods('extprop'):
		print "=================================="
		print "===      Spike Probability     ==="
		spprop = 0
		for etx in xrange(methods['external']['count']):
			lidx = int( np.floor(methods['external']['start']+methods['external']['interval']*etx) )
			ridx = int( np.floor(lidx + methods['external']['interval']*methods['extprop']) )
			spprop += float( np.sum(spbins[lidx:ridx]) )
		spprop /= methods['external']['count']*methods["ncell"]
		methods['extprop-results'] = spprop
		print " > Spike group probability      :",spprop
		print "==================================\n"
			

	if checkinmethods("peakDetec") or checkinmethods("R2") or checkinmethods("SAC") or checkinmethods('N2NHI') or checkinmethods('N2NHI-netISI') or checkinmethods('spike2net-dist'):
		print "=================================="
		print "===         Peak Detector      ==="
		print "==================================\n"
		kernel = np.arange(-methods["gkernel"][1],methods["gkernel"][1],1.)
		kernel = np.exp(kernel**2/methods["gkernel"][0]/(-methods["gkernel"][0]))
		module = np.convolve(spbins,kernel)
		module = module[kernel.size/2:1-kernel.size/2]
		#spbinmax = (np.diff(np.sign(np.diff(module))) < 0).nonzero()[0] + 1
		#spbinmin = (np.diff(np.sign(np.diff(module))) > 0).nonzero()[0] + 1
		spbinmark = []
		for idx in (np.diff(np.sign(np.diff(module))) < 0).nonzero()[0] + 1:
			spbinmark.append([idx,1])
		for idx in (np.diff(np.sign(np.diff(module))) > 0).nonzero()[0] + 1:
			spbinmark.append([idx,-1])
		peakmark  = []
		spc,ccnt = 0.,0.
		if checkinmethods("SPC-std"):
			spccun = []
		if len(spbinmark) > 2:
			spbinmark.sort()
			spbinmark = np.array(spbinmark)
			for mx in np.where( spbinmark[:,1] > 0 )[0]:
				if mx <= 2 or mx >= (spbinmark.shape[0] -2):continue
				if spbinmark[mx-1][1] > 0 or spbinmark[mx+1][1] > 0 or spbinmark[mx][1] < 0:continue
				peakmark.append(spbinmark[mx])
				ccnt    += 1
				curespc = np.sum(spbins[spbinmark[mx-1][0]:spbinmark[mx+1][0]])
				if checkinmethods("SPC-std"):spccun.append(curespc)
				spc += curespc
		else:
			spbinmark = None
		if ccnt > 0:
			spc /= ccnt
		
	
	if checkinmethods("SAC") and not spbinmark is None:
		print "=================================="
		print "===    Spikes Autocorrelation  ==="
		#sac_v_max = len( peakmark )- 2
		#sac_vector = np.zeros( (methods['ncell'],sac_v_max) )
		#sac_peaks  = [ (lp[0],mp[0],rp[0]) for lp,mp,rp in zip(peakmark[:-2],peakmark[1:-1],peakmark[2:]) ]
		#for pi,(lp,mp,rp) in enumerate(sac_peaks):
			#Pl = (lp - mp)/2 + mp
			#Pr = (rp - mp)/2 + mp
			#for ni, n in enumerate(neurons):
				#sac_vector[ni,pi] = 0 if len( np.where( (n.spks>Pl)*(n.spks<Pr) )[0] ) == 0 else 1
				##sac_vector[ni,pi] = -1 if len( np.where( (n.spks>Pl)*(n.spks<Pr) )[0] ) == 0 else 1
		#sac_c_max = methods["SAC"] if type(methods["SAC"]) is int else sac_v_max/2-1
		#sac_events = []
		#for sac_ac_idx in range(1,sac_c_max):
			#sac_ac = np.sum(sac_vector[:,:-sac_ac_idx]*sac_vector[:,sac_ac_idx:],axis=1)
			#sac_events.append( (float(np.amax(sac_ac))/float(sac_vector.shape[1]-sac_ac_idx-1),sac_ac_idx,np.where(sac_ac==int(np.amax(sac_ac)) )[0] ) )
		#sac_events.sort()
		#for sac_ev in  sac_events:
			#print sac_ev
			
			
		##DB>>
		#np.savetxt("sac_ac.txt",sac_vector, fmt='%-1d')
		#with open("sac_ev.txt","w") as fd:
			#for sac_ev in  sac_events:
				#fd.write("{}\n".format(sac_ev) )
			
		##exit(0)
		##<<DB 
					
		print "==================================\n"
	if methods['tracetail'] == 'R2':
		costR2p = float( methods['cont-R2'] )
		contR2t = np.arange(0,costR2p,1.)
		sinkern = np.sin(np.pi*2.*contR2t/costR2p)
		coskern = np.cos(np.pi*2.*contR2t/costR2p)
		cntkern = np.ones(contR2t.shape)
		sinkern = np.convolve(spbins,sinkern)
		coskern = np.convolve(spbins,coskern)
		cntkern = np.convolve(spbins,cntkern)
		contR2  = (coskern/cntkern)**2+(sinkern/cntkern)**2
		contR2  = contR2[contR2t.shape[0]/2:1-contR2t.shape[0]/2]
		contSPC = cntkern[contR2t.shape[0]/2:1-contR2t.shape[0]/2]/methods["ncell"]
		if checkinmethods("cont-R2-smooth"):
			if type(methods["cont-R2-smooth"]) is float:
				costR2p = methods["cont-R2-smooth"]
			else:
				costR2p *= 3.
			contR2t = np.arange(0,costR2p*2,1.)
			contR2t = np.exp(-((contR2t-costR2p)/costR2p)**2)
			contR2  = np.convolve(contR2, contR2t)/np.sum(contR2t)
			contR2  = contR2[contR2t.shape[0]/2:1-contR2t.shape[0]/2]
			contSPC = np.convolve(contSPC, contR2t)/np.sum(contR2t)
			contSPC = contSPC[contR2t.shape[0]/2:1-contR2t.shape[0]/2]
		print "=================================="
		print "===        Continues R2        ==="
		print " > Preiod (Frequency)            :",methods['cont-R2'],"(ms) /",1000./methods['cont-R2'],"(Hz)"
		if checkinmethods("cont-R2-smooth"):
			print " > Smooth kernel                 :",costR2p,"(ms)"
		print "==================================\n"
		
		
	if checkinmethods('jitter-rec'):
		print "=================================="
		print "===       Jitter Detector      ==="
		print "==================================\n"
		jkernel = np.arange(-methods["gkernel"][1],methods["gkernel"][1],methods['timestep'])
		jkernel = np.exp(jkernel**2/methods["gkernel"][0]/(-methods["gkernel"][0]))
		jmodule = np.convolve(jallspikes,jkernel)
		jmodule = jmodule[jkernel.size/2:1-jkernel.size/2]
		jpeaks = []
		#for idx in (np.diff(np.sign(np.diff(jmodule))) < 0).nonzero()[0] + 1:
			#jpeaks.append(float(idx))
		#for il,ic,ir in zip(jpeaks[:-2],jpeaks[1:-1],jpeaks[2:]):
			#for ik in jmodule[(il+ic)/2:ic]:
			#???*methods['timestep']

	##R2
	##Per
	if checkinmethods("R2"):
		methods["R2-results"] = {}
		if ccnt > 0:
			methods["R2-results"]['spc'] = spc
		else:
			methods["R2-results"]['spc'] = None
		print "=================================="
		print "===             R2             ==="
		X,Y,Rcnt,netpermean,netpercnt,netfrqmean=0.,0.,0.,0.0,0.0,0.0
		phydist  = []
		if not spbinmark is None:
			for mx in np.where( spbinmark[:,1] > 0 )[0]:
				#if mx >= (spbinmark.shape[0]/2 - 3):continue
				if mx >= (spbinmark.shape[0] - 3):continue
				if spbinmark[mx+1][1] > 0 or spbinmark[mx+2][1] < 0 or spbinmark[mx][1] < 0:continue
				Pnet = float(spbinmark[mx+2][0] - spbinmark[mx][0])
				netpermean += Pnet
				if Pnet > 0.:
					netfrqmean += 1000./Pnet
				netpercnt  += 1.
#				for n,i in map(None,spbins[spbinmark[mx][0]:spbinmark[mx+2][0]],xrange(spbinmark[mx+2][0] - spbinmark[mx][0])):
				for i,n in enumerate(spbins[spbinmark[mx][0]:spbinmark[mx+2][0]]):
					phyX = np.cos(np.pi*2.*float(i)/Pnet)
					phyY = np.sin(np.pi*2.*float(i)/Pnet)
					X += n*phyX
					Y += n*phyY
					Rcnt += n
					if methods['sycleprop']:
						#phydist.append( (360.*np.arctan2(phyY,phyX)/2/np.pi,n) )
						#phydist.append( (np.arctan2(phyY,phyX),n) )
						phydist.append( (np.pi*2.*float(i)/Pnet,n) )
		if Rcnt > 0.:
			R2 = (X/Rcnt)**2+(Y/Rcnt)**2
			methods["R2-results"]["R2"] = R2
		else:
			methods["R2-results"]["R2"] = None
		if netpercnt > 1.:
			netpermean /= ( netpercnt - 1)
			methods["R2-results"]["netPmean"] = netpermean
			netfrqmean /= ( netpercnt - 1)
			methods["R2-results"]["netFmean"] = netfrqmean
		else:
			methods["R2-results"]["netPmean"] = None
			methods["R2-results"]["netFmean"] = None
		if not(methods["R2-results"]["netFmean"] is None or methods["nrnFmean"] is None):
			nrnfr = np.array(nrnfr)
			methods["R2-results"]["mean_Fr/Fnet"] = np.mean(nrnfr/methods["R2-results"]["netFmean"])
			methods["R2-results"]["stdr_Fr/Fnet"] = np.std(nrnfr/methods["R2-results"]["netFmean"])
		else:
			methods["R2-results"]["mean_Fr/Fnet"] = methods["R2-results"]["stdr_Fr/Fnet"] = None
		print "  > R2       =           ",methods["R2-results"]["R2"]
		print "  > SPC      =           ",methods["R2-results"]['spc']
		print "  > netPmean =           ",methods["R2-results"]["netPmean"]
		print "  > netFmean =           ",methods["R2-results"]["netFmean"]
		if not(methods["R2-results"]["netFmean"] is None or methods["nrnFmean"] is None):
			print "  > Fsr/Fnet =           ",methods["R2-results"]["mean_Fr/Fnet"], "+-",methods["R2-results"]["stdr_Fr/Fnet"]
		print "==================================\n"
		if checkinmethods('sycleprop'):
			phydist = np.array(phydist)
			phydist[:,1] /= np.sum(phydist[:,1])
			phyhist,phyhistbins = np.histogram(phydist[:,0], bins=37, weights=phydist[:,1],range=(-np.pi/36,2.*np.pi+np.pi/36))
			methods['sycleprop-results'] = { 'histogram':phyhist, 'bins-bounders':phyhistbins }


	if 10 < methods["netISI"] < 3000 or checkinmethods('N2NHI-netISI'):
		print "=================================="
		print "===          NET ISI           ==="
		print "==================================\n"
		#netisi	= np.zeros(methods["netISI"])
		#lock = threading.RLock()
		#def calcnetisi(ns):
			#global netisi, lock
			#scans	= np.zeros(methods["ncell"],dtype=int)
			#localnetisi = np.zeros(methods["netISI"])
			#for n in ns:
				#for sp in n.spks:
					#for (idx,m) in map(None,xrange(methods["ncell"]),neurons):
						#if m.spks.size < 2 : continue
						#while m.spks[scans[idx]] <= sp and scans[idx] < m.spks.size - 1 : scans[idx] += 1
						#if m.spks[scans[idx]] <= sp : continue
						#if m == n and m.spks[scans[idx]] - sp < 1e-6 : continue
						#aisi = m.spks[scans[idx]] - sp
						#if int(round(aisi)) >= methods["netISI"] : continue
						#localnetisi[ int(round(aisi)) ] += 1
			#with lock:
				#netisi += localnetisi
		#pids = [ threading.Thread(target=calcnetisi, args=(neurons[x::methods['corefunc']],)) for x in xrange(methods['corefunc']) ]
		#for pidx in pids:
			#pidx.start()
			##print pidx, "starts"
		#for pidx in pids:
			#pidx.join()
			##print pidx,"finishs"
		#methods['netISI-results'] = netisi
		netisi   = []
		netsp    = [ np.array(n.spks) for n in neurons ]
		#lock = threading.RLock()
		#def calcnetisi(ns):
			#global netisi, lock
			#localnetisi = [ np.amin(x-netsp[n][np.where(netsp[n] <= x)],initial=100000) for n in xrange(methods['ncell']) for nisi in ns for x in nisi  ]
			#with lock:
				#netisi += localnetisi
		#pids = [ threading.Thread(target=calcnetisi, args=(netsp[x::methods['corefunc']],)) for x in xrange(methods['corefunc']) ]
		#for pidx in pids: pidx.start()
		#for pidx in pids: pidx.join()
		netisi   = np.array([ np.amin(x-netsp[n][np.where(netsp[n] <= x)],initial=100000) for n in xrange(methods['ncell']) for nisi in netsp for x in nisi  ])
		if checkinmethods('N2NHI-netISI'):
			cg_netisi,_ = np.histogram(np.array(netisi), range=(0,300), bins=100)
			cg_netisi   = cg_netisi.astype(float)
			##DB>>
			#print "CG NetISI:",cg_netisi
			#print "CG NetBIN:",_
			##<<DB
		netisi,_ = np.histogram(np.array(netisi), range=(0,methods["netISI"]), bins=int(round(methods["netISI"]/1.) ))
		netisi   = netisi.astype(float)
		methods['netISI-results'] = netisi

	if checkinmethods('N2NHI') or checkinmethods('N2NHI-netISI'):
		def peakdetector(isi,bins):
			#kernel = np.arange(-25,25,1.)
			#kernel = np.exp(kernel**2/(-9))
			#module = np.convolve(isi,kernel)
			#module = module[kernel.size/2:1-kernel.size/2]
			#return [
				#idx for idx in (np.diff(np.sign(np.diff(module))) < 0).nonzero()[0] + 1
			#]
			return [
				idx for idx in (np.diff(np.sign(np.diff(isi))) < 0).nonzero()[0] + 1
			]
		
		cg_nrnisi,x = np.histogram(np.array(cg_nrnisi), range=(0,300), bins=100)
		cg_nrnisi   = cg_nrnisi.astype(float)
		cg_nrnbin   = (x[1:]+x[:-1])/2.
		##DB>>
		#print "CG NrnISI:",cg_nrnisi
		#print "CG NrnBIN:",_
		##<<DB
		if checkinmethods('N2NHI'):
			print "=================================="
			print "===  Carmen's Cluster Indices  ==="
			Pnet = methods["R2-results"]["netPmean"]
			if Pnet is None:
				ccc_clsidx = [ None for clsidx in xrange(10) ]
			else:
				ccc_clsidx = [ int(round(Pnet*(clsidx+1)))/3 for clsidx in xrange(10) if int(round(Pnet*(clsidx+1)))/3+1 < cg_nrnisi.shape[0] ]
				ccc_clsidx = [ sum(cg_nrnisi[clsidx-1:clsidx+1])/sum(cg_nrnisi) for clsidx in ccc_clsidx ]
			methods['N2NHI']=ccc_clsidx
			for clsidx,clsize in enumerate(ccc_clsidx):
				print " > Harmonic %02d                    :"%clsidx,clsize
			print "==================================\n"
		if checkinmethods('N2NHI-netISI'):
			clidx = peakdetector(cg_netisi,None)
			print "=================================="
			print "===    RTH's Cluster Indices   ==="
			rth_clsidx = [ cg_nrnisi[xidx]/sum(cg_nrnisi) for xidx in clidx if xidx < cg_nrnisi.shape[0] ]
			methods['N2NHI-netISI'] = rth_clsidx
			for clsidx,clsize in enumerate(rth_clsidx):
				print " > Harmonic %02d                    :"%clsidx,clsize
			print "==================================\n"

		#methods["netISI"] = 300
		#methods["nrnISI"] = 300
	#if 
		#methods["nrnISI"] = 300
		
	if checkinmethods('spike2net-dist'):
		if not spbinmark is None:
			piskd = np.zeros(200)
			for mx in np.where( spbinmark[:,1] > 0 )[0]:
				if mx >= (spbinmark.shape[0] - 3):continue
				#if spbinmark[mx+1][1] > 0 or spbinmark[mx+2][1] < 0 or spbinmark[mx][1] < 0:continue
				LPnet = float(spbinmark[mx  ][0] - spbinmark[mx-2][0])
				RPnet = float(spbinmark[mx+2][0] - spbinmark[mx  ][0])
				for i,n in enumerate( spbins[spbinmark[mx-2][0]:spbinmark[mx][0]] ):
					binID = int( round( float(i)*100./LPnet ) )
					piskd[binID] +=  n
				for i,n in enumerate( spbins[spbinmark[mx][0]:spbinmark[mx+2][0]] ):
					binID = int( round( float(i)*100./RPnet+100 ) )
					if binID>=200: continue
					piskd[binID] +=  n
			methods['spike2net-dist-result'] = piskd.tolist()
		else:
			methods['spike2net-dist-result'] = None
		
	if checkinmethods("T&S") or checkinmethods('lastspktrg'):
		print "=================================="
		print "===           T & S            ==="
		print "==================================\n"
		allspikes,activeneurons = [],0.
		for n in neurons:
			allspikes += list(n.spks)
			if n.spks.size !=0 :activeneurons += 1.
		allspikes.sort()
		allspikes = np.array(allspikes)
		TaSisi = allspikes[1:]-allspikes[:-1]
		if checkinmethods('lastspktrg'):
			lastspktrg = int( np.mean(allspikes) > methods['tstop']/4. )
			methods['lastspktrg-results'] = lastspktrg
		
		del allspikes
		if checkinmethods("T&S"):
			if bool(lastspktrg):
				mean1TaSisi = np.mean(TaSisi)
				TaSindex	= (np.sqrt(np.mean(TaSisi**2) - mean1TaSisi**2)/mean1TaSisi - 1.)/np.sqrt(activeneurons) 
				methods['T&S-results'] = TaSindex
			else:
				methods['T&S-results'] = None
		
	if checkinmethods("Delay-stat"):
		print "=================================="
		print "===     Delays distribution    ==="
		delays = np.array([ x[0].delay for x in connections])
		mdly, sdly,mxdly,Mxdly = np.mean(delays), np.std(delays), np.min(delays), np.max(delays)
		if not type(methods["Delay-stat"]) is tuple:
			methods["Delay-stat"] = (0., 15., 500)
		dlyhist,dlybins = np.histogram(delays, bins=methods["Delay-stat"][2], normed=True, range=methods["Delay-stat"][0:2] )
		dlyhist /= np.sum(dlyhist)
		print "  > Delays mean  =           ",mdly
		print "  > Delays stdev =           ",sdly
		print "  > Delays CV    =           ",sdly/mdly
		print "  > Delays min   =           ",mxdly
		print "  > Delays max   =           ",Mxdly
		methods['Delay-stat-results'] = {
			'mean': mdly, 'stdev': sdly, 'min': mxdly, 'max': Mxdly, 'histogram':dlyhist, 'bins-bounders':dlybins
		}
		print "==================================\n"
		
	if checkinmethods('Gtot-dist') :
		print "=================================="
		print "===   G-total  distribution    ==="
		#gsk = [ n.gsynscale for n in neurons ]
		gsk = [ n.gtotal for n in neurons ]
		mgto, sgto,mxgto,Mxgto = np.mean(gsk), np.std(gsk), np.min(gsk), np.max(gsk)
		if not type(methods["Gtot-dist"]) is tuple:
			methods["Gtot-dist"] = (0,0.06,200)
		#gskhist,gskbins = np.histogram(gsk, bins=methods["ncell"]/25, normed=True, range=[0,Mxgto])#/10)#, normed=True)
		gskhist,gskbins = np.histogram(gsk, bins=methods["Gtot-dist"][2], normed=True, range=methods["Gtot-dist"][0:2] )#/10)#, normed=True)
		gskhist /= np.sum(gskhist)
		print "  > Total Syn. Cond mean  =  ",mgto
		print "  > Total Syn. Cond stdev =  ",sgto
		print "  > Total Syn. Cond CV    =  ",sgto/mgto
		print "  > Total Syn. Cond min   =  ",mxgto
		print "  > Total Syn. Cond max   =  ",Mxgto
		methods['Gtot-dist-results'] = { 
			'mean': mgto, 'stdev': sgto, 'min': mxgto, 'max': Mxgto,
			'histogram':gskhist, 'bins-bounders':gskbins 
		}
		if checkinmethods('Gtot-rec') :
			methods['Gtot-rec'] = gsk
		else:
			del gsk
		print "==================================\n"
		
	if checkinmethods('Gtot-stat'):
		print "=================================="
		print "===     G-total Statistics     ==="
		agtot = np.array([ n.gtotal/n.concnt for n in neurons ])
		#DB>>
		#print [n.gtotal for n in neurons ]
		#print [n.concnt for n in neurons ]
		#exit(0)
		#<<DB
		mgtot = np.mean(agtot)
		sgtot = np.std(agtot)
		print "  > mean   gtotal norm    =  ",mgtot
		print "  > stderr gtotal norm    =  ",sgtot
		print "  > CV     gtotal norm    =  ",sgtot/mgtot
		print "==================================\n"
		methods['Gtot-stat-results'] = { 'mean':mgtot, 'stdev':sgtot, 'CV':sgtot/mgtot }

	
	if checkinmethods('2cintercon'):
		print "=================================="
		print "===  2 clusters connectivity   ==="
		tims = methods['tstop']*3./4.
		if pcnt != 0:
			halfpnet = pmean/pcnt/2.
			clslst = []
			print "  >  Searching for clusters       "
			rarr = np.array(neurons[0].spks)
			tims = rarr[ np.where( rarr > tims ) ]
			del rarr
			tims = tims[0] + halfpnet/2.
			print "  >  Time to search              :",tims
			for idx, n in enumerate(neurons):
				getlest = np.array(n.spks)
				getlest = getlest[ np.where( getlest>tims )]
				if getlest.shape[0] < 1: continue
				if getlest[0] > tims+halfpnet:
					clslst.append(True)
				else:
					clslst.append(False)
			print "  >  Searching for connectivity index"
			WithinA, WithinB, CrossAB, CrossBA = 0, 0, 0, 0
			countA, countB = 0, 0
			fullstat = checkinmethods('2clrs-stat')
			if fullstat:
				within,cross = np.zeros(methods['ncell']),np.zeros(methods['ncell'])
			for idx, cnt in enumerate(OUTList):
				if clslst[idx]:
					#Cluster A
					countA += 1
					for c in cnt:
						if clslst[c]:
							WithinA += 1
							if fullstat : within[idx] += 1
						else:
							CrossAB += 1
							if fullstat : cross[idx] += 1
				else:
					#cluster B
					countB += 1
					for c in cnt:
						if clslst[c]:
							CrossBA += 1
							if fullstat : cross[idx] += 1
						else:
							WithinB += 1
							if fullstat : within[idx] += 1
			print "  >  Cells in the Cluster A      :",countA
			print "  >  Cells in the Cluster B      :",countB
			print "  >  Within Cluster A            :",WithinA
			print "  >  Within Cluster B            :",WithinB
			print "  >  From A to B                 :",CrossAB
			print "  >  From A to B                 :",CrossBA
			print "  >  Total Within Both Clusters  :",WithinA + WithinB
			print "  >  Total Between Both Clusters :",CrossAB + CrossBA
			print "  >  Ratio Between to Within     :",float(CrossAB + CrossBA)/float(WithinA + WithinB)
			if fullstat:
				print "  >  FUUL STATISTICS "
				#print "  >  ",
				#for idx,(win,btwn) in enumerate(zip(within,cross)):
					#print "{}:{}".format(win,btwn),
					#if not bool((idx+1)%6): print "\n  >  ",
				#print "  >  RATIOS "
				#print "  >  ",
				#for idx,(win,btwn) in enumerate(zip(within,cross)):
					#print "{}".format(btwn/win),
					#if not bool((idx+1)%6): print "\n  >  ",
				#print "  >  "
				withinA = within[ np.where( np.array(clslst) ) ]
				crossA = cross[ np.where( np.array(clslst) ) ]
				print "  >  Cluster A within mean,stdev :",np.mean(withinA), np.std(withinA)
				print "  >  Cluster A to B mean,stdev   :",np.mean(crossA), np.std(crossA)
				withinB = within[ np.where(  (1-1*np.array(clslst).astype(int)).astype(bool) )]
				crossB = cross[ np.where(  (1-1*np.array(clslst).astype(int)).astype(bool) )]
				print "  >  Cluster B within mean,stdev :",np.mean(withinB), np.std(withinB)
				print "  >  Cluster B to A mean,stdev   :",np.mean(crossB), np.std(crossB)
			methods['2cintercon-results'] = {
				"cells-in-A":countA, "cells-in-B":countB,
				'connections-in-A':WithinA, 'connections-in-B':WithinB,
				"connections-A2B":CrossAB, "connections-B2A":CrossBA,
				'total in A&B':WithinA + WithinB,
				'total between A&B':CrossAB + CrossBA,
				'total ratio between/in':float(CrossAB + CrossBA)/float(WithinA + WithinB)
			}
		else:
			print "  >  Pnet isn't defined...,       "
			methods['2cintercon-results'] = None
		print "==================================\n"
	
	if checkinmethods('CtrlISI'):
		print "=================================="
		print "=== Controlled ISI calculation ==="
		if type(methods['CtrlISI']) is not dict:
			methods['CtrlISI'] = {'bin'   : 5.,'max'   : 120.,}
		if not checkinmethods('CtrlISI/bin'): methods['CtrlISI']['bin'] =   5.
		if not checkinmethods('CtrlISI/max'): methods['CtrlISI']['max'] = 120.
		xbin,xmax = methods['CtrlISI']['bin'], methods['CtrlISI']['max']
		CtrINT = np.linspace(0,xmax,xmax/xbin)+xbin/2
		CtrISI = np.zeros(int(round(xmax/xbin)))
		for n in neurons:
			for isi in n.spks[1:]-n.spks[:-1]:
				isiid = int( round(isi/xbin) )
				if isiid < CtrISI.shape[0]:
					CtrISI[isiid] += 1
		CtrISI = np.column_stack((CtrINT,CtrISI))
		methods['CtrlISI']['result'] = CtrISI
	print "==================================\n"
	
	if checkinmethods('T1vsT2/spikerate'):
		T1cnt,T2cnt = 0, 0
		T1spk,T2spk = 0, 0
		for n in neurons:
			if n.type21 == 1:
				T1cnt += 1
				T1spk += n.spks.shape[0]
			elif n.type21 == 2:
				T2cnt += 1
				T2spk += n.spks.shape[0]
		methods['T1vsT2']['spikerate'] = {
			'T1' : float(T1spk)/float(T1cnt)/methods['tstop']*1000. if T1cnt != 0 else 0.,
			'T2' : float(T2spk)/float(T2cnt)/methods['tstop']*1000. if T2cnt != 0 else 0.
		}
		print "=================================="
		print "===   Type I vs. Type2 Rate    ==="
		print " > Type I  rate                 :",methods['T1vsT2']['spikerate']['T1'],"[Hz]"
		print " > Type II rate                 :",methods['T1vsT2']['spikerate']['T2'],"[Hz]"
		print "==================================\n"

	if checkinmethods('get-steadystate'):
		if type(methods['get-steadystate']) is str:
			ssthr=+30.
			ssfilename = methods['get-steadystate']
		elif type(methods['get-steadystate']) is float or type(methods['get-steadystate']) is int:
			ssthr = float(methods['get-steadystate'])
			if type(methods["neuron"]["Vinit"]) is str:
				ssfilename =  methods["neuron"]["Vinit"]+"-ss.dat"
			else:
				ssfilename = 'get-steadystate.dat'
		else:
			ssthr=+30.
			ssfilename = 'get-steadystate.dat'
		print "=================================="
		print "===     Write Steady State     ==="
		print "  >  Threshold                   :",ssthr
		print "  >  Output File                 :",ssfilename
		ssvec = np.array(neurons[0].volt)
		ssmsk = np.where( (t>methods['tstop']*4./5.)*(ssvec >= ssthr) )[0]
		if ssmsk.shape[0] == 0:
			print "Error Cannot Get Voltage above Threshold!"
			with open(ssfilename,"w") as fd:
				fd.write("None")
				for n in neurons[1:]:
					fd.write(" None")
				fd.write("\n")
		else:
			ssmsk = int(ssmsk[0])
			with open(ssfilename,"w") as fd:
				fd.write("%g"%ssvec[ssmsk])
				for n in neurons[1:]:
					fd.write(" %g"%(np.array(n.volt)[ssmsk]))
				fd.write("\n")
		print "==================================\n"
	#EN
	#p.set_title("Mean individual Period = %g, Sychrony(Entropy) = %g(%g)"%(pmean/pcnt,1./(1.+en),en))

	if methods['tracetail'] == 'p2eLFP' or checkinmethods('PAC-Canolty06') or checkinmethods('PAC-Tort10') or checkinmethods('PAC-VS'):
		print "=================================="
		print "===       Calculating LFP      ==="
		lfp=np.zeros(module.shape[0])
		x1,x2=0.,0,
		for i,s in enumerate(module):
			x1 = s+x1-x1/2.
			x2 = s+x2-x2/5.
			lfp[i]=x2-x1
		#DB>>
		#with open("lfp.csv","w") as fd:
			#for x1,x2 in enumerate(lfp):
				#fd.write("{},{}\n".format(float(x1),x2) )
		#exit(0)
		#<<DB
		if checkinmethods("p2eLFP/BPF"):
			print "===       Band-filter LFP      ==="
			from scipy.signal import butter, lfilter, freqz
			nyq = 0.5 * 1000.				#SAMPLING EVERY 1 ms
			if type(methods["p2eLFP"]["BPF"]) is float or type(methods["p2eLFP"]["BPF"]) is int:
				normal_cutoff = float(methods["p2eLFP"]["BPF"]) / nyq		#lowpass 100Hz
				b, a = butter(5, normal_cutoff, btype='low', analog=False)
				lfp = lfilter(b, a, lfp)
				print " > Filtered low-pass            : ",float(methods["p2eLFP"]["BPF"]),"Hz"
			elif (type(methods["p2eLFP"]["BPF"]) is tuple or type(methods["p2eLFP"]["BPF"]) is list) and len(methods["p2eLFP"]["BPF"]) == 2:
				normal_cutoff = float(methods["p2eLFP"]["BPF"][0]) / nyq	#lowpass 100Hz
				b, a = butter(5, normal_cutoff, btype='low', analog=False)
				normal_cuton  = float(methods["p2eLFP"]["BPF"][1]) / nyq	#highpass 25Hz
				c, d = butter(5, normal_cuton , btype='high', analog=False)
				lfp = lfilter(b, a, lfilter(c, d, lfp) )
				print " > Filtered  low-pass             : ",float(methods["p2eLFP"]["BPF"][0]),"Hz"
				print " > Filtered high-pass             : ",float(methods["p2eLFP"]["BPF"][1]),"Hz"
		print "==================================\n"
	
	if checkinmethods('PAC-Canolty06'):
		from scipy.signal import butter, lfilter, freqz
		from scipy.signal import hilbert
		from scipy.stats import norm
		print "=================================="
		print "===       Calculating PAC      ==="
		print "===      Canolty et al 2006    ==="
		if type(methods['PAC-Canolty06']) is not dict: methods['PAC-Canolty06']={}
		if not 'BPF' in methods['PAC-Canolty06']:
			methods['PAC-Canolty06']['BPF'] = 100.,25. # 25 - 100 Hz Gamma range
		nyq = 0.5 * 1000.				#SAMPLING EVERY 1 ms
		normal_cutoff = float(methods["PAC-Canolty06"]["BPF"][0]) / nyq	#lowpass 100Hz
		normal_cuton  = float(methods["PAC-Canolty06"]["BPF"][1]) / nyq	#highpass 25Hz
		b, a = butter(5, np.array([normal_cuton,normal_cutoff]), btype='bandpass', analog=False)
		vlfp = lfilter(b, a, lfp) 
		print " > Filtered  low-pass             : ",float(methods["PAC-Canolty06"]["BPF"][0]),"Hz"
		print " > Filtered high-pass             : ",float(methods["PAC-Canolty06"]["BPF"][1]),"Hz"
		singen = neurons[0].sin if checkinmethods('/sinmod') else neurons[0].sing
		modstart,modstop,modper =\
			singen.tstart,singen.tstop,singen.per
		modtime   = np.arange(lfp.shape[0]) + methods['cliptrn']
		#modsignal = - (1.-np.cos(np.pi*2./modper*(modtime-modstart)))/2.
		#modsignal[np.where(modtime<modstart)] = 0.
		#modsignal[np.where(modtime>modstop )] = 0.
		#phase		= np.angle( hilbert(modsignal) )
		phase       = (np.pi*2./modper*(modtime-modstart))%(np.pi*2.)
		amplitude	= np.abs(  hilbert( vlfp )  )
		z           = amplitude*np.exp(1.j*phase)
		m_raw       = np.mean(z)
		surrogate_m = np.array([\
			np.abs(\
				np.mean(\
					np.hstack((amplitude[s:],amplitude[:s]))*np.exp(1.j*phase)\
				)\
			)\
			for s in np.random.randint(amplitude.shape[0],size=200) ])
		surrogate_mean,surrogate_std = norm.fit(surrogate_m)
		m_norm_length = (np.abs(m_raw)-surrogate_mean)/surrogate_std;
		m_norm_phase  = np.angle(m_raw);
		#m_norm        = m_norm_length*np.exp(1j*m_norm_phase);
		methods['PAC-Canolty06']['length'] = np.abs(m_norm_length),
		methods['PAC-Canolty06']['phase' ] = m_norm_phase
		print " > Normalized score length        : ",np.abs(m_norm_length)
		print " > Normalized score phase         : ",m_norm_phase
		print "==================================\n"
		#MATLAB code>>
		#From Canolty et al 2006 - supporting material
		#	srate=2003;                %% sampling rate used in this study, in Hz 
		#	numpoints=length(x);       %% number of sample points in raw signal
		#	numsurrogate=200;          %% number of surrogate values to compare to actual value 
		#	minskip=srate;             %% time lag must be at least this big 
		#	maxskip=numpoints-srate;   %% time lag must be smaller than this
		#	skip=ceil(numpoints.*rand(numsurrogate*2,1)); 
		#	skip(find(skip>maxskip))=[]; 
		#	skip(find(skip<minskip))=[]; 
		#	skip=skip(1:numsurrogate,1);
		#	surrogate_m=zeros(numsurrogate,1);
		#	%% HG analytic amplitude time series, uses eegfilt.m from EEGLAB toolbox
		#	%% (http://www.sccn.ucsd.edu/eeglab/)
		#	amplitude=abs(hilbert(eegfilt(x,srate,80,150)));
		#	%% theta analytic phase time series, uses EEGLAB toolbox
		#	phase=angle(hilbert(eegfilt(x,srate,4,8)));
		#	%% complex-valued composite signal
		#	z=amplitude.*exp(i*phase); 
		#	%% mean of z over time, prenormalized value 
		#	m_raw=mean(z);  
		#	%% compute surrogate values
		#	for s=1:numsurrogate
		#		surrogate_amplitude=[amplitude(skip(s):end) amplitude(1:skip(s)-1)];
		#		surrogate_m(s)=abs(mean(surrogate_amplitude.*exp(i*phase)));
		#		disp(numsurrogate-s)
		#	end 
		#	%% fit gaussian to surrogate data, uses normfit.m from MATLAB Statistics toolbox [
		#	surrogate_mean,surrogate_std]=normfit(surrogate_m);
		#	%% normalize length using surrogate data (z-score)
		#	m_norm_length=(abs(m_raw)-surrogate_mean)/surrogate_std;
		#	m_norm_phase=angle(m_raw);
		#	m_norm=m_norm_length*exp(i*m_norm_phase); 
		#<<<<<<<<
		#<<DB
	if checkinmethods('PAC-Tort10'):
		from scipy.signal import butter, lfilter, freqz
		from scipy.signal import hilbert
		from scipy.stats import norm
		print "=================================="
		print "===     PAC-modulation depth   ==="
		print "===        Tort et al 2010     ==="
		if type(methods['PAC-Tort10']) is not dict: methods['PAC-Tort10']={}
		if not 'BPF' in methods['PAC-Tort10']:
			methods['PAC-Tort10']['BPF'] = 100.,25. # 25 - 100 Hz Gamma range
		nyq = 0.5 * 1000.				#SAMPLING EVERY 1 ms
		normal_cutoff = float(methods["PAC-Tort10"]["BPF"][0]) / nyq	#lowpass 100Hz
		normal_cuton  = float(methods["PAC-Tort10"]["BPF"][1]) / nyq	#highpass 25Hz
		b, a = butter(5, np.array([normal_cuton,normal_cutoff]), btype='bandpass', analog=False)
		vlfp = lfilter(b, a, lfp) 
		print " > Filtered  low-pass             : ",float(methods["PAC-Tort10"]["BPF"][0]),"Hz"
		print " > Filtered high-pass             : ",float(methods["PAC-Tort10"]["BPF"][1]),"Hz"
		singen = neurons[0].sin if checkinmethods('/sinmod') else neurons[0].sing
		modstart,modstop,modper =\
			singen.tstart,singen.tstop,singen.per
		modtime   = np.arange(lfp.shape[0]) + methods['cliptrn']
		#modsignal = np.sin(np.pi*2./modper*(modtime-modstart))
		#modsignal[np.where(modtime<modstart)] = 0.
		#modsignal[np.where(modtime>modstop )] = 0.
		#phase		= np.angle( hilbert(modsignal) )
		phase       = (np.pi*2./modper*(modtime-modstart))%(np.pi*2.)
		amplitude	= np.abs(  hilbert( vlfp )  )
		#z           = amplitude*np.exp(1.j*phase)
		#m           = np.abs( np.mean(z) )
		
		hist,bed = np.histogram(phase, weights=amplitude, bins=36)
		hist = hist.astype(float)/np.sum( hist.astype(float) )
		bct = (bed[:-1]+bed[1:])/2.
		methods['PAC-Tort10']["MI"] = (np.log2(float(hist.shape[0])) + np.sum( hist* np.log2(hist) ) )/np.log2(float(hist.shape[0]))
		print " > Modulation index               : ",methods['PAC-Tort10']["MI"]
		print "==================================\n"
		##DB>>
		#plt.subplot(nplot+3,sharex=p)
		#plt.plot(modtime,vlfp[tproc:tprin.size+tproc])
		#plt.show()
		#exit(0)
		##<<DB
	if checkinmethods('PAC-VS'):
		from scipy.signal import butter, lfilter, freqz
		from scipy.signal import hilbert
		from scipy.stats import norm
		print "=================================="
		print "===    PAC - Vector Strength   ==="
		#print "===        Tort et al 2010     ==="
		if type(methods['PAC-VS']) is not dict: methods['PAC-VS']={}
		if not 'BPF' in methods['PAC-VS']:
			methods['PAC-VS']['BPF'] = 100.,25. # 25 - 100 Hz Gamma range
		nyq = 0.5 * 1000.				#SAMPLING EVERY 1 ms
		normal_cutoff = float(methods["PAC-VS"]["BPF"][0]) / nyq	#lowpass 100Hz
		normal_cuton  = float(methods["PAC-VS"]["BPF"][1]) / nyq	#highpass 25Hz
		b, a = butter(5, np.array([normal_cuton,normal_cutoff]), btype='bandpass', analog=False)
		vlfp = lfilter(b, a, lfp) 
		print " > Filtered  low-pass             : ",float(methods["PAC-VS"]["BPF"][0]),"Hz"
		print " > Filtered high-pass             : ",float(methods["PAC-VS"]["BPF"][1]),"Hz"
		singen = neurons[0].sin if checkinmethods('/sinmod') else neurons[0].sing
		modstart,modstop,modper =\
			singen.tstart,singen.tstop,singen.per
		modtime   = np.arange(lfp.shape[0]) + methods['cliptrn']
		#modsignal = np.sin(np.pi*2./modper*(modtime-modstart))
		#modsignal[np.where(modtime<modstart)] = 0.
		#modsignal[np.where(modtime>modstop )] = 0.
		#phase		= np.angle( hilbert(modsignal) )
		phase       = (np.pi*2./modper*(modtime-modstart))%(np.pi*2.)
		amplitude	= np.abs(  hilbert( vlfp )  )
		z           = np.mean( amplitude*np.exp(1.j*phase) )
		m           = np.abs( z )
		methods['PAC-VS']['length'] = np.abs( z )
		methods['PAC-VS']['phase']  = np.angle( z )
		print " > Vector length                  : ",methods['PAC-VS']['length']
		print " > Vector  phase                  : ",methods['PAC-VS']['phase']
		print "==================================\n"


	##R2
	if checkinmethods('gui'):
		title = methods["MainFigTitle"] if checkinmethods("MainFigTitle") else ""
		title += "Mean individual Period = %s"%("NONE" if pcnt == 0 else "%g"%(pmean/pcnt))
		if checkinmethods('popfr'):
			title += 'Mean FR =%g'%popfr
		if checkinmethods("R2"):
			if Rcnt > 0 :
				title += r", $R^2$ = %g, Mean network Period = %g, Spike per cycle = %g"%(R2,netpermean,spc)
			else:
				title += ", *Fail to estimate network period*"
		elif checkinmethods("peakDetec"):
			title += ", Spike per cycle = %g. "%(spc)
		if checkinmethods('T&S'):
			title += ", TaS = %g"%TaSindex
		if checkinmethods('lastspktrg'):
			title += ", LST = %g"%lastspktrg
		p.set_title(title)

		
		plt.subplot(nplot+2,sharex=p)
		if checkinmethods('cliptrn'):
			nppoints = np.arange(methods['tv'][0]+methods['cliptrn'],methods['tv'][1],1.0)
			plt.bar(nppoints,spbins[:methods['tv'][1]-methods['cliptrn']],0.5,color="k")
			hight = spbins[:methods['tv'][1]-methods['cliptrn']].max()
			if (checkinmethods("peakDetec") or checkinmethods("R2")) and not spbinmark is None :
				for mark in spbinmark:
					if mark[0]+methods['cliptrn'] < methods['tv'][0] or mark[0]+methods['cliptrn'] > methods['tv'][1]: continue
					if mark[1] > 0:
						plt.plot([mark[0]+methods['cliptrn'],mark[0]+methods['cliptrn']],[0,hight],"r--")
					else:
						plt.plot([mark[0]+methods['cliptrn'],mark[0]+methods['cliptrn']],[0,hight],"b--")
		else:
			nppoints = np.arange(methods['tv'][0],methods['tv'][1],1.0)
			plt.bar(nppoints,spbins[int(methods['tv'][0]):int(methods['tv'][1])],0.5,color="k")
			hight = spbins[int(methods['tv'][0]):int(methods['tv'][1])].max()
			if (checkinmethods("peakDetec") or checkinmethods("R2")) and not spbinmark is None :
				for mark in spbinmark:
					if mark[0] < methods['tv'][0] or mark[0] > methods['tv'][1]: continue
					if mark[1] > 0:
						plt.plot([mark[0],mark[0]],[0,hight],"r--")
					else:
						plt.plot([mark[0],mark[0]],[0,hight],"b--")
#			plt.plot(nppoints,module[methods['tv'][0]:methods['tv'][1]]/np.sum(kernel),"k--")
#			plt.plot(nppoints,module[methods['tv'][0]:methods['tv'][1]],"k--")
			##DB>>
			#print peakmark
			#plt.plot(np.array(peakmark)[:,0],np.array(peakmark)[:,1],"k^",ms=20)
			##<<DB
		plt.ylabel("Rate (ms$^{-1}$)", fontsize=16)

	if checkinmethods('gui'):
		if methods['tracetail'] == 'mean total current' or methods['tracetail'] == 'MTI':
			meancur = meancur / float(-methods["ncell"])
		elif methods['tracetail'] == 'mean total synaptic current' or methods['tracetail'] == 'MTSI':
			meancur = meancur / float(-methods["ncell"])
		elif methods['tracetail'] == 'mean total conductance' or methods['tracetail'] == 'MTG':
			meancur = meancur / float(methods["ncell"])
	
	if checkinmethods('gui'):
		
		plt.subplot(nplot+3,sharex=p)
			
		if methods['tracetail'] == "R2":
			plt.ylabel(r"$R^2$/SPC", fontsize=16)
			xvcrv, = plt.plot(np.arange(0,float(contR2.shape[0])),contR2,'k-',lw=2)
			xvcrv, = plt.plot(np.arange(0,float(contSPC.shape[0])),contSPC,"r--",lw=2)
			plt.ylim(0.,1.)
			#DB>>
			#print contR2
			#print contSPC
			#<<DB
		elif methods['tracetail'] == 'total current' or methods['tracetail'] == 'TI' or methods['tracetail'] == 'mean total current' or methods['tracetail'] == 'MTI'\
		  or methods['tracetail'] == 'total synaptic current' or methods['tracetail'] == 'TSI' or methods['tracetail'] == 'mean total synaptic current' or methods['tracetail'] == 'MTSI':
			plt.ylabel("Current (nA)", fontsize=16)
			plt.plot(tprin,-meancur[tproc:tprin.size+tproc])
			plt.plot([tprin[0],tprin[-1]],[0.,0.],"k--")
			#plt.plot([tprin[0],tprin[-1]],[-methods["neuron"]["Iapp"],-methods["neuron"]["Iapp"]],"r--")
			#print np.amax(meancur[tproc+tprin.size/2:tprin.size+tproc]),methods["neuron"]["Iapp"]
		elif methods['tracetail'] == 'total conductance' or methods['tracetail'] == 'mean total conductance' or methods['tracetail'] == 'TG' or methods['tracetail'] == 'MTG':
			plt.ylabel("Total Conductance (nS)", fontsize=16)
			plt.plot(tprin,meancur[tproc:tprin.size+tproc]*1e5)
		elif methods['tracetail'] == 'firing rate' and ( methods["peakDetec"] or methods["R2"] ):
			plt.ylabel("Firing Rate (ms$^{-1}$)", fontsize=16)
			tvl,tvr = int( round(methods['tv'][0]) ), int( round(methods['tv'][1]) )
			plt.plot(nppoints,module[tvl:tvr]/np.sum(kernel),"k--")
			hight = np.max(module[tvl:tvr]/np.sum(kernel))
			if not spbinmark is None :
				for mark in spbinmark:
					if mark[0] < methods['tv'][0] or mark[0] > methods['tv'][1]: continue
					if mark[1] > 0:
						plt.plot([mark[0],mark[0]],[0,hight],"k--")
		elif methods['tracetail'] == 'conductance':
			plt.ylabel("Conductance (nS)", fontsize=16)
			xvcrv, = plt.plot(tprin,np.array(neurons[vindex].isyng)[tproc:tprin.size+tproc]*1e5,'k-',lw=2)
		elif methods['tracetail'] == 'current':
			plt.ylabel(r"Current ($\mu$A)", fontsize=16)
			xvcrv, = plt.plot(tprin,np.array(neurons[vindex].isyni)[tproc:tproc+tprin.size]*1e5,'k-',lw=2)
		elif methods['tracetail'] == 'LFP':
			plt.ylabel(r"LFP", fontsize=16)
			xvcrv, = plt.plot(np.arange(0,float(module.shape[0])),module,'k-',lw=2)
		elif methods['tracetail'] == 'p2eLFP':
			vlfp = lfp[int(round(methods["tv"][0])):int(round(methods["tv"][1]))]
			plt.plot(np.arange(vlfp.shape[0])+methods["tv"][0],vlfp,'k-',lw=2)
			if not checkinmethods("p2eLFP/BPF"):
				plt.ylim(ymin=0)
			if checkinmethods("p2eLFP_max"):
				plt.ylim(ymax=methods["p2eLFP_max"])
		
		if checkinmethods('simvar'):
			simvarrec = np.array(simvarrec)
			plt.subplot(nplot+4,sharex=p)
			#plt.ylabel(methods['simvar']['var'], fontsize=16)
			plt.ylim(min([methods['simvar']["a0"],methods['simvar']["a1"]]),max([methods['simvar']["a0"],methods['simvar']["a1"]]))
			plt.plot(tprin,simvarrec[tproc:tprin.size+tproc])
		if checkinmethods('spectrogram'):
			plt.subplot(nplot+(5 if checkinmethods('simvar') else 4),sharex=p)
			populationcurrent
			#NFFT = 131072       # the length of the windowing segments
			NFFT = 65535       # the length of the windowing segments
			Fs = int(1000.0/methods['timestep'])  # the sampling frequency
			from scipy.signal import spectrogram
			f, tf, Sxx = spectrogram(populationcurrent, fs=Fs, nperseg=NFFT,noverlap=NFFT*1020/1024,window='hanning')
			Sxx = Sxx[np.where(f<100),:][0]
			f   = f[np.where(f<100)]
			Sxx = Sxx[np.where(f>20),:][0]
			f   = f[np.where(f>20)]
			print "T SHAPE",tf.shape
			print "T      ",tf
			print "F SHAPE",f.shape
			print "F      ",f
			print "S SHAPE",Sxx.shape
			print "s      ",Sxx
			plt.pcolormesh(tf*1e3, f, Sxx)

		plt.xlabel("time (ms)", fontsize=16)


	
	if (checkinmethods("netFFT") or checkinmethods("nrnFFT")) and checkinmethods('gui'):
		plt.figure(2)
		if checkinmethods("netFFT") and checkinmethods("nrnFFT"):
			pl=plt.subplot(211)
		elif checkinmethods("netFFT"):
			pl=plt.subplot(111)
		if checkinmethods("netFFT"):
			plt.title( (methods["MainFigTitle"] if checkinmethods("MainFigTitle") else "")+ "Network spectrum")
			plt.bar(specX[1:],fft[1:pnum],0.75,color="k",edgecolor="k")
		if checkinmethods("netFFT") and checkinmethods("nrnFFT"):
			plt.subplot(212,sharex=pl)
		elif checkinmethods("nrnFFT"):
			plt.subplot(111)
		if checkinmethods("nrnFFT"):
			plt.title((methods["MainFigTitle"] if checkinmethods("MainFigTitle") else "")+"Neuronal spectrum")
			plt.bar(specX[1:],specN[1:],0.75,color="k",edgecolor="k")

	#plt.subplot(313,sharex=p)
	#specX =np.arange(0.0,methods['tstop']+h.dt,h.dt)
	#specX *= 1000.0/methods['tstop']/h.dt
	#pnum = specX.size/2
	#plt.title("Voltage spectrum")
	#plt.plot(specX[1:pnum],specV[1:pnum])
	#plt.xlim(0,200)
	
	if 10 < methods["netISI"] <= 3000 and sum(netisi) > 0: netisi /= sum(netisi)
	if 10 < methods["nrnISI"] <= 3000 and sum(isi) > 0: isi /= sum(isi)
	if (10 < methods["netISI"] <= 3000 or 10 < methods["nrnISI"] <= 3000) and methods['gui']:
		plt.figure(3)
		if 10 < methods["netISI"] <= 3000 and 10 < methods["nrnISI"] <= 3000:
			pl=plt.subplot(211)
		elif 10 < methods["netISI"] <= 3000 :
			plt.subplot(111)
		if 10 < methods["netISI"] <= 3000: 
			plt.title((methods["MainFigTitle"] if checkinmethods("MainFigTitle") else "")+"Network ISI")
			plt.ylabel("Normalized number of interspike intervals", fontsize=16)
			plt.bar(np.arange(methods["netISI"]),netisi,0.75,color='k')
			plt.xlim(0,methods["netISI"])
		if 10 < methods["netISI"] <= 3000 and 10 < methods["nrnISI"] <= 3000:
			plt.subplot(212)#,sharex=pl)
		elif 10 < methods["nrnISI"] <= 3000:
			plt.subplot(111)
		if 10 < methods["nrnISI"] <= 3000:
			plt.ylabel("Normalized number of interspike intervals", fontsize=16)
			plt.title((methods["MainFigTitle"] if checkinmethods("MainFigTitle") else "")+"Neuronal ISI")
			plt.bar(np.arange(methods["nrnISI"]),isi,0.75,color='k')
			plt.xlim(0,methods["nrnISI"])
			plt.xlabel("Interspike intervals (ms)", fontsize=16)
	
	if ( checkinmethods('traceView') or checkinmethods('pop-pp-view') )and checkinmethods('gui'):
		if not checkinmethods('PhaseLims'):
			methods['PhaseLims'] = [ (-76.,40.),(0.,1.) ]
		elif not (type(methods['PhaseLims']) is list or type(methods['PhaseLims']) is tuple ) or len(methods['PhaseLims']) != 2:
			print "/PhaseLims should be a list of two tuples for x and y coordinats.\n Not list given\n"
			exit(1) 
		elif not (type(methods['PhaseLims'][0]) is list or type(methods['PhaseLims'][0]) is tuple ) or len(methods['PhaseLims'][0]) != 2:
			print "/PhaseLims should be a list of two tuples for x and y coordinats\n First coordinat isn'a a tuple"
			exit(1) 
		elif not (type(methods['PhaseLims'][1]) is list or type(methods['PhaseLims'][1]) is tuple ) or len(methods['PhaseLims'][1]) != 2:
			print "/PhaseLims should be a list of two tuples for x and y coordinats\n Second coordinat isn'a a tuple"
			exit(1) 
	if checkinmethods('traceView') and checkinmethods('gui'):
#>>>>>>>			
		def numsptk(postidx,idxrange):
			prespikes = np.array([])
			trange=t[idxrange]
			sptk = np.zeros(trange.size)
			for nidx in OUTList[postidx]:
				sptime = np.array(neurons[nidx].spks)
				sptime = sptime[ np.where( (sptime > trange[0]) * (sptime < trange[-1]) ) ]
				prespikes = np.append(prespikes,sptime)
			
			prespikes = np.sort(prespikes)
			#print prespikes
			accumulator = 0
			for tm in trange:
				mp = np.where(prespikes < tm)[0]
				sptk[np.where( trange == tm )] = mp.size
			return sptk
			
		def getprespikes(postidx,tl,tr):
			postspk = []
			for nidx in OUTList[postidx]:
				for nspk in neurons[nidx].spks[ np.where( (neurons[nidx].spks >= tl)*(neurons[nidx].spks < tr) ) ]:
					postspk.append([nspk,nidx] )
				
			return np.array( postspk )

		def zoolyupdate(imax):
			zoolyclickevent.spikesymbol = "."
			zoolyclickevent.imax = imax
			onclick1.lines[imax].set_linewidth(4)
			onclick1.lines[imax].set_ls("--")
			zooly.canvas.draw()
			
			zoolyclickevent.v = np.array(neurons[imax].volt)
			zoolyclickevent.u = np.array(neurons[imax].svar)
			zoolyclickevent.g = np.array(neurons[imax].isyng)
			zoolyclickevent.i = np.array(neurons[imax].inoise)
			if hasattr(neurons[imax], "iandnoise"):
				zoolyclickevent.i += np.array(neurons[imax].iandnoise)
			if checkinmethods('sinmod'):
				tstart = methods["sinmod"]["tstart"] if checkinmethods('sinmod/tstart') else 200.
				tstop  = methods["sinmod"]["tstop" ] if checkinmethods('sinmod/tstop' ) else 2200.
				per    = methods["sinmod"]["per"   ] if checkinmethods('sinmod/per'   ) else 2000.
				amp    = methods["sinmod"]["amp"   ] if checkinmethods('sinmod/amp'   ) else 7e-7
				bias   = methods["sinmod"]["bias"  ] if checkinmethods('sinmod/bias'  ) else 0.
				ppa_simmod =  bias - amp/2*(1.-np.cos(np.pi*2./per*(tprin-tstart))) 
				ppa_simmod[np.where(tprin<tstart)] = 0.
				ppa_simmod[np.where(tprin>tstop )] = 0.
				zoolyclickevent.i += ppa_simmod
			if checkinmethods('singmod'):
				gmodtstart = methods["singmod"]["tstart"] if checkinmethods('singmod/tstart') else 200.
				gmodtstop  = methods["singmod"]["tstop" ] if checkinmethods('singmod/tstop' ) else 2200.
				gmodper    = methods["singmod"]["per"   ] if checkinmethods('singmod/per'   ) else 2000.
				gmodmax    = methods["singmod"]["gmax"  ] if checkinmethods('singmod/gmax'  ) else 7e-7
				gmodbias   = methods["singmod"]["bias"  ] if checkinmethods('singmod/bias'  ) else 0.
				zoolyclickevent.gmodE  = methods["singmod"]["E"     ] if checkinmethods('singmod/E'     ) else -75.
				zoolyclickevent.gmod =  gmodbias - gmodmax/2*(1.-np.cos(np.pi*2./gmodper*(tprin-gmodtstart))) 
				zoolyclickevent.gmod[np.where(tprin<gmodtstart)] = 0.
				zoolyclickevent.gmod[np.where(tprin>gmodtstop )] = 0.
			else:
				zoolyclickevent.gmod  = np.zeros(zoolyclickevent.v.shape[0])
				zoolyclickevent.gmodE = 0.
				
			zoolyclickevent.rst = getprespikes(imax,onclick1.tl, onclick1.tr)
			moddyupdate(idx)
			
		def zoolyclickevent(event):
			if not hasattr(onclick1,"lines"): return
			et = event.xdata
			ev = event.ydata
			idx = np.where( np.abs(t-et)<h.dt)[0][0]
			#DB>>
			#print idx, et,ev
			#<<DB
			vmax = abs(neurons[0].volt.x[idx] - ev)
			zoolyclickevent.imax = 0
#			for ind,n in map(None,xrange(methods["ncell"]),neurons):
			for ind,n in enumerate(neurons):
				onclick1.lines[ind].set_linewidth(1)
				onclick1.lines[ind].set_ls("-")
				if vmax > abs(n.volt.x[idx] - ev) :
					vmax = abs(n.volt.x[idx] - ev)
					zoolyclickevent.imax = ind
				#print vmax,n.volt.x[idx],ev

		def moddyclickevent(event):
			et = event.xdata
			idx = np.where( np.abs(t-et)<h.dt)[0][0]
			moddyupdate(idx)


		def moddyupdate(idx):
			if not hasattr(moddyupdate,"tail"):
				moddyupdate.tail = False
			if moddyupdate.tail:
				ridx = np.where( onclick1.idx == idx)[0][0]+1
				moddyupdate.tailsize = 300
				lidx = ridx-moddyupdate.tailsize
				if lidx < 0 :
					lidx = 0
			moddyupdate.idx = idx
			vmin,vmax = methods["PhaseLims"][0]
			nmin,nmax = methods["PhaseLims"][1]
			n0c,v0c,v0n,thc,thn,tpy = getnulls(vindex,vmin,vmax,zoolyclickevent.g[idx], zoolyclickevent.i[idx],neurons[zoolyclickevent.imax].innp.mean,zoolyclickevent.gmod[idx],zoolyclickevent.gmodE)
			#DB>>
			#print "\n\nDB>> THC=",thc,"THN=",thn,tpy != 2 and not thc is None,"\n\n"
			#<<DB
			moddyupdate.rst  = getprespikes(zoolyclickevent.imax,onclick1.tl, onclick1.tr)
			###!!!!zoolyclickevent.lines

			if not hasattr(moddyupdate,"lines"):
				
				dsynmax = np.max(zoolyclickevent.g[onclick1.idx]) if not checkinmethods("tV-synmax") else methods["tV-synmax"]
				#DB>>
				print  "\n\n----\nDB! dsynmax",dsynmax,"\n----\n\n"
				#<<DB
				moddyupdate.lines = [
					faxi.plot([],[], "k"+zoolyclickevent.spikesymbol,ms=9,lw=5)[0]\
						if zoolyclickevent.rst.shape[0] == 0 else \
						faxi.plot(zoolyclickevent.rst[:,0],zoolyclickevent.rst[:,1],"k"+zoolyclickevent.spikesymbol,ms=9,lw=5)[0],
					vaxi.plot(tprin[onclick1.idx],zoolyclickevent.v[onclick1.idx],"k-")[0],
					uaxi.plot(tprin[onclick1.idx],zoolyclickevent.u[onclick1.idx],"k-")[0],
					gaxi.plot(tprin[onclick1.idx],zoolyclickevent.g[onclick1.idx],"k-")[0],
					iaxi.plot(tprin[onclick1.idx],zoolyclickevent.i[onclick1.idx],"k-")[0],
					naxi.scatter(zoolyclickevent.v[onclick1.idx],zoolyclickevent.u[onclick1.idx],\
						c=zoolyclickevent.g[onclick1.idx]/dsynmax,cmap=cmap,vmin=0., vmax=1.,linewidths=0)\
						if checkinmethods("color-phase") else\
					naxi.plot(zoolyclickevent.v[onclick1.idx],zoolyclickevent.u[onclick1.idx],"k-")[0],
					naxi.plot(n0c[:,0],n0c[:,1],"r:",label="n\'=0")[0],
					naxi.plot(v0c[:,0],v0c[:,1],"b-.",label="v\'=0",lw=2)[0],
					#None if checkinmethods("non-isnt-vnul") else naxi.plot(v0n[:,0],v0n[:,1],"b.",mfc="b",mec="b",ms=9)[0],
					None if checkinmethods("non-isnt-vnul") else naxi.plot(v0n[:,0],v0n[:,1],"b--",mfc="b",mec="b",ms=9)[0],
					naxi.plot([zoolyclickevent.v[idx]],[zoolyclickevent.u[idx]],"r.",mfc="r",mec="r",ms=24)[0],
					naxi.plot([],[],"k--",lw=3)[0] if tpy != 2 or thc is None else naxi.plot(thc[:,0],thc[:,1],"k--",lw=3)[0],
					naxi.plot([],[],"k-.",lw=3)[0] if tpy != 2 or thn is None else naxi.plot(thn[:,0],thn[:,1],"k-.",lw=3)[0],
				]
				#try:
					#with open("nulls/threshould-JR.pkl",'rb') as fd:
						#thx = pickle.load(fd)
						#moddyupdate.lines.append(naxi.plot(thx[:,0],thx[:,1],"g--",label="dv/dn=1"),)
				#except:
					#pass
				naxi.legend(loc=0)
			else:
				if zoolyclickevent.rst.shape[0] == 0 :
					moddyupdate.lines[0].set_xdata([])
					moddyupdate.lines[0].set_ydata([])
				else:
					moddyupdate.lines[0].set_xdata(zoolyclickevent.rst[:,0])
					moddyupdate.lines[0].set_ydata(zoolyclickevent.rst[:,1])
				moddyupdate.lines[1].set_xdata(tprin[onclick1.idx])
				moddyupdate.lines[1].set_ydata(zoolyclickevent.v[onclick1.idx])
				moddyupdate.lines[2].set_xdata(tprin[onclick1.idx])
				moddyupdate.lines[2].set_ydata(zoolyclickevent.u[onclick1.idx])
				moddyupdate.lines[3].set_xdata(tprin[onclick1.idx])
				moddyupdate.lines[3].set_ydata(zoolyclickevent.g[onclick1.idx])
				moddyupdate.lines[4].set_xdata(tprin[onclick1.idx])
				moddyupdate.lines[4].set_ydata(zoolyclickevent.i[onclick1.idx])
				##
				if checkinmethods("color-phase"):pass
				elif moddyupdate.tail:
					moddyupdate.lines[5].set_xdata(zoolyclickevent.v[onclick1.idx[lidx:ridx]])
					moddyupdate.lines[5].set_ydata(zoolyclickevent.u[onclick1.idx[lidx:ridx]])
				else:
					moddyupdate.lines[5].set_xdata(zoolyclickevent.v[onclick1.idx])
					moddyupdate.lines[5].set_ydata(zoolyclickevent.u[onclick1.idx])
				moddyupdate.lines[6].set_xdata(n0c[:,0])
				moddyupdate.lines[6].set_ydata(n0c[:,1])
				moddyupdate.lines[7].set_xdata(v0c[:,0])
				moddyupdate.lines[7].set_ydata(v0c[:,1])
				if not checkinmethods("non-isnt-vnul"):
					moddyupdate.lines[8].set_xdata(v0n[:,0])
					moddyupdate.lines[8].set_ydata(v0n[:,1])
				moddyupdate.lines[9].set_xdata([zoolyclickevent.v[idx]])
				moddyupdate.lines[9].set_ydata([zoolyclickevent.u[idx]])
				if tpy == 2 and not thc is None:
					moddyupdate.lines[10].set_xdata(thc[:,0])
					moddyupdate.lines[10].set_ydata(thc[:,1])
				else:
					moddyupdate.lines[10].set_xdata([])
					moddyupdate.lines[10].set_ydata([])
				if tpy == 2 and not thn is None:
					moddyupdate.lines[11].set_xdata(thn[:,0])
					moddyupdate.lines[11].set_ydata(thn[:,1])
				else:
					moddyupdate.lines[11].set_xdata([])
					moddyupdate.lines[11].set_ydata([])
			faxi.set_ylim(0,methods["ncell"])
			vaxi.set_ylim(-85.,40.)
			uaxi.set_ylim(0.,1.)
			if not checkinmethods("tV-synmax"):
				gaxi.set_ylim(0.,zoolyclickevent.g[onclick1.idx].max())
			else:
				gaxi.set_ylim(0.,methods["tV-synmax"])
			#iaxi.set_ylim(zoolyclickevent.i[onclick1.idx].min(),zoolyclickevent.i[onclick1.idx].max())
			faxi.set_xlim(onclick1.tl, onclick1.tr)
			naxi.set_xlim(vmin,vmax)
			naxi.set_ylim(nmin,nmax)
			if not hasattr(moddyupdate,"markers"):
				moddyupdate.markers= [
					faxi.plot([t[idx],t[idx]],[0,methods['ncell']],"r--")[0],
					vaxi.plot([t[idx],t[idx]],[vmin,vmax],"r--")[0],
					uaxi.plot([t[idx],t[idx]],[0,1],"r--")[0],
					gaxi.plot([t[idx],t[idx]],[0,zoolyclickevent.g[onclick1.idx].max()],"r--")[0],
					iaxi.plot([t[idx],t[idx]],iaxi.get_ylim(),"r--")[0]
				]

			else:
				for line in moddyupdate.markers:
					line.set_xdata([t[idx],t[idx]])
			moddy.canvas.draw()

			
		def zoolykeyevent(event):
			if not hasattr(zoolyclickevent,"lines"): return
			if event.key == "K":
				v,u,g,i = (
					np.array(neurons[zoolyclickevent.imax].volt),
					np.array(neurons[zoolyclickevent.imax].svar),
					np.array(neurons[zoolyclickevent.imax].isyng),
					np.array(neurons[zoolyclickevent.imax].inoise)
					)
				sptk = numsptk(zoolyclickevent.imax,onclick1.idx)
				rst = getprespikes(zoolyclickevent.imax,onclick1.tl, onclick1.tr)
				moddyupdate.lines.append(faxi.plot(zoolyclickevent.rst[:,0],zoolyclickevent.rst[:,1],zoolyclickevent.spikesymbol,ms=7,lw=5)[0])
				moddyupdate.lines.append(vaxi.plot(tprin[onclick1.idx],v[onclick1.idx])[0])
				moddyupdate.lines.append(uaxi.plot(tprin[onclick1.idx],u[onclick1.idx])[0])
				moddyupdate.lines.append(gaxi.plot(tprin[onclick1.idx],g[onclick1.idx])[0])
				moddyupdate.lines.append(naxi.plot(v[onclick1.idx],u[onclick1.idx])[0])
				moddyupdate.lines.append(iaxi.plot(tprin[onclick1.idx],i[onclick1.idx])[0])
				#zoolyclickevent.lines.append(saxi.plot(tprin[onclick1.idx],sptk)[0])
			elif event.key == "X":
				for lin in moddyupdate.lines:
					lin.remove()
				del zoolyclickevent.lines
			moddy.canvas.draw()	

		def moddykeyevent(event):
			if event.key == "K" or event.key == "X":
				zoolykeyevent(event)
			elif event.key == "left":
				moddyupdate.idx -= 1
				if moddyupdate.idx < onclick1.idx[0] :
					moddyupdate.idx = onclick1.idx[0]
				moddyupdate(moddyupdate.idx)
			elif event.key == "right":
				moddyupdate.idx += 1
				if moddyupdate.idx > onclick1.idx[-1] :
					moddyupdate.idx = onclick1.idx[-1]
				moddyupdate(moddyupdate.idx)
			elif event.key == "pageup":
				moddyupdate.idx -= 10
				if moddyupdate.idx < onclick1.idx[0] :
					moddyupdate.idx = onclick1.idx[0]
				moddyupdate(moddyupdate.idx)
			elif event.key == "pagedown":
				moddyupdate.idx += 10
				if moddyupdate.idx > onclick1.idx[-1] :
					moddyupdate.idx = onclick1.idx[-1]
				moddyupdate(moddyupdate.idx)
			elif event.key == "home":
				moddyupdate.idx = onclick1.idx[0]
				moddyupdate(moddyupdate.idx)
			elif event.key == "end":
				moddyupdate.idx = onclick1.idx[-1]
				moddyupdate(moddyupdate.idx)
			elif event.key == "T":
				moddyupdate.tail = not moddyupdate.tail
				moddyupdate(moddyupdate.idx)
			elif event.key == "M":
				ridx = np.where( onclick1.idx == moddyupdate.idx)[0][0]+1
				nmax = len(onclick1.idx[ridx::5])
				moddy.set_size_inches(18.5, 10.5, forward=True)
				timestamp = time.strftime("%Y%m%d%H%M%S")
				moviedir = methods["movie-dir"] if checkinmethods("movie-dir") else "/home/rth/tmp/"
				print "=================================="
				print "===        Making MOVIE        ==="
				print "  > Time Stamp                 : ",timestamp
				print "  > Fraim step (mc)            : ",5. * methods['timestep']
				print "  > Tail length (ms)           : ",float(moddyupdate.tailsize) * methods['timestep']
				print "  > Movie Dir                  : ",moviedir
				for ndx,idx in enumerate(onclick1.idx[ridx::5]):
					moddyupdate(idx)
					#moddy.savefig("/home/rth/media/rth-storage-old/rth/tmp/%s-%04d.jpg"%(timestamp,ndx))
					moddy.savefig("%s/%s-%04d.jpg"%(moviedir,timestamp,ndx))
					sys.stderr.write("  > Frame:%04d of %04d         : Saved\r"%(ndx,nmax))
				print "\n==================================\n"
			elif event.key == "S":
				if hasattr(moddykeyevent,"spx"):
					spx.remove()
					del spx
				else:
					with open("separatrix/separatrix.pkl",'rb') as fd:
						spx = pickle.load(fd)
					#spx = np.genfromtxt("quasithresh.dat")[:,1:]
					for fx in np.linspace(0.0,0.27,6):
						naxi.fill_between(spx[:,0],spx[:,1]+fx, spx[:,1]-fx, facecolor='grey', alpha=0.3-fx*0.3/0.27)
					
					spx, = naxi.plot(spx[:,0],spx[:,1],"k--")
				moddy.canvas.draw()	
			#elif event.key == "J":
				#if hasattr(moddykeyevent,"thx"):
					#thx.remove()
					#del thx
				#else:
					##with open("nulls/threshould.pkl",'rb') as fd:
					#with open("nulls/threshould-JR.pkl",'rb') as fd:
						#thx = pickle.load(fd)
					#print thx
					#spx, = naxi.plot(thx[:,0],thx[:,1],"g--")
				#moddy.canvas.draw()	
			elif event.key == "D":
				if hasattr(moddykeyevent,"d10p"):
					for x in moddykeyevent.d10p: x.remove()
					del moddykeyevent.d10p
				else:
					vdata = np.array(neurons[zoolyclickevent.imax].volt)
					udata = np.array(neurons[zoolyclickevent.imax].svar)
					gdata = np.array(neurons[zoolyclickevent.imax].isyng)
					moddykeyevent.d10p = []
					maxg = methods['maxg'] if checkinmethods('maxg') else np.max(gdata)*0.1 
					print "MAXG:",maxg
					d10idx = np.where(gdata<maxg)[0]
					d10idx = [ idx0 for idx0, idx1 in zip(d10idx[:-1],d10idx[1:]) if idx1 != idx0+1 and idx0>=onclick1.idx[0] and idx0<=onclick1.idx[-1]]+\
							 [ idx1 for idx0, idx1 in zip(d10idx[:-1],d10idx[1:]) if idx0 != idx1-1 and idx1>=onclick1.idx[0] and idx1<=onclick1.idx[-1]]
					d10idx.sort()
					moddykeyevent.d10p.append( gaxi.plot(tprin[d10idx], gdata[d10idx], "kX",ms=12)[0] )
					moddykeyevent.d10p.append( naxi.plot(vdata[d10idx], udata[d10idx], "kX",ms=12)[0] )
					
					d10idx = (np.diff(np.sign(np.diff(gdata))) < 0).nonzero()[0] + 1
					d10idx = [ idx for idx in d10idx if idx>=onclick1.idx[0] and idx<=onclick1.idx[-1]]
					d10idx.sort()
					moddykeyevent.d10p.append( gaxi.plot(tprin[d10idx], gdata[d10idx], "rX",ms=12)[0] )
					moddykeyevent.d10p.append( naxi.plot(vdata[d10idx], udata[d10idx], "rX",ms=12)[0] )
				moddy.canvas.draw()
			elif  event.key == "G":
				print np.max(np.array(neurons[zoolyclickevent.imax].isyng) )
			else:
				print event.key
				

#<<<<<<<		
		zooly = plt.figure(4 )
		zooly.canvas.mpl_connect('button_press_event', zoolyclickevent)
		zooly.canvas.mpl_connect('key_press_event', zoolykeyevent)
		moddy = plt.figure(5,figsize=(16,7) )
		faxi = plt.subplot2grid((6,10),(0,0),colspan=4,rowspan=2)
		faxi.set_ylabel("Presynaptic spikes",fontsize=12)
		vaxi = plt.subplot2grid((6,10),(2,0),colspan=4,sharex=faxi)
		vaxi.set_ylabel("V[mV]",fontsize=12)
		uaxi = plt.subplot2grid((6,10),(3,0),colspan=4,sharex=faxi)
		uaxi.set_ylabel('n',fontsize=12)
		gaxi = plt.subplot2grid((6,10),(4,0),colspan=4,sharex=faxi)
		gaxi.set_ylabel(r"$g_{syn} [uS]$",fontsize=12)
		iaxi = plt.subplot2grid((6,10),(5,0),colspan=4,sharex=faxi)
		iaxi.set_ylabel(r"$I_{noise} [nA]$",fontsize=12)
		#saxi = plt.subplot2grid((6,10),(5,0),colspan=4,sharex=faxi)
		naxi = plt.subplot2grid((6,10),(0,5),colspan=6,rowspan=6)
		naxi.set_ylabel('n',fontsize=12)
		naxi.set_xlabel("V[mV]",fontsize=12)
		moddy.canvas.mpl_connect('key_press_event', moddykeyevent)# zoolykeyevent)
		moddy.canvas.mpl_connect('button_press_event', moddyclickevent)


	if checkinmethods('GPcurve') and checkinmethods('gui'):
		plt.figure(7)
		f  = np.array([ [n.gsynscale,n.spks.size]       for n in neurons])
		#f = np.sort(f, axis=0)
		plt.plot(f[:,0] ,f[:,1],"k+")
			
	if checkinmethods('sycleprop') and checkinmethods('gui'):
		plt.figure(8)
		polarax = plt.subplot(111, polar=True)
		#bars = polarax.bar(phydist[:,1], phydist[:,0], width=0.25, bottom=0.0)
		#np.histogram(phydist[:,0], bins=180, weights=phydist[:,1])
		#polarax.hist(phydist[:,0], bins=36, weights=phydist[:,1])
		polarax.bar(phyhistbins[:-1],phyhist,width=phyhistbins[1]-phyhistbins[0],bottom=0)
		#DB>>
		plt.figure(9)
		plt.bar(phyhistbins[:-1],phyhist,width=phyhistbins[1]-phyhistbins[0],bottom=0)
		#<<DB
	if checkinmethods('Gtot-dist') and checkinmethods('gui'):
		plt.figure(10)
		plt.bar(gskbins[:-1]+(gskbins[1]+gskbins[0])/2.,gskhist,width=(gskbins[1]-gskbins[0]),color="k")
		#plt.hist(gsk,bins=methods["ncell"]/50)
		plt.title((methods["MainFigTitle"] if checkinmethods("MainFigTitle") else "")+"mean total synaptic conductance={}(uS), stdr total synaptic conductance={}(uS)".format(mgto,sgto) )
		plt.ylabel("Probability")
		plt.xlabel("Toaol conductance (uS)")

	if checkinmethods("Delay-stat") and checkinmethods('gui'):
		plt.figure(11)
		plt.bar((dlybins[1:]+dlybins[:-1])/2.,dlyhist,width=dlybins[1]-dlybins[0],color="k")
		plt.title((methods["MainFigTitle"] if checkinmethods("MainFigTitle") else "")+"mean delay={}(ms), stdr delay={}(ms)".format(mdly,sdly))
		plt.ylabel("Probability")
		plt.xlabel("delay(ms)")
	
	if checkinmethods('T1vsT2/spikerate') and checkinmethods('gui'):
		plt.figure(12)
		plt.bar([1,2],[methods['T1vsT2']['spikerate']['T1'],methods['T1vsT2']['spikerate']['T2']],width=0.1,color="k")
		plt.title((methods["MainFigTitle"] if checkinmethods("MainFigTitle") else "")+"Mean firing rate")
		plt.ylabel("Firing rate spike/sec")
		plt.xlabel("Type")

	if checkinmethods('pop-pp-view') and checkinmethods('gui'):
###>>>		
		def getVfield(vmap,nmap,gsyn,iapp,gstim = 0.,estim = 0.):
			type21,gna,ena,gk,ek,gl,el,n0,sn,t0,st,v12,sl,a,b,s,es = getType21(0)
			## N-component
			dndt = lambda v,n:(n0 + sn/(1.+np.exp(-(v-v12)/sl )) - n)/( t0 + st/(1.+np.exp( (v+40.)/20.)) ) 
			## V-component
			dvdt = lambda v,n:(-iapp- gstim*(v-estim) - gsyn*(v- es))*1e-3 /s+gl*(el-v)+gna*(1./(1.+np.exp(-(v+40.)/9.5)))**3*(b*n+a)*(ena-v)+gk*n**4*(ek-v)	
			return np.array([ v         for v in vmap for n in nmap]),np.array([ n         for v in vmap for n in nmap]),\
			       np.array([ dvdt(v,n) for v in vmap for n in nmap]),np.array([ dndt(v,n) for v in vmap for n in nmap])
###------------------------------------
		def PopPPview_update(idx):
			ppp_ppp = np.array(
				[ (n.volt.x[idx],n.svar.x[idx]) for n in neurons]
			)
			if checkinmethods('pop-pp-view-color'):
				ppp_pfpp.set_offsets(ppp_ppp)
			else:
				ppp_pfpp.set_xdata(ppp_ppp[:,0])
				ppp_pfpp.set_ydata(ppp_ppp[:,1])
			vmin,vmax = methods["PhaseLims"][0]
			nmin,nmax = methods["PhaseLims"][1]
			if checkinmethods('pop-pp-vector-field'):
				_,_,vfl.U,vfl.V = getVfield(vec_v,vec_n,\
					float(ppp_av_gsyn[idx]),
					float(ppp_av_mean+ppp_simmod[idx-tproc]),\
					float(ppp_simgmod[idx-tproc]), ppp_gmodE)
				vfl.U /= (vmax-vmin)
				vfl.V /= (nmax-nmin)
				if checkinmethods('pop-pp-minmax'):
					_,_,vflmin.U,vflmin.V = getVfield(vec_v,vec_n,float(ppp_av_gsyn[ppp_idx]),float(ppp_av_min+ppp_simmod[ppp_idx-tproc]),ppp_simgmod[ppp_idx-tproc], ppp_gmodE)
					vflmin.U /= vmax-vmin
					vflmin.V /= nmax-nmin
				 	_,_,vflmax.U,vflmax.V = getVfield(vec_v,vec_n,float(ppp_av_gsyn[ppp_idx]),float(ppp_av_max+ppp_simmod[ppp_idx-tproc]),ppp_simgmod[ppp_idx-tproc], ppp_gmodE)
					vflmax.U /= vmax-vmin
					vflmax.V /= nmax-nmin
				

			
			n0c,v0c,v0n,thc,thn,type21 = getnulls(0,vmin,vmax,\
				float(ppp_av_gsyn[idx]),\
				float(ppp_av_mean+ppp_simmod[idx-tproc]),\
				float(ppp_av_mean+ppp_simmod[idx-tproc]),\
				float(ppp_simgmod[idx-tproc]), ppp_gmodE)
			ppp_sax_m.set_xdata(t[idx])
			ppp_vax_m.set_xdata(t[idx])
			ppp_nax_m.set_xdata(t[idx])
			ppp_gax_m.set_xdata(t[idx])
			ppp_pfn0.set_xdata(n0c[:,0])
			ppp_pfn0.set_ydata(n0c[:,1])
#			ppp_pfv0.set_xdata(v0c[:,0])
#			ppp_pfv0.set_ydata(v0c[:,1])
			ppp_pfvN.set_xdata(v0n[:,0])
			ppp_pfvN.set_ydata(v0n[:,1])
#			ppp_pfth0.set_xdata(thc[:,0])
#			ppp_pfth0.set_ydata(thc[:,1])
			ppp_pfthi.set_xdata(thn[:,0] if type21 == 2 else [])
			ppp_pfthi.set_ydata(thn[:,1] if type21 == 2 else [])
			if checkinmethods('pop-pp-minmax'):
				_,_,v0n,_,thn,_ = getnulls(0,vmin,vmax,float(ppp_av_gsyn[ppp_idx]),float(ppp_av_max+ppp_simmod[ppp_idx-tproc]),float(ppp_av_max+ppp_simmod[ppp_idx-tproc]),ppp_simgmod[ppp_idx-tproc], ppp_gmodE)	
				ppp_pfvNmax.set_xdata(v0n[:,0])
				ppp_pfvNmax.set_ydata(v0n[:,1])
				ppp_pfthi_max.set_xdata(thn[:,0] if type21 == 2 else [])
				ppp_pfthi_max.set_ydata(thn[:,1] if type21 == 2 else [])
				_,_,v0n,_,thn,_ = getnulls(0,vmin,vmax,float(ppp_av_gsyn[ppp_idx]),float(ppp_av_min+ppp_simmod[ppp_idx-tproc]),float(ppp_av_min+ppp_simmod[ppp_idx-tproc]),ppp_simgmod[ppp_idx-tproc], ppp_gmodE)
				ppp_pfvNmin.set_xdata(v0n[:,0])
				ppp_pfvNmin.set_ydata(v0n[:,1])
				ppp_pfthi_min.set_xdata(thn[:,0] if type21 == 2 else [])
				ppp_pfthi_min.set_ydata(thn[:,1] if type21 == 2 else [])


			PopPPview.canvas.draw()
		
		def PopPPview_keyevent(event):
			global ppp_idx, ppp_lines, ppp_av_kpos_cnt
			if event.key == "K":
				n0c,v0c,v0n,thc,thn,type21 = getnulls(0,vmin,vmax,float(ppp_av_gsyn[ppp_idx]),float(ppp_av_mean+ppp_simmod[ppp_idx-tproc]),float(ppp_av_mean+ppp_simmod[ppp_idx-tproc]) ,ppp_simgmod[ppp_idx-tproc], ppp_gmodE)
				ppp_lines.append( ppp_ppax.plot(v0n[:,0],v0n[:,1],"-",ms=3)[0] )
				if type21 == 2:
					ppp_lines.append(ppp_ppax.plot(thn[:,0],thn[:,1],"-",lw=1)[0])
			if event.key == "X":
				for lin in ppp_lines.append:
					lin.remove()
				ppp_lines.append = []
			elif event.key == "left":      ppp_idx -= 1
			elif event.key == "right":     ppp_idx += 1
			elif event.key == "pageup":    ppp_idx -= 10
			elif event.key == "pagedown":  ppp_idx += 10
			elif event.key == "home":      ppp_idx  = tproc
			elif event.key == "end":       ppp_idx  = tproc+tprin.size
			elif event.key == "<"  :
				ppp_av_kpos_cnt -= 1
				if ppp_av_kpos_cnt < 0:ppp_av_kpos_cnt = ppp_av_kpos.shape[0]-1
				ppp_idx  = ppp_av_kpos[ppp_av_kpos_cnt]
			elif event.key == ">"  :
				ppp_av_kpos_cnt += 1
				if ppp_av_kpos_cnt > ppp_av_kpos.shape[0]-1:ppp_av_kpos_cnt = 0
				ppp_idx  = ppp_av_kpos[ppp_av_kpos_cnt]
			elif event.key == "M":
				#PopPPview.set_size_inches(18.5, 10.5, forward=True)
				nmax = tprin.size/5
				timestamp = time.strftime("%Y%m%d%H%M%S")
				moviedir = methods["movie-dir"] if checkinmethods("movie-dir") else "/home/rth/tmp/"
				print "=================================="
				print "===        Making MOVIE        ==="
				print "  > Time Stamp                 : ",timestamp
				print "  > Fraim step (mc)            : ",5. * methods['timestep']
				#print "  > Tail length (ms)           : ",float(moddyupdate.tailsize) * methods['timestep']
				print "  > Movie Dir                  : ",moviedir
				

				for ndx,idx in enumerate(range(tproc,tproc+tprin.size,5)):
					PopPPview_update(idx)
					#PopPPview.savefig("/home/rth/tmp/movies/%s-%04d.jpg"%(timestamp,ndx))
					#PopPPview.savefig("/home/rth/media/rth-media/rth-media/tmp/movie/%s-%04d.jpg"%(timestamp,ndx))
					PopPPview.savefig("%s/%s-%04d.jpg"%(moviedir,timestamp,ndx))
					sys.stderr.write("  > Frame:%04d of %04d         : Saved\r"%(ndx,nmax))
				os.system("ffmpeg -r 10 -f image2  -i %s/%s-%%04d.jpg -vcodec libx264 -crf 25  -pix_fmt yuv420p %s/%s.mp4"%(moviedir,timestamp,moviedir,timestamp) )
				os.system("rm %s/%s-*.jpg"%(moviedir,timestamp) )
				print "\n==================================\n"
				
			if ppp_idx < tproc:            ppp_idx  = tproc
			if ppp_idx > tproc+tprin.size: ppp_idx  = tproc+tprin.size
			PopPPview_update(ppp_idx)
			
		def PopPPview_clickevent(event):
			global ppp_idx
			et = event.xdata
			ppp_idx = np.where( np.abs(t-et)<h.dt)[0][0]
			PopPPview_update(ppp_idx)
		
		PopPPview=plt.figure(13,figsize=(16,7) )
		PopPPview.suptitle(methods["MainFigTitle"] if checkinmethods("MainFigTitle") else "")
		ppp_lines = []
		ppp_idx = tprin.size/2+tproc
		ppp_sax = plt.subplot2grid((4,3),(0,0) )#,colspan=4,rowspan=2)
		ppp_sax.set_ylabel("neuron",fontsize=12)

		if checkinmethods('singmod'):
			gmodtstart = methods["singmod"]["tstart"] if checkinmethods('singmod/tstart') else 200.
			gmodtstop  = methods["singmod"]["tstop" ] if checkinmethods('singmod/tstop' ) else 2200.
			gmodper    = methods["singmod"]["per"   ] if checkinmethods('singmod/per'   ) else 2000.
			gmodmax    = methods["singmod"]["gmax"  ] if checkinmethods('singmod/gmax'  ) else 7e-7
			gmodbias   = methods["singmod"]["bias"  ] if checkinmethods('singmod/bias'  ) else 0.
			ppp_gmodE  = methods["singmod"]["E"     ] if checkinmethods('singmod/E'     ) else -75.
			ppp_simgmod =  gmodbias - gmodmax/2*(1.-np.cos(np.pi*2./gmodper*(tprin-gmodtstart))) 
			ppp_simgmod[np.where(tprin<gmodtstart)] = 0.
			ppp_simgmod[np.where(tprin>gmodtstop )] = 0.
		else:
			ppp_simgmod = np.zeros(tprin.shape)
			ppp_gmodE   = 0.

		if checkinmethods('sinmod'):
			tstart = methods["sinmod"]["tstart"] if checkinmethods('sinmod/tstart') else 200.
			tstop  = methods["sinmod"]["tstop" ] if checkinmethods('sinmod/tstop' ) else 2200.
			per    = methods["sinmod"]["per"   ] if checkinmethods('sinmod/per'   ) else 2000.
			amp    = methods["sinmod"]["amp"   ] if checkinmethods('sinmod/amp'   ) else 7e-7
			bias   = methods["sinmod"]["bias"  ] if checkinmethods('sinmod/bias'  ) else 0.
			ppp_simmod =  bias - amp/2*(1.-np.cos(np.pi*2./per*(tprin-tstart))) 
			ppp_simmod[np.where(tprin<tstart)] = 0.
			ppp_simmod[np.where(tprin>tstop )] = 0.
		else:
			ppp_simmod = np.zeros(tprin.shape)

		vmin,vmax = methods["PhaseLims"][0]
		nmin,nmax = methods["PhaseLims"][1]
		
		ppp_av_mean = np.mean( [ n.innp.mean for n in neurons ] )
		if checkinmethods('pop-pp-minmax'):
			ppp_av_min  = np.max ( [ n.innp.mean for n in neurons ] )
			ppp_av_max  = np.min ( [ n.innp.mean for n in neurons ] )
			
		ppp_av_gsyn = np.array(neurons[0].isyng)
		for n in neurons[1:]:
			ppp_av_gsyn += np.array(n.isyng)
		ppp_av_gsyn /= methods['ncell']
		
		if checkinmethods('pop-pp-view-color'):
			ppp_color = np.array([ -n.innp.mean for n in neurons ])
			ppp_color /= np.amax(ppp_color)
			if not checkinmethods('pop-pp-view-nrnsize'): methods['pop-pp-view-nrnsize'] = 52

		if checkinmethods('pop-pp-view-color'):
			ppp_s_color = np.array([ -neurons[nindex[int(n)][1]].innp.mean for n in rast[:,1] ])
			ppp_s_color /= np.amax(ppp_s_color)
			if not checkinmethods('pop-pp-view-spsize'): methods['pop-pp-view-spsize'] = 3

		if checkinmethods('pop-pp-view-color'):
			ppp_sax.scatter(rast[:,0],rast[:,1], s=methods['pop-pp-view-spsize'],c=ppp_s_color, cmap=matplotlib.cm.get_cmap('rainbow'))#,ms=9)
		else:
			ppp_sax.plot(rast[:,0],rast[:,1],"k"+methods['rstmark'],mew=0.75,ms=methods['rstmarksize'])

		ppp_sax_m, = ppp_sax.plot([t[ppp_idx],t[ppp_idx]],[0,methods['ncell']],"r--")
		ppp_vax = plt.subplot2grid((4,3),(1,0),sharex=ppp_sax) 
		ppp_vax.set_ylabel("V[mV]",fontsize=12)
		xmin,xmax = methods["PhaseLims"][0]
		for n in neurons:
			ppp_vax.plot(tprin,np.array(n.volt)[tproc:tprin.size+tproc],'-',c="#000000",lw=0.1)
			xm, xM = np.min(np.array(n.volt)[tproc:tprin.size+tproc]),np.max(np.array(n.volt)[tproc:tprin.size+tproc])
			if xmin > xm : xmin = xm
			if xmax < xM : xmax = xM
		ppp_vax_m, = ppp_vax.plot([t[ppp_idx],t[ppp_idx]],[xmin,xmax],"r--")

		ppp_nax = plt.subplot2grid((4,3),(2,0),sharex=ppp_sax)
		ppp_nax.set_ylabel('n',fontsize=12)
		xmin,xmax = methods["PhaseLims"][1]
		for n in neurons:
			ppp_nax.plot(tprin,np.array(n.svar)[tproc:tprin.size+tproc],'-',c="#000000",lw=0.1)
			xm, xM = np.min(np.array(n.svar)[tproc:tprin.size+tproc]),np.max(np.array(n.svar)[tproc:tprin.size+tproc])
			if xmin > xm : xmin = xm
			if xmax < xM : xmax = xM
		ppp_nax_m, = ppp_nax.plot([t[ppp_idx],t[ppp_idx]],[xmin,xmax],"r--")

		ppp_gax = plt.subplot2grid((4,3),(3,0),sharex=ppp_sax)
		ppp_gax.set_ylabel(r"mean $g_{syn} [uS]$",fontsize=12)

		ppp_av_kpos  = [idx for idx in (np.diff(np.sign(np.diff(ppp_av_gsyn))) < 0).nonzero()[0] + 1]
		ppp_av_kpos += [idx for idx in (np.diff(np.sign(np.diff(ppp_av_gsyn))) > 0).nonzero()[0] + 1]
		ppp_av_kpos.sort()
		ppp_av_kpos += \
			[ np.argmin( abs(ppp_av_gsyn[l:r]-(ppp_av_gsyn[l]+ppp_av_gsyn[r])/2.) )+l\
				for l,r in zip(ppp_av_kpos[:-1],ppp_av_kpos[1:]) ]
		ppp_av_kpos.sort()
		ppp_av_kpos = np.array(ppp_av_kpos)
		ppp_av_kpos = ppp_av_kpos[np.where( (ppp_av_kpos>=tproc)*( ppp_av_kpos<=(tprin.size+tproc) ) ) ]
		ppp_av_kpos_cnt = 0
		plt.plot(tprin[ppp_av_kpos-tproc], np.zeros(ppp_av_kpos.shape[0]), "^")

		ppp_gax.plot(tprin,np.array(ppp_av_gsyn)[tproc:tprin.size+tproc],"k-")
		xmin,xgmax = np.min(ppp_av_gsyn),np.max(ppp_av_gsyn)
		ppp_gax_m, = ppp_gax.plot([t[ppp_idx],t[ppp_idx]],[xmin,xgmax],"r--")
		
		ppp_ppax = plt.subplot2grid((4,3),(0,1),colspan=2,rowspan=5)
		ppp_ppax.set_ylabel('n',fontsize=12)
		ppp_ppax.set_xlabel("V[mV]",fontsize=12)
		
		ppp_ppp = np.array(
			[ (n.volt.x[ppp_idx],n.svar.x[ppp_idx]) for n in neurons]
		)
		if checkinmethods('pop-pp-vector-field'):
			if   type(methods['pop-pp-vector-field']) is tuple and len(methods['pop-pp-vector-field']) == 6:
				flvmin,flvmax = methods['pop-pp-vector-field'][2],methods['pop-pp-vector-field'][3]
				flnmin,flnmax = methods['pop-pp-vector-field'][4],methods['pop-pp-vector-field'][5]
				methods['pop-pp-vector-field'] = methods['pop-pp-vector-field'][0],methods['pop-pp-vector-field'][1]
			elif type(methods['pop-pp-vector-field']) is tuple and len(methods['pop-pp-vector-field']) == 4:
				fldv,fldn = methods['pop-pp-vector-field'][2],methods['pop-pp-vector-field'][3]
				methods['pop-pp-vector-field'] = methods['pop-pp-vector-field'][0],methods['pop-pp-vector-field'][1]
				flvmin,flvmax = vmin,vmin+(vmax-vmin)*fldv
				flnmin,flnmax = nmin,nmin+(nmax-nmin)*fldn
			elif type(methods['pop-pp-vector-field']) is not tuple or len(methods['pop-pp-vector-field']) !=2:
				fldv,fldn = 1., 1.
				methods['pop-pp-vector-field'] = 20,10
				flvmin,flvmax = vmin,vmax
				flnmin,flnmax = nmin,nmax
			else:
				fldv,fldn = 1., 1.
				flvmin,flvmax = vmin,vmax
				flnmin,flnmax = nmin,nmax
				
			vec_v = np.linspace(flvmin,flvmax,methods['pop-pp-vector-field'][0])
			vec_n = np.linspace(flnmin,flnmax,methods['pop-pp-vector-field'][1])

			map_v,map_n,vf,nf = getVfield(vec_v,vec_n,float(ppp_av_gsyn[ppp_idx]),float(ppp_av_mean+ppp_simmod[ppp_idx-tproc]),ppp_simgmod[ppp_idx-tproc], ppp_gmodE)
			vf /= vmax-vmin
			nf /= nmax-nmin
			
			if checkinmethods('pop-pp-minmax'):
				_,_,vfmin,nfmin = getVfield(vec_v,vec_n,float(ppp_av_gsyn[ppp_idx]),float(ppp_av_min+ppp_simmod[ppp_idx-tproc]),ppp_simgmod[ppp_idx-tproc], ppp_gmodE)
				vfmin /= vmax-vmin
				nfmin /= nmax-nmin
				
				_,_,vfmax,nfmax = getVfield(vec_v,vec_n,float(ppp_av_gsyn[ppp_idx]),float(ppp_av_max+ppp_simmod[ppp_idx-tproc]),ppp_simgmod[ppp_idx-tproc], ppp_gmodE)
				vfmax /= vmax-vmin
				nfmax /= nmax-nmin
				vflmax = ppp_ppax.quiver(map_v,map_n,vfmax,nfmax,color='r',scale=3.,headwidth=1.7,width=0.005*0.5)
				vfl = ppp_ppax.quiver(map_v,map_n,vf,nf,scale=3.,headwidth=1.7,width=0.005*0.5)
				vflmin = ppp_ppax.quiver(map_v,map_n,vfmin,nfmin,color='b',scale=3.,headwidth=1.7,width=0.005*0.5)
			else:
				vfl = ppp_ppax.quiver(map_v,map_n,vf,nf,scale=3.)
			
			
		if checkinmethods('pop-pp-view-color'):
			ppp_pfpp = ppp_ppax.scatter(ppp_ppp[:,0],ppp_ppp[:,1], s=methods['pop-pp-view-nrnsize'],c=ppp_color, cmap=matplotlib.cm.get_cmap('rainbow'))#,ms=9)
		else:
			ppp_pfpp, = ppp_ppax.plot(ppp_ppp[:,0],ppp_ppp[:,1],"ro",ms=9)
		if checkinmethods('sinmod'):
			ppp_gax.plot(tprin,(bias-ppp_simmod/amp)*xgmax,"b--")
		if checkinmethods('singmod'):
			ppp_gax.plot(tprin,(gmodbias-ppp_simgmod/gmodmax)*xgmax,"r--")
			
		
		n0c,v0c,v0n,thc,thn,type21 = getnulls(0,vmin,vmax,float(ppp_av_gsyn[ppp_idx]),float(ppp_av_mean+ppp_simmod[ppp_idx-tproc]),float(ppp_av_mean+ppp_simmod[ppp_idx-tproc]),ppp_simgmod[ppp_idx-tproc], ppp_gmodE)
		ppp_pfn0, = ppp_ppax.plot(n0c[:,0],n0c[:,1],"g-",lw=2)
		#ppp_pfv0, = ppp_ppax.plot(v0c[:,0],v0c[:,1],"c--",lw=1)
		ppp_pfvN, = ppp_ppax.plot(v0n[:,0],v0n[:,1],"k-",ms=3)
		#ppp_pfth0 = ppp_ppax.plot(thc[:,0],thc[:,1],"c-.",lw=1)[0] if type21 == 2 else  ppp_ppax.plot([],[],"c-.",lw=1)[0]
		ppp_pfthi = ppp_ppax.plot(thn[:,0],thn[:,1],"k-.",lw=1)[0] if type21 == 2 else ppp_ppax.plot([],[],"c-.",lw=1)[0]
		if checkinmethods('pop-pp-minmax'):
			_,_,v0n,_,thn,_ = getnulls(0,vmin,vmax,float(ppp_av_gsyn[ppp_idx]),float(ppp_av_min+ppp_simmod[ppp_idx-tproc]),float(ppp_av_min+ppp_simmod[ppp_idx-tproc]),ppp_simgmod[ppp_idx-tproc], ppp_gmodE)	
			ppp_pfvNmin , = ppp_ppax.plot(v0n[:,0],v0n[:,1],"b-",ms=3)
			ppp_pfthi_min = ppp_ppax.plot(thn[:,0],thn[:,1],"b-.",lw=1)[0] if type21 == 2 else ppp_ppax.plot([],[],"b-.",lw=2)[0]
			_,_,v0n,_,thn,_ = getnulls(0,vmin,vmax,float(ppp_av_gsyn[ppp_idx]),float(ppp_av_max+ppp_simmod[ppp_idx-tproc]),float(ppp_av_max+ppp_simmod[ppp_idx-tproc]),ppp_simgmod[ppp_idx-tproc], ppp_gmodE)	
			ppp_pfvNmax , = ppp_ppax.plot(v0n[:,0],v0n[:,1],"r-",ms=3)
			ppp_pfthi_max = ppp_ppax.plot(thn[:,0],thn[:,1],"r-.",lw=1)[0] if type21 == 2 else ppp_ppax.plot([],[],"r-.",lw=2)[0]
		ppp_ppax.set_ylim(nmin,nmax)
		
		ppp_ppax.set_xlim(vmin,vmax)
		PopPPview.canvas.mpl_connect('key_press_event',    PopPPview_keyevent)# zoolykeyevent)
		PopPPview.canvas.mpl_connect('button_press_event', PopPPview_clickevent)

	if checkinmethods("ttFFT"):
		if methods['tracetail'] == 'total current' or methods['tracetail'] == 'TI' or methods['tracetail'] == 'mean total current' or methods['tracetail'] == 'MTI'\
		  or methods['tracetail'] == 'total synaptic current' or methods['tracetail'] == 'TSI' or methods['tracetail'] == 'mean total synaptic current' or methods['tracetail'] == 'MTSI':
			ttfft_dt = np.mean(tprin[1:]-tprin[:-1])
			ttfft_sr = meancur[tproc:tprin.size+tproc]
			ttfft_dr = (tprin[-1] - tprin[0])
			ttspecX	= np.arange(ttfft_sr.shape[0], dtype=float)
			ttspecX	*= 1000.0/ttfft_dr
			ttpnum 	= int(200.*ttfft_dr/1000.0)
			ttspecX	= ttspecX[:ttpnum]
			ttfft = np.abs( np.fft.fft(ttfft_sr)*1.0/ttfft_dr )**2
			if checkinmethods('gui'):
				ttftpf = plt.figure(12)
				plt.bar(ttspecX[1:],ttfft[1:ttpnum],0.75,color="k",edgecolor="k")
				plt.ylabel("Power ($nA^2$)", fontsize=16)
				plt.xlabel("Frequency (Hz)", fontsize=16)
				if methods['tracetail'] == 'total current' or methods['tracetail'] == 'TI':
					plt.title((methods["MainFigTitle"] if checkinmethods("MainFigTitle") else "")+"Power Spectrum of Total Current", fontsize=16)
				elif methods['tracetail'] == 'mean total current' or methods['tracetail'] == 'MTI':
					plt.title((methods["MainFigTitle"] if checkinmethods("MainFigTitle") else "")+"Power Spectrum of Mean Total Current", fontsize=16)
				elif methods['tracetail'] == 'total synaptic current' or methods['tracetail'] == 'TSI':
					plt.title((methods["MainFigTitle"] if checkinmethods("MainFigTitle") else "")+"Power Spectrum of Total Synaptric Current", fontsize=16)
				elif methods['tracetail'] == 'mean total synaptic current' or methods['tracetail'] == 'MTSI':
					plt.title((methods["MainFigTitle"] if checkinmethods("MainFigTitle") else "")+"Power Spectrum of Mean Total Synaptric Current", fontsize=16)
		elif methods['tracetail'] == 'total conductance' or methods['tracetail'] == 'mean total conductance' or methods['tracetail'] == 'TG' or methods['tracetail'] == 'MTG':
			ttfft_dt = np.mean(tprin[1:]-tprin[:-1])
			ttfft_sr = meancur[tproc:tprin.size+tproc]
			ttfft_dr = (tprin[-1] - tprin[0])
			ttspecX	= np.arange(ttfft_sr.shape[0], dtype=float)
			ttspecX	*= 1000.0/ttfft_dr
			ttpnum 	= int(200.*ttfft_dr/1000.0)
			ttspecX	= ttspecX[:pnum]
			ttfft = np.abs( np.fft.fft(ttfft_sr)*1.0/ttfft_dr )**2			
			if checkinmethods('gui'):
				ttftpf = plt.figure(12)
				plt.bar(ttspecX[1:],ttfft[1:pnum],0.75,color="k",edgecolor="k")
				plt.ylabel("Power ($\mu S^2$)", fontsize=16)
				plt.xlabel("Frequency (Hz)", fontsize=16)
				if   methods['tracetail'] == 'total conductance' or methods['tracetail'] == 'TG':
					plt.title((methods["MainFigTitle"] if checkinmethods("MainFigTitle") else "")+"Power Spectrum of Total Conductance", fontsize=16)
				elif methods['tracetail'] == 'mean total conductance' or methods['tracetail'] == 'MTG':
					plt.title((methods["MainFigTitle"] if checkinmethods("MainFigTitle") else "")+"Power Spectrum of Mean Total Conductance", fontsize=16)
		elif methods['tracetail'] == 'p2eLFP':			
			ttfft_lfp = lfp[int(round(methods['cliptrn'])):]
			ttfft_dt = 1.
			ttfft_dr = float(ttfft_lfp.shape[0])
			ttspecX	= np.arange(ttfft_lfp.shape[0], dtype=float)
			ttspecX	*= 1000.0/ttfft_dr
			ttpnum 	= int(200.*ttfft_dr/1000.0)
			ttspecX	= ttspecX[:pnum]
			ttfft = np.abs( np.fft.fft(ttfft_lfp)*1.0/ttfft_dr )**2	
			if checkinmethods('gui'):
				ttftpf = plt.figure(12)
				plt.bar(ttspecX[1:],ttfft[1:pnum],0.75,color="k",edgecolor="k")
				plt.ylabel("Power ($\mu S^2$)", fontsize=16)
				plt.xlabel("Frequency (Hz)", fontsize=16)
				plt.title((methods["MainFigTitle"] if checkinmethods("MainFigTitle") else "")+"Power Spectrum of Generated LFP", fontsize=16)
		methods["ttFFT-Retults"]={
			'Freq'  : ttspecX[1:],
			'Power' : ttfft[1:pnum]
			}
	if ( checkinmethods('N2NHI') or checkinmethods('N2NHI-netISI') ) and checkinmethods('gui'): 
		plt.figure(19)
		if checkinmethods('N2NHI') and checkinmethods('N2NHI-netISI'):
			clsAx = plt.subplot(121)
			plt.bar(np.arange(len(ccc_clsidx))+1,np.array(ccc_clsidx),0.75,color="k",edgecolor="k")
			plt.subplot(122, sharex=clsAx,sharey=clsAx)
			plt.bar(np.arange(len(rth_clsidx))+1,np.array(rth_clsidx),0.75,color="k",edgecolor="k")
			plt.xlim(1,10)
		elif checkinmethods('N2NHI'):
			plt.bar(np.arange(len(ccc_clsidx))+1,np.array(ccc_clsidx),0.75,color="k",edgecolor="k")
			plt.xlim(1,10)
		elif checkinmethods('N2NHI-netISI'):
			plt.bar(np.arange(len(rth_clsidx))+1,np.array(rth_clsidx),0.75,color="k",edgecolor="k")
			plt.xlim(1,10)
		plt.figure(20)
		plt.bar(cg_nrnbin,cg_nrnisi,0.75,color="k",edgecolor="k")

	if checkinmethods('vpop-view') and checkinmethods('gui'):
		import matplotlib as mpl
		from mpl_toolkits.mplot3d import Axes3D
		vpop_view=plt.figure(21)
		ax = vpop_view.add_subplot(111, projection='3d')
		for nidx in xrange(methods['ncell']-1,-1,-1):
			ax.plot(tprin,float(nidx)*np.ones(tprin.shape[0]),np.array(neurons[nindex[nidx][1]].volt)[tproc:tproc+tprin.size],"-")
	if checkinmethods('CtrlISI') and checkinmethods('gui'):
		plt.figure(22)
		plt.bar(CtrISI[:,0],CtrISI[:,1],methods["CtrlISI"]["bin"],color="k",edgecolor="k")
		plt.xlim(0,methods["CtrlISI"]["bin"]+methods["CtrlISI"]["max"])
		plt.xlabel("ISI (ms)")
		plt.ylabel("Number of spikes")
		if Pnet is not None:
			pnets = np.arange(methods["CtrlISI"]["max"]/Pnet)*Pnet+methods["CtrlISI"]["bin"]/2
			pnets = pnets[1:]
			xmax = np.amax(CtrISI[:,1])
			plt.plot(pnets,np.ones(pnets.shape[0])*xmax*1.1,"kv")
			plt.text(methods["CtrlISI"]["max"]/2,xmax/2,r"$F_{{nrn}}/F_{{net}}$"+"\n{}".format(Pnet/methods["nrnPmean"]),fontsize=24)
	if checkinmethods('spike2net-dist') and checkinmethods('gui') and checkinmethods('spike2net-dist-result'):
		plt.figure(23)
		peakd = np.array(methods['spike2net-dist-result'])
		settd = np.arange(-100.,100)
		plt.bar(settd,peakd,0.75,color="k",edgecolor="k")
		plt.ylabel("Number of spikes")
		plt.xlabel(r"Phase of network period ( \% )")
	if checkinmethods('nrnFRhist') and checkinmethods('gui'):
		plt.figure(24)
		if type(methods['nrnFRhist']) is dict:
			#filter for names:
			hparam = {}
			for n in methods['nrnFRhist']:
				if n in ["bins","range","normed","weights","density"]:
					hparam[n]=methods['nrnFRhist'][n]
				elif n == "part" and type(methods['nrnFRhist'][n]) is bool and methods['nrnFRhist'][n]:
					nrnfr /= netfrqmean
			nrnfrhist,nrnfrages=np.histogram(nrnfr,**hparam)
			if "xnorm" in methods['nrnFRhist'] and methods['nrnFRhist']["xnorm"]:
				nrnfrhist = nrnfrhist.astype(float)/float(np.sum(nrnfrhist))
			nrnfrages = (nrnfrages[1:]+nrnfrages[:-1])/2.
			plt.bar(nrnfrages,nrnfrhist,width=(nrnfrages[1]-nrnfrages[0])*0.9,edgecolor='k',facecolor='k')
			if "ymax" in methods['nrnFRhist']:
				plt.ylim(ymax=methods['nrnFRhist']["ymax"])
			if "ymin" in methods['nrnFRhist']:
				plt.ylim(ymin=methods['nrnFRhist']["ymin"])
		else:
			nrnfrhist,nrnfrages=np.histogram(nrnfr)
			nrnfrages = (nrnfrages[1:]+nrnfrages[:-1])/2.
			plt.bar(nrnfrages,nrnfrhist,width=(nrnfrages[1]-nrnfrages[0])*0.9,edgecolor='k',facecolor='k')
	if checkinmethods('PAC-Tort10') and checkinmethods('gui'):
		# or checkinmethods('PAC-Canolty06') or checkinmethods('PAC-Tort10') or checkinmethods('PAC-VS')
		plt.figure(25)
		plt.bar(bct,hist,color='k')
		
			
	if checkinmethods('git'):
		os.system("git commit -a &")
	if checkinmethods('beep'):
		os.system("beep &")
	if checkinmethods('corelog'):
		def writetree(tree,fd,prefix):
			for name in tree:
				if type(tree[name]) is dict:
					writetree(tree[name],fd,prefix+name+'/')
				elif isinstance(tree[name],np.ndarray):
					fd.write(":{}={}".format(prefix+name,tree[name].tolist()))
				else:
					#DB>>
					#print prefix,name,tree[name]
					#<<DB
					fd.write(":{}={}".format(prefix+name,tree[name]))
		with open(methods['corelog']+".simdb","a") as fd:
			fcntl.lockf(fd, fcntl.LOCK_EX)
			now = datetime.now()
			fd.write("SIMDB/TIMESTAMP=(%04d,%02d,%02d,%02d,%02d,%02d)"%(now.year, now.month, now.day, now.hour, now.minute, now.second) )
			writetree(methods,fd,"/")
			fd.write("\n")
			#fcntl.lockf(fd, fcntl.LOCK_UN)
	if methods['gui']:
		if methods['gif']:
			plt.savefig(methods['gif'])
		else:
			plt.show()
	if not methods['noexit']:
		sys.exit(0)

Loading data, please wait...