/************************************************************ 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 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 */