// calcisilag.hoc
// Olfactory bulb network model: calculate inter-spike interval
//                             : and lag time statistics
// Andrew Davison, The Babraham Institute, 2000.

/* The following procedures are defined for writing results to file
 * fileroot is the filename root. A suffix will be added to
 * this, e.g. fileroot.synch for print_si().
 * i,j are the mitral cell indices
 *
 * print_smooth_hist(variance,fileroot)
 * print_gran_smooth_hist(variance,fileroot)
 * print_spiketimes(i,j,fileroot)
 * print_raster(fileroot)
 * print_gran_raster(fileroot)
 * print_isis(i,j,fileroot)
 * print_isi_stats(fileroot)
 * print_lags(i,j,fileroot) 	
 * print_si(fileroot)
 */

// Variables used in this file
strdef tstr
objref work, work2, outputarray, isi, lags, hist, tmpspike
work 	= new Vector()
work2 	= new Vector()
outputarray = new Matrix(nmitx,nmity)
isi 	= new Vector()
lags 	= new Vector()
tmpspike = new Vector()

// Procedures for processing spike times -------------------------------

proc calc_isis() { local i,j,k,n // 3 args - indices of mitral cell, transient time
  if ($1 > nmitx || $2 > nmity) {
    print "Sorry - index out of range. Please try again."
    return
  }
  i = int($1)
  j = int($2)
  isi.resize(0)
  n = mit[i][j].spiketimes.size()
  if (n > 1) {
    for k = 1,n-1 {
      if (mit[i][j].spiketimes.x[k-1] > $3) {
        isi.append(mit[i][j].spiketimes.x[k]-mit[i][j].spiketimes.x[k-1])
      }
    }
  }
}

proc calc_gran_isis() { local i,j,k,n // 3 args - indices of granule cell, transient time
  if ($1 > ngranx || $2 > ngrany) {
    print "Sorry - index out of range. Please try again."
    return
  }
  i = int($1)
  j = int($2)
  isi.resize(0)
  n = gran[i][j].spiketimes.size()
  if (n > 1) {
    for k = 1,n-1 {
      if (gran[i][j].spiketimes.x[k-1] > $3) {
        isi.append(gran[i][j].spiketimes.x[k]-gran[i][j].spiketimes.x[k-1])
      }
    }
  }
}

func minisi() { local i,j,min	// find shortest mean ISI
  min = 1e6
  for i = 0, nmitx-1 {
    for j = 0, nmity-1 {
      calc_isis(i,j,ttrans)
      if (isi.size() > 0) {
        if (isi.mean() < min) { min = isi.mean() }
      }
    }
  }
  return min
}
      
  
proc rate_array() { local i,j
  for i = 0,nmitx-1 {
    for j = 0, nmity-1 {
      calc_isis(i,j,ttrans)
      if (isi.size() > 0) {
        outputarray.x[i][j] = 1000/isi.mean()
      } else {
        outputarray.x[i][j] = 0
      }
    }
  }  
  print "Max: ",arraymax(outputarray)
  outputarray.muls(1/arraymax(outputarray))
} 

proc calc_lags() { local i1,j1,i2,j2,k,minidx,min // 5 args - indices of mitral cells, transient time
    if ($1 > nmitx || $2 > nmity || $3 > nmitx || $4 > nmity) {
    print "Sorry - index out of range. Please try again."
    return
  }
  i1 = int($1)
  j1 = int($2)
  i2 = int($3)
  j2 = int($4)
  lags.resize(0)
  // for each spiketime in cell 1, find closest spike in cell 2
  // Note: first and last spikes ignored since can't calculate previous ISI
  if (mit[i2][j2].spiketimes.size > 0) {
    for k = 1,mit[i1][j1].spiketimes.size()-2 {
      if (mit[i1][j1].spiketimes.x[k] > $5) {
        work = mit[i2][j2].spiketimes.c.add(-mit[i1][j1].spiketimes.x[k])
        minidx = work.c.abs.min_ind()
        min = work.x[minidx]
        isiprev = mit[i1][j1].spiketimes.x[k-1]-mit[i1][j1].spiketimes.x[k]
        isinext = mit[i1][j1].spiketimes.x[k+1]-mit[i1][j1].spiketimes.x[k]
        if (min > isiprev/2 && min < isinext/2) {
          lags.append(min)
        }
      }
    }
  }
}

proc calc_phase_lags() { local i1,j1,i2,j2,k,minidx,min // 5 args - indices of mitral cells, transient time
    if ($1 > nmitx || $2 > nmity || $3 > nmitx || $4 > nmity) {
    print "Sorry - index out of range. Please try again."
    return
  }
  i1 = int($1)
  j1 = int($2)
  i2 = int($3)
  j2 = int($4)
  lags.resize(0)
  // for each spiketime in cell 1, find closest spike in cell 2
  // Note: first and last spikes ignored since can't calculate previous ISI
  if (mit[i2][j2].spiketimes.size > 0) {
    for k = 1,mit[i1][j1].spiketimes.size()-2 {
      if (mit[i1][j1].spiketimes.x[k] > $5) {
        work = mit[i2][j2].spiketimes.c.add(-mit[i1][j1].spiketimes.x[k])
        minidx = work.c.abs.min_ind()
        min = work.x[minidx]
        isiprev = mit[i1][j1].spiketimes.x[k-1]-mit[i1][j1].spiketimes.x[k]
        isinext = mit[i1][j1].spiketimes.x[k+1]-mit[i1][j1].spiketimes.x[k]
        if (min > isiprev/2 && min < isinext/2) {
          if (min < 0) {
            lags.append(min/isiprev)
          } else {
            lags.append(min/isinext)
          }
        }
      }
    }
  }
}

proc calc_gran_lags() { local i1,j1,i2,j2,k,minidx,min // 5 args - indices of granule cells, transient time
    if ($1 > ngranx || $2 > ngrany || $3 > ngranx || $4 > ngrany) {
    print "Sorry - index out of range. Please try again."
    return
  }
  i1 = int($1)
  j1 = int($2)
  i2 = int($3)
  j2 = int($4)
  lags.resize(0)
  // for each spiketime in cell 1, find closest spike in cell 2
  // Note: first and last spikes ignored since can't calculate previous ISI
  if (gran[i2][j2].spiketimes.size > 0) {
    for k = 1,gran[i1][j1].spiketimes.size()-2 {
      if (gran[i1][j1].spiketimes.x[k] > $5) {
        work = gran[i2][j2].spiketimes.c.add(-gran[i1][j1].spiketimes.x[k])
        minidx = work.c.abs.min_ind()
        min = work.x[minidx]
        isiprev = gran[i1][j1].spiketimes.x[k]-gran[i1][j1].spiketimes.x[k-1]
        isinext = gran[i1][j1].spiketimes.x[k]-gran[i1][j1].spiketimes.x[k+1]
        if (min < isiprev/2 && min > isinext/2) {
          lags.append(min)
        }
      }
    }
  }
}

proc time_hist() { // 1 arg - time step
  work.resize(0)
  for i = 0, nmitx-1 {
    for j = 0, nmity-1 {
      work.append(mit[i][j].spiketimes)
    }
  }
  hist = work.histogram(0,tstop,$1)
  hist.printf("%d\n")
}

func synch_index() { local i1,j1,i2,j2,n
  synchindex = 0
  n = 0
  for i1 = 0, nmitx-1 {
    for j1 = 0, nmity-1 {
      if (mit[i1][j1].spiketimes.size() > 0) {
        for i2 = 0, nmitx-1 {
          for j2 = 0, nmity-1 {
            if (i1 != i2 || j1 != j2) {
              calc_phase_lags(i1,j1,i2,j2,ttrans)
              n += lags.size()
              synchindex += lags.reduce("abs",0)
            }
          }
        }
      }
    }
  }
  if (n > 0) {
    return synchindex/n
  } else {
    return 1e6
  }
}


func phaselock_index() { local n,i1,j1,i2,j2
  synchindex = 0
  n = 0
  for i1 = 0, nmitx-1 {
    for j1 = 0, nmity-1 {
      if (mit[i1][j1].spiketimes.size() > 0) {
        for i2 = 0, nmitx-1 {
          for j2 = 0, nmity-1 {
            if (i1 != i2 || j1 != j2) {
              calc_phase_lags(i1,j1,i2,j2,ttrans)
              if (lags.size() > 1) {
                synchindex += lags.var()
                n += 1
              }
            }
          }
        }
      }
    }
  }
  if (n > 0) {
    synchindex = sqrt(synchindex/n)
    return synchindex
  } else {
    return 1e6
  }
}

// Procedures for writing out data --------------------------------

proc calc_smooth_hist() { // 1 arg - variance
  work.resize(0)
  for i = 0, nmitx-1 {
    for j = 0, nmity-1 {
      work.append(mit[i][j].spiketimes)
    }
  }
  hist = work.sumgauss(0,tstop,1,$1)
}

proc print_smooth_hist() { // 2 args - variance, filename root
  calc_smooth_hist($1)
  sprint(filename,"%s.smhist",$s2)
  outfile.wopen(filename)
  outfile.printf("# Mitral cell smoothed histogram\n")
  hist.printf(outfile,"%8.3f\n")
  outfile.close()
/*
  work.resize(0)
  hist.remove(0,ttrans)
  work.spctrm(hist)
  sprint(filename,"%s.pow",$s2)
  outfile.wopen(filename)
  outfile.printf("# Power spectrum of Mitral cell smoothed histogram\n")
  work.printf(outfile,"%9.5f\n")
  outfile.close()
  */
}

proc calc_gran_smooth_hist() { // 1 arg - variance
  work.resize(0)
  for i = 0, ngranx-1 {
    for j = 0, ngrany-1 {
      work.append(gran[i][j].spiketimes)
    }
  }
  hist = work.sumgauss(0,tstop,1,$1)
}

proc print_gran_smooth_hist() { // 2 args - variance, filename root
  calc_gran_smooth_hist($1)
  sprint(filename,"%s.gran.smhist",$s2)
  outfile.wopen(filename)
  outfile.printf("# Granule cell smoothed histogram\n")
  hist.printf(outfile,"%8.3f\n")
  outfile.close()
  //work.resize(0)
  //hist.remove(0,ttrans)
  //work.spctrm(hist)
  //sprint(filename,"%s.gran.pow",$s2)
  //outfile.wopen(filename)
  //outfile.printf("# Power spectrum of Granule cell smoothed histogram\n")
  //work.printf(outfile,"%9.5f\n")
  //outfile.close()
}

proc print_gran_hist() { // 2 args - binsize, filename root
  work.resize(0)
  pnm.pc.context("post_raster(ngranx, ngrany, grangid)\n")
  for i = 0, ngranx-1 {
    for j = 0, ngrany-1 {
	take_raster(i, j, grangid)
      work.append(tmpspike)
    }
  }
  hist = work.histogram(0,tstop,$1)
  sprint(filename,"%s.gran.hist",$s2)
  outfile.wopen(filename)
  outfile.printf("# Granule cell unsmoothed histogram\n")
  hist.printf(outfile,"%8.3f\n")
  outfile.close()
}

proc print_spiketimes() { // 3 args - indices of mitral cell plus filename root
  if (numarg() == 3) {
    sprint(filename,"%s_%d_%d.isi",$s3,$1,$2)
    outfile.wopen(filename)
    outfile.printf("# Spiketimes for mitral cell [%d][%d]",$1,$2)
    mit[$1][$2].spiketimes.printf(outfile,"%10.3f")
  }
  if (numarg() == 1) {
    sprint(filename,"%s.isi",$s1)
    outfile.wopen(filename)
    for i = 0, nmitx-1 {
      for j = 0, nmity-1 {
        outfile.printf("[%d][%d]",i,j)
        mit[i][j].spiketimes.printf(outfile,"%10.3f")
      }
    }
  }
  outfile.close()
}

proc post_raster() {local i, j
  for i = 0, $1-1 {
    for j = 0, $2-1 {
	if (pnm.pc.gid_exists($o3.x[i][j])) {
		sprint(tstr, "%d %d spikes", i, j)		
		if ($o3 == mitgid) {
			pnm.pc.post(tstr, mit[i][j].spiketimes)
		}else{
			pnm.pc.post(tstr, gran[i][j].spiketimes)
		}
	}
    }
  }
}

proc take_raster() {
	if (pnm.pc.gid_exists($o3.x[$1][$2])) {
		if ($o3 == mitgid) {
			tmpspike = mit[$1][$2].spiketimes.c
		}else{
			tmpspike = gran[$1][$2].spiketimes.c
		}
	}else{
		sprint(tstr, "%d %d spikes", $1, $2)
		pnm.pc.take(tstr, tmpspike)
	}
}

proc print_raster() { // 1 arg - filename root

  sprint(filename,"%s.ras",$s1)
  outfile.wopen(filename)
  outfile.printf("# Mitral cell raster plot\n")
  pnm.pc.context("post_raster(nmitx, nmity, mitgid)\n")
  for i = 0, nmitx-1 {
    for j = 0, nmity-1 {
	take_raster(i, j, mitgid)
      for k = 0, tmpspike.size()-1 {
        outfile.printf("%d %d %d %10.3f\n",i,j,i*nmity+j,tmpspike.x[k])
      }
    }
  }
  outfile.close()
}

proc print_gran_raster() { // 1 arg - filename root
  sprint(filename,"%s.gran.ras",$s1)
  outfile.wopen(filename)
  outfile.printf("# Granule cell raster plot\n")
  pnm.pc.context("post_raster(ngranx, ngrany, grangid)\n")
  for i = 0, ngranx-1 {
    for j = 0, ngrany-1 {
	take_raster(i, j, grangid)
      for k = 0, tmpspike.size()-1 {
        outfile.printf("%d %d %d %10.3f\n",i,j,i*ngrany+j,tmpspike.x[k])
      }
    }
  }
  outfile.close()
}

proc print_isis() { // 3 args - indices of mitral cell plus filename root
  calc_isis($1,$2,ttrans)
  sprint(filename,"%s_%d_%d.isi",$s3,$1,$2)
  outfile.wopen(filename)
  outfile.printf("# Interspike intervals for mitral cell [%d][%d]",$1,$2)
  isi.printf(outfile,"%10.3f")
  outfile.close()
}

proc print_isi_stats() { // 1 arg - filename root
  sprint(filename,"%s.stats",$s1)
  outfile.wopen(filename)
  outfile.printf("#Interspike interval statistics for mitral cells\n")
  outfile.printf("# i  j   n   mean   median   stdev  \n")
  for i = 0, nmitx-1 {
    for j = 0, nmity-1 {
      calc_isis(i,j,ttrans)
      outfile.printf("%3d%3d%4d",i,j,isi.size())
      if (isi.size() > 0) {
        outfile.printf("%8.2f%8.2f",isi.mean(),isi.median())
        if (isi.size() > 1) {
          outfile.printf("%8.2f\n",isi.stdev())
        } else {
          printf("\n")
        }
      } else { outfile.printf("\n") }
    }
  }
  outfile.printf("#Interspike interval statistics for granule cells\n")
  outfile.printf("# i  j   n   mean   median   stdev  \n")
  for i = 0, ngranx-1 {
    for j = 0, ngrany-1 {
      calc_gran_isis(i,j,ttrans)
      outfile.printf("%3d%3d%4d",i,j,isi.size())
      if (isi.size() > 0) {
        outfile.printf("%8.2f%8.2f",isi.mean(),isi.median())
        if (isi.size() > 1) {
          outfile.printf("%8.2f\n",isi.stdev())
        }
      } else { outfile.printf("\n") }
    }
  }
  outfile.close()
}

proc print_lags() { local i,j // 3 args - indices of mitral cell + filename root
  sprint(filename,"%s_%d_%d.lags",$s3,$1,$2)
  outfile.wopen(filename)
  outfile.printf("# Lag times for mitral cell [%d][%d]\n",$1,$2)
  for i = 0, nmitx-1 {
    for j = 0, nmity-1 {
      calc_lags($1,$2,i,j,ttrans)
      outfile.printf("[%d,%d]",i,j)
      lags.printf(outfile,"%10.3f")
    }
  }
  outfile.close()
}

proc print_si() { // 1 arg - fileroot
  print "Calculating synchronization indices"
  sprint(filename,"%s.synch",$s1)
  outfile.wopen(filename)
  outfile.printf("Synchronization index: %10.3f\n",synch_index())
  outfile.printf("Phase-locking index:   %10.3f\n",phaselock_index())
  outfile.close()
}