#from mpi4py import MPI
from neuron import h, run
from math import ceil
#h.nrnmpi_init()


import sys
import numpy as np
#import pylab as pl
import os
affix = sys.argv[1]


v_init = -60
RUN = 0
GABARUN=0
GABA_ARRAY=0
PRINT = 1
IofV = 0
FGI = 0
RANDOM=0
PULSE = 0
STEADY = 1
GNAVARY = 0
SYNAPTIC = 0
SYNSTEPS = 0
SYNPULSE=0
DINGRUN = 0
DOCUMENT = 1
VOLTSTEPRUN=0
STEP2=0
STEPLEN=1000
OFFSET = 500
HOLD=0
VSTEP=0
JAC = 0
ARRAY = 0
NTHREADS=64
FAST=0
GRAPH =0
USEOLD=0
VHOLD=-60
KV4PULSE=0
NMDAFUNC=0
UNIFORM=1
COSTIM=0
GABA=0
PRINTALL=0
RASTER=0
BIFURCATION =0
VERBOSE=0

READVALS = 1
ALIGN=0
ALIAT=5000
custom_name = None

MORPHO='template_simple.hoc'
#MORPHO='160216_3.hoc'


rec_list = ['gka_hhb','q1_hhb','q2_hhb','ik_erg','cai']
argdict = {'maxi':100,'mini':10,'freq':4} # this is probably overkill
from paramsdict import p_init as pstable
from paramsdict import argdict

axis =9

EXTERN = 'array_parms.py'
affix2=''
for things in sys.argv:
	if len(things.split('='))>1:
		blah = things.split('=')[0]
		try:
			val = float(things.split('=')[1])
		except:
			val=0
		if blah in pstable:
			if blah == 'both':
				pstable['gnahh']=pstable['gnahh']*val
				pstable['gkhh']=pstable['gkhh']*val
			else:
				temp = things.split('=')
				try:
					exec('pstable[\'%s\'] = %s' % (temp[0], temp[1]))
				except:
					exec('pstable[\'%s\'] = \'%s\'' % (temp[0], temp[1]))
				
				#print 'pstable[\'%s\'] = %s' % (things.split('=')[0], things.split('=')[1])
			if VERBOSE or temp[0] in ['nmdaamp','gabaamp','idel','idur','regtype','number']: # control might need to be here, but had to kludge sort without it
				affix2+='_'
				affix2+=things.split('=')[0]
				try:
					affix2+='%.2e' % float(things.split('=')[1])
				except:
					affix2+=things.split('=')[1]
			continue
		else:
			try: 
				exec(things)
			except:
				exec('%s =\'%s\'' % (things.split('=')[0],things.split('=')[1]))
			if things.split('=')[0] in ['VSTEP','STEP2','STEPLEN','OFFSET','SYNAPTIC','COSTIM','UNIFORM','NMDAFUNC']:
				temp = things.split('=')
				affix2+='_%s%s' % (temp[0],temp[1])

if custom_name !=None:
	affix2+='_%s' % str(custom_name)
	
				
print(pstable['syntype'])

if DOCUMENT:
	path = os.getcwd() + '/' + affix+ '/'
	try:
		os.mkdir(path)
	except:
		pass
	#os.system('cp blockhh3.mod %s.' %path)
	#os.system('cp NaMarkov.mod %s.' %path)
	os.system('cp paramsdict.py %s.' %path)
	
	
	

h.load_file("nrngui.hoc")
#h.load_file("template_SNc.hoc") # for rebound

h.load_file('%s' % MORPHO)

if MORPHO!= 'template_simple.hoc':
	fp = open(MORPHO,'r')
	line = fp.readline()
	temp = line.split()
	TYPE=temp[1]
else:
	TYPE = 'SIMPLE'


from init_cell_depblock import init, to_shreds_you_say #depblock only has somatic assignments for single compartment
import new_geomnseg as geom

import Morpho_Func as MF

print(TYPE)
FIXED_DISTRO = None



def poisson_gaba_network(target,dur,amp=1e-3,dist=None,start=0,freq=10,ncells=1):
	# create a synapse and an event vector
	
	stim = h.NetStim()
	synapse = h.Exp2Syn(target(0.5))
	synapse.tau1=1
	synapse.tau2=6
	synapse.e = -70
	
	ncobj=h.NetCon(None, synapse)
	ncobj.delay=0 
	ncobj.weight[0] = amp
	
	
	if dist==None:
		spikes = []
		nspikes = 100
		for i in range(ncells):
			isi = np.random.poisson(freq,nspikes)
			cs = np.cumsum(isi)
			for things in cs:
				spikes.append(things)
		print(spikes)
		#spikes.sort() # arange spikes from all input cells in temporal order
		# might not be required
		event_vec = h.Vector(spikes)
		#event_vec.add(start)
		event_vec.sort()
	else:
		try:
			event_vec = dist
			event_vec.mul(1) # crude type check
		except:
			print('returning an empty event vector because')
			print(str(dist))
			print(' is not a neuron vector')
			event_vec = h.Vector()
	
	return ncobj, event_vec
	

class dopa:
	def __init__(self,i, parms, args=[0,1]):
		h.pop_section()
		self.p1 = args[0]
		self.p2 = args[1]
		self.internal = parms
		#print 'make'
		exec('self.nrn = h.%s()' % TYPE)
		#self.nrn = (h.SIMPLE())
		#self.nrn = MF.add_AIS(self.nrn,35,100,0.9)
		
		#print self.ais.L, self.ais.diam
		
		#to_shreds_you_say(self.nrn,parms['chop'])# create morphology template
		self.dt = 0.05 # 20khz
		#print 'init'
		init(self.nrn,self.internal) # from initcell - using dict parameters
		
		geom.geom_nseg(self.nrn)
		
		if len(args)>3:
			if args[2] in parms:
				exec('self.internal[\'%s\'] = %s' % (args[2], args[3]))
			if args[2] == 'gnmda':
				for s in self.nrn.all:
					s.insert('nmda')
					s.gnmdabar_nmda = self.internal['gnmda']
					s.cMg_nmda = self.internal['cmg']
					s.cafrac_nmda = parms['cafrac']
		#temp = np.arange(1000,4000,self.dt)
		self.srec = h.Vector()# spike times

		self.ntc = h.NetCon(self.nrn.soma(0.5)._ref_v,None,sec=self.nrn.soma)
		self.ntc.threshold = -20
		self.ntc.record(self.srec)
		#self.vrec.record(self.nrn.soma(0.5)._ref_v,self.tvec)#"""
		#print 'areas'
		for s in self.nrn.all:
			#s.insert('erg')
			#s.gergbar_erg = 20e-6
			if TYPE=='SIMPLE':
				s.L = parms['size']
				s.diam = 5
				if parms['gcat']> 0: # normally, gcat is not somatic
					s.insert('catchan')
					s.gcatbar_catchan = 100e-6*parms['gcat']
			s.dist_NaMark = parms['slow']
			#s.insert('erg')
			#if s == self.nrn.soma:
			#	s.L = self.internal['size']
			s.dist_NaMark = parms['slow']
			s.v = -50 # -40 starts in db
			s.cai = 1e-9
			s.hshift_NaMark = self.internal['na_hshift']
			s.nslope_hhb = self.internal['nslope']
			s.nshift_hhb = self.internal['nshift']
			s.kchip_hhb = self.internal['kchip']
			s.fchip_hhb = self.internal['fchip']
			if h.ismembrane('kca',sec=s):
				s.tausk_kca = 5
				s.km_kca *=1.0
			if h.ismembrane('calmark', sec=s):
				s.gcalbar_calmark *=0.0
				s.moff_calmark = -30
				s.mslope_calmark = 5.0
			if h.ismembrane('calhh', sec=s):
				#print 'derp'
				s.gcalbar_calhh *=1
				s.mhalf_calhh = -30
				s.mslope_calhh = 5.0
				s.pf_calhh = 0.6
				#s.coup_calmark = 0
			if h.ismembrane('canchan'):
				s.gcanbar_canchan = 50e-6
			#s.oerg_ergkin =0.5
			#s.vshift_erg = 0
			if h.ismembrane("hcn",sec=s):
				s.scale_hcn =3.0
				s.mhalf_hcn = -75.0
			if h.ismembrane('cabalstore',sec=s): # convert to instant buffering
				#s.MitoBuffer_cabalthin = 0.03
				for seg in s:
					seg.TotalBuffer_cabalstore = 0.03
					seg.SCALE_cabalstore = 1.0
					seg.shellfrac_cabalstore = min(0.1/seg.diam,1.0)
					seg.tog_cabalstore = 0 # syntax for calcium release via IP3 is present but not used. 
					#seg.
					seg.DCa_cabalstore = self.internal['dca'] # calcium diffusion constant (radial)
			#s.ashift_hhb = 0
			if h.ismembrane('catchan',sec=s):
				s.phi_h_catchan = 1.5
				s.shift_catchan = 0
				s.buff_catchan=parms['catbuff']
			if h.ismembrane('nmda'):
				s.cafrac_nmda = parms['cafrac']
		
		
		self.areas = []
		for s in self.nrn.all:
			#s.ashift_hhb = -10.0
			if s.nseg > 1:
				locarea = 0
				for seg in s:
					locarea += h.area(seg.x,sec=s)
			else:
				locarea = h.area(0.5,sec=s)
			self.areas.append(locarea)
	def __del__(self):
		if VERBOSE:
			print('killing', self.__class__.__name__)
		else:
			pass


def sine_func(t,freq=pstable['freq'],maxi=40,mini=10):
	if t==0:
		print(freq, maxi, mini)
	temp = (1.0+np.sin(2*np.pi*1e-3*freq*t))/2.0
	temp *= (maxi-mini)
	temp += mini
	return temp
	
def cos_func(t,freq=1,maxi=60,mini=2):
	if t==0:
		print(freq, maxi, mini)
	temp = (1.0+np.cos(2*np.pi*1e-3*freq*t))/2.0
	temp *= (maxi-mini)
	temp += mini
	return temp
	
def window_func(t,on=pstable['won'],off=pstable['woff'],fon=pstable['fon'],foff=pstable['foff']):
	if t%int(on+off)<on:
		temp=fon
	else:
		temp=foff
	return temp
	
def randomized_window_func(t,on=100,off=150,fmax=100,fmin=2,foff=0):
	if t%int(on+off)<on:
		temp=np.random.random()*(fmax-fmin)+fmin
	else:
		temp=foff
	return temp

#cell = dopa(0,[1,2],pstable)
#nc, stim = gen_netstim(sine_func,cell.nrn.soma,syn_type='gabaa',nspikes=1,duration=1e4)

def synpulse_func():
	
	cv = h.CVode()

	cell = dopa(0,pstable)
	rec_vectors = [h.Vector(),h.Vector()]


			
	iclamp = h.IClamp(cell.nrn.soma(0.5), sec=cell.nrn.soma)
	iclamp.dur = 1e9
	iclamp.amp = 0

	
	dt = 0.1
	rec_vectors = []
	for j in range(2):
		rec_vectors.append(h.Vector())
	i=2
	rec_vectors[0].record(iclamp,h._ref_t,dt,sec=cell.nrn.soma) # the first argument is just a pointer reference to the correct cell/compartment
	rec_vectors[1].record(iclamp,cell.nrn.soma(0.5)._ref_v,dt,sec=cell.nrn.soma)
	#print rec_list
	for things in rec_list:
		try:
			rec_vectors.append(h.Vector())
			exec('rec_vectors[%d].record(iclamp,cell.nrn.soma(0.5)._ref_%s,dt,sec=cell.nrn.soma)' % (i,things) ) in locals()
			i+=1
		except:
			print('failed to record %s'  % things	)
	
	syn_type = pstable['syntype']
	if syn_type not in ['nmda','gabab']:
		synapse = h.Exp2Syn(cell.nrn.soma(0.5))
		synapse.tau1=1
		synapse.tau2=6
		if syn_type=='gabaa':
			synapse.e = -70
		if syn_type=='ampa':
			synapse.e = 0
			synapse.tau2=3
	if syn_type == 'nmda':
		synapse = h.Exp2NMDA(target(0.5))
	if syn_type == 'gabab':
		#synapse = h.GABAb(target(0.5)) # this doesn't have right syntax
		synapse = h.Exp2Syn(target(0.5))
		synapse.tau1=50
		synapse.tau2=1000
		synapse.e = -90
			
			
	#ncobj.delay=2

	# make a netstim with appropriate parameters
	ns = h.NetStim()
	nc = h.NetCon(ns,synapse)
	
	if syn_type=='gabaa':
		nc.weight[0] = pstable['gabaamp']
	if syn_type=='nmda':
		nc.weight[0] = pstable['nmdaamp']
	if syn_type=='gabab':
		nc.weight[0] = pstable['gababamp']	
	
	if pstable['freq'] > 0:
		ns.interval = 1000.0/float(pstable['freq'])
		ns.number = pstable['number']
		ns.noise = pstable['noise']
		ns.start = pstable['idel']
	else:
		ns.interval = 1e9
		ns.number = 0
		ns.noise = 0
	# have that netstim be the source of a netcon
	
	# set the conductance of the gaba synapse with pstable
	
	# play
	
	cv.active(1)
	h.t=0
	h.finitialize()
	cv.solve(2*pstable['idel']+ns.number*ns.interval)
	
	fp= open('pulsedgaba_%s%s.dat' % (affix,affix2),'w')
	
	size = len(rec_vectors[0])
	nstates = len(rec_vectors)
	
	for i in range(size):
		for j in range(nstates):
			fp.write('%e  ' % rec_vectors[j][i])
		fp.write('\n' )
	fp.close()
	
	
	


if SYNPULSE:
	synpulse_func()

		

if GABARUN:
	from synapse_funcs import gaba_run, gen_netstim
	if not NMDAFUNC:
		gaba_run(rec_list,dopa,pstable=pstable,func=sine_func,funcargs={},args=[affix,affix2],ncells=int(pstable['ncells']),stype=pstable['syntype'])
	elif COSTIM:
		gaba_run(rec_list,dopa,pstable=pstable,func=window_func,funcargs={},args=[affix,affix2],ncells=int(pstable['ncells']),stype=pstable['syntype'],nmdafunc=window_func,nmdaargs={},nmdacells=1)	
	else:
		gaba_run(rec_list,dopa,pstable=pstable,func=window_func,funcargs={'fon':pstable['fon']},args=[affix,affix2],ncells=int(pstable['ncells']),stype=pstable['syntype'],nmdafunc=window_func,nmdaargs={'on':20000,'off':0,'fon':pstable['foff'],'foff':0},nmdacells=1) # yes foff	

def run_func(STEPLEN):
	cell = dopa(0,pstable)
	blah = h.IClamp(cell.nrn.soma(0.5))
	blah2 = h.IClamp(cell.nrn.soma(0.5))
	blah2.delay = 2*pstable['idel'] + 3*pstable['idur']#-4000
	blah2.amp = pstable['iamp']#- pstable['basal']
	blah2.dur = STEPLEN#pstable['idur']
	blah.delay = pstable['idel']
	blah.amp = pstable['iamp']#0e-3 # max at 3e-3
	blah.dur = pstable['idur']
    
	if STEADY or SYNAPTIC:
		dtime = pstable['idel']+STEPLEN-OFFSET
	else:
		dtime = 16000+STEPLEN-OFFSET

	if VSTEP:
		blah3 = h.SEClamp(cell.nrn.soma(0.5))
		blah3.amp1 = STEP2
		blah3.rs = 1e9 # turn off clamp
		blah3.dur1 = 1e9

	else:
		blah3 = h.IClamp(cell.nrn.soma(0.5))
		blah3.delay = dtime
		blah3.amp = 1e-3*STEP2
		blah3.dur = 200
	
	
	ampmax = pstable['iamp']
	ramp = []
	basal = pstable['basal']
	nc = h.NetCon(cell.nrn.soma(0.5)._ref_v,None,sec=cell.nrn.soma)
	nc.threshold = -10
	freqvec = h.Vector()
	nc.record(freqvec)
	times = np.arange(0,blah.dur+blah.delay,1)
	tvec = h.Vector(times)
	if not GNAVARY:
		for i in range(pstable['idel']):
			ramp.append(basal)
		for i in range(int(blah.dur)):
			if i < (blah.dur)/2.0:
				ramp.append(basal+2.0*ampmax*(i)/float(blah.dur)) # @ i = blah.dur/2 this is basal + ampmax = 75
				#if i%100 == 0:
				#	print i, basal+2*ampmax*(i)/float(blah.dur)
			else:
				ramp.append(basal+2.0*ampmax*(blah.dur-i)/float(blah.dur))
		for i in range(pstable['idel']):
			ramp.append(basal)

		blah.dur = blah.dur+2*blah.delay
		blah.delay = 2000*STEADY
		if not STEADY:
			ampvec= h.Vector(ramp)
			ampvec.play(blah._ref_amp, 1)	
		else:
			blah.amp = 0
			blah2.delay = pstable['idel']
	else:
		gnabar = cell.nrn.soma.gnabar_NaMark

		for i in range(int(blah.delay)):
			ramp.append(0*gnabar)
		for i in range(int(blah.dur)):
			if i < (blah.dur)/2.0:
				ramp.append(gnabar*(2*i/float(blah.dur)))
			else:
				ramp.append(2*gnabar*(blah.dur-i)/float(blah.dur))
		for i in range(pstable['idel']):
			ramp.append(0*gnabar)
		ampvec = h.Vector(ramp)
		for s in cell.nrn.all:
			s.dist_NaMark = 0
			ampvec.play(s(0.5)._ref_gnabar_NaMark,1)
	
	if SYNAPTIC:
		#delete all the synaptic stuff
		blah = None
		blah2 = None
		#blah3 = None  # keep this one
		gnmdamax = pstable['gnmda']
		step = []
		for i in range(pstable['idel']):#+pstable['idur']):
			step.append(0)
		
		for i in range(STEPLEN):#pstable['idur']):
			step.append(gnmdamax)
			
		for i in range(pstable['idel']):
			step.append(0)
		
		ampvec = h.Vector(step)
		#add an NMDA channel to cell
		for s in cell.nrn.all:
			s.insert('nmda')
			s.cafrac_nmda = pstable['cafrac']
			s.cMg_nmda = pstable['cmg']
			#s.tog_nmda = 0 # set in variable conductance mode
			ampvec.play(s(0.5)._ref_gnmdabar_nmda,1)
			
			
		# create a square pulse waveform
		
		# scale waveform by NMDA conductance
	
	
	if GABA:
		ncgaba, gaba_events = poisson_gaba_network(cell.nrn.soma,2000,start=1000,freq=10,ncells=1,amp=pstable['gabaamp'])

	if PRINT:
		
		#h.t = 0
		cv = h.CVode()
		cv.active(True)
		if JAC:
			cv.jacobian(1)
		#fp = open('pacing_%s.dat' %affix,'w')
		loc = cell.nrn.soma
		#while h.t < max(pstable['idel']-1600,500):
		#	h.fadvance()
		#while h.t < max(pstable['idel']-500,500):
		#	h.fadvance()
		#	fp.write('%e  %e  0  %e  %e  %e  %e  %e\n' %(h.t, loc.v, loc.o1_NaMark, loc.c1_NaMark, loc.c2_NaMark, loc.i1_NaMark,loc.i2_NaMark))
		#fp.close()
		h.finitialize()
		if not STEADY:
			if not SYNAPTIC:
				fp = open('ramp_%s_%f.dat' % (affix+affix2,pstable['iamp']),'w')
				while h.t < pstable['idel']-1000:
					h.fadvance()
			else:
				fp = open('nmda_%s_%.1f.dat' % (affix+affix2,1e6*pstable['gnmda']),'w')
				if HOLD:
					vc = h.SEClamp(cell.nrn.soma(0.5))
					vc.dur1 = pstable['idel']
					vc.dur2 = 0
					vc.dur3 = 0
					vc.amp1 = VHOLD
				h.finitialize()
		else:
			fp = open('steady_%s%s.dat' % (affix,affix2),'w')
			if HOLD:
				vc = h.SEClamp(cell.nrn.soma(0.5))
				vc.dur1 = pstable['idel']
				vc.dur2 = 0
				vc.dur3 = 0
				vc.amp1 = VHOLD
			h.finitialize()
			
		if VSTEP:
			def vfunc(val):
				blah3.rs = val
				return
			#cv.event(dtime,vfunc(1e-4),blah3,None)
			#cv.event(dtime+200,vfunc(1e9),blah3,None)
		vold = loc.v
		told = float(h.t)-0.02
		dvdtold = 0
		offset = 0
		FaA = 0

		while h.t < 3*pstable['idel']+3*pstable['idur']+STEPLEN:
			h.fadvance()
			if h.t > 2*pstable['idel']+STEPLEN and (SYNAPTIC or STEADY):
				break
			if dtime+3000 > h.t > dtime and VSTEP:
				blah3.rs=1e-3
			elif VSTEP:
				blah3.rs=1e9
			else:
				pass
				################### yeet this into its own function #################
			if JAC and (STEADY or SYNAPTIC) and h.t > pstable['idel']+STEPLEN - 100:
				pstates = h.Vector()
				cv.states(pstates)
				import neuron_jacobian as nj
				jac = h.Matrix()
				nj.get_jacobian(cv,jac,relstep=1e-6)
				temp = h.Vector()
				jac_p = []
				sz = int(jac.ncol())
				#print jac, jac.ncol()
				for i in range(sz):
					jac.getcol(i,temp)
					jac_p.append(temp.to_python())
					a = h.ref('')
					cv.statename(i,a)
					j=0
					for things in temp.x:
						b = h.ref('')
						cv.statename(j,b)
						print(a, b, things)
						j+=1
						
				jac_n = np.array(jac_p)
				jac_n.resize([sz,sz])
				
				
				from scipy import linalg
				a,b = linalg.eig(jac_n)
				print(a)
				quit()
				##################### end yeet #######################################
			if h.t - told > 0:
				dvdt = (loc.v-vold)/(float(h.t)-told)
			else:
				dvdt = dvdtold
			dvdtold = dvdt
			vold = loc.v
			told = float(h.t)
			rat = pstable['ratio']

			if not SYNAPTIC:
				if ALIGN and not FaA: # align to fixed point on AP to compare AHP, shape between models.
					if h.t < ALIAT:
						continue
					elif h.t > ALIAT and not FaA and loc.v <-40:
						continue
					elif ALIGN and loc.v > -40 and not FaA: # FaA = First spike after alignment offset
						offset = h.t
						FaA = 1		
				time = h.t-offset	
				#if (2000<h.t < 4000) or not STEADY:
				gkainf = pow(loc.pinf_hhb,3)*loc.qinf_hhb
				gka = pow(loc.p_hhb,3)*loc.q_hhb
				if h.t%1000 < h.dt:
					print(h.t)
				fp.write('%e  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e  ' %(time, loc.v, blah.amp, loc.o1_NaMark, loc.c1_NaMark, loc.c2_NaMark, loc.i1_NaMark,rat*loc.i2_NaMark,loc.i1_NaMark + loc.i2_NaMark,dvdt,1.0/3.14159*1e6*(loc.ina+loc.ik+loc.ica+loc.ihcn_hcn)-1.0/3.14159*1e8*blah2.amp/h.area(0.5,loc)))
				fp.write('%e  %e  %e  %e  %e  %e  %e  %e' %(loc.ina_NaMark, loc.ikhh_hhb, loc.ika_hhb, loc.ik_kca, loc.ik, gka, gkainf, loc.ika_hhb+loc.ina_NaMark+loc.ik_leak+loc.ina_leak))
			else:
				#if (2000<h.t < 4000):
					fp.write('%e  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e  ' %(h.t, loc.v, loc.inmda_nmda+loc.ica_nmda, loc.o1_NaMark, loc.c1_NaMark, loc.c2_NaMark, loc.i1_NaMark,rat*loc.i2_NaMark, loc.i1_NaMark + loc.i2_NaMark, loc.m_kca,dvdt,1.0/3.14159*1e6*(loc.ina+loc.ik+loc.ica+loc.ihcn_hcn+loc.inmda_nmda)))
			fp.write('\n')
		fp.close()
		interval = -1
		minterval = 1e9
		if not SYNAPTIC:
			for t in freqvec:
	
				if interval < 0:
					try:
						interval = freqvec[1] - freqvec[0]
						told = freqvec[0]
						continue
					except:
						print('no spikes')
						break
				else:
					interval = t - told
					if t < 6000 and pstable['idur'] > 1000 and not STEADY:
						print(1000/(t-told), ampvec[int(t+told)/2])
					if interval < minterval:
						minterval = interval
					told = t
		#print 1000.0/minterval

	else:
		h.finitialize()
		h.t=0
		from neuron import gui
		return cell
	
	
	
	return 0
	
if RUN:
	cell = run_func(STEPLEN)
	#print cell.nrn.soma.gkhhbar_hhb

def dingrun_fun():
	cells = [dopa(0,pstable),dopa(0,pstable)]
	i=0
	for cell in cells:
		for s in cell.nrn.all:
			#s.L = pstable['size']
			s.dist_NaMark = pstable['slow']
			if i:
				s.gnabar_NaMark = 0
				#s.hshift_NaMark = -pstable['na_shift']"""
		i+=1
	blah0 = h.VClamp(cells[0].nrn.soma(0.5))
	blah0.dur[0] = 1e9
	blah1 = h.VClamp(cells[1].nrn.soma(0.5))
	blah1.dur[0] = 1e9
	
	recvecs = [h.Vector(),h.Vector()] 
	recvecs[0].record(blah0._ref_i,0.1)
	recvecs[1].record(blah1._ref_i,0.1)
	
	volpy = []
	
	for j in range(10000):
		volpy.append(-70)
	
	for i in range(10):
		for j in range(1000):
			volpy.append(-70)
		for j in range(30):
			volpy.append(0)
			
	for j in range(10000):
		volpy.append(-70)
		
	volh = h.Vector(volpy)
	volh.play(blah0._ref_amp[0],0.1)
	volh.play(blah1._ref_amp[0],0.1)
	cv = h.CVode()
	h.finitialize()
	while h.t < 0.1*len(volh):
		h.fadvance()
		#print h.t, cells[0].nrn.soma.v, blah0.i, blah1.i, cells[0].nrn.soma.i2_NaMark
		
	recvecs[0].sub(recvecs[1])
	
	for i in range(len(volpy)):
		print(i*0.1, recvecs[0].x[i]/recvecs[0].min(), volpy[i])#"""
	return 0

if DINGRUN:
	dingrun_fun()

if ARRAY:
	from array_funcs import array_func
	array_func(EXTERN,dopa,pstable,args=[affix,affix2])	

if SYNSTEPS:
	from array_funcs import synsteps
	synsteps(STEPLEN)

if GABA_ARRAY:
	from array_funcs import gaba_array
	if NMDAFUNC:
		if COSTIM:
			gaba_array(EXTERN,dopa,window_func,UNIFORM=UNIFORM,pstable=pstable,args=[affix,affix2],funcargs={},ncells=int(pstable['ncells']),syn_type=pstable['syntype'],nmdafunc=window_func,nmdaargs={},nmdacells=int(pstable['ncells']))
		else:
			
			#print pstable['fon']
			#quit()
			gaba_array(EXTERN,dopa,window_func,UNIFORM=UNIFORM,pstable=pstable,args=[affix,affix2],funcargs={'fon':pstable['fon']},ncells=int(pstable['ncells']),syn_type=pstable['syntype'],nmdafunc=window_func,nmdaargs={'on':20000,'off':0,'fon':pstable['foff'],'foff':0},nmdacells=1) # foff for balanced state
	else:
		gaba_array(EXTERN,dopa,window_func,pstable=pstable,args=[affix,affix2],funcargs={},ncells=int(pstable['ncells']),syn_type=pstable['syntype'])
	
if RANDOM:
	from array_funcs import random_cells
	random_cells(STEPLEN)
	
def clean_vec(vector,thresh,last=0):
	vlen = len(vector)
	for i in range(vlen-2):
		j=i+1
		if vector[j] != 0:
			delta = abs((vector[j+1]-vector[j])/vector[j])
		else:
			delta= 1
		#print delta, vector[j],vector[j+1],vector[j-1]
		if delta > thresh:
			if not last:
				vector.x[j] = vector[j+1]/2.0+vector[j-1]/2.0
			elif j>10:
				vector.x[j+1] = vector[j]


def kv4_pulse():
	#does Kv4 experiments from Tang et al 2013 - except with difference current on full cell.
	actsteps = np.arange(-80,35.1,5)
	timesteps = [5,10,25,50,100,200,500,1000]
	

	
	aln = len(actsteps)
	tln = len(timesteps)
	
	ahold = -100
	thold = -80
	
	timecond = -40
	timetest = -40
	pstable['gnahh'] = 0
	#pstable['gcal'] = 0
	#pstable['gkhh'] = 0
	
	warm = 3000 # warm up time prior to pulse
	condpulse = 1000 # for voltage response
	spikepulse = 1000
	cells = []
	clamps = []
	recvectors = []

	h.dt = 0.1 # evaluation dt
	dt = 0.1 # recording dt
	start = 0.2
	t_record = np.arange(start,1000,dt) # generic t_record vector - offset

	atvec = h.Vector(t_record)
	atvec.add(warm)
	atvecs = []
	toffs = []
	for i in range(tln):	
		atvecs.append(h.Vector(t_record))
		t_off = timesteps[i]#1000.0*pow(2.0,timesteps[i])
		toffs.append(t_off)
		atvecs[i].add(warm+2000+t_off)
		#print atvecs[i].min(), atvecs[i].max()
	
	
	i=0 # index of current cell cause we are doing both runs at once for some reason
	if PRINTALL:
		vrecvectors=[]
	for j in range(2*aln):
		h.pop_section() # kludge fix to section stack overflow

		cells.append(dopa(0,pstable))

		if MORPHO=='template_simple.hoc':
			cells[i].nrn.soma.L = 0.1 # small cells are fast cells
		clamps.append(h.SEClamp(cells[i].nrn.soma(0.5)))
		clamps[i].dur1 = warm
		clamps[i].dur2 = spikepulse
		clamps[i].dur3 = warm
		clamps[i].amp1 = ahold
		if j >= aln:
			clamps[i].amp2 = actsteps[j-aln]
			for s in cells[i].nrn.all:
				s.gkabar_hhb = 0
		else:
			clamps[i].amp2 = actsteps[j]
		clamps[i].amp3 = ahold
		recvectors.append(h.Vector())
		#recvectors[i].record(clamps[i]._ref_i,atvec)
		if PRINTALL:
			vrecvectors.append(h.Vector())
			recvectors[i].record(clamps[i]._ref_i,dt)
			vrecvectors[i].record(cells[i].nrn.soma(0.5)._ref_v,dt)
		else:
			recvectors[i].record(clamps[i]._ref_i,atvec)
		i+=1
	vplay_vecs = []
	tplay_vecs = []
	
	
	for j in range(2*tln):
		cells.append(dopa(0,pstable))
		clamps.append(h.SEClamp(cells[i].nrn.soma(0.5)))

		if MORPHO=='template_simple.hoc':
			cells[i].nrn.soma.L = 0.1 # small cells are fast cells
		#prestep
		clamps[i].dur1 = warm #toffs[j%tln]
		#clamps[i].dur2 = timesteps[j%tln]
		clamps[i].dur2 = 500
		clamps[i].dur3 = 1e9 # record time
		if j < tln:
			tempt = [warm+500,warm+1000,warm+1000+timesteps[j%tln]]
			tempv = [timetest,thold,timetest]
			vplay_vecs.append(h.Vector(tempv))
			tplay_vecs.append(h.Vector(tempt))
		clamps[i].amp1 = timecond # -40
		clamps[i].amp2 = thold # -80
		clamps[i].amp3 = timetest # -40
		
		vplay_vecs[j%tln].play(clamps[i]._ref_amp3,tplay_vecs[j%tln])
		
		
		
		
		
		if j >= tln:
			for s in cells[i].nrn.all:
				s.gkabar_hhb = 0
			
		recvectors.append(h.Vector())
		if PRINTALL:
			vrecvectors.append(h.Vector())
			recvectors[i].record(clamps[i]._ref_i,dt)
			vrecvectors[i].record(cells[i].nrn.soma(0.5)._ref_v,dt)
		else:
			recvectors[i].record(clamps[i]._ref_i,atvec)
		
		#print( clamps[i].dur2+clamps[i].dur1, atvecs[j%tln].min())
		i+=1
		

	pc = h.ParallelContext()
	pc.nthread(NTHREADS)
	h.finitialize()
	h.t = 0
	TSTOP = max(warm+spikepulse +2000+ max(toffs),2*warm+spikepulse)
	while h.t < TSTOP:
		if h.t%100 <h.dt:
			print( h.t)
		h.fadvance()
	#pc.psolve(TSTOP) # this does non-desirable things
	
	for i in range(aln):
		recvectors[i].sub(recvectors[i+aln])
		clean_vec(recvectors[i],0.1) # difference curren
	

	#for i in range(tln):
	#	recvectors[i+2*aln].sub(recvectors[i+tln+2*aln])
		#clean_vec(recvectors[i+2*aln],0.3,last=1)
		#clean_vec(recvectors[i+2*aln],0.1)
		#clean_vec(recvectors[i+2*aln],0.1)
	
	#pc.nthread(1)
	dln = len(recvectors[0]) # should be int(spikepulse/dt)
	maxmcur = 0
	
	from scipy.optimize import curve_fit
	
	def biexp(t,a,b,c,d):
		temp = d+a*pow(1.0-np.exp(-t/b),3)*np.exp(-t/c)
		return temp
		
	bounds = ([0,0.1,0.1,-1],[100,2000,2000,1])
	initial = [0.003,1,1,0]

	fp=open('Kv4_activation_%s%s.dat' % (affix,affix2),'w')
	time_constants = []

	for i in range(aln):
		temp = recvectors[i].to_python() # these are the ones that contain difference current
		to_fit = np.array(temp[0:len(temp)]) # (msteps[i]-50.0) # chop off ends that may have RC effects
		temp = None # free memory (not really needed tbh)
		
		xdata = np.arange(0,len(to_fit)) # signs are correct
		for k in range(len(xdata)):
			xdata[k] *= dt # time
			xdata[k] += start
		
		#try:
			#popt, cov = curve_fit(biexp,xdata,to_fit,p0=initial,bounds=bounds)
			#print msteps[i], popt
		#except:
		#	pass
		#time_constants.append([popt[0],popt[1],popt[2]])
	tdata = [dt*j for j in range(dln)]
	ydata = recvectors[aln-1]
	try:
		popt, cov = curve_fit(biexp,tdata,ydata,p0=initial,bounds=bounds)
		#print(popt[0], popt[1], popt[2])
	except:
		print('curvefit failed')
		quit()
		
	for j in range(dln):	# print out curves	
		#if mt1 > j or j > mt2:
		#	continue
		fp.write('%e  ' % (dt*j))	# this is inefficient
		for i in range(aln):
			fp.write('%e ' % (recvectors[i][j]))#, biexp(dt*j,time_constants[i][0],time_constants[i][1],time_constants[i][2])))\
		fp.write('%e  ' %biexp(dt*j,popt[0],popt[1],popt[2],popt[3]))
		fp.write('\n')
	
	fp.close()	



	tpeak = None
	locmini = []
	for i in range(aln):
		if actsteps[i] == -90:
			continue
		recvectors[i].div(actsteps[i]+90) # convert to conductance
		temp = recvectors[i].max(4,len(recvectors[i])-4) # conductance is positive definite
		locmini.append(temp)
		if actsteps[i] == timetest and temp != 0:
			tpeak = temp
		if temp > maxmcur:
			maxmcur = temp
	for i in range(aln):
			print(actsteps[i], locmini[i]/maxmcur)
	
	
	
	if PRINTALL:
		for i in range(tln):
			fp = open('recov_%s%s_%e.dat'%(affix,affix2,timesteps[i]),'w')
			for j in range(len(recvectors[2*aln])):
				if j < 100:
					continue
				t = j*dt*1e-3
				cur = recvectors[2*aln+i].x[j]
				v = vrecvectors[2*aln+i].x[j]*1e-3
				fp.write('%e  %e  %e\n' % (t,cur,v))
			
			
			fp.close()
	else:
		fp=open('kv4_tshift_%s%s.dat' %( affix,affix2),'w')
		
	
	
		
		fp.write('#ramp  ') # commented outline for xmgrace
		for i in range(tln):
			recvectors[i+2*aln].div(timetest+90.0) # normalize
			#if tpeak != None:
				#recvectors[i+2*aln].div(tpeak)
			fp.write('%d  ' % toffs[i])
		
	
		print( tln, aln, len(recvectors))
		for j in range(len(recvectors[2*aln])):
			
			fp.write('%e  ' % (start+dt*j))
			for i in range(tln):
				fp.write('%e  ' % recvectors[i+2*aln][j])
				if j==0:
				  print( toffs[i%tln], recvectors[i+2*aln].max())
			fp.write('\n')
	
	
		fp.close()
	
	
if KV4PULSE:
	kv4_pulse()


def pulse_func(): # simulates activation/ (fast) inactivation experiments on nucleated patch
	msteps = np.arange(-120,41,5)
	hsteps = np.arange(-120,41,5)
	mln = len(msteps)
	hln = len(hsteps)
	hhold = -100
	mhold = -100
	htest = 0
	warm = 1000
	condpulse = 50
	spikepulse = 5
	cells = []
	clamps = []
	recvectors = []
	i=0
	h.dt = 0.02
	mtrecord = np.arange(0,warm+condpulse+spikepulse+100, 0.01)
	mtvec = h.Vector(mtrecord)
	htvec = h.Vector(mtrecord)
	#print mtvec.x[0]
	
	mt1 = warm
	mt2 = warm+spikepulse
	
	ht1 = warm+condpulse
	ht2 = warm+condpulse+spikepulse
	
	mt1*=100
	mt2*=100
	ht1*=100
	ht2*=100
	
	mt1+=5
	mt2-=5
	
	ht1+=5
	ht2-=5
	
	
	for j in range(2*mln):
		h.pop_section() # kludge fix to section stack overflow
		cells.append(dopa(0,pstable))
		cells[i].nrn.soma.L = 0.1 # required to have space clamp
		clamps.append(h.SEClamp(cells[i].nrn.soma(0.5)))
		clamps[i].dur1 = warm
		clamps[i].dur2 = spikepulse
		clamps[i].dur3 = condpulse
		clamps[i].amp1 = mhold
		if j >= mln:
			clamps[i].amp2 = msteps[j-mln]
			for s in cells[i].nrn.all:
				s.gnabar_NaMark = 0
				if h.ismembrane('FastNaMark'):
					s.gnabar_FastNaMark=0
		else:
			clamps[i].amp2 = msteps[j]
		clamps[i].amp3 = mhold
		recvectors.append(h.Vector())
		recvectors[i].record(clamps[i]._ref_i,mtvec)
		i+=1

	for j in range(2*hln):
		cells.append(dopa(0,pstable))
		clamps.append(h.SEClamp(cells[i].nrn.soma(0.5)))
		cells[i].nrn.soma.L = 0.1
		clamps[i].dur1 = warm
		clamps[i].dur2 = condpulse
		clamps[i].dur3 = spikepulse
		clamps[i].amp1 = hhold
		if j >= hln:
			clamps[i].amp2 = hsteps[j-hln]
			for s in cells[i].nrn.all:
				s.gnabar_NaMark = 0
				if h.ismembrane('FastNaMark'):
					s.gnabar_FastNaMark=0
		else:
			clamps[i].amp2 = hsteps[j]
		clamps[i].amp3 = htest
		recvectors.append(h.Vector())
		recvectors[i].record(clamps[i]._ref_i,htvec)
		i+=1
		

	pc = h.ParallelContext()
	pc.nthread(NTHREADS)
	#cv = h.CVode()
	h.finitialize()
	h.t = 0
	TSTOP = warm + condpulse+spikepulse
	while h.t < warm + condpulse+spikepulse+100:
		h.fadvance()
	#pc.psolve(TSTOP) # this does non-desirable things
	
	for i in range(mln):
		recvectors[i].sub(recvectors[i+mln])
	

	for i in range(hln):
		recvectors[i+2*mln].sub(recvectors[i+hln+2*mln])
	
	#pc.nthread(1)
	dln = len(recvectors[0])
	maxmcur = 0
	
	from scipy.optimize import curve_fit
	
	def biexp(t,a,b,c):
		temp = a*pow(1.0-np.exp(-b*t),3)*np.exp(-c*t)
		return temp
		
	bounds = ([-2,0,0],[0,1e9,1e9])
	initial = [-1e-3,1,1]

	fp=open('mdiff_%s.dat' % affix,'w')
	time_constants = []

	for i in range(mln):
		temp = recvectors[i].to_python()
		#print temp[mt1-4:mt2]
		#quit()
		to_fit = np.array(temp[mt1:mt2]) # (msteps[i]-50.0)
		temp = None
		
		xdata = np.arange(4.0,len(to_fit)+4.0)
		xdata *= 0.01 # time
		
		try:
			popt, cov = curve_fit(biexp,xdata,to_fit,p0=initial,bounds=bounds)
			#print msteps[i], popt
		except:
			pass
		time_constants.append([popt[0],popt[1],popt[2]])

	for j in range(dln):

		
		
		if mt1 > j or j > mt2:
			continue
		#if abs(j-100*warm) <20 or abs(j-100*(warm+spikepulse))<20: # ignore data within 0.2 ms?
		#	continue
		fp.write('%e  ' % (0.01*j))	
		for i in range(mln):
			#if i==1 and j>warm*100:
				#print j, recvectors[i].x[j], recvectors[i+mln].x[j]
			fp.write('%e  %e  ' % (recvectors[i][j], biexp(0.01*j-1000,time_constants[i][0],time_constants[i][1],time_constants[i][2])))
		fp.write('\n')
	fp.close()
	for i in range(mln):
		temp = recvectors[i].min(mt1,mt2)/(msteps[i]-50.0)
		if abs(temp) > abs(maxmcur) and i > 0:
			maxmcur = temp
	for i in range(mln):
		if mln != hln:
			print( msteps[i], abs(recvectors[i].min(mt1,mt2)/(maxmcur*(msteps[i]-50.0))) )
	
	fp=open('hdiff_%s.dat' % affix,'w')
	maxhcur = 0
	for i in range(hln):
		temp =recvectors[i+2*mln].min(ht1,ht2)
		if abs(temp) > abs(maxhcur):
			maxhcur = temp
	for i in range(hln):
		if mln != hln:
			print( hsteps[i], recvectors[i+2*mln].min(int((condpulse)*50),len(recvectors[i+2*mln])-1)/maxhcur)
	for j in range(len(recvectors[2*mln])):
		if j < mt1 or j > mt2:
			continue
		fp.write('%e  ' % (0.01*j))
		for i in range(hln):
			fp.write('%e  ' % recvectors[i+2*mln][j])
		fp.write('\n')
	fp.close()
	if hln==mln:
		for i in range(hln):
			print(hsteps[i], recvectors[i].min(mt1,mt2), recvectors[i+2*mln].min(ht1,ht2), abs(recvectors[i].min(mt1,mt2)/(maxmcur*(msteps[i]-50.0))),recvectors[i+2*mln].min(ht1,ht2)/maxhcur, 1.0/time_constants[i][1], 1.0/time_constants[i][2])
	# for each pair find difference and max (according to abs) current
	
if PULSE:
	pulse_func()


def voltsteprun():
	nsteps = pstable['number']
	von = -40
	voff = pstable['ecl']
	
	cv = h.CVode()

	cell = dopa(0,pstable)
	rec_vectors = [h.Vector(),h.Vector()]


			
	vclamp = h.SEClamp(cell.nrn.soma(0.5), sec=cell.nrn.soma)
	vclamp.dur1=1e9
	#vclamp.rs=1e-6

	
	dt = 0.1
	rec_vectors = []
	for j in range(3):
		rec_vectors.append(h.Vector())
	
	i=3
	rec_vectors[0].record(vclamp,h._ref_t,dt,sec=cell.nrn.soma) # the first argument is just a pointer reference to the correct cell/compartment
	rec_vectors[1].record(vclamp,cell.nrn.soma(0.5)._ref_v,dt,sec=cell.nrn.soma)
	rec_vectors[2].record(vclamp._ref_i,dt)
	#print rec_list
	for things in rec_list:
		try:
			rec_vectors.append(h.Vector())
			exec('rec_vectors[%d].record(vclamp,cell.nrn.soma(0.5)._ref_%s,dt,sec=cell.nrn.soma)' % (i,things) ) in locals()
			i+=1
		except:
			print ('failed to record %s'  % things	)
	
	
	vclamp_array = [von,voff]
	tclamp_array = [0,4000]
	t = 4000
	
	for i in range(nsteps):
		vclamp_array.append(von)
		vclamp_array.append(voff)
		t+=pstable['won']		
		tclamp_array.append(t)
		t+=pstable['woff']
		tclamp_array.append(t)
	
	vclamp_array.append(von)
	tclamp_array.append(t)	
	t+=2000

	vclamp_array.append(von)
	tclamp_array.append(t)


	vclamp_vector = h.Vector(vclamp_array)
	tclamp_vector = h.Vector(tclamp_array)
	
	vclamp_vector.play(vclamp._ref_amp1,tclamp_vector)
	
	
	h.t=0
	h.finitialize()
	cv.solve(t)
	
	fp= open('vsteps_%s%s.dat' % (affix,affix2),'w')
	
	size = len(rec_vectors[0])
	nstates = len(rec_vectors)
	
	for i in range(size):
		for j in range(nstates):
			fp.write('%e  ' % rec_vectors[j][i])
		fp.write('\n' )
	fp.close()
	
	
	
if VOLTSTEPRUN:
	voltsteprun()
		
	

if FGI:	
	Iarray1 = np.arange(-8e-3,5e-3,1e-4)
	Iarray2 = np.arange(6e-3,101e-3,1e-3)
	Iarray = np.concatenate([Iarray1,Iarray2])

	iln = len(Iarray)
	cv = h.CVode()
	pc = h.ParallelContext()
	pc.nthread(64)
	pcnt = int(pc.nthread())
	nspikes = []
	freq = []

	nbatches = int(ceil(len(Iarray)/float(pcnt)))
	#print nbatches
	breakcon = 0
	anit = 0
	prespikes = 0
	WARMUP = 2000
	RUNTIME = 8000
	for i in range(nbatches):
		if breakcon:
			break
		ncs = []
		rvs = []
		cells = []
		clamps = []
		for j in range(pcnt):
			h.pop_section() # kludge fix to section stack overflow
			# without this neuron thinks all cells sections are part of an ever expanding subtree
			cells.append(dopa(0,pstable))
			clamps.append(h.IClamp(cells[j].nrn.soma(0.5)))
			clamps[j].dur = 1e9
			clamps[j].delay = WARMUP # 0 for 'from rest'
			for s in cells[j].nrn.all:
				#if s in cells[j].nrn.somatic:
				#	s.L = pstable['size']
				#else:
				#	s.L = 1e-6
				s.dist_NaMark = pstable['slow']
				s.v = -45 #-60
			try:
				clamps[j].amp = Iarray[j+pcnt*i]
			except: # klude to kill cells out range
				#print 'nope'
				cells[j] = None
				clamps[j] = None
				break
			ncs.append(h.NetCon(cells[j].nrn.soma(0.5)._ref_v,None,sec=cells[j].nrn.soma))
			ncs[j].threshold = -20
			rvs.append(h.Vector())
			ncs[j].record(rvs[j])
		
		 
		h.finitialize()
		h.frecord_init()
		h.t = 0
		#while h.t < RUNTIME:
		#	h.fadvance()
		pc.psolve(RUNTIME+WARMUP)
		#print 'batch %d of %d done' % (i, nbatches)
		for j in range(len(cells)):
			if j == len(rvs):
				break
			if len(rvs[j].x) < 3:
				tfreq1 = 0
				tfreq2 = 0
			elif rvs[j].x[-2] < WARMUP:
				tfreq1 = 0
				tfreq2 = 0
			else:
				tfreq1 = 1000.0/(rvs[j].x[-1] - rvs[j].x[-2])
				tfreq2 = 1000.0/(rvs[j].x[-2] - rvs[j].x[-3])
			
			freq.append(tfreq1)
			if not anit: # only need to check number of spike in warmup once as it is identical
				for k in range(len(rvs[j])):
					if rvs[j].x[k] > WARMUP:
						anit = 1
						prespikes = k
					break
			tnum = len(rvs[j]) -  prespikes# number of spikes after square wave starts
			nspikes.append(tnum)
			if tfreq1 > 0:
				if 2*(rvs[j].x[-1] - rvs[j].x[-2]) < RUNTIME-rvs[j].x[-1] or (rvs[j].x[-1] - rvs[j].x[-2]) > 5*(rvs[j].x[1]-rvs[j].x[0]):
					#print (rvs[j].x[-1] - rvs[j].x[-2]),  RUNTIME-rvs[j].x[-1], rvs[j].x[1] - rvs[j].x[0], rvs[j].x[-1] - rvs[j].x[-2]
					continue

			print( clamps[j].amp, )
			print( max(tfreq1,tfreq2), min(tfreq1,tfreq2), tnum, )
			if len(rvs[j].x) > 0:
				print( rvs[j].x[-1])
			else:
				print( -1)
			

		# get number of spikes
		# get amplitude of last spike? not yet
		# get time of last spike
		# if 6000 - tlast > 2*t2ndtolast - tlast assume block
	

			
	

if IofV:
	v_array = np.arange(-85,45,1)
	import nullcline_funcs as nf
    
	f, vout, iout, cells = nf.IV_curve(dopa,args=(0,pstable),volt_bounds=[-70,-30],voltstep=0.1)
    
    
	ln = len(vout)
    
	try:
		fp = open('%s/IofV_%s_%s.dat' % (affix, affix,affix2),'w')
	except:
		fp = open('IofV_%s_%s.dat' % (affix,affix2),'w')
		
	for i in range(ln):
		iout[i] -= pstable['iamp']
		fp.write('%e  %e  %e\n' % (vout[i],iout[i],f(vout[i])))
		
	fp.close()
			

	
if BIFURCATION:
	import nullcline_funcs as nf
	

	if SYNAPTIC:
		iarray = np.arange(0,200e-6, 1e-6)
	else:
		iarray = np.arange(-20e-3, 50e-3, 1.0e-4)
	
	cells = []
	
	
	# create IV curve or equivalent
	if SYNAPTIC:
		if USEOLD and (os.path.exists('%s/nmda_v_nullcline_%s.dat' %(affix,affix)) or os.path.exists('nmda_v_nullcline_%s.dat' %(affix))):
			gvpairs = []
			try:
				fp = open('%s/nmda_v_nullcline_%s.dat' %(affix,affix),'r')
			except:
				fp = open('nmda_v_nullcline_%s.dat' %(affix),'r')
			
			for lines in fp:
				temp = lines.split()
				g = float(temp[0])
				v = float(temp[1])
				print( g, v)
				gvpairs.append([g,v])
			
			fp.close()		
			
		else: # this part is time consuming so only do it once per problem if possible
			gvpairs = []
			for conductance in iarray:
				if cells ==[]:
					test, cells = nf.v_root(dopa, args=(0,pstable),volt_bounds=[-50,-20],voltstep=1)
				else:
					for cell in cells:
						for s in cell.nrn.all:
							s.gnmdabar_nmda = conductance
					test, cells = nf.v_root(dopa, cells=cells, args=(0,pstable),volt_bounds=[-50,-20],voltstep=1)
				print( conductance, test[0])
				gvpairs.append([conductance,test[0]])
			
			fp = open('nmda_v_nullcline_%s.dat' % affix,'w')
			for pairs in gvpairs:

				fp.write('%e  %e\n' % (pairs[0],pairs[1]))
			fp.close()

	else:
		f, varray, curarray, cells = nf.IV_curve(dopa,args=(0,pstable))
	    	    
		#print cells
		vb = [-80,-20]
		fs, extrema, iex, cells = nf.VI_curves(dopa,args=(0,pstable),volt_bounds=vb)
		# last bit allows for sampling over a uniform current step
	
	
	ncells = len(cells)
	#print ncells
	
	if not SYNAPTIC:
		stim = []

		
	kicks = []
	vecs = []
	i = 0
	for current in iarray:
		if i >= ncells:
			cells.append(dopa(0,pstable))
			ncells+=1
		vecs.append(h.Vector())
		vecs[i].record(cells[i].nrn.soma(0.5)._ref_v, 0.1)
		kicks.append(h.IClamp(cells[i].nrn.soma(0.5)))
		kicks[i].amp = -50e-3
		kicks[i].dur = 50
		kicks[i].delay = 50
		if SYNAPTIC:
			for s in cells[i].nrn.all:
				if not h.ismembrane('nmda',sec=s):
					s.insert('nmda')
				s.gnmdabar_nmda = current
				s.cMg_nmda = pstable['cmg']
				s.cafrac_nmda=pstable['cafrac']
				s.v = -60
				kicks[i] = None
				kicks[i] = h.SEClamp(cells[i].nrn.soma(0.5))
				kicks[i].amp1 = -50
				kicks[i].dur1 = 1000
				kicks[i].dur2 = 0
				kicks[i].dur3 = 0
		else: # previous currents SHOULD only exist in function
			stim.append(h.IClamp(cells[i].nrn.soma(0.5)))
			stim[i].amp = current
			if current <= 0:
				kicks[i] = None
				kicks[i] = h.SEClamp(cells[i].nrn.soma(0.5))
				kicks[i].amp1 = -50
				kicks[i].dur1 = 1000
				kicks[i].dur2=0
				kicks[i].dur3=0
				kicks[i].rs = 1e-3
			if current > 0:
				kicks[i].amp -= current
			stim[i].delay = 0
			stim[i].dur = 1e9
		i+=1
	
	while len(cells) > len(iarray):
		cells[-1] = None
		cells.pop()
		ncells-=1
			
	pc = h.ParallelContext()
	pc.nthread(64)
	
	h.finitialize()

	pc.psolve(STEPLEN)
	
	
	#free up some un-needed currents
	
	del kicks
	if not SYNAPTIC:
		del stim
	
	print( affix, affix2)
	fp = open('envelope_%s.dat' %(affix+affix2) , 'w')
	
	print(len(cells), ncells, len(vecs))
	
	for i in range(ncells):
		if STEPLEN > 2000:
			vecs[i].remove(0,10*(STEPLEN-2000)) # remove all but last second of data
		lmin = vecs[i].min()
		lmax = vecs[i].max()
		fp.write('%e  %e  %e\n' % (iarray[i],lmin,lmax))

	
	del vecs
	fp.close()


	fp = open('IV_%s.dat' %(affix+affix2) , 'w')
	
	
	pc.nthread(1)
	del pc
	del cells
	#cv = h.CVode()
	
	
	if SYNAPTIC:
		for pairs in gvpairs:
			max_eig, comp, ncomp, eigs = nf.stability_check(dopa,pairs[1],args=(0,pstable)) # this may be broken now
			print( pairs[0], pairs[1], max_eig, comp, ncomp)
			fp.write('%e  %e  %e  %e  %d\n' % (pairs[0],pairs[1],max_eig,comp, ncomp))
	else:
		iex.sort()
		for i in iarray:
			if 3>len(iex)>1:
				fstemp = []
				if i < iex[1]:
					fstemp.append(fs[0])
				if i > iex[0]:
					fstemp.append(fs[2])
				if iex[0] < i <iex[1]:
					fstemp.append(fs[1])
				
			else:
				fstemp = fs
			for f in fstemp:
				try:
					v = f(i)
					if v < vb[0] or v > vb[1]:
						continue
				except:
					continue
				max_eig, comp, ncomp, eigs = nf.stability_check(dopa,v ,args=(0,pstable),offset=i)
				print( v, i, max_eig, comp, ncomp)
				fp.write('%e  %e  %e  %e  %d\n' % (v,i,max_eig,comp,ncomp))
	
	#for i in range(len(varray)):
		#fp.write('%e  %e\n' %(varray[i],curarray[i]))
	fp.close()
	

if RUN and PRINT and GRAPH:
	if STEADY:
		os.system('xmgrace steady_%s_*.dat' % (affix))
	elif SYNAPTIC:
		os.system('xmgrace nmda_%s_*.dat' % (affix))
	else:
		os.system('xmgrace ramp_%s_*.dat' % (affix))
		
if GABARUN and GRAPH:
	os.system('xmgrace GABA_*%s*.dat' % (affix))

if KV4PULSE and GRAPH:
	os.system('xmgrace -nxy Kv4_*%s*.dat' % (affix))

if GABA_ARRAY and GRAPH:
	os.system('python heat_plot.py stats_*%s%s.dat %d' % (affix,affix2,axis))
	
if SYNPULSE and GRAPH:
	os.system('xmgrace pulsedgaba_%s*.dat' % (affix))

if DOCUMENT and False:
	try:
		os.system('mv *%s*.dat %s/.' % (affix, affix))
	except:
		pass

if not RUN or PRINT:

	quit()
else:
	h.load_file('depses.ses')