{load_file("nrngui.hoc")}  // load the GUI and standard run libraries

//////////////////////////////////
// Step 1: Define the cell classes
//////////////////////////////////

{load_file("cellid.hoc")}
{load_file("ranstream.hoc")}  // to give each cell its own sequence generator

//////////////////////////////////////////////////////////////
// Steps 2 and 3 are to create the cells and connect the cells
//////////////////////////////////////////////////////////////

NCELL = 20  // total number of cells in the ring network
C_E = 3  // # of excitatory connections received by each cell
         // 2 gives more sustained activity!
connect_random_low_start_ = 1  // low seed for mcell_ran4_init()

objref cells, nclist  // will be Lists that hold all network cell
    // and NetCon instances, respectively
objref ranlist  // for RandomStreams, one per cell

proc mknet() {
  mkcells($1)  // create the cells
  connectcells()  // connect them together
}

// creates the cells and appends them to a List called cells
// argument is the number of cells to be created
proc mkcells() {local i,j  localobj cell
  cells = new List()
  ranlist = new List()
  for i=0, $1-1 {
    cell = new B_BallStick()
    for j=0, cell.synlist.count-1 {
      cell.synlist.o(j).cid = i
      cell.synlist.o(j).sid = j
    }
    cells.append(cell)
    ranlist.append(new RandomStream(i))
  }
}

// connects the cells
// appends the NetCons to a List called nclist
// each target will receive exactly C_E unique non-self random connections
proc connectcells() {local i, nsyn, r  localobj src, syn, nc, rs, u
  // initialize the pseudorandom number generator
  mcell_ran4_init(connect_random_low_start_)
  u = new Vector(NCELL)  // for sampling without replacement
  nclist = new List()
  for i=0, cells.count-1 {
    // target synapse is synlist.object(0) on cells.object(i)
    syn = cells.object(i).synlist.object(0)
    rs = ranlist.object(i)  // the corresponding RandomStream
    rs.start()
    rs.r.discunif(0, NCELL-1)  // return integer in range 0..NCELL-1
    u.fill(0)  // u.x[i]==1 means spike source i has already been chosen
    nsyn = 0
    while (nsyn < C_E) {
      r = rs.repick()
      // no self-connection, & only one connection from any source
      if (r != i) if (u.x[r] == 0) {
        // set up connection from source to target
        src = cells.object(r)
        nc = src.connect2target(syn)
        nclist.append(nc)
        nc.delay = 1
        nc.weight = 0.01
        u.x[r] = 1
        nsyn += 1
      }
    }
  }
}

mknet(NCELL)  // go ahead and create the net!

//////////////////////////
// Report net architecture
//////////////////////////

proc tracenet() { local i  localobj src, tgt
  printf("source\ttarget\tsynapse\n")
  for i = 0, nclist.count-1 {
    src = nclist.o(i).precell
    tgt = nclist.o(i).syn
    printf("%d\t%d\t%d\n", src.synlist.o(0).cid, tgt.cid, tgt.sid)
  }
}

tracenet()

//////////////////////////////////////////////////
// Instrumentation, i.e. stimulation and recording
//////////////////////////////////////////////////

// stim will be an artificial spiking cell that generates a "spike" event
// that is delivered to the first cell in the net by ncstim
// in order to initiate network spiking.
// We won't bother including this "external stimulus source" or its NetCon
// in the network's lists of cells or NetCons.
objref stim, ncstim
proc mkstim() {
  stim = new NetStim()
  stim.number = 1
  stim.start = 0
  ncstim = new NetCon(stim, cells.object(0).synlist.object(0))
  ncstim.delay = 0
  ncstim.weight = 0.01
}

mkstim()

objref tvec, idvec  // will be Vectors that record all spike times (tvec)
        // and the corresponding id numbers of the cells that spiked (idvec)
proc spikerecord() {local i  localobj nc, nil
  tvec = new Vector()
  idvec = new Vector()
// the next line causes problems if there are more NetCons than cells
//  for i=0, nclist.count-1 {
  for i=0, cells.count-1 {
    nc = cells.object(i).connect2target(nil)
    nc.record(tvec, idvec, i)
    // the Vector will continue to record spike times
    // even after the NetCon has been destroyed
  }
}

spikerecord()

/////////////////////
// Simulation control
/////////////////////

tstop = 100
run()

////////////////////////////
// Report simulation results
////////////////////////////

proc spikeout() { local i
  printf("\ntime\t cell\n")
  for i=0, tvec.size-1 {
    printf("%g\t %d\n", tvec.x[i], idvec.x[i])
  }
}

spikeout()

quit()