// $Id: shufmua.hoc,v 1.16 2012/01/11 19:39:35 samn Exp $ 

ncells = 300
nsec = 30
binsz = 5
refrac = 2.5
jitterdt = 0
declare("reg",0)
declare("noisespks", 0)

objref vs[ncells],vh[ncells],vmua,rdm,vy,vx
vmua=new Vector()
rdm=new Random()
vy=new Vector()
vx=new Vector()

sz =  (1e3/binsz) * nsec
rate  = 1 // rate in Hz

sampr = 1e3 / binsz

objref nqp[300]

SPECTY=0
PRESM=0

//* checkrefrac(vector,refrac)
proc checkrefrac () { local i localobj v1
  v1=$o1
  for i=1,v1.size-1 {
    if(v1.x(i)-v1.x(i-1)<$2) v1.x(i) = v1.x(i) + $2
  }
}

//* myshuf(vec,nshuffles,rdm)
proc myshuf () { local idx,i,j,k,n localobj v1,rdm,vm,vr,vr2
  v1=$o1 n=$2 rdm=$o3
  vm=new Vector(v1.size)
  vm.resize(0)
  for i=0,v1.size-1 if(v1.x(i)>0) vm.append(i)
  vr=new Vector(n)
  vr2=new Vector(n)
  rdm.discunif(0,vm.size-1)
  vr.setrand(rdm)
  rdm.discunif(0,v1.size-1)
  vr2.setrand(rdm)
  for i=0,n-1 {
    idx = vm.x(vr.x(i))
    j = vr2.x(i)
    k = v1.x(idx)
    v1.x(idx) = v1.x(j)
    v1.x(j) = k
  }
}

//* applyjitter(vec,rdm,dt)
proc applyjitter () { local i,jdt localobj v1,rdm,vj
  v1=$o1 rdm=$o2 jdt=$3
  vj=new Vector(v1.size)
  rdm.uniform(-jdt,jdt)
  vj.setrand(rdm)
  v1.add(vj)
  for i=0,v1.size-1 if(v1.x(i)<0) v1.x(i)=0// make sure no neg #s
  v1.sort()
}

//* addnoise(numspikes,rdm)
proc addnoise () { local i,ns localobj vec
  ns=$1 rdm.uniform(0,nsec*1e3)
  vec=new Vector(ns)
  for i=0,ncells-1 {
    vec.setrand(rdm)
    vs[i].append(vec)
    vs[i].sort()
  }
}

//* initcells
proc initcells () { local tt,isi,i,nshuf localobj vr
  nspks =  rate * nsec
  rdm.ACG(1234*nshuf)  
  for i=0,ncells-1 if(vs[i]==nil) vs[i]=new Vector() else vs[i].resize(0)
  if(reg) {
    isi = 1e3 / rate
    tt = isi
    for i=0,nspks-1 {
      vs[0].append(tt)
      tt += isi
    }
  } else {
    vs[0].resize(nspks)
    rdm.uniform(0,nsec*1e3)
    vs[0].setrand(rdm)
    vs[0].sort()
  }
  checkrefrac(vs[0],refrac)
  for i=1,ncells-1 {
    vs[i].copy(vs[0])
    applyjitter(vs[i],rdm,jitterdt)
    checkrefrac(vs[i],refrac)
  }
  if(noisespks) addnoise(noisespks,rdm)
}

//* histcells - make spike counts per time for each cell
proc histcells () { local i
  for i=0,ncells-1 {
    if(vh[i]==nil) vh[i]=new Vector() 
    vh[i].hist(vs[i],0,sz,binsz)
  }
}

//* shufhist(nshuf)
proc shufhist () { local i,nshuf
  nshuf=$1
  for i=1,ncells-1 {
    vh[i].copy(vh[0])
    myshuf(vh[i],nshuf,rdm)
    if(refrac) checkrefrac(vh[i],refrac)
  }
}

//* mkmua
proc mkmua () { local i
  vmua.resize(vh[0].size())
  vmua.fill(0)
  for i=0,ncells-1 vmua.add(vh[i])
  vmua.sub(vmua.mean())
}

//* plotrast
proc plotrast () { local i
  vrsz(0,vx,vy)
  for i=0,ncells-1 {
    for j=0,vs[i].size-1 {
      vx.append(vs[i].x(j))
      vy.append(i)
    }
  }
  vy.mark(g,vx,"O",2,1)
  g.exec_menu("View = plot")
}

//* setjitter(jitterdt)
proc setjitter () {
  jitterdt = $1
  initcells()
  histcells()
  mkmua()
}

//* jittertest(maxjitter,jitterinc)
proc jittertest () { local i
  for(jitterdt=0;jitterdt<=$1;jitterdt+=$2) {
    print "jitterdt is " , jitterdt
    setjitter(jitterdt)
    {nqsdel(nqp[i]) nqp[i]=getspecnq(vmua,sampr,SPECTY,PRESM)}
    // nqp.gr("pow","f",0,1,1)
    i += 1
  }
}

proc nqpg () {
  nqp[$1].gr("pow","f",0,1,1)
  g.exec_menu("View = plot")
}

//* main

setjitter(0)
gg()