main.hoc

{load_file("nrngui.hoc")}  // loads the gui and standard run libraries
{load_file("netparmpi.hoc")} // parallelnetmanager definitions

objref nil
strdef modelName, datasetPrefix, datasetName, datasetPath, resultsPath

// parameters must be the name of a file specified on the command line
// as follows: nrniv -c "strdef parameters" -c "parameters=\"parameters.hoc\""
{load_file(parameters)}

{cvode.active(use_cvode)}
{cvode.use_local_dt(use_local_dt)}
{cvode.cache_efficient(use_cache_efficient)}

{xopen("lib.hoc")}

// path to the dataset directory
sprint(datasetPath, "%s/%s", datasetPrefix, datasetName)

//nrn_load_dll("mechanisms/x86_64/.libs/libnrnmech.so")

// used to calculate model construction times and run time
mkcellstime = 0
connectcellstime = 0
connectgjstime = 0

// used for executing command strings
strdef cmd

tstop = mytstop

// path for various datafiles (used for multiple files)
strdef datafilePath

// randomstream template
{load_file("./templates/ranstream.hoc")}
// stimulus cell template
{load_file("./templates/StimCell.hoc")}

// used for recording membrane potentials
objref indicesVrecord
if (!(name_declared("vrecordFraction"))) {
    vrecordFraction = 0
}
if (vrecordFraction > 1) {
    vrecordFraction = 1
}
if (vrecordFraction < 0) {
    vrecordFraction = 0
}
if (!(name_declared("vrecordLimit"))) {
    vrecordLimit = 100
}
if (vrecordLimit < 0) {
    vrecordLimit = 0
}

// Builds a list of indices of cells to record from
proc mkvrecord() { local nCellTypes, numCells, numRec, i, j, k, offset, n, p, fraction  localobj pnm, cellTypes, cellType, indices, rec_indices, rs

    pnm = $o1    
    rec_indices = $o2
    cellTypes = $o3
    fraction = $4

    nCellTypes = cellTypes.count()
    
    for k=0, nCellTypes-1 { 
        
        rs = new RandomStream(k, 0, 0)
        
        cellType = cellTypes.o(k)
        numCells = cellType.getPropertyScalar("numCells")
        numCellsPerHost = int(numCells / pnm.pc.nhost)
        offset   = cellType.getPropertyScalar("offset")

        if ((numCellsPerHost < 1) ) {
            if (pnm.myid < numCells) { 
              numRec = 1
            } else {
              numRec = 0
            }
        } else {
            numRec = int(fraction * numCellsPerHost)
            if (numRec < 1) {
               numRec = 1
            }
        }
        if ((vrecordLimit > 0) && (numRec > vrecordLimit)) {
            numRec = vrecordLimit
        }

        if (verbose > 0) {
            printf ("cell type %s: numCells = %d numRec = %d\n", cellType.getPropertyString("typeName").s, numCells, numRec)
        }
        
        if (cellType.hasProperty("indices")) {
            indices = cellType.getPropertyObject("indices")
        } else {
            indices = new Vector()
            indices.resize(numCells)
            for i = 0, numCells-1 {
                indices.x[i] = i+offset
            }
        }

        if (numRec > 1) {
          rs.r.discunif(0,numCellsPerHost-1)
          for i=0, numRec-1 { 
              n = rs.repick()
              j = pnm.myid + n * pnm.pc.nhost
              rec_indices.prepend(new Value(0,indices.x[j]))
          }
        } else if (numRec == 1) {
          rec_indices.prepend(new Value(0,indices.x[pnm.myid]))
       }
    }
    
}



// used to keep track of artificial stimulus cells
objref stimTypes
stimTypes = new List()
objref ncstimlist

// reads the cell types file and populates celltypes with the information
sprint (datafilePath, "%s/stim.dat", datasetPath)
sprint(cmd, "test -e %s", datafilePath) // check if the path exist
// status = system(cmd)
status = 0
if (status == 0) {
   loadStimulationInfo (stimTypes, datafilePath, datasetPath)
   if (verbose > 0) {
      printf ("%s: loaded stimulation info\n", modelName)
   }
}



// list with cell type information
objref cellTypes
cellTypes = new List()

// reads the cell types file and populates celltypes with the information
sprint (datafilePath, "%s/celltypes.dat", datasetPath)
loadCellCategoryInfo (cellTypes, datafilePath, datasetPath)

if (verbose > 0) {
    printf ("%s: loaded cell info\n", modelName)
}

// used for synapse information
objref synapseTypes
synapseTypes = new List()
// reads the synapse types file and populates synapseTypes with the information
loadSynapseInfo (synapseTypes, datasetPath)

if (verbose > 0) {
    printf ("%s: loaded synapse info\n", modelName)
}

// total number of cells in the network
ncells = 0
maxindex = 0
strdef templatename
objref celltype, indices
for i=0, cellTypes.count()-1 { 
    celltype = cellTypes.o(i)
    if (celltype.hasProperty("indices")) {
       indices = celltype.getPropertyObject("indices")
       maxindex = indices.max()
       if (ncells < maxindex) {
          ncells = maxindex+1
       }
    } else {
       ncells = ncells + celltype.getPropertyScalar("numCells") + 1
    }
}
if (verbose > 0) {
    printf ("%s: ncells = %d\n", modelName, ncells)
}


// list with connectivity information
objref connectivityTypes
connectivityTypes = new List()

// reads the connectivity types file and populates connectivitytypes with the information
sprint (datafilePath, "%s/connectivity.dat", datasetPath)
loadConnectivityInfo (connectivityTypes, datafilePath)
if (verbose > 0) {
    printf ("connectivityTypes.count = %d\n", connectivityTypes.count())
}

// list with gap junction information
objref gjTypes
gjTypes = new List()

// reads the connectivity types file and populates connectivitytypes with the information
sprint (datafilePath, "%s/gapjunctions.dat", datasetPath)
loadGapJunctionInfo (gjTypes, datafilePath)

// pnm.cells is a list that holds all instances of network cells that
// exist on this host

// pnm.nclist holds all netcon instances that exist on this host

objref pnm

pnm = new ParallelNetManager(ncells)

ncells = 0
for i=0, cellTypes.count()-1 { 
    celltype = cellTypes.o(i)
    if (celltype.hasProperty("indices")) {
       indices = celltype.getPropertyObject("indices")
       for j = 0, indices.size()-1 {
          pnm.set_gid2node(indices.x[j], j % pnm.pc.nhost)
       }
       maxindex = indices.max()
       if (ncells < maxindex) {
          ncells = maxindex+1
       }
    } else {
       for j = ncells, ncells+celltype.getPropertyScalar("numCells")-1 {
          pnm.set_gid2node(j, j % pnm.pc.nhost)
       }
       ncells = ncells + celltype.getPropertyScalar("numCells") + 1
    }
}



// connectivity list local to the current mpi process
objref conlist
if (makegraph > 0) {
    conlist = new List()
}

// gap junction connectivity list local to the current mpi process
objref gjconlist
if (makegraph > 0) {
    gjconlist = new List()
}

// gap junction list local to the current mpi process
objref gjlist
gjlist = new List()

/* mkcells and connectcells create the cells and connect the cells */

objref cell

// creates the cells and appends them to a list called cells
// argument is the number of cells to be created
proc mkcells() { local nCellTypes, numCells, i, j, p, gid, hasSynapseProps, templateType  localobj f, cellTypes, cellType, datadirIndex, indices, synProps, synRiseTime, synDecayTime, synErev, synDend, synComp, rsdend, rscomp

    
    strdef templateName, typeName, datadir, templatePath
    
    cellTypes = $o1
    nCellTypes = cellTypes.count()
    
    for k=0, nCellTypes-1 { 
        
        
        cellType = cellTypes.o(k)
        
        typeName = cellType.getPropertyString("typeName").s
        templateName = cellType.getPropertyString("templateName").s
        numCells = cellType.getPropertyScalar("numCells")
        offset   = cellType.getPropertyScalar("offset")
        templateType   = cellType.getPropertyScalar("templateType")

        if (cellType.hasProperty("indices")) {
            indices = cellType.getPropertyObject("indices")
        } else {
            indices = new Vector()
            indices.resize(numCells)
            for i = 0, numCells-1 {
                indices.x[i] = i+offset
            }
            cellType.addPropertyObject("indices", indices)
        }
        

        hasSynapseProps = hasSynapseProperties(synapseTypes,typeName)
        if (hasSynapseProps) {
           synProps = getSynapseProperties(synapseTypes,typeName)
           synRiseTime = synProps.getPropertyObject("rise time")
           synDecayTime = synProps.getPropertyObject("decay time")
           synErev = synProps.getPropertyObject("reversal potential")
           synDend = synProps.getPropertyObject("dendrite")
           synComp = synProps.getPropertyObject("compartment")
       }
       
       
       sprint(templatePath, "./templates/%s.hoc", templateName)
       sprint(cmd, "load_file(\"%s\")", templatePath)
       execute1(cmd)

       for i=0, numCells-1 { 
           gid = indices.x(i)

           if (pnm.gid_exists(gid)) {
               if (verbose > 0) {
                   printf ("host %d: mkcells: i = %d gid = %d\n", pnm.myid, i, gid)
               }
               // used for choosing dendrite and compartment for synapses
               rsdend = new RandomStream(gid, k, 1)
               rscomp = new RandomStream(gid, k, 2)
               sprint(cmd, "cell = new %s(%d, %d, %d, \"%s/%s/\")", templateName, i, gid, gid, datasetPath, typeName)
               execute(cmd)
               if (hasSynapseProps) {
                   mksyn(cell,synRiseTime,synDecayTime,synErev,synDend,synComp,rsdend,rscomp,typeName)
                   if (i == 0) { printCellSections (typeName,cell) }
               }
	       pnm.register_cell(gid, cell)
               pnm.spike_record(gid)
            }
        }
    }
    
}
    
// connects the cells
// appends the NetCon objects to nclist
proc connectcells() {local i, k, status, srcid, destid, srcid1, destid1, del, dist, velocity, wdType, order, presynapticSize, postsynapticSize, stimsrc localobj presynapticType, postsynapticType, cellTypes, connectivityType, connectivityTypes, nc, edges, m, f, stim, synType, w, rs
    
    strdef cmd, connectivityPath, connectivityName, presynapticPopulation, postsynapticPopulation
    connectivityTypes = $o1
    cellTypes = $o2
    nConnectivityTypes = connectivityTypes.count()
    if (verbose > 0) {
       printf ("nConnectivityTypes = %d\n", nConnectivityTypes)
   }
   
   for k=0, nConnectivityTypes-1 { 
       connectivityType = connectivityTypes.o(k)
       order = connectivityType.getPropertyScalar("order")
       synType = connectivityType.getPropertyObject("synType")
       wdType = connectivityType.getPropertyScalar("wdType")
       connectivityName = connectivityType.getPropertyString("name").s
       assert(synType.is_undefined == 0, "synType is undefined")
       if (verbose >= 0) {
           printf ("connectivityName = %s\n", connectivityName)
           printf ("wdType = %d\n", wdType)
	   print "synType.is_number = ", synType.is_number
	   print "synType.is_list = ", synType.is_list
           if (synType.is_number) {
   	      printf ("synType.n = %d\n", synType.n)
           }
	   printf ("pnm.nclist.count() = %d\n", pnm.nclist.count())
       }
       presynapticPopulation = connectivityType.getPropertyString("presynapticPopulation").s
       presynapticType = cellTypes.o(findCellType(cellTypes,presynapticPopulation))
       presynapticSize = presynapticType.getPropertyScalar("numCells")
       presynapticOffset = presynapticType.getPropertyScalar("offset")
       postsynapticPopulation = connectivityType.getPropertyString("postsynapticPopulation").s
       postsynapticType = cellTypes.o(findCellType(cellTypes,postsynapticPopulation))
       postsynapticSize = postsynapticType.getPropertyScalar("numCells")
       postsynapticOffset = postsynapticType.getPropertyScalar("offset")
       if (verbose >= 0) {
           printf ("presynaptic population = %s\n", presynapticPopulation)
           printf ("presynaptic size = %d\n", presynapticSize)
           printf ("postsynaptic population = %s\n", postsynapticPopulation)
           printf ("postsynaptic size = %d\n", postsynapticSize)
       }
       sprint(connectivityPath, "%s/%s.dat", datasetPath, connectivityName)
       sprint(cmd, "test -e %s", connectivityPath) // check if the path exist
       // status = system(cmd)
       // on bg-q we can't run system commands on compute node
       status = 0
       if (status == 0) {
           f = new File(connectivityPath)
           pernode = 0
       } else {
           // Try connectivity file local to this MPI node
           sprint(connectivityPath, "%s/%s.%d.dat", datasetPath, connectivityName, pnm.myid+1)
           f = new File(connectivityPath)
           pernode = 1
       }
       f.ropen()
       m = new Matrix()
       if (!f.eof()) {
         m.scanf(f)
       }
       f.close()
       
       if (m.nrow > 1) {
           if (wdType == 1) {
               for i=0, m.nrow-1 {
                   srcid = m.x[i][0]
                   destid = m.x[i][1]
                   if (order == 2) {
                       // relative numbering of pre-synaptic cells, absolute numbering of post synaptic cells
                       assert (srcid <= presynapticSize, "source cell id is greater than the size of the presynaptic population")
                       srcid1 = srcid + presynapticOffset
                       destid1 = destid
                   } else if (order == 1) {
                       // absolute numbering of pre/post synaptic cells
                       srcid1 = srcid
                       destid1 = destid
                   } else if (order == 0) {
                       // relative numbering of pre/post synaptic cells -- add the respective offsets
                       assert (srcid <= presynapticSize, "source cell id is greater than the size of the presynaptic population")
                       assert (destid <= postsynapticSize, "destination cell id is greater than the size of the postsynaptic population")
                       srcid1 = srcid + presynapticOffset
                       destid1 = destid + postsynapticOffset
                   } else {
                       execerror ("unknown connectivity order type")
                   }
                   if (!pnm.gid_exists(destid1)) {
                       if (!pernode) { 
                           continue 
                       } else {
                           printf("host %d: connectivityName = %s wdType = %d destid1 = %d\n", pnm.myid, connectivityName, wdType, destid1)
                           execerror ("destination id mismatch in per-node connectivity file")
                       }
                   }
                   w = m.x[i][2]
                   del = m.x[i][3]
	           nc_appendsyn(pnm, srcid1, destid1, synType, w, del)
                   if (verbose > 0) {
                       printf("host %d: connection src = %d dest = %d w = %g del = %g\n", pnm.myid, srcid1, destid1, w, del)
                   }
	           if (makegraph > 0) {
                       edges = new Vector(2)
                       edges.x[0] = srcid1
                       edges.x[1] = destid1
                       conlist.prepend(edges)
	           }
               }
               m.resize(1,1)
           } else if (wdType == 2) {
               w = connectivityType.getPropertyObject("standardWeight")
               velocity = connectivityType.getPropertyScalar("standardVelocity")
               for i=0, m.nrow-1 {
                   srcid = m.x[i][0]
                   destid = m.x[i][1]
                   
                   if (order == 2) {
                       // relative numbering of pre-synaptic cells, absolute numbering of post synaptic cells
                       assert (srcid <= presynapticSize, "source cell id is greater than the size of the presynaptic population")
                       srcid1 = srcid + presynapticOffset
                       destid1 = destid
                   } else if (order == 1) {
                       // absolute numbering of pre/post synaptic cells
                       srcid1 = srcid
                       destid1 = destid
                   } else if (order == 0) {
                       // relative numbering of pre/post synaptic cells -- add the respective offsets
                       assert (srcid <= presynapticSize, "source cell id is greater than the size of the presynaptic population")
                       assert (destid <= postsynapticSize, "destination cell id is greater than the size of the postsynaptic population")
                       srcid1 = srcid + presynapticOffset
                       destid1 = destid + postsynapticOffset
                   } else {
                       execerror ("unknown connectivity order type")
                   }
                   
                   if (!pnm.gid_exists(destid1)) {
                       if (!pernode) { 
                           continue 
                       } else {
                           printf("host %d: connectivityName = %s wdType = %d destid1 = %d\n", pnm.myid, connectivityName, wdType, destid1)
                           execerror ("destination id mismatch in per-node connectivity file")
                       }
                   }
                   dist = m.x[i][2] + m.x[i][3]
                   del = (dist / velocity) + 1.0
                   if (verbose > 0) {
                       printf("host %d: connection src = %d dest = %d del = %g\n", pnm.myid, srcid1, destid1, del)
                   }
	           nc_appendsyn(pnm, srcid1, destid1, synType, w, del)
	           if (makegraph > 0) {
                       edges = new Vector(2)
                       edges.x[0] = srcid1
                       edges.x[1] = destid1
                       conlist.prepend(edges)
	           }
               }
               m.resize(1,1)
           } else if (wdType == 3) {
               // connectivity of the form:
               // src dest distance layer section node/point
               w = connectivityType.getPropertyObject("weightHistogram")
               velocity = connectivityType.getPropertyScalar("standardVelocity")
               for i=0, m.nrow-1 {
                   srcid   = m.x[i][0]
                   destid  = m.x[i][1]
                   dist    = m.x[i][2]
                   layer   = m.x[i][3]
                   section = m.x[i][4]
                   node    = m.x[i][5]
                   
                   if (order == 2) {
                       // relative numbering of pre-synaptic cells, absolute numbering of post synaptic cells
                       assert (srcid <= presynapticSize, "source cell id is greater than the size of the presynaptic population")
                       srcid1 = srcid + presynapticOffset
                       destid1 = destid
                   } else if (order == 1) {
                       // absolute numbering of pre/post synaptic cells
                       srcid1 = srcid
                       destid1 = destid
                   } else if (order == 0) {
                       // relative numbering of pre/post synaptic cells -- add the respective offsets
                       assert (srcid <= presynapticSize, "source cell id is greater than the size of the presynaptic population")
                       assert (destid <= postsynapticSize, "destination cell id is greater than the size of the postsynaptic population")
                       srcid1 = srcid + presynapticOffset
                       destid1 = destid + postsynapticOffset
                   } else {
                       execerror ("unknown connectivity order type")
                   }
                   
                   if (!pnm.gid_exists(destid1)) {
                       if (!pernode) { 
                           continue 
                       } else {
                           printf("host %d: connectivityName, wdType = %d destid1 = %d\n", pnm.myid, connectivityName, wdType, destid1)
                           execerror ("destination id mismatch in per-node connectivity file")
                       }
                   }
                   del = (dist / velocity) + 1.0
                   if (verbose > 0) {
                       printf("host %d: connection src = %d dest = %d del = %g\n", pnm.myid, srcid1, destid1, del)
                   }
	           nc_appendsyn_lsn(pnm, srcid1, destid1, synType, w, del, layer, section, node)
	           if (makegraph > 0) {
                       edges = new Vector(2)
                       edges.x[0] = srcid1
                       edges.x[1] = destid1
                       conlist.prepend(edges)
	           }
               }
               m.resize(1,1)
           } else if (wdType == 4) {
               w = connectivityType.getPropertyObject("standardWeight")
               velocity = connectivityType.getPropertyScalar("standardVelocity")
               for i=0, m.nrow-1 {
                   srcid = m.x[i][0]
                   destid = m.x[i][1]
                   
                   if (order == 2) {
                       // relative numbering of pre-synaptic cells, absolute numbering of post synaptic cells
                       assert (srcid <= presynapticSize, "source cell id is greater than the size of the presynaptic population")
                       srcid1 = srcid + presynapticOffset
                       destid1 = destid
                   } else if (order == 1) {
                       // absolute numbering of pre/post synaptic cells
                       srcid1 = srcid
                       destid1 = destid
                   } else if (order == 0) {
                       // relative numbering of pre/post synaptic cells -- add the respective offsets
                       assert (srcid <= presynapticSize, "source cell id is greater than the size of the presynaptic population")
                       assert (destid <= postsynapticSize, "destination cell id is greater than the size of the postsynaptic population")
                       srcid1 = srcid + presynapticOffset
                       destid1 = destid + postsynapticOffset
                   } else {
                       execerror ("unknown connectivity order type")
                   }
                   
                   if (!pnm.gid_exists(destid1)) {
                       if (!pernode) { 
                           continue 
                       } else {
                           printf("host %d: connectivityName = %s wdType = %d destid1 = %d\n", pnm.myid, connectivityName, wdType, destid1)
                           execerror ("destination id mismatch in per-node connectivity file")
                       }
                   }                   
		   dist = m.x[i][2]
                   del = (dist / velocity) + 1.0
                   if (verbose > 0) {
                       printf("host %d: connection src = %d dest = %d del = %g\n", pnm.myid, srcid1, destid1, del)
                   }
	           nc_appendsyn(pnm, srcid1, destid1, synType, w, del)
	           if (makegraph > 0) {
                       edges = new Vector(2)
                       edges.x[0] = srcid1
                       edges.x[1] = destid1
                       conlist.prepend(edges)
	           }
               }
               m.resize(1,1)
           }
       }
   }
}



proc mkgap() { local gid, sgid, dgid, w localobj cell, gj, gjlist
    
    gjlist = $o1
    gid    = $2
    branch = $3
    sec    = $4
    ggid   = $5
    w      = $6
    
    cell = pnm.pc.gid2cell(gid)
    
    cell.dendrites[branch][sec] { 
        gj = new ggap(0.5) 
        pnm.pc.target_var(gj, &gj.vgap, ggid)
        pnm.pc.source_var(&v(0.5), ggid)
    }
    
    gjlist.append(gj)
    gj.g = w
    
    if (verbose > 0) {
        printf ("host %d: gap junction: gid = %d branch = %d sec = %d coupling = %g ggid = %d\n", pnm.myid, gid, branch, sec, w, ggid)
    }
    
}

// creates the gap junctions between cells
proc connectgjs() {local i, k, w, ggid, srcid, destid, srcid1, destid1, srcbranch, destbranch, srcsec, destsec, order, srcSize, destSize localobj srcType, destType, cellTypes, gjType, gjTypes, edges, m, f
    
    strdef gjPath, gjName, srcPopulation, destPopulation
    gjTypes = $o1
    cellTypes = $o2
    nGJTypes = gjTypes.count()
    ggid = 2*ncells
    
    for k=0, nGJTypes-1 { 
        gjType = gjTypes.o(k)
        order  = gjType.getPropertyScalar("order")
        gjName = gjType.getPropertyString("name").s
        if (verbose > 0) {
            printf ("gjName = %s\n", gjName)
        }
        srcPopulation  = gjType.getPropertyString("srcPopulation").s
        destPopulation = gjType.getPropertyString("destPopulation").s
        if (verbose > 0) {
            printf ("src population = %s\n", srcPopulation)
            printf ("dest population = %s\n", destPopulation)
        }
        srcType    = cellTypes.o(findCellType(cellTypes,srcPopulation))
        destType   = cellTypes.o(findCellType(cellTypes,destPopulation))
        srcSize    = srcType.getPropertyScalar("numCells")
        destSize   = destType.getPropertyScalar("numCells")
        srcOffset  = srcType.getPropertyScalar("offset")
        destOffset = destType.getPropertyScalar("offset")
        if (verbose > 0) {
            printf ("src size = %d\n", srcSize)
            printf ("dest size = %d\n", destSize)
        }
        sprint(gjPath, "%s/%s.dat", datasetPath, gjName)
        print "gjPath = ", gjPath
        m = new Matrix()
        f = new File(gjPath)
        f.ropen()
        m.scanf(f)
        f.close()
        for i=0, m.nrow-1 {
            srcid = m.x[i][0]
            destid = m.x[i][1]
            srcbranch = m.x[i][2]
            srcsec = m.x[i][3]
            destbranch = m.x[i][4]
            destsec = m.x[i][5]
            w = m.x[i][6]
            
            if (order == 2) {
                // relative numbering of pre-synaptic cells, absolute numbering of post synaptic cells
                assert (srcid <= presynapticSize, "source cell id is greater than the size of the presynaptic population")
                srcid1 = srcid + presynapticOffset
                destid1 = destid
            } else if (order == 1) {
                // absolute numbering of pre/post synaptic cells
                srcid1 = srcid
                destid1 = destid
            } else if (order == 0) {
                // relative numbering of pre/post synaptic cells -- add the respective offsets
                assert (srcid <= presynapticSize, "source cell id is greater than the size of the presynaptic population")
                assert (destid <= postsynapticSize, "destination cell id is greater than the size of the postsynaptic population")
                srcid1 = srcid + presynapticOffset
                destid1 = destid + postsynapticOffset
            } else {
                   execerror ("unknown connectivity order type")
            }
            
	    if (pnm.gid_exists(srcid1)) { 
		mkgap(gjlist, srcid1, srcbranch, srcsec, ggid, w)
                edges = new Vector(2)
                edges.x[0] = srcid1
                edges.x[1] = destid1
                if (makegraph > 0) {
                  gjconlist.prepend(edges)
                }
	    }
            
	    if (pnm.gid_exists(destid1)) {
		mkgap(gjlist, destid1, destbranch, destsec, ggid+1, w)
	    }
            ggid = ggid+2
        }
    }
}


proc mknet() { strdef fname
    startsw()
    mkcells(cellTypes)  // create the cells
    mkcellstime = stopsw()
    if (pnm.myid == 0) {
        printf ("*** Cells created in %g seconds\n", mkcellstime)
    }
    if (vrecordFraction > 0) {
        indicesVrecord = new List()
        mkvrecord(pnm,indicesVrecord,cellTypes,vrecordFraction)
    }
    startsw()
    connectcells(connectivityTypes,cellTypes)  // connect them together
    connectcellstime = stopsw()
    if (pnm.myid == 0) {
        printf ("*** Synaptic connections created in %g seconds\n", connectcellstime)
    }
    startsw()
    connectgjs(gjTypes,cellTypes)  // creates the gap junctions
    connectgjstime = stopsw()
    if (pnm.myid == 0) {
        printf ("*** Gap junctions created\n")
    }
    sprint(fname, "%s/%s_%d.dat", resultsPath, modelName, pnm.myid)
}


mknet(ncells)

// Configuration of artificial spiking cells

objref nc, stimcells
stimcells = new List()
ncstimlist = new List()

proc mkstim() { local i, k, nCellTypes, numStimTypes, gid, stimgid, weight, noise, number, start localobj stim, stimulationParameters, cellType, cellTypes, stimType, stimTypes, edges
    
    strdef templateName, typeName
    
    datasetPath = $s1
    cellTypes = $o2
    nCellTypes = cellTypes.count()
    stimTypes = $o3
    numStimTypes = stimTypes.count()
    
    // Used to keep track of artificial stimulus cells
    typeName = "stim"
    sprint (cmd, "objref %slist, nc%slist", typeName, typeName)
    execute(cmd)
    sprint (cmd, "%slist = new List()", typeName)
    execute(cmd)
    sprint (cmd, "nc%slist = new List()", typeName)
    execute(cmd)
    
    stimgid = maxindex+1+pnm.myid*numStimTypes+1
    
    if (verbose >= 0) {
        printf ("host %d: mkstim: maxindex = %d numStimTypes = %d stimgid = %d\n", pnm.myid, maxindex, numStimTypes, stimgid)
    }
    
    for i=0, numStimTypes-1 { 
        
        stimType = stimTypes.o(i)
        
        indices = stimType.getPropertyObject("indices")
        weight = stimType.getPropertyScalar("weight")
        noise = stimType.getPropertyScalar("noise")
        number = stimType.getPropertyScalar("number")
        start = stimType.getPropertyScalar("start")
        rate = stimType.getPropertyScalar("rate")
        
        if (verbose > 0) {
            printf ("mkstim:  i = %d stimgid = %d weight = %g noise = %g number = %g start = %g rate = %g\n", i, stimgid, weight, noise, number, start, rate)
        }
        
        stim = new StimCell(stimgid,0,0,rate,number,start,noise)
        stimcells.append(stim)
        
        //pnm.register_cell(stimgid, stim)
        //pnm.spike_record(stimgid) 
        
        for j=0, indices.size()-1 {
            
            gid = indices.x(j)
            if (pnm.gid_exists(gid)) {
                nc = stim.connect2target(pnm.pc.gid2obj(gid).syn_stim, nc)
                nc.weight = weight
                nc.delay = 0.1
                ncstimlist.append(nc)
                
	        if (makegraph > 0) {
                    edges = new Vector(2)
                    edges.x[0] = stimgid
                    edges.x[1] = gid
                    conlist.prepend(edges)
	        }	
            }
        }
        stimgid = stimgid+1
        
    }
}

if (stimTypes.count() > 0) {
    mkstim(datasetPath, cellTypes, stimTypes)
}
if (pnm.myid == 0) {
    printf ("*** Stimulus sources initialized\n")
}

if (makegraph > 0) {
        graphout(fname,conlist)
        conlist = nil
}


connectivityTypes = nil
cellTypes = nil

/* Simulation control */


{pnm.pc.setup_transfer()}
{pnm.pc.set_maxstep(10)}

print "dt = ", dt

objref fi_status
if (pnm.myid == 0) {
  fi_status = new FInitializeHandler("simstatus()")
}

objref fi_vrecord
if (vrecordFraction > 0) {
   fi_vrecord = new FInitializeHandler("vrecord(pnm)")
}

stdinit()
if (pnm.myid == 0) {
    printf ("*** Initialization completed\n")
}

proc prun() { localobj po
  if(dump_coreneuron_model) {
     pnm.pc.nrnbbcore_write("coredat")
  }
  if (use_coreneuron || use_gpu) {
    if (!nrnpython("from neuron import coreneuron")) {
      printf("Python not available to run CoreNEURON\n")
      return
    }
    po = new PythonObject()
    po.coreneuron.enable = 1
    if(use_gpu) {
      po.coreneuron.gpu = 1
    }
    if(pnm.pc.id == 0) {
      printf("nrncore_arg: |%s|\n", po.coreneuron.nrncore_arg(tstop))
    }
  } else {
    tsav = startsw()
  }
  pnm.pc.psolve(tstop)
  if (!use_coreneuron && pnm.pc.id == 0) {
    printf("psolve time %g\n", startsw() - tsav)
  }

  if (pnm.myid == 0) {
    printf ("*** Simulation completed\n")
  }
}

prun()

/* Write simulation results */

strdef spikeout_fname
sprint(spikeout_fname, "%s/%s/%s_spikeout_%d.dat", outdir, resultsPath, modelName, pnm.myid)
spikeout(spikeout_fname,pnm.spikevec,pnm.idvec)

strdef vrecord_fname
if (vrecordFraction > 0) {
  sprint(vrecord_fname, "%s/%s/%s_vrecord_%d.dat", outdir, resultsPath, modelName, pnm.myid)
  vrecordout(vrecord_fname,indicesVrecord)
}

comptime = pnm.pc.step_time
avgcomp  = pnm.pc.allreduce(comptime, 1)/pnm.pc.nhost
maxcomp  = pnm.pc.allreduce(comptime, 2)

if (pnm.myid == 0) {
    printf("Execution time summary for host 0:\n")
    printf("  created cells in %g seconds\n", mkcellstime)
    printf("  connected cells in %g seconds\n", connectcellstime)
    printf("  created gap junctions in %g seconds\n", connectgjstime)
    printf("  ran simulation in %g seconds\n", comptime)
    if (maxcomp > 0) {
        if (pnm.myid == 0) { printf("  load_balance = %g\n", avgcomp/maxcomp) }
    }
}

if (pnm.myid == 0) {
    printf ("*** Simulation results written\n")
}

pnm.pc.barrier()

if (pnm.myid == 0) {
    printf ("*** Ending simulation\n")
}