"""
izh.py

Python/NEURON GUI for the different celltypes of Izhikevich neuron (versions from 2 publications). 

* 2003 Izhikevich artificial neuron model from 
EM Izhikevich "Simple Model of Spiking Neurons"
IEEE Transactions On Neural Networks, Vol. 14, No. 6, November 2003 pp 1569-1572

* 2007 Izhikevich artificial neuron model from 
EM Izhikevich (2007)  "Dynamical systems in neuroscience" MIT Press

Cell types available from Izhikevich, 2007 book:
    1. RS - Layer 5 regular spiking pyramidal cell (fig 8.12 from 2007 book)
    2. IB - Layer 5 intrinsically bursting cell (fig 8.19 from 2007 book)
    3. CH - Cat primary visual cortex chattering cell (fig8.23 from 2007 book)
    4. LTS - Rat barrel cortex Low-threshold  spiking interneuron (fig8.25 from 2007 book)
    5. FS - Rat visual cortex layer 5 fast-spiking interneuron (fig8.27 from 2007 book)
    6. TC - Cat dorsal LGN thalamocortical (TC) cell (fig8.31 from 2007 book)
    7. RTN - Rat reticular thalamic nucleus (RTN) cell  (fig8.32 from 2007 book)

Implementation by: Salvador Dura-Bernal, Cliff Kerr, Bill Lytton 
(salvadordura@gmail.com; cliffk@neurosim.downstate.edu; billl@neurosim.downstate.edu)
"""

# adapted from /u/billl/nrniv/sync/izh.hoc
import os, sys, collections
import numpy as np
from neuron import h, gui
h.load_file('stdrun.hoc')
import izhi2007Wrapper as izh07
import __main__
py = __main__
h.tstop=500
h.cvode_active(0)
h.dt=0.1
izh, cell07 = None, None # must be declared here since needs to be imported elsewhere

type2003 = collections.OrderedDict([
  #                                 a         b     c         d    vviv      tstop
 ('regular spiking (RS)'        , (0.02   ,  0.2 , -65.0 ,   8.0 , -63.0 ,   150.0)) ,
 ('intrinsically bursting (IB)' , (0.02   ,  0.2 , -55.0 ,   4.0 , -70.0 ,   150.0)) ,
 ('chattering (CH)'             , (0.02   ,  0.2 , -50.0 ,   2.0 , -70.0 ,   150.0)) ,
 ('fast spiking (FS)'           , (0.1    ,  0.2 , -65.0 ,   2.0 , -70.0 ,   150.0)) ,
 ('thalamo-cortical (TC)'       , (0.02   , 0.25,  -65.0 , 0.05 , -63.0 ,   150.0)) ,
 ('thalamo-cortical burst (TC)' , (0.02   , 0.25,  -65.0 , 0.05 , -87.0 ,   150.0)) ,
 ('resonator (RZ)'              , (0.1   ,  0.26 , -65.0 ,   2.0 , -70.0 ,   100.0)) ,
 ('low-threshold spiking (LTS)' , (0.02   , 0.25 , -65.0 ,   2.0 , -63.0 ,   250.0))])

type2004 = collections.OrderedDict([
  #                                 a         b     c         d    vviv      tstop
 ('tonic spiking'               , (0.02   ,  0.2 , -65.0 ,   6.0 , -70.0 ,   100.0)) ,
 ('mixed mode'                  , (0.02   ,  0.2 , -55.0 ,   4.0 , -70.0 ,   160.0)) ,
 ('spike latency'               , (0.02   ,  0.2 , -65.0 ,   6.0 , -70.0 ,   100.0)) ,
 ('rebound spike'               , (0.03   , 0.25 , -60.0 ,   4.0 , -64.0 ,   200.0)) ,
 ('Depolarizing afterpotential' , (1.0    ,  0.2 , -60.0 , -21.0 , -70.0 ,    50.0)) ,
 ('phasic spiking'              , (0.02   , 0.25 , -65.0 ,   6.0 , -64.0 ,   200.0)) ,
 ('spike frequency adaptation'  , (0.01   ,  0.2 , -65.0 ,   8.0 , -70.0 ,    85.0)) ,
 ('subthreshold oscillations'   , (0.05   , 0.26 , -60.0 ,   0.0 , -62.0 ,   200.0)) ,
 ('rebound burst'               , (0.03   , 0.25 , -52.0 ,   0.0 , -64.0 ,   200.0)) ,
 ('accomodation'                , (0.02   ,  1.0 , -55.0 ,   4.0 , -65.0 ,   400.0)) ,
 ('tonic bursting'              , (0.02   ,  0.2 , -50.0 ,   2.0 , -70.0 ,   220.0)) ,
 ('Class 1'                     , (0.02   , -0.1 , -55.0 ,   6.0 , -60.0 ,   300.0)) ,
 ('resonator'                   , (0.1    , 0.26 , -60.0 ,  -1.0 , -62.0 ,   400.0)) ,
 ('threshold variability'       , (0.03   , 0.25 , -60.0 ,   4.0 , -64.0 ,   100.0)) ,
 ('inhibition-induced spiking'  , (-0.02  , -1.0 , -60.0 ,   8.0 , -63.8 ,   350.0)) ,
 ('phasic bursting'             , (0.02   , 0.25 , -55.0 ,  0.05 , -64.0 ,   200.0)) ,
 ('Class 2'                     , (0.2    , 0.26 , -65.0 ,   0.0 , -64.0 ,   300.0)) ,
 ('integrator'                  , (0.02   , -0.1 , -55.0 ,   6.0 , -60.0 ,   100.0)) ,
 ('bistability'                 , (0.1    , 0.26 , -60.0 ,   0.0 , -61.0 ,   300.0)) ,
 ('inhibition-induced bursting' , (-0.026 , -1.0 , -45.0 ,  -2.0 , -63.8 , 350.0))])

choices = collections.OrderedDict([
  ('2003 PP model' ,  (lambda: h.Izhi2003a(0.5,sec=cell03),  lambda: izh._ref_V,      type2003)),
  ('2003 Sec model',  (lambda: h.Izhi2003b(0.5,sec=cell03), lambda: cell03(0.5)._ref_v, type2003)),
  ('2004 PP model' ,  (lambda: h.Izhi2003a(0.5,sec=cell03),  lambda: izh._ref_V,      type2004)),
  ('2004 Sec model',  (lambda: h.Izhi2003b(0.5,sec=cell03), lambda: cell03(0.5)._ref_v, type2004)),
  ('2007 PP model' ,  (lambda: izh07.IzhiCell(host=izh07.dummy), lambda: izh._ref_V, izh07.type2007)),
  ('2007 Sec model' , (lambda: izh07.IzhiCell(), lambda: cell07.sec(0.5)._ref_v, izh07.type2007))])

ch=list(choices.keys())
def newmodel (ty=None) : 
  "2003,2004 was the orig model; 2007 is the redesign; look at global izhtype if no "
  return izhtype.find('2007') > -1 if ty is None else ty.find('2007') > -1

#* setup the cell
izhtype='2004 PP model'
def cellset ():
  global cell07, cell03, izh, vref, uvvset, fih, izhtype
  if newmodel():
    cell07 =  choices[izhtype][0]()
    izh = cell07.izh
    def uvvset () : pass
  else: 
    cell03 = h.Section(name="cell2003") # this cell will be used for 2003/4; different cell created in izhi2007Wrapper for those
    izh =  choices[izhtype][0]()
    def uvvset () : vref[0], izh.u = vviv, vviv*izh.b
    cell03.L, cell03.diam = 6.37, 5 # empirically tuned -- cell size only used for Izh1
  fih = [h.FInitializeHandler(uvvset), h.FInitializeHandler(0,Isend)]
  vref = choices[izhtype][1]() # can define this afterwards even though used in uvvset above
  # h('objref izh'); h.izh = izh # if need to access from hoc

#* parameters for different cell types
playvec, playtvec = [h.Vector() for x in range(2)]

# initialization routines
name, params = None, None

def p (nm, pm=None) :
  global name, vviv, params, vvset
  import izhi2007Figs as iz07fig
  if pm is None : pm = choices[izhtype][2][nm]
  name, params = nm, pm
  if newmodel():
    izh.C, izh.k, izh.vr, izh.vt, izh.vpeak, izh.a, izh.b, izh.c, izh.d, izh.celltype = params
    h.tstop=1000
  else: 
    izh.a, izh.b, izh.c, izh.d, vviv, h.tstop = params
  g.size(0,h.tstop,-100,50)
  try: 
    if newmodel():
      graphx() # interviews graphics
      iz07fig.recorder(cell07, choices[izhtype][1]()) # vectors to draw under matplotlib
      iz07fig.test1(cell07, nm, izhtype)
    else:
      iz07fig.closeFig()
      graphx()
      playinit()
      h.run()
  except: print((sys.exc_info()[0],' :',sys.exc_info()[1]))

def ivwrap (func, label=''):
  wrapper = h.VBox()
  wrapper.intercept(1)
  func()
  wrapper.intercept(0)
  wrapper.map(label)
  return wrapper

def graphx ():
  g.erase_all()
  g.addvar("v", choices[izhtype][1](), 2,2)
  g.addvar("u", izh._ref_u, 3,1)
  g.addvar("Iin", izh._ref_Iin if newmodel() else izh._ref_Iin, 4,2)
  try: g.addvar("gsyn", izh._ref_gsyn, 1, 1)
  except: pass

I0=I1=T1=0
def playinit () : 
  global I0,I1,T1
  try: izh.f, izh.g= 5, 140 # standard params: V'=0.04*V^2 + 5*V + 140 - u + Iin
  except: pass
  bub.label[0] = '%s'%(name)
  if name=='Depolarizing afterpotential': bub.label[0] = "%s -- REPEATED SPIKING"%(bub.label[0])
  if name=='accomodation': bub.label[0] = "%s -- NOT IMPLEMENTED (different functional form;see izh.mod)"%(bub.label[0])
  if name=='inhibition-induced bursting': bub.label[0] = "%s -- NOT IMPLEMENTED (convergence problems)"%(bub.label[0])
  g.label(0.1,0.9,bub.label[0])
  print(bub.label[0])
  playvec.play_remove()
  playtvec.resize(0); playvec.resize(0)
  if name=='Class 1' :
    T1=30
    playtvec.append(0,T1,h.tstop)
    playvec.append(0,0,0.075*(h.tstop-T1))
  elif name=='Class 2' : #  (H) Class 2 exc.
    T1=30
    playtvec.append(0,T1,h.tstop)
    playvec.append(-0.5, -0.5,-0.05+0.015*(h.tstop-T1))
  elif name=='accomodation' : #  (R) accomodation
    playtvec.append(0, 200,    200.001, 300,     312.5, 312.501, h.tstop)
    playvec.append( 0, 200/25, 0      , 0      , 4    , 0      , 0)
  if name in ['Class 1', 'Class 2', 'accomodation'] : playvec.play(izh._ref_Iin, playtvec, 1)
  if name in ['Class 1', 'integrator'] :
    try: izh.f, izh.g = 4.1, 108 # don't exist in all the models
    except: pass

def synon () : 
  "Turn on a synapse"
  global ns, nc
  ns = h.NetStim()
  nc = h.NetCon(ns,izh,0,1,10)
  ns.start, ns.interval, ns.number = 10, 10, 10
  nc.weight[0] = 2
  izh.taug = 3

#* box of buttons
class Bubox :

  def __init__ (self, type, li) :
    self.izhtype = type
    vbox, hbox, hbox1 = h.VBox(), h.HBox(), h.HBox()
    self.vbox = vbox
    lil = len(li)
    self.cols, self.rows = {20:(4,5), 8:(4,2), 9:(3,3)}[lil]
    self.label=h.ref('================================================================================')
    vbox.intercept(1)
    h.xpanel("")
    h.xvarlabel(self.label)
    if newmodel(self.izhtype):
      h.xlabel("V' = (k*(V-vr)*(V-vt) - u + Iin)/C     if (V>vpeak) V=c  [reset]")
      h.xlabel("u' = a*(b*(V-vr) - u)                  if (V>vpeak) u=u+d")
    else: 
      h.xlabel("v' = 0.04*v*v + f*v + g - u + Iin;     if (v>thresh) v=c [reset]")
      h.xlabel("u' = a*(b*v - u);                    if (v>thresh) u=u+d")
    h.xpanel()
    hbox1.intercept(1)
    h.xpanel(""); h.xbutton("RUN",h.run); h.xpanel()
    self.xvalue('I0','I0')
    self.xvalue('I1','I1')
    self.xvalue('T1','T1')
    hbox1.intercept(0); hbox1.map("")
    hbox.intercept(1)
    for ii,(k,v) in enumerate(li.items()):
      if ii%self.rows==0: h.xpanel("")
      h.xbutton(k, (lambda f, arg1, arg2: lambda: f(arg1,arg2))(p, k, v)) # alternative is to use functools.partial
      if ii%self.rows==self.rows-1: h.xpanel()
    hbox.intercept(0); hbox.map("")
    vbox.intercept(0); vbox.map("Spike patterns")
    self.label[0]=""

  def pr (): pass

  def xvalue (self,name,var,obj=py,runner=pr):
    h.xpanel("")
    h.xvalue(name,(obj, var),0,runner)
    h.xpanel()

  def xpvalue (self,name,ptr,runner=pr):
    "Doesn't work currently"
    h.xpanel("")
    h.xpvalue(name,ptr,1,runner)
    h.xpanel()

  def transpose (self,x) : return int(x/self.rows) + x%self.rows*self.cols
# end class Bubox

# current injections for specific models
def Isend () :
  global T1,I0,I1
  if I0!=0 or I1!=0:
    Iin = I0
    Isend1(T1,I1)
    return
  T1=h.tstop/10
  if not newmodel(): izh.Iin=0
  if name=='tonic spiking':   #  (A) tonic spiking
    Isend1(T1,14)
  elif name=='phasic spiking':  #  (B) phasic spiking
    T1=20
    Isend1(T1,0.5)
  elif name=='tonic bursting':  #  (C) tonic bursting
    T1=22
    Isend1(T1,15)
  elif name=='phasic bursting':  #  (D) phasic bursting
    T1=20
    Isend1(T1,0.6)
  elif name=='mixed mode':  #  (E) mixed mode
    Isend1(T1,10)
  elif name=='spike frequency adaptation':  #  (F) spike freq. adapt
    Isend1(T1,30)
  elif name=='Class 1':  #  (G) Class 1 exc. -- playvec
    pass
  elif name=='Class 2':  #  (H) Class 2 exc. -- playvec
    pass
  elif name=='spike latency':  #  (izh.Iin) spike latency
    Isend1(T1,7.04)
    Isend1(T1+3,0.0)
  elif name=='subthreshold oscillations':  #  (J) subthresh. osc.
    Isend1(T1,2)
    Isend1(T1+5,0)
  elif name=='resonator':  #  (K) resonator
    T2, T3 = T1+20, 0.7*h.tstop
    T4 = T3+40
    Isend1(T1,0.65) ;   Isend1(T2,0.65)  ;  Isend1(T3,0.65)  ;  Isend1(T4,0.65)
    Isend1(T1+4,0.) ;   Isend1(T2+4,0.)  ;  Isend1(T3+4,0.)  ;  Isend1(T4+4,0.)
  elif name=='integrator':  #  (L) integrator
    T1, T3 = h.tstop/11, 0.7*h.tstop
    T2, T4 = T1+5, T3+10
    Isend1(T1,9)    ; Isend1(T2,9)    ; Isend1(T3,9)    ; Isend1(T4,9)
    Isend1(T1+2,0.) ; Isend1(T2+2,0.) ; Isend1(T3+2,0.) ; Isend1(T4+4,0.)
  elif name=='rebound spike':  #  (M) rebound spike
    T1=20
    Isend1(T1,-15)
    Isend1(T1+5,0)
  elif name=='rebound burst':  #  (N) rebound burst
    T1=20
    Isend1(T1,-15)
    Isend1(T1+5,0)
  elif name=='threshold variability':  #  (O) thresh. variability
    T1, T2, T3 =10, 70, 80
    Isend1(T1,1)    ; Isend1(T2,-6)   ; Isend1(T3,1)
    Isend1(T1+5,0.) ; Isend1(T2+5,0.) ; Isend1(T3+5,0.)
  elif name=='bistability':  #  (P) bistability
    T1, T2, izh.Iin = h.tstop/8, 216, 0.24
    Isend1(T1,1.24) ;   Isend1(T2,1.24)
    Isend1(T1+5,0.24);  Isend1(T2+5,0.24)
  elif name=='Depolarizing afterpotential':  #  (Q) DAP depolarizing afterpotential
    T1 = 10
    Isend1(T1-1,20)
    Isend1(T1+1,0)
  elif name=='accomodation':  #  (R) accomodation -- playvec
    pass
  elif name=='inhibition-induced spiking':  #  (S) inhibition induced spiking
    izh.Iin=80
    Isend1(50,75)
    Isend1(250,80)
  elif name=='inhibition-induced bursting':  #  (T) inhibition induced bursting
    izh.Iin=80
    Isend1(50,80) # Isend1(50,75) -- will crash simulator
    Isend1(250,80)
  elif name=='regular spiking (RS)':  #  regular spiking (RS)
    Isend1(T1,14)
  elif name=='intrinsically bursting (IB)':  #  intrinsically bursting (IB)
    Isend1(T1,11)
  elif name=='chattering (CH)':  #  chattering (CH)
    Isend1(T1,10)
  elif name=='fast spiking (FS)':  #  fast spiking (FS)
    Isend1(T1,10)
  elif name=='thalamo-cortical (TC)':  #  thalamo-cortical (TC)
    Isend1(2*T1,1.5)
  elif name=='thalamo-cortical burst (TC)':  #  thalamo-cortical burst (TC)
    Isend1(0,-25)
    Isend1(3*T1,0)
  elif name=='resonator (RZ)':  #  resonator (RZ)
    Isend1(0,-2)
    Isend1(T1,-0.5)
    Isend1(T1+50,10)
    Isend1(T1+55,-0.5)
  elif name=='low-threshold spiking (LTS)':  #  low-threshold spiking (LTS)
    Isend1(T1,10)    
  elif name == 'TC_burst': # thalamo-cortical burst (TC) (2007)
    Isend1(0,-1200) 
    Isend1(120,110)
  elif name == 'RTN_burst': # reticular thalamic nucleus burst (TC) (2007)
    Isend1(0,-350) 
    Isend1(120,90)

def Isend1 (tm, Iin) :
  def my_event():
    izh.Iin = Iin
    h.CVode().re_init()
  h.cvode.event(tm, my_event)

# izhstim() sets up a single stim into izh cell
# effect easily seen by running "Class 1"
def izhstim () :
  stim=h.NetStim(0.5)
  stim.number = stim.start = 1
  nc = h.NetCon(stim,izh)
  nc.delay = 2
  nc.weight = 0.1
  izh.erev = -5

#* plotting & printing
g, nmenu, bub = None, None, None
def isinstanceh (objref,objtype) : return objref.hname().startswith(objtype.hname()[:-2])

def winup (izht=izhtype):
  global bub, g, nmenu, izhtype
  izhtype = izht  # swap in the new one
  cellset()
  if g is None: 
    g=h.Graph(0)
    h.graphList[0].append(g)
  if g.view_count()<1: 
    g.view(-0.1*h.tstop,-90,1.2*h.tstop,150,300,200,400,200)
    g.size(0,h.tstop,-80,40)
  if not bub is None: bub.vbox.unmap()
  bub = Bubox(izhtype,choices[izhtype][2])
  bub.label[0] = izhtype
  if not nmenu is None: nmenu.unmap()
  nmenu = ivwrap(lambda: h.nrnpointmenu(izh), izh.hname())

def chwin ():
  "Launch windows from model list"
  h.xpanel("Izhikevich models")
  # outer lambda returns inner lambda so as to pass arg to winup() -- the innermost routine
  for c in ch:
    h.xbutton(c, (lambda f, arg1: lambda: f(arg1))(winup,c)) 
  h.xpanel()

def vtvec(vv): return np.linspace(0, len(vv)*h.dt, len(vv), endpoint=True)

if __name__ == '__main__': chwin()
