/************************************************************
For more information, consult ModelDoc.pdf
************************************************************/
loadstart = startsw() // record the start time of the set up
/***********************************************************************************************
I. LOAD LIBRARIES
***********************************************************************************************/
{load_file("nrngui.hoc")} // Standard definitions - NEURON library file
{load_file("netparmpi.hoc")} // Contains the template that defines the properties of
// the ParallelNetManager class, which is used to set
// up a network that runs on parallel processors
{load_file("./setupfiles/ranstream.hoc")} // Contains the template that defines a RandomStream
// class used to produce random numbers
// for the cell noise (what type of noise?)
{load_file("./setupfiles/CellCategoryInfo.hoc")} // Contains the template that defines a
// CellCategoryInfo class used to store
// celltype-specific parameters
{load_file("./setupfiles/SynStore.hoc")} // Contains the template that defines a
{load_file("./setupfiles/defaultvar.hoc")} // Contains the proc definition for default_var proc
{load_file("./setupfiles/parameters.hoc")} // Loads in operational and model parameters that can
// be changed at command line
{load_file("./setupfiles/set_other_parameters.hoc")}// Loads in operational and model parameters
// that can't be changed at command line
//{load_file("./setupfiles/check_for_files.hoc")} // Checks that the necessary chosen files exist
// (datasets, connectivity, and stimulation)
default_var("gidOI", 0) // Run Comments // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("cellindOI", 0) // Run Comments // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("spkflag", 0) // Run Comments // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("stimflag", 0) // Run Comments // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("resultsfolder", "00001") // Run Comments // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("runname", "None") // Run Comments // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("origRunName", "CutSma_040_11") // Run Comments // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("celltype2use", "axoaxoniccell") // Run Comments // -c 'strdef RunName ' 'RunName = "Woo"'
celltypeOI=cellindOI
//gidOI=0
i=0
ij=0
gid=0
/*strdef origRunName
origRunName = "CutSma_040_11"
PrintConnDetails=2*/
/***********************************************************************************************
II. SET MODEL SIZE, CELL DEFINITIONS
***********************************************************************************************/
celsius=34
{load_file("./setupfiles/load_cell_category_info.hoc")} // Reads the 'cells2include.hoc' file and
// loads the information into one
// 'CellCategoryInfo' object for each cell
// type (bio or art cells?). Info includes
// number of cells, gid ranges, type name
strdef path2ConnData
sprint(path2ConnData, "../networkclamp_results/%s/%s/conndata.dat", origRunName, resultsfolder)
strdef path2SynData
sprint(path2SynData, "../networkclamp_results/%s/%s/syndata.dat", origRunName, resultsfolder)
{load_file("./setupfiles/load_cell_conns.hoc")} // Load in the cell connectivity info
{load_file("./setupfiles/load_cell_syns.hoc")} // Load in the cell connectivity info
strdef tempFileStr // Define a string reference to store the name of the
// current cell template file
proc loadCellTemplates(){local i // Proc to load the template that defines (each) cell class
//for i=0, numCellTypes-1 { // Iterate over each cell type in cells2include (and art cells?)
i = celltypeOI
sprint(tempFileStr,"./cells/class_%s.hoc",celltype2use) // Concatenate the
//sprint(tempFileStr,"./cells/class_%s.hoc",cellType[i].technicalType) // Concatenate the
// path and file
load_file(tempFileStr) // Load the file with the template that defines the class
// for each cell type
//}
if (stimflag==1) {
load_file("./cells/class_sintrain.hoc") // Load the file with the template that defines the class
//tstop = 100
//SimDuration = tstop
print "Using PhasicData, SimDuration = ", SimDuration
} else {
load_file("./cells/class_ppvec.hoc") // Load the file with the template that defines the class
}
}
loadCellTemplates() // Run the newly defined proc
proc calcNetSize(){local i // Calculate the final network size (after any cell death)
//cellType[0].numCells = cellType[numCellTypes-1].cellEndGid - cellType[0].cellEndGid + 4*2 // just added now for spontburst
//cellType[0].updateGidRange(0)
totalCells = 0 // Initialize totalCells (which counts the number of 'real'
// cells) so we can add to it iteratively in the 'for' loop
ncell = cellType[0].numCells // Initialize ncell (which counts all 'real' and 'artificial'
// cells) so we can add to it iteratively in the 'for' loop
for i=1,numCellTypes-1 { // Run the following code for 'real' cell types only - need a different way of singling out real cells?
if (cellType[i].layerflag==0) { // For cell types in layer 0, which are susceptible to death
cellType[i].numCells = int(cellType[i].numCells * ((100-PercentCellDeath)/100))
// Calculate the number of cells surviving after cell loss
if (cellType[i].numCells == 0) {cellType[i].numCells = 1} // If all cells of one type
// are killed, let 1 live
}
cellType[i].updateGidRange(cellType[i-1].cellEndGid+1) // Update the gid range for each
// cell type
totalCells = totalCells + cellType[i].numCells // Update the total number of cells
// after sclerosis, not including
// artificial cells
ncell = ncell + cellType[i].numCells // Update the total number of cells
// after sclerosis, including
// artificial cells
}
}
calcNetSize()
proc calcBinSize(){local NumGCells
for i=0, numCellTypes-1 { // Using the specified dimensions of the network (in um) and
// the total number of cells of each type, set the number
// of bins in X, Y, Z dimensions such that the cells will be
// evenly distributed throughout the allotted volume
// just changed this so even the stim cells will be allotted, as now we have some
// stimulation protocols that incorporate stim cell position
cellType[i].setBins(scLongitudinalLength,scTransverseLength,LayerVector.x[cellType[i].layerflag])
// For the z length, use the height of the layer in which the
// cell somata are found for this cell type
}
}
calcBinSize()
/***********************************************************************************************
III.SET UP PARALLEL CAPABILITY AND WRITE OUT RUN RECEIPT
***********************************************************************************************/
objref pnm, pc, nc, nil
proc parallelizer() {
pnm = new ParallelNetManager(ncell) // Set up a parallel net manager for all the cells
pc = pnm.pc
pnm.round_robin() // Incorporate all processors - cells 0 through ncell-1
// are distributed throughout the hosts
// (cell 0 goes to host 0, cell 1 to host 1, etc)
}
parallelizer()
iterator pcitr() {local i2, startgid // Create iterator for use as a standard 'for' loop
// throughout given # cells usage:
// for pcitr(&i1, &i2, &gid, it_start, it_end) {do stuff}
// it_start and it_end let you define range over
// which to iterate
// i1 is the index of the cell on the cell list for that host
// i2 is the index of that cell for that cell type on that host
numcycles = int($4/pc.nhost)
extra = $4%pc.nhost
addcycle=0
if (extra>pc.id) {addcycle=1}
startgid=(numcycles+addcycle)*pc.nhost+pc.id
i1 = numcycles+addcycle // the index into the cell # on this host.
i2 = 0 // the index of the cell in that cell type's list on that host
if (startgid<=$5) {
for (i3=startgid; i3 <= $5; i3 += pc.nhost) { // Just iterate through the cells on
// this host(this simple statement
// iterates through all the cells on
// this host and only these cells because
// the roundrobin call made earlier dealt
// the cells among the processors in an
// orderly manner (like a deck of cards)
$&1 = i1
$&2 = i2
$&3 = i3
// if (i3 == $5) {iterator_statement}
iterator_statement
i1 += 1
i2 += 1
}
}
}
iterator fastitr() {local i2, startgid // Create iterator for use as a standard 'for' loop
// throughout given # cells usage:
// for fastitr(&i1, &i2, &gid, it_start, it_end) {do stuff}
// it_start and it_end let you define range over
// which to iterate
// i1 is the index of the cell on the cell list for that host
// i2 is the index of that cell for that cell type on that host
numcycles = int($4/pc.nhost)
extra = $4%pc.nhost
addcycle=0
if (extra>pc.id) {addcycle=1}
startgid=(numcycles+addcycle)*pc.nhost+pc.id
i1 = numcycles+addcycle // the index into the cell # on this host.
i2 = 0 // the index of the cell in that cell type's list on that host
if (startgid<=$5) {
for (i3=startgid; i3 <= $5; i3 += pc.nhost) { // Just iterate through the cells on
// this host(this simple statement
// iterates through all the cells on
// this host and only these cells because
// the roundrobin call made earlier dealt
// the cells among the processors in an
// orderly manner (like a deck of cards)
$&1 = i1
$&2 = i2
$&3 = i3
if (i3 == $5) {iterator_statement}
// iterator_statement
i1 += 1
i2 += 1
}
}
}
objref strobj
strobj = new StringFunctions()
strdef direx
if (strcmp(UID,"0")==0 && pc.id==0) {
type = unix_mac_pc() // 1 if unix, 2 if mac, 3 if mswin, or 4 if mac osx darwin
if (type==3) {
{system("cscript //NoLogo setupfiles/uuid.vbs", direx)} // pc
} else {
{system("uuidgen", direx)} // unix or mac
strobj.left(direx, strobj.len(direx)-1)
}
UID = direx
}
strdef cmdstr, cmd
//{load_file("./setupfiles/save_run_info.hoc")}
loadtime = startsw() - loadstart // Calculate the set up time (now - recorded start time) in seconds
if (pc.id == 0) {printf("\nTIME HOST 0: %g seconds (set up)\n************\n", loadtime)}
createstart = startsw() // Record the start time of the cell creation
/***********************************************************************************************
IV. CREATE, UNIQUELY ID, AND POSITION CELLS
***********************************************************************************************/
objref memfile, topfile //, pccc
//{pccc = new ParallelContext()}
{zzz=0}
if (pc.id==0) {
{sprint(cmd,"./networkclamp_results/%s/%s/memory.dat", runname, resultsfolder)}
{memfile = new File(cmd)}
{memfile.wopen()}
{sprint(cmd,"./networkclamp_results/%s/%s/topoutput.dat", runname, resultsfolder)}
{topfile = new File(cmd)}
{topfile.wopen()}
}
strdef direx1, direx2, TopCommand
func mallinfo() {local m // this function is a wrapper for the system function mallinfo
m = nrn_mallinfo(0)
if (pc.id == 0) {
if (PrintTerminal>1) {printf("Memory (host 0) - %s: %ld. Since last call: %ld\n", $s2, m, m - $1)}
memfile.printf("%s\t%f\t%ld\t%ld\n", $s2, startsw()-loadstart, m, m - $1) // Print the spike time and spiking cell gid
sprint(TopCommand,"top -p `pgrep %s | head -n20 | tr \"\\n\" \",\" | sed 's/,$//'` -n1 -b | tail -n2 | head -n1", TopProc)
//print "TopCommand = ", TopCommand
if (strcmp(TopProc,"")!=0) { system(TopCommand, direx1)}
//system("top -p `pgrep special | head -n20 | tr \"\\n\" \",\" | sed 's/,$//'` -n1 -b | tail -n2 | head -n1",direx1)
//system("top -p `pgrep nrniv | head -n20 | tr \"\\n\" \",\" | sed 's/,$//'` -n1 -b | tail -n2 | head -n1",direx2)
//if (strcmp(direx1,direx2)<0) {
// topfile.printf("%s\t%s\n", $s2, direx2) // Print the spike time and spiking cell gid
//} else {
topfile.printf("%s\t%s\n", $s2, direx1) // Print the spike time and spiking cell gid
//}
}
return m
}
{zzz = mallinfo(zzz, "Prior to creating cells")}
objref cells, ransynlist, ranstimlist, raninitlist
cells = new List()
ransynlist = new List()
ranstimlist = new List()
raninitlist = new List()
{load_file("./setupfiles/clamp/netclamp_create_cells_pos.hoc")} // Creates each cell on its assigned host
// and sets its position using the algorithm
// defined above
strdef cmd
createtime = startsw() - createstart // Calculate time taken to create the cells
if (pc.id == 0) {printf("\nTIME HOST 0: %g seconds (created cells)\n************\n", createtime)}
connectstart = startsw() // Grab start time of cell connection
oldtimeout = pc.timeout(0)
zzz = mallinfo(zzz, "After creating cells")
/***********************************************************************************************
V. CONNECT THE MODEL CELLS AND CONNECT THE STIMULATION CELLS TO SOME MODEL CELLS
***********************************************************************************************/
celsius=34
{load_file("./setupfiles/nc_append_functions.hoc")} // Defines nc_append and nc_appendo, which
// are procs (?) needed to create the netcon
// object associated with each connection
// (synapse)
{pc.barrier()}
{load_file("./setupfiles/clamp/list_randfaststim_connections.hoc")} // Loads in a particular definition of the connectCells function,
// depending on the value of the Connectivity string
zzz = mallinfo(zzz, "After connecting cells")
{pc.barrier()}
if (stimflag==0) {
{load_file("./setupfiles/clamp/netclamp_stimulation.hoc")} // Loads in a particular stimulation protocol, including the spike
// vectors for the artificial cells and their connections to the
// 'real' cells (where do we define the # of art cells and create them?)
} else {
if (stimflag==1) {
{load_file("./setupfiles/clamp/netclamp_oscillation_stimulation.hoc")}
}
}
zzz = mallinfo(zzz, "After defining stimulation")
{load_file("./setupfiles/write_conns_functions.hoc")} // Defines two funcs (procs?): printNumConFile()
// writes the number of synapses between every
// combination of pre- and post-synaptic cell type
// and tracenet() lists the presynaptic cell type,
// postsynaptic cell type, and synapse type for
// every connection
if (PrintConnSummary==1) {printNumConFile()} // Write "numcons.dat": a QUICK, SMALL summary file of the # of
// connections by pre and post cell types
/*
if (PrintConnDetails==2) {allCellConnsParallel()} // Write "connections.dat": a VERY SLOW, LARGE file of all cell
// connections (pre and post gids, synapse type, host on which
// (target?) cell exists)
if (PrintConnDetails==1) {allCellConnsParallel()} // Write the file of all cell connections (pre and post gids, synapse type,
// host on which cell exists), "connections.dat"
if ((PrintConnDetails==0) && (PrintVoltage==2)) {recordedCellConnsParallel()} // If detailed connections are not being written out, at least write them out for cells that are being traced
*/
connecttime = startsw() - connectstart // calculate time taken to connect the cells
if (pc.id == 0) {printf("\nTIME HOST 0: %g seconds (connected cells)\n************\n", connecttime)}
/***********************************************************************************************
VI. INITIALIZE AND RUN NETWORK, OUTPUT RESULT FILES
***********************************************************************************************/
objref mytrace, myvrec, cell
objref lfpkmatrix
lfpkmatrix = new Vector()
cell = pc.gid2cell(gidOI)
/*********************
inserted stuff right below here:
************/
/*
Computes xyz coords of nodes in a model cell whose topology & geometry are defined by pt3d data.
Code by Ted Carnevale.
*/
proc interpxyz() { local ii, nn, kk, xr, nsegs localobj xx, yy, zz, xint, yint, zint, ll, range
// get the data for the section
nn = $1
nsegs = $2
xx = $o3
yy = $o4
zz = $o5
ll = $o6
xint = $o7
yint = $o8
zint = $o9
// to use Vector class's .interpolate()
// must first scale the independent variable
// i.e. normalize length along centroid
ll.div(ll.x[nn-1])
// initialize the destination "independent" vector
range = new Vector(nsegs+2)
range.indgen(1/nsegs)
range.sub(1/(2*nsegs))
range.x[0]=0
range.x[nsegs+1]=1
// length contains the normalized distances of the pt3d points
// along the centroid of the section. These are spaced at
// irregular intervals.
// range contains the normalized distances of the nodes along the
// centroid of the section. These are spaced at regular intervals.
// Ready to interpolate.
xint.interpolate(range, ll, xx)
yint.interpolate(range, ll, yy)
zint.interpolate(range, ll, zz)
}
proc setup_npole_lfp() { local i, j, k, m, n, r, h, s, avglfp, typei, celltype, ex, ey, ez, sx, sy, sz, lfptype, ii, nn localobj myr, dists, lst, xx, yy, zz, ll, xint, yint, zint strdef ks
rho = 333.0 // extracellular resistivity, [ohm cm]
fdst = 0.1 // percent of distant cells to include in the computation
ex = mypoint.x[0] //$1
ey = mypoint.x[1] //$2
ez = mypoint.x[2] //$3
//printf ("host %d: entering setup_npole_lfp\n", pc.id)
// Vector for storing longitudinal and perpendicular distances
// between recording electrode and compartments
dists = new Vector(2)
m = 0
n = 0
// Relative to the recording electrode position
forsec cell.all {
insert extracellular
n = n + nseg
}
lfpkmatrix = new Vector(n+1)
//printf ("host %d: npole_lfp: lfpkmatrix.nrow = %d lfpkmatrix.ncol = %d\n", pc.id, lfpkmatrix.nrow, lfpkmatrix.ncol)
// Iterates over each compartment of the cell
j = 0
forsec cell.all {
if (ismembrane("extracellular")) {
nn = n3d()
xx = new Vector(nn)
yy = new Vector(nn)
zz = new Vector(nn)
ll = new Vector(nn)
for ii = 0,nn-1 {
xx.x[ii] = x3d(ii)
yy.x[ii] = y3d(ii)
zz.x[ii] = z3d(ii)
ll.x[ii] = arc3d(ii)
}
xint = new Vector(nseg+2)
yint = new Vector(nseg+2)
zint = new Vector(nseg+2)
interpxyz(nn,nseg,xx,yy,zz,ll,xint,yint,zint)
j = 0
for (x,0) {
sx = xint.x[j]
sy = yint.x[j]
sz = zint.x[j]
// l = L/nseg is compartment length
// r is the perpendicular distance from the electrode to a line through the compartment
// h is longitudinal distance along this line from the electrode to one end of the compartment
// s = l + h is longitudinal distance to the other end of the compartment
l = L/nseg
r = sqrt((ex-sx)*(ex-sx) + (ey-sy)*(ey-sy) + (ez-sz)*(ez-sz))
h = l/2
s = l + h
k = 0.0001 * area(x) * (rho / (4 * PI * l)) * abs(log((sqrt(h*h + r*r) - h) / (sqrt(s*s + r*r) - s)))
sprint(ks,"%g",k)
if ((strcmp(ks,"nan") == 0) || (strcmp(ks,"-nan") == 0)) {
k = 0
}
lfpkmatrix.x[j] = k
j = j + 1
}
}
}
}
setup_npole_lfp()
/************** done inserting *********/
func npole_lfp() { local vlfp
vlfp = 0 // Initialize the LFP variable
forsec cell.all {
if (ismembrane("extracellular")) {
j = 0
for (x,0) {
vlfp = vlfp + (i_membrane(x) * lfpkmatrix.x[j])
j = j + 1
}
}
} // vlfp = cell.apical[cell.NumApical-1].v(0.5) - cell.basal[0].v(0.5)
return vlfp
}
objref lfplist, lfptracelist, vec
lfplist = new List()
lfptracelist = new List()
lfp_dt = 0.5 // sampling interval for LFP
objref mytraceidxlist, lfptracelist
proc sample_lfp() { local i, avglfp, startgid, endgid, gid, celltype localobj vec, lst
avglfp = npole_lfp(MaxEDist) // x, y, z position of electrode in microns
if (pc.id == 0) {
vec = new Vector()
vec.append(t, avglfp)
lfplist.append(vec.c)
}
cvode.event(t + lfp_dt, "sample_lfp()") // execute sample_lfp lfp_dt ms from now
}
strdef filestr
sprint(filestr, "./networkclamp_results/%s/%s/allsyns.dat", origRunName, resultsfolder)
objref fAll, fSyn
fAll = new File(filestr)
fAll.wopen()
sprint(filestr, "./networkclamp_results/%s/%s/synlist.dat", origRunName, resultsfolder)
fSyn = new File(filestr)
fSyn.wopen()
strdef cmdstr, tempFileStr
objref newSecRef
if (pc.gid_exists(gidOI)) {
cell = pc.gid2cell(gidOI)
print "cell is located at: ", cell.x, ", ", cell.y, ", ", cell.z
mytrace = new Vector((tstop-tstart)/dt)
mytrace.record(&cell.soma.v(0.5))
ist = 0
ien = 0
myi = 0
access cell.soma[0]
{distance()}
for precellType = 0, numCellTypes-1 {
ist = myi
for r=0, cellType[celltypeOI].SynList[precellType].count()-1 {
sprint(cmdstr,"newSecRef = cell.%s",cellType[celltypeOI].SynList[precellType].object(r).SecRefStr)
execute(cmdstr) // sets newSecRef
forsec newSecRef {
for (x,0) {
//y=0
execute(cellType[celltypeOI].SynList[precellType].object(r).CondStr)
if (y==1) {
fSyn.printf("%s %s %g %s\n", cellType[celltypeOI].cellType_string, cellType[precellType].cellType_string, myi, secname())
myi = myi + 1
}
}
}
}
ien = myi - 1
if (ien>=ist) {fAll.printf("%s %s %g %g\n", cellType[celltypeOI].cellType_string, cellType[precellType].cellType_string, ist, ien)}
}
fAll.close()
fSyn.close()
/*forsec cell.all {
insert extracellular
insert xtra
}
myvrec = new Vector((tstop-tstart)/dt)
myvrec.record(&lfp_xtra)*/
} else {
print "couldnt make cell ", gidOI
}
/*
load_file("clamp/interpxyz.hoc") // only interpolates sections that have extracellular
load_file("clamp/setpointers.hoc") // automatically calls grindaway() in interpxyz.hoc
load_file("clamp/calcrxc.hoc") // computes transfer r between segments and recording electrodes
*/
proc init() { local dtsav, temp, secsav, secondordersav // Redefine the proc that initializes the
// simulation (why?)
dtsav = dt // Save the desired dt value to reset it after temporarily
// changing dt to run a quick presimulation to allow the
// network to reach steady state before we start 'recording'
secondordersav = secondorder // Save the desired secondorder value to reset it after
// temporarily changing secondorder to run a quick presimulation
// to allow the network to reach steady state before we start
// 'recording' its results
finitialize(v_init) // Call finitialize from within the custom init proc, just as the default
// init proc does. Note that finitialize will call the INITIAL block for
// all mechanisms and point processes inserted in the sections and set the
// initial voltage to v_init for all sections
t = -200 // Set the start time for (pre) simulation; -500 ms to allow the network to
// reach steady state before t = 0, when the real simulation begins
dt= 10 // Set dt to a large value so that the presimulation runs quickly
secondorder = 0 // Set secondorder to 0 to set the solver to the default fully implicit backward
// euler for numerical integration (see NEURON ref)
temp= cvode.active() // Check whether cvode, a type of solver (?) is on
if (temp!=0) {cvode.active(0)} // If cvode is on, turn it off while running the presimulation
while(t<-100) { fadvance() if (PrintTerminal>2) {print t}} // Run the presimulation from t = -500
// to t = -100 (why not 0?) to let the
// network and all its components reach
// steady state. Integrate all section
// equations over the interval dt,
// increment t by dt, and repeat until
// t at -100
if (temp!=0) {cvode.active(1)} // If cvode was on and then turned off, turn it back on now
t = tstart // Set t to the start time of the simulation
dt = dtsav // Reset dt to the specified value for the simulation
secondorder = secondordersav // Reset secondorder to the specified value for the simulation
if (cvode.active()){
cvode.re_init() // If cvode is active, initialize the integrator
} else {
fcurrent() // If cvode is not active, make all assigned variables
// (currents, conductances, etc) consistent with the
// values of the states
}
frecord_init() // see email from ted - fadvance() increments the recorder, so we need to fix the index it ends up at
}
{load_file("./setupfiles/write_results_functions.hoc")} // Defines the procs spikeout() and timeout()
// which write out the spikeraster (spikeout)
// and the time taken to run each section of
// the code (timeout)
use_cache_efficient=1
get_spike_hist=0
use_bin_queue=0
use_spike_compression=0
if (use_spike_compression==1) {
maxstepval = 2.5
} else {
maxstepval = 10
}
proc rrun(){ // Run the network simulation and write out the results
local_minimum_delay = pc.set_maxstep(maxstepval) // Set every machine's max step size to minimum delay of
// all netcons created on pc using pc.gid_connect, but
// not larger than 10
if (pc.id==0 && use_spike_compression==1) {print "Host ", pc.id, " has local_minimum_delay=", local_minimum_delay}
stdinit() // Call the init fcn (which is redefined in this code) and
// then make other standard calls (to stop watch, etc)
runstart = startsw() // grab start time of the simulation
pc.psolve(tstop) // Equivalent to calling cvode.solve(tstop) but for parallel NEURON;
// solve will be broken into steps determined by the result of
// set_maxstep
runtime = startsw() - runstart // Calculate runtime of simulation
// Print a time summary message to screen
writestart = startsw()
comptime = pc.step_time
avgcomp = pc.allreduce(comptime, 1)/pc.nhost
maxcomp = pc.allreduce(comptime, 2)
if (maxcomp>0) {
if (pc.id == 0) { printf("load_balance = %.3f\n", avgcomp/maxcomp)}
if (pc.id == 0) { printf("exchange_time = %.2f ms\n", runtime - maxcomp) }
} else {
if (pc.id == 0) { printf("no load balance info available\nno spike exchange info available\n")}
}
/*
spikeoutfast()
TimeOutSerial() // Write out a file of run times for each code section
if (PrintVoltage>0) {voltageout()} // Write out any voltages that were recorded during simulation
if (PrintConnDetails==1) {allCellConnsParallel()} // Write the file of all cell connections (pre and post gids, synapse type,
// host on which cell exists), "connections.dat"
if ((PrintConnDetails==0) && (PrintVoltage==1)) {recordedCellConnsParallel()} // If detailed connections are not being written out, at least write them out for cells that are being traced
SummaryOut(runtime) // Write out a file with the total number of cells, spikes, and connections (including contributions from artificial cells)
*/
writetime = startsw() - writestart // Calculate write time of program
totaltime = startsw() - loadstart
allottedtime = JobHours*3600
if (pc.id == 0) {printf("****\nTIME SUMMARY for host 0\nset up in %g seconds\ncreated cells in %g seconds\nconnected cells in %g seconds\nran simulation in %g seconds\nwrote results in %g seconds\nTOTAL TIME = %g seconds\nALLOTTED TIME = %g seconds\n************\nThis run is called: %s %s\n************\n", loadtime, createtime, connecttime, runtime, writetime, totaltime, allottedtime, origRunName, resultsfolder)}
}
objref fih
// ComputeNpoleLFP=1 ComputeDipoleLFP=0
if (ComputeNpoleLFP > 0) {
// execute sample_lfp() at t = 0,
// right after the mechanism INITIAL blocks have been executed
fih = new FInitializeHandler("sample_lfp()")
}
objref fihw
fihw = new FInitializeHandler(2, "midbal()")
walltime = startsw()
strdef cmd, cmdo
newtstop = tstop
warningflag=0
proc midbal() {local wt, thisstep
wt = startsw()
if (t>0) {
thisstep = wt - walltime
simleft = tstop - t
compleft = JobHours*3600 - (startsw() - loadstart)
//print "t=",t,", tstop=",tstop,", simleft=",simleft,", jobtime=",JobHours*3600,", timeused=",startsw() - loadstart,", compleft=", compleft
if (warningflag==0 && (simleft/StepBy*thisstep+EstWriteTime)>compleft && pc.id == 0) {
newtstop = int((compleft-EstWriteTime)/thisstep)*StepBy + t
print "Not enough time to complete ", tstop, " ms simulation, simulation will likely stop around ", newtstop, " ms"
warningflag=1
}
if (pc.id == 0) { printf("%g ms interval at t=%g ms was %g s\n", StepBy, t, thisstep) }
if ((thisstep+EstWriteTime)>compleft && t<tstop) { // not enough time for another step, end now
{pc.barrier()}
sprint(cmd,"cat results/%s/runreceipt.txt | sed -e 's/^SimDuration = [^ ]*/SimDuration = %g;/' > x ; mv x results/%s/runreceipt.txt", RunName, t, RunName)
if (pc.id == 0) { {system(cmd,cmdo)} }
if (pc.id == 0) { print "simulation stopped early at t=", t, " ms"}
tstop = t
stoprun=1
}
}
walltime = wt
cvode.event(t+StepBy, "midbal(0)")
}
{cvode.cache_efficient(use_cache_efficient)} // always double check that this addition does not affect the spikeraster (via pointers in mod files, etc)
if (use_bin_queue==1) {
use_fixed_step_bin_queue = 1 // boolean
use_self_queue = 0 // boolean - this one may not be helpful for me, i think it's best for large numbers of artificial cells that receive large numbers of inputs
{mode = cvode.queue_mode(use_fixed_step_bin_queue, use_self_queue)}
}
if (use_spike_compression==1) {
nspike = 3 // compress spiketimes or not
gid_compress = 0 //only works if fewer than 256 cells on each proc
{nspike = pc.spike_compress(nspike, gid_compress)}
}
objref spkhist
if (get_spike_hist==1) {
spkhist = new Vector(pc.nhost)
if (pc.id==0) {
pc.max_histogram(spkhist)
}
}
zzz = mallinfo(zzz, "Before running simulation")
{rrun()} // Run the network simulation (in proc rrun)
zzz = mallinfo(zzz, "After running simulation")
if (pc.id==0) {
{memfile.close()}
{topfile.close()}
}
objref f4
strdef cmd
if (pc.id==0 && get_spike_hist==1) {
sprint(cmd,"./results/%s/spkhist.dat", RunName)
f4 = new File(cmd)
f4.wopen()
// Open for appending to file
for i=0, pc.nhost-1 {
f4.printf("%d\t%d\n", i, spkhist.x[i]) // Print the spike time and spiking cell gid
}
f4.close()
}
strdef outfile
objref f
if (pc.gid_exists(0)) { // If cell exists on this machine
sprint(outfile, "./networkclamp_results/%s/%s/mytrace_%d_soma.dat", runname, resultsfolder, gidOI)
//ssprint(outfile, "./results/%s/mytrace_%d_soma.dat", RunName, gidOI)
f = new File(outfile)
f.wopen()
f.printf("t\tv\n")
for i=0, (tstop-tstart)/dt-1 {
f.printf("%g\t%g\n", i*dt, mytrace.x[i])
}
f.close()
/*
sprint(outfile, "../networkclamp_results/%s/%s/myvrec_%d_soma.dat", RunName, resultsfolder, gidOI)
f = new File(outfile)
f.wopen()
f.printf("t\tv\n")
for i=0, (tstop-tstart)/dt-1 {
f.printf("%g\t%f\n", i*dt, myvrec.x[i])
}
f.close()*/
}
proc SpikeOutSerial() {local i, rank localobj f // Write out a spike raster (cell, spike time)
pc.barrier() // Wait for all ranks to get to this point
sprint(cmd,"./networkclamp_results/%s/%s/spikeraster.dat", runname, resultsfolder)
f = new File(cmd)
if (pc.id == 0) { // Write header to file 1 time only
f.wopen()
f.close()
}
for rank = 0, pc.nhost-1 { // For each processor, allow processor to append to file the spike times of its cells
if (rank == pc.id) { // Ensure that each processor runs once
f.aopen() // Open for appending to file
for i=0, pnm.idvec.size-1 {
f.printf("%.3f %d\n", pnm.spikevec.x[i], pnm.idvec.x[i]) // Print the spike time and spiking cell gid
}
f.close()
}
pc.barrier()
}
}
SpikeOutSerial()
proc allCellConnsSerial() {local i, rank, srcid localobj tgt, f // Write out the connections list, including synapse type
pc.barrier() // Wait for all ranks to get to this point
sprint(cmd,"./networkclamp_results/%s/%s/connections.dat", runname, resultsfolder)
f = new File(cmd)
if (pc.id == 0) { // Write header to file 1 time only
f.wopen()
f.printf("source\ttarget\tsynapse\n")
f.close()
}
for rank = 0, pc.nhost-1 { // For each processor, allow processor to append to file its connection info
if (rank == pc.id) { // Ensure that each processor runs once
f.aopen() // Open for appending to file
for i=0, nclist.count -1 { // For each connection in the list
srcid = nclist.o(i).srcgid() // Get the gid of the source cell
tgt = nclist.o(i).syn // Get a reference to the target cell
f.printf("%d\t%d\t%d\n", srcid, tgt.cid, tgt.sid) // Print source gid, target gid, synapse id
}
f.close()
}
pc.barrier()
}
}
allCellConnsSerial()
proc writeCells() {
pc.barrier() // Wait for all ranks to get to this point
sprint(cmd,"./networkclamp_results/%s/%s/celltype.dat", runname, resultsfolder)
f = new File(cmd)
if (pc.id == 0) { // Write header to file 1 time only
f.wopen()
f.printf("celltype\ttypeIndex\trangeStart\trangeEnd\n")
for i=0, numCellTypes-1 {
f.printf("%s\t%s\t%d\t%d\t%d\n", cellType[i].cellType_string, cellType[i].technicalType, i, cellType[i].cellStartGid, cellType[i].cellEndGid)
}
f.close()
}
}
writeCells()
proc writeLFP() { localobj f, lfpvec, lfptrace
sprint(cmd,"./networkclamp_results/%s/%s/lfp.dat", runname, resultsfolder)
f = new File(cmd)
f.wopen()
// Open for appending to file
for i=0, lfplist.count()-1 {
lfpvec = lfplist.o(i)
f.printf("%g\t%g\n", lfpvec.x[0], lfpvec.x[1]) // Prints time and average LFP value
}
f.close()
}
writeLFP()
{pc.barrier}
{quit()}
/*
// Old quit method
{pc.runworker()} // Everything after this line is executed only by the host node
// The NEURON website describes this as "The master host returns immediately. Worker hosts
// start an infinite loop of requesting tasks for execution."
{pc.done()} // Sends the quit message to the worker processors, which then quit NEURON
quit() // Sends the quit message to the host processor, which then quits NEURON
*/