// Template for NetLayer class.
// Intended to be used together with the LayerConn and ObjectArray classes.
// A NetLayer is an array of cells (either artificial cells such as the built-in
// IntFire4, or compartmental models encapsulated in a template) all of the same
// type. `Layer' is used as a generic term intended to include layers, columns,
// nuclei, etc., of cells.This class is intended to make management of large
// neuronal networks easier. The class provides convenient procedures for
// setting parameters of all cells in the layer at the same time. The LayerConn
// class can be used to easily specify synaptic connections between two
// NetLayers.

// SYNTAX:
// nl = new NetLayer(dim,n,"celltype")
// nl = new NetLayer(dim,n,"celltype",arg[,label])
// nl = new NetLayer(dim,n,"celltype",paramvec[,label])
//
// dim is the number of dimensions (1-3) of the cell array.
// n is the number of cells in each dimension. At the moment, all dimensions
//   will have the same size, i.e. only square and cubic arrays are possible.
//   Allowing n to be different for each dimension would add considerably to the
//   complexity of the code, and it didn't seem worthwhile. However, this could
//   be added if needed.
// "celltype" is the name of the artificial cell (e.g. IntFire4) or the name in
//   the begintemplate declaration (e.g. MyPyramidalCell).
// arg is a numerical or string argument, if the cells are created with a single argument
// paramvec is a vector containing the parameter values for the cell if there
// are two or more (otherwise use arg), e.g. if
//   celltype = "IntFire1" it is a 2-element vector containing the values of tau
//   and refrac.

// cell
//   Address an individual cell in the array, e.g. nl.cell[3][2][6], like x in
//   the Vector and Matrix classes.

// set("parametername",value)
//   Sets the value of parametername to value for every cell in the array. 
//   e.g. nl.set("tau",20.0)

// tset("parametername",&valueObj)
//   'Topographic' set. Sets the value of parametername to valueObj.x[...] 
//   for every cell i in the array. valueObj can be a Vector or a Matrix.

// rset("parametername",randomobj)
//   'Random' set. Sets the value of parametername to a value taken from
//   the randomobj Random object.

// call("procname","arguments")
//   Calls the procedure procname(arguments) for every cell in the array.
//   e.g. nl.call("set_background","0.1") if the cell template defines the
//   procedure set_background()

// tcall("procname","objarr")
//   `Topographic' call. Calls the procedure procname() for every cell in the 
//   array. The argument to the procedure depends on the coordinates of the
//   cell. objarr is a Vector, a Matrix or an ObjectArray of the same dimension
//   and size as the NetLayer.
//   e.g. nl.tcall("memb_init",matobj) 
//   calls nl.cell[i][j].memb_init(matobj.x[i][j]) for all i,j
//   and sets the initial membrane potential for each cell to the values stored
//   in the Matrix matobj.

// print_spikes(fileobj)
//   Assuming the cell object records spiketimes in a Vector spiketimes, this
//   procedure prints spike times to file in the two-column format
//   "spiketime cell_id" where cell_id is the index of the cell counting along
//   rows and down columns (and the extension of that for 3-D). This allows
//   easy plotting of a `raster' plot of spiketimes, with one line for each
//   cell.

// mean_spike_count()
//   Assuming the cell object records spiketimes in a Vector spiketimes, this
//   function returns the mean number of spikes per neuron

// dimensions()
//   Returns the number of dimensions of the NetLayer.

// size()
//   Returns the size of each dimension. If the NetLayer class is extended to
//   allow different sizes in different dimensions, this function could take an
//   argument specifying which dimension to return the size for.

// label
//   Allows a label for the LayerConn to be set or retreived. It doesn't do
//   anything at the moment, but could be used in a GUI.

// celltype
//   A string containing the type of cell used in the array. Changing this
//   string will not change the cell type, however.

// random_init(randobj)
//   Sets initial membrane potentials for all the
//   cells in the array to random values.

// Andrew Davison, UNIC, CNRS, August 2003-February 2006


begintemplate NetLayer
  public cell, tset, set, call, tcall, print_spikes, dimensions
  public size, plot_spikes, gui, label, celltype, rset
  public mean_spike_count, random_init, print_v
  
  create cell_holder  
  objref cell, this, cellParams
  strdef command, label, celltype, fmt, arg
  
  proc init() { local i,j,k,a4  // =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
    // $1:  Number of dimensions (max 3)
    // $2:  Size of each dimension (all dimensions same size)
    // $s3: Name of NetworkCell
    // $o4: Vector containing list of cell parameters
    
    create cell_holder  // for holding artificial cells/point processes
    access cell_holder
    m = $1  // should catch error if 1>m>3
    n = $2
    celltype = $s3
    if (numarg() > 3) {
      a4 = argtype(4)
      if (a4==0) { // number
	sprint(arg,"%g",$4)
      } else if (a4==2) { // string
	sprint(arg,"\"%s\"",$s4) // need to add quotes around
	//arg = $s4
      } else if (a4==1) { // object
	cellParams = $o4.c()  // copy here so the same vector can be used
                              // to initialize other layers, but the
                              // parameters can then be changed independently
        arg = "cellParams"
      } else { // error
	print "Error in NetLayer creation. Incorrect 4th argument"
	arg = ""
      }
      if (numarg() > 4) {
	label = $s5
      }
    } else {
      arg = ""           // if the `cell' is a NetStim or similar
    }
    sprint(command,"objref cell")
    for i = 0, m-1 {
      sprint(command,"%s[%d]",command,n)
    }
    execute1(command,this) // declare the objref array.
    
    // Create the cells.
    for i = 0,n-1 {
      if (m > 1) {
	for j = 0,n-1 {
	  if (m > 2) {
	    for k = 0,n-1 {
	      sprint(command,"cell[%d][%d][%d] = new %s(%s)",i,j,k,celltype,arg)
	      execute1(command,this)
	    }
	  } else {
	    sprint(command,"cell[%d][%d] = new %s(%s)",i,j,celltype,arg)
	    execute1(command,this)
	  }
	}
      } else {
	sprint(command,"cell[%d] = new %s(%s)",i,celltype,arg)
	execute1(command,this)
      }
    }
  }
  
  func dimensions() {  // =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
    // number of dimensions (1-3)
    return m
  }
  
  func size() {  // =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
    // size of each dimension
    return n
  }
  
  proc set() { local i,j,k // =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
    // for each cell in the array, set a parameter (1st argument) to a value
    // (2nd argument).
    for i = 0, n-1 {
      if (m > 1) {
	for j = 0,n-1 {
	  if (m > 2) {
	    for k = 0,n-1 {
	      sprint(command,"%s.cell[%d][%d][%d].%s = %g",this,i,j,k,$s1,$2)
	      execute1(command)
	    }
	  } else {
	    sprint(command,"%s.cell[%d][%d].%s = %g",this,i,j,$s1,$2)
	    execute1(command)
	  }
	}
      } else {
	sprint(command,"%s.cell[%d].%s = %g",this,i,$s1,$2)
	execute1(command)
      }
    }
  }
  
  proc tset() { local i,j,k // =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
    // for each cell in the array, set a parameter (1st argument) to a value
    // (2nd argument).
    if (m > 2) {
      print "Cannot be used for 3-D Layers"
      return
    }
    for i = 0, n-1 {
      if (m > 1) {
	for j = 0,n-1 {
	  sprint(command,"%s.cell[%d][%d].%s = %s.x[%d][%d]",this,i,j,$s1,$s2,i,j)
	  execute1(command)
	}
      } else {
	sprint(command,"%s.cell[%d].%s = %s.x[%d]",this,i,$s1,$s2,i)
	execute1(command)
      }
    }
  }
  
  proc rset() { local i,j,k // =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
    // for each cell in the array, set a parameter (1st argument) to a value
    // from the Random object given as the 2nd argument. An optional third
    // argument can specify a condition on the current value, e.g. "> 0"
    for i = 0, n-1 {
      if (m > 1) {
	for j = 0,n-1 {
	  if (m > 2) {
	    for k = 0,n-1 {
	      if (numarg() > 2) {
		sprint(command,"if (%s.cell[%d][%d][%d].%s %s) { %s.cell[%d][%d][%d].%s = %g }",this,i,j,k,$s1,$s3,this,i,j,k,$s1,$o2.repick())
	      } else {
		sprint(command,"%s.cell[%d][%d][%d].%s = %g",this,i,j,k,$s1,$o2.repick())
	      }
	      execute1(command)
	    }
	  } else {
	    if (numarg() > 2) {
	      sprint(command,"if (%s.cell[%d][%d].%s %s) { %s.cell[%d][%d].%s = %g }",this,i,j,$s1,$s3,this,i,j,$s1,$o2.repick())
	    } else {
	      sprint(command,"%s.cell[%d][%d].%s = %g",this,i,j,$s1,$o2.repick())
	    }
	    execute1(command)
	  }
	}
      } else {
	if (numarg() > 2) {
	  sprint(command,"if (%s.cell[%d].%s %s) { %s.cell[%d].%s = %g }",this,i,$s1,$s3,this,i,$s1,$o2.repick())
	} else {
	  sprint(command,"%s.cell[%d].%s = %g",this,i,$s1,$o2.repick())
	}
	execute1(command)
      }
    }
  }
  
  proc call() { local i,j,k  // =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
    // for each cell in the array, call a procedure (name is 1st argument) with
    // arguments given by the 2nd argument string.
    for i = 0, n-1 {
      if (m > 1) {
	for j = 0,n-1 {
	  if (m > 2) {
	    for k = 0,n-1 {
	      sprint(command,"err = %s.cell[%d][%d][%d].%s(%s)",this,i,j,k,$s1,$s2)
	      execute1(command)
	    }
	  } else {
	    sprint(command,"err = %s.cell[%d][%d].%s(%s)",this,i,j,$s1,$s2)
	    execute1(command)
	  }
	}
      } else {
	sprint(command,"err = %s.cell[%d].%s(%s)",this,i,$s1,$s2)
	execute1(command)
      }
    }
  }
  
  proc tcall() { local i,j,k // =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
    // t for topographic.
    // for each cell in the array, call a procedure (name is 1st argument) with
    // argument given by the value of the element of the 2nd argument (Vector,
    // Matrix or ObjectArray) at the corresponding location.
    for i = 0, n-1 {
      if (m > 1) {
	for j = 0,n-1 {
	  if (m > 2) {
	    for k = 0,n-1 {
	      sprint(command,"err = %s.cell[%d][%d][%d].%s(%s.x[%d][%d][%d])",this,i,j,k,$s1,$s2,i,j,k)
	      execute1(command)
	    }
	  } else {
	    sprint(command,"err = %s.cell[%d][%d].%s(%s.x[%d][%d])",this,i,j,$s1,$s2,i,j)
	    execute1(command)
	  }
	}
      } else {
	sprint(command,"err = %s.cell[%d].%s(%s.x[%d])",this,i,$s1,$s2,i)
	execute1(command)
      }
    }
  }
  
  proc print_spikes() { local i,j,k  // =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
    // Write spike times to file. Assumes that the cell template defines a 
    // Vector spiketimes that is used to record spikes.
    for i = 0, n-1 {
      if (m > 1) {
	for j = 0,n-1 {
	  if (m > 2) {
	    for k = 0,n-1 {
	      sprint(fmt,"%s\t%d\t%d\t%d\t%d\n","%.2f",i*n*n+j*n+k,i,j,k)
	      this.cell[i][j][k].spiketimes.printf($o1,fmt)
	    }
	  } else {
	    sprint(fmt,"%s\t%d\t%d\t%d\n","%.2f",i*n+j,i,j)
	    this.cell[i][j].spiketimes.printf($o1,fmt)
	  }
	}
      } else {
	sprint(fmt,"%s\t%d\n","%.2f",i)
	this.cell[i].spiketimes.printf($o1,fmt)
      }
    }
  }
  
  func mean_spike_count() { local i,j,k,nspikes
    // return the average number of spikes per cell. Assumes that the cell template defines a 
    // Vector spiketimes that is used to record spikes.
    nspikes = 0
    for i = 0, n-1 {
      if (m > 1) {
	for j = 0,n-1 {
	  if (m > 2) {
	    for k = 0,n-1 {
	      nspikes += this.cell[i][j][k].spiketimes.size()
	    }
	  } else {
	    nspikes += this.cell[i][j].spiketimes.size()
	  }
	}
      } else {
	nspikes += this.cell[i].spiketimes.size()
      }
    }
    return nspikes/(n^m)
  }
  
  proc print_v() { local i,j,k  // =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
    // Write membrane potential traces to file. Assumes that the cell
    // template defines a Vector vrecord that is used to record membrane
    // potential
    for i = 0, n-1 {
      if (m > 1) {
	for j = 0,n-1 {
	  if (m > 2) {
	    for k = 0,n-1 {
	      sprint(fmt,"%s\t%d\t%d\t%d\t%d\n","%.3g",i*n*n+j*n+k,i,j,k)
	      this.cell[i][j][k].vrecord.printf($o1,fmt)
	    }
	  } else {
	    sprint(fmt,"%s\t%d\t%d\t%d\n","%.3g",i*n+j,i,j)
	    this.cell[i][j].vrecord.printf($o1,fmt)
	  }
	}
      } else {
	sprint(fmt,"%s\t%d\n","%.3g",i)
	this.cell[i].vrecord.printf($o1,fmt)
      }
    }
  }
  
  proc random_init() { local r
    // randomly set cell membrane initial values 
    // $o1 = random object
    // $2  = lower bound of uniform distribution
    // $3  = upper bound of uniform distribution
    r = $o1.uniform($2,$3)
    for i = 0, n-1 {
      if (m > 1) {
	for j = 0,n-1 {
	  if (m > 2) {
	    for k = 0,n-1 {
	      r = $o1.repick()
	      sprint(command,"err = %s.cell[%d][%d][%d].memb_init(%f)",this,i,j,k,r)
	      execute1(command)
	    }
	  } else {
	    r = $o1.repick()
	    sprint(command,"err = %s.cell[%d][%d].memb_init(%f)",this,i,j,r)
	    execute1(command)
	  }
	}
      } else {
	r = $o1.repick()
	sprint(command,"err = %s.cell[%d].memb_init(%f)",this,i,r)
	execute1(command)
      }
    }
  }
  
  proc plot_spikes() { local i,j,k,l
    // $o1 = Graph
    // $s2 = mark symbol, e.g. "o", "|"
    // $3  = mark size
    // $4  = mark color
    // $5 =  mark brush
    for i = 0, n-1 {
      if (m > 1) {
	for j = 0,n-1 {
	  if (m > 2) {
	    for k = 0,n-1 {
	      for l = 0, this.cell[i][j][k].spiketimes.size()-1 {
		$o1.mark(this.cell[i][j][k].spiketimes.x[l],i*n*n+j*n+k,$s2,$3,$4,$5)
	      }
	    }
	  } else {
	    for l = 0, this.cell[i][j].spiketimes.size()-1 {
	      $o1.mark(this.cell[i][j].spiketimes.x[l],i*n+j,$s2,$3,$4,$5)
	    }
	  }
	}
      } else {
	for l = 0, this.cell[i].spiketimes.size()-1 {
	  $o1.mark(this.cell[i].spiketimes.x[l],i,$s2,$3,$4,$5)
	}
      }
    }
  }
  
  // TO DO
  //proc gui() {}    // provide a graphical interface to the NetLayer object.
  //proc save() {}   // for session files
  //proc load() {}
  
  // N.B. if a NetLayer is removed, are all the constituent cells
  // removed automatically, or do I need to clean up with unref() ?
  
endtemplate NetLayer