# plot.py # Mark Rowan, School of Computer Science, University of Birmingham, UK # April 2012 # Use NEURON / neurosim's grvec hoc routines to read saved spike data # and then to allow manipulation and plotting in a Python+pyplot environment # Run interactively, where is a directory containing spks, .spks, aux and vars files: # python plot.py # or in batch: # python plot.py [options] [plots] # Options and plot types are listed at the bottom of this file, and can be # viewed by running interactively (see above). A good default is 'all [noinhib]'. # Each option and plot type should be separated by a colon, e.g: # python plot.py detail 4: activity 1-15: info noinhib ############## DEFINE ROUTINES ############### # loadspks(filepath) # takes name of file to be read and calls grvec read procedure to load file def loadspks(filepath): filename = filepath + "/spks" if os.path.exists(filename): h.grv_.read_vfile(filename) # open file using grvec print "There are %d segments in this file" % numsegs() else: print "ERROR: No spks file found at path %s!" % filename sys.exit() # numsegs() # returns number of data segments in the file def numsegs(): return int(h.grv_.llist.count()) # loadvars(filepath) # Reads variables saved in the vars file into memory def loadvars(filepath): filename = filepath + "/vars" if os.path.exists(filename): # Read the file into 'auxvars' import imp f = open(filename) auxvars = imp.load_source('auxvars', '', f) # Read parameters (e.g. numcells) from aux file f.close() else: print "ERROR: No vars file found at path %s!\n" % filename sys.exit() return auxvars # getauxdata(filepath, cellID) # Returns a list of all aux data items for a given cell ID, read via random access from the given aux file # Note: does not sanity-check that requested cell ID < numcells! # Note: all cells must be present and listed in order at each save step, with one 't = ...' line separator # e.g. cell 0 must always be on line 1, numcells + 1 + 1, numcells * 2 + 1 + 1 def getauxdata(filepath, cell): global detail # Declare that we want to use the global 'detail' value filename = filepath + "/aux" datalist = [] # Create empty list if os.path.exists(filename): lineno = 2 + cell # Start at line 2 (after "t = ..." line) curline = "-1" while curline != "": # linecache.getline returns "" on read errors, e.g. eof curline = linecache.getline(filename, lineno) lineno += (auxvars.numcells + 1) * detail # +1 accounts for the "t = ..." lines between segments splitline = curline.split(",") # Check that this is correct cell before adding to list; exit otherwise! if len(splitline) > 1: if int(splitline[CELLID]) == cell: datalist.append(splitline) else: print "ERROR: Expected cell %d, found cell %s at line %d!" % (cell, splitline[CELLID], lineno) sys.exit() else: print "ERROR: No aux file found at path %s!\n" % filename sys.exit() return datalist; # getIcellIDs(filepath) # Returns a list of all cell IDs for I-cells (so they can be explicitly plotted / ignored) def getIcellIDs(filepath): cell_list = [] for i in range(auxvars.numcells): celldata = getauxdata(filepath, i) # Read data celltype = int(celldata[0][TYPE]) if (h.strm(h.CTYP.o(celltype).s,"^I")): # Cell is I type cell_list.append(int(celldata[0][CELLID])) return cell_list # getEcellIDs(filepath) # Returns a list of all cell IDs for E-cells (so they can be explicitly plotted / ignored) def getEcellIDs(filepath): cell_list = [] for i in range(auxvars.numcells): celldata = getauxdata(filepath, i) # Read data celltype = int(celldata[0][TYPE]) if (h.strm(h.CTYP.o(celltype).s,"^I")) == False: # Cell is E type cell_list.append(int(celldata[0][CELLID])) return cell_list # read(index) # takes index of file segment to be read (1 to numsegs) # returns vec and tvec spiking data in a 2-D array # data[0] is the spike time vector # data[1] is the vector of cell IDs which spike at times in data[0] def read(index): print "Reading grvec data" data = np.zeros( (1, 1) ) # Preassign zero-length data variable # Check we're not requesting an index beyond the extent of the file if index > numsegs(): print "Index must be <= %d" % numsegs() else: h.grv_.rv_readvec(index, h.tvec, h.vec) # Read segment 'index' into vec/tvec # Convert h.tvec and h.vec into Python / numpy arrays. data = np.zeros( (2, h.vec.size()) ) # Preallocate for speed for i in range(int(h.vec.size())): data[0,i] = h.tvec.x[i] data[1,i] = h.vec.x[i] #print data.shape return data # Depending on the length of the time series of scale and activity graphs, # return an x-axis delimited in seconds, hours, or days. # Returns two values: a list with the x axis points, and a string for the units def makexaxis(start, numpoints, interval, detail): # Default to s units = 0.001 unitstring = "seconds" # If > 0.8e7 ms (8000s), use 'hours' if start + numpoints * interval * detail > 0.8e7: units = 0.001 / 60 / 60 unitstring = "hours" # If > 0.8e8 ms (80 000s, 22 hrs), use 'days' if start + numpoints * interval * detail > 0.8e8: units = 0.001 / 60 / 60 / 24 unitstring = "days" return np.arange(float(start*units), float((start + numpoints * interval * detail)*units), float((interval * detail)*units)), unitstring ############# DEFINE DIFFERENT PLOT TYPES ################# # raster # basic spike raster (cell index vs time) def raster(filepath, rasterfrom, rasterto): global titleson # Declare that we want to use the global 'titleson' variable rasterfrom = float(rasterfrom) rasterto = float(rasterto) if rasterto == rasterfrom: # By default (when rasterto = 0), draw a separate plot for each data segment. # If rasterfrom != rasterto, draw only a single plot for the data between rasterfrom and rasterto # Create save directory if not already present if (not os.path.isdir("%s/rasters" % filepath)): os.mkdir("%s/rasters" % filepath) for i in range(numsegs()): print i data = read(i) # Read data for segment i pyplot.plot(data[0], data[1], '.', markersize=1, color='black') pyplot.xlabel("Time (ms)") pyplot.ylabel("Cell ID") if titleson: pyplot.suptitle(filepath) pyplot.savefig("%s/rasters/%d.png" % (filepath, i), dpi=300, format="png") # Save at 300 DPI as 'filepath/raster/segment.png' pyplot.hold(False) # Don't add next segment's data to current figure pyplot.clf() # Clear figure for next plot else: segmentstart = int(rasterfrom / auxvars.buffertime) segmentend = int(rasterto / auxvars.buffertime) print "%d - %d" % (segmentstart, segmentend) plotx = [] ploty = [] for i in range (segmentstart, segmentend + 1): print "Reading segment %d" % i data = read(i) # Read data for segment i for j in range(len(data[0])): if data[0,j] >= rasterfrom and data[0,j] <= rasterto: # If data time is within requested bounds, add to plotdata plotx.append(data[0,j]) ploty.append(data[1,j]) print "Plotting..." pyplot.plot(plotx, ploty, '.', markersize=1, color='black') pyplot.xlabel("Time (ms)") pyplot.ylabel("Cell ID") if titleson: pyplot.suptitle(filepath) pyplot.savefig("%s/raster.png" % filepath, dpi=300, format="png") # Save at 300 DPI pyplot.clf() # Clear figure for next plot # info # Information contribution of each neuron in the network, per data segment # Using mlabwrap to call Crumiller et al.'s MATLAB Fourier information calculation code def info(filepath, cells, noinhib): # Spike data is saved in alternating 'unique' and 'repeat' stimulation blocks # Crumiller's code requires a N-by-M matrix (N = number of trials, M = number of cells), # in which each element is a vector containing the spike times of cell M in trial N. # This procedure takes the first n saved data segments (up to the auxvars.deletionstart time) # and packs all the unique/repeat trials' spike times into separate CSV files print "Processing spike data into CSV for information contribution calculations..." stimlength = auxvars.infotriallength # Length of 'repeat' or 'unique' stimulation block (ms) buffertime = auxvars.buffertime # Length of whole data segment (ms) ignoretime = auxvars.GetInfoDur - (auxvars.numinfotrials * stimlength * 2) # How much of the start of the data is not supposed to be part of the info trials? if cells == "": # Make a list of cells whose data we want to write out # (Either we requested all cells, or all-but-I cells) if noinhib: cells = getEcellIDs(filepath) else: cells = range(auxvars.numcells) if numsegs() * buffertime < auxvars.GetInfoDur: print "Not enough spike segments to obtain pre-deletion information (needs at least %d)" % int(auxvars.GetInfoDur / buffertime) else: uniquesfile = open("%s/info-uniques.csv" % filepath, 'w') repeatsfile = open("%s/info-repeats.csv" % filepath, 'w') for i in range(int( ignoretime/buffertime ), int( auxvars.GetInfoDur / buffertime)): print "Scanning segment %d (up to %d) to add to CSV file" % (i, int(auxvars.GetInfoDur/buffertime)-1) # For each data segment in the file data = read(i) # Read NEURON data for segment i for index in range(len(data[0])): spiketime = data[0, index] # What time was this cell spike? cellID = data[1, index] # Which cell just fired? if cellID in cells: trialnum = int((spiketime-ignoretime) / (stimlength * 2)) # Obtain trial number. # Force int, as in Python 2 and 3 '/' operator means one or other of integer or float div... # Find out whether we are in 1st stimlength ms (repeats), or last stimlength ms (uniques) uniquestrial = ((spiketime % float(stimlength * 2)) / float(stimlength * 2)) > 0.5 # Tests (stimlength = 4000): # spiketime = 0 => trialnum = 0, uniquestrial = false # spiketime = 3999 => trialnum = 0, uniquestrial = false # spiketime = 4000 => trialnum = 0, uniquestrial = true # spiketime = 7999 => trialnum = 0, uniquestrial = true # spiketime = 8000 => trialnum = 1, uniquestrial = false # spiketime = 11999 => trialnum = 1, uniquestrial = false # spiketime = 12000 => trialnum = 1, uniquestrial = true # spiketime = 15999 => trialnum = 1, uniquestrial = true # spiketime = 16000 => trialnum = 2, uniquestrial = false # Remove offset so that each unique/repeat block has spikes beginning at 0ms # Also div by 1000 to convert ms to s spiketime = (spiketime % stimlength) / 1000 # If spike was recorded in second half of stimlength, then write to file # (we ignore first-half spikes from each trial, so the network can safely # transition from the previous trial and settle into a stable state). if spiketime >= ((stimlength/1000) / float(2)): # Dump line to file if uniquestrial: uniquesfile.write("%d,%d,%f\n" % (cellID,trialnum,spiketime)) else: repeatsfile.write("%d,%d,%f\n" % (cellID,trialnum,spiketime)) uniquesfile.close() repeatsfile.close() # Now find out how long each cell lived for, and write to a third CSV file print "Finding cell death times..." celldeathtimes = open("%s/info-celldeaths.csv" % filepath, 'w') # Store time-of-death for each cell for cellID in cells: # For each cell in the requested list celldata = getauxdata(filepath, cellID) # Get a list of all entries for this particular cell # Find out at which point this cell died (if at all), and its scalefactor at this time for entrynum in range(len(celldata)): if int(celldata[entrynum][DEAD]) or entrynum >= len(celldata)-1: # If cell has just died, or we've reached the end of the data... break # Stop processing cellscale = float(celldata[entrynum][SCALE]) # Get scalefactor deathtime = auxvars.t_start + (entrynum * auxvars.recording_interval) # Time of death celldeathtimes.write("%d,%d,%f\n" % (cellID,deathtime,cellscale)) celldeathtimes.close() print "Done" # scale # Scalefactors over time, read from the aux file # 'cells' is a list of cells which should be plotted (or leave blank to plot all) def scale(filepath, cells, noinhib): print "Plotting scale factors..." drawlegend = False if len(cells) < 8 and len(cells) > 1: drawlegend = True # Only draw a legend when examining 2-7 specific cells plotavg = False if cells == "": plotavg = True # Make a list of cells whose data we want to write out # (Either we requested all cells, or all-but-I cells) if noinhib: cells = getEcellIDs(filepath) else: cells = range(auxvars.numcells) celldata = getauxdata(filepath, 0) # Get data for an arbitrary cell, to obtain number of entries for pre-allocating array numentries = len(celldata) yaxis = np.zeros([len(cells), numentries]) # Collect cell data for plotting average values #xaxis = range(auxvars.t_start, (auxvars.t_start + numentries * auxvars.recording_interval * detail), auxvars.recording_interval * detail) # Create x-axis xaxis, xaxisunits = makexaxis(auxvars.t_start, numentries, auxvars.recording_interval, detail) for cellnum in range(len(cells)): # For each cell in the requested list celldata = getauxdata(filepath, int(cells[cellnum])) # Get a list of all entries for this particular cell # If the cell is not dead, add the data point to the y-axis for entrynum in range(len(celldata)): if int(celldata[entrynum][DEAD]): yaxis[cellnum, entrynum] = -1 else: yaxis[cellnum, entrynum] = float(celldata[entrynum][SCALE]) if plotavg: yavg = np.zeros(numentries) stdy = np.zeros(numentries) for i in range(numentries): yavg[i] = np.mean(yaxis[ yaxis[:,i]!=-1, i ] ) # Get avg scale for each non-dead time point stdy[i] = np.std(yaxis[ yaxis[:,i]!=-1, i ] ) # Get std dev for each non-dead time point pyplot.errorbar(xaxis, yavg, yerr=stdy, ecolor='grey', linestyle='-', marker='.', markersize=1.0) else: pyplot.plot(xaxis, yaxis) # Plot each cell individually # Set axes limits to appropriate values pyplot.ylim(ymin=0) pyplot.xlim(xmax=xaxis[len(xaxis)-1], xmin=0) # Cut graph off at end of data, rather than leaving a margin until the next 'significant point' # Draw labels pyplot.xlabel("Time (%s)" % xaxisunits) pyplot.ylabel("Scale factor") if titleson: pyplot.suptitle(filepath) # Draw legend if only a small number of cells if drawlegend: pyplot.legend(cells, loc='upper right') # Draw legend # Save at 300 DPI as 'filepath/scale.pdf' pyplot.savefig("%s/scale.pdf" % filepath, dpi=300, format="pdf") # Save processed plot data to a file for later analysis np.savez("%s/scale" % filepath, x=xaxis, y=yavg, err=stdy) pyplot.clf() # Clear figure for next plot print "Done" # activity # Activity values over time, read from the aux file # 'cells' is a list of cells which should be plotted (or leave blank to plot all) def activity(filepath, cells, noinhib): global titleson # Declare that we want to use the global 'titleson' variable print "Plotting activity values..." drawlegend = False if len(cells) < 8 and len(cells) > 1: drawlegend = True # Only draw a legend when examining 2-7 specific cells plotavg = False if cells == "": plotavg = True # Make a list of cells whose data we want to write out # (Either we requested all cells, or all-but-I cells) if noinhib: cells = getEcellIDs(filepath) else: cells = range(auxvars.numcells) celldata = getauxdata(filepath, 0) # Get data for an arbitrary cell, to obtain number of entries for pre-allocating array numentries = len(celldata) yaxis = np.zeros([len(cells), numentries]) # Collect cell data for plotting average values #xaxis = range(auxvars.t_start, (auxvars.t_start + numentries * auxvars.recording_interval * detail), auxvars.recording_interval * detail) # Create x-axis xaxis, xaxisunits = makexaxis(auxvars.t_start, numentries, auxvars.recording_interval, detail) for cellnum in range(len(cells)): # For each cell in the requested list celldata = getauxdata(filepath, int(cells[cellnum])) # Get a list of all entries for this particular cell # If the cell is not dead, add the data point to the y-axis for entrynum in range(len(celldata)): if int(celldata[entrynum][DEAD]): yaxis[cellnum, entrynum] = -1 else: yaxis[cellnum, entrynum] = float(celldata[entrynum][ACTIVITY]) * 1000 if plotavg: yavg = np.zeros(numentries) stdy = np.zeros(numentries) for i in range(numentries): yavg[i] = np.mean(yaxis[ yaxis[:,i]!=-1, i ] ) # Get avg activity for each non-dead time point stdy[i] = np.std(yaxis[ yaxis[:,i]!=-1, i ] ) # Get std dev for each non-dead time point pyplot.errorbar(xaxis, yavg, yerr=stdy, ecolor='grey', linestyle='-', marker='.', markersize=1.0) else: pyplot.plot(xaxis, yaxis) # Plot each cell individually # Set axes limits to appropriate values pyplot.ylim(ymin=0) pyplot.xlim(xmax=xaxis[len(xaxis)-1], xmin=0) # Cut graph off at end of data, rather than leaving a margin until the next 'significant point' # Draw labels pyplot.xlabel("Time (%s)" % xaxisunits) pyplot.ylabel("Activity (Hz)") if titleson: pyplot.suptitle(filepath) # Draw legend if only a small number of cells if drawlegend: pyplot.legend(cells, loc='upper right') # Draw legend # Save at 300 DPI as 'filepath/activity.pdf' pyplot.savefig("%s/activity.pdf" % filepath, dpi=300, format="pdf") # Save processed plot data to a file for later analysis np.savez("%s/activity" % filepath, x=xaxis, y=yavg, err=stdy) pyplot.clf() # Clear figure for next plot print "Done" # deletionscale # Plots relative scaling factors of cells (compared to mean) at time-of-deletion, vs total time def deletionscale(filepath, noinhib): global titleson # Declare that we want to use the global 'titleson' variable print "Plotting deletion scale values..." if noinhib: cells = getEcellIDs(filepath) else: cells = range(auxvars.numcells) celldata = getauxdata(filepath, 0) # Get data for an arbitrary cell, to obtain number of entries for pre-allocating array numentries = len(celldata) yaxis = [] # Individual cell scale values at time x ymeanscale = [] # Mean scale values over all cells at time x # Get xaxisunits # But throw away the actual xaxis vector, as we don't want fixed-size intervals xaxis, xaxisunits = makexaxis(auxvars.t_start, numentries, auxvars.recording_interval, detail) xaxis = [] # Get mean scale of all the population at each entry in the x-axis print "Finding mean scale values..." # Go through all data once for every entrynum and calculate mean population scale at that point allscales = np.zeros([len(cells), numentries]) # Collect cell data for plotting average values for cellnum in range(len(cells)): # For each cell in the requested list celldata = getauxdata(filepath, int(cells[cellnum])) # Get cell data # If the cell is not dead, add the data point to the y-axis for entrynum in range(len(celldata)): if int(celldata[entrynum][DEAD]): allscales[cellnum, entrynum] = -1 else: allscales[cellnum, entrynum] = float(celldata[entrynum][SCALE]) meanscales = np.zeros(numentries) for i in range(numentries): meanscales[i] = np.mean(allscales[ allscales[:,i]!=-1, i ] ) # Get avg scale for each time point # Get individual scale factor for each cell at time of deletion print "Finding individual scale values..." for cellnum in range(len(cells)): # For each cell in the requested list celldata = getauxdata(filepath, int(cells[cellnum])) # Get a list of all entries for this particular cell # Find out at which point this cell died (if at all), and its scalefactor at this time for entrynum in range(len(celldata)): if int(celldata[entrynum][DEAD]) or entrynum >= len(celldata)-1: # If cell has just died, or we've reached the end of the data... cellscale = float(celldata[entrynum][SCALE]) # Get scalefactor yaxis.append(cellscale) # Add scalefactor to y-axis ymeanscale.append(meanscales[entrynum]) # Obtain mean scalefactor and add to y xaxis.append(auxvars.t_start + (entrynum * auxvars.recording_interval)) # Add time to x break # Stop processing # Data is in arbitrary order (x-axis time depends on cell ID, rather than x-axis incrementing in time) # This makes it impossible to draw a line from 'start' to 'finish'. So sort the x,y,mean(y) data first. sortedmeanscales = sorted(zip(xaxis, ymeanscale)) sortedscales = sorted(zip(xaxis, yaxis)) # Now unzip xaxis = [point[0] for point in sortedscales] yaxis = [point[1] for point in sortedscales] ymeanscale = [point[1] for point in sortedmeanscales] # Draw plots pyplot.plot(xaxis, ymeanscale, '', linestyle='-', color='r') # Plot line for mean scale pyplot.plot(xaxis, yaxis, '.', linestyle='', color='b') # Plot each cell's individual scale # Set axes limits to appropriate values pyplot.ylim(ymin=0) pyplot.xlim(xmax=xaxis[len(xaxis)-1], xmin=0) # Cut graph off at end of data, rather than leaving a margin until the next 'significant point' # Draw labels pyplot.xlabel("Time (%s)" % xaxisunits) pyplot.ylabel("Scale factor") if titleson: pyplot.suptitle(filepath) # Save at 300 DPI as 'filepath/deletionscale.pdf' pyplot.savefig("%s/deletionscale.pdf" % filepath, dpi=300, format="pdf") # Save processed plot data to a file for later analysis np.savez("%s/deletionscale" % filepath, x=xaxis, y=yaxis) pyplot.clf() # Clear figure for next plot print "Done" # getmeanpopscale # Helper function for deletionscale # Returns mean scalefactor of the currently-alive population at a given time def getmeanpopscale(cells,entrynum): totalscale = 0 livingcells = 0 for cellnum in range(len(cells)): celldata = getauxdata(filepath, int(cells[cellnum])) # Obtain each cell's data if not int(celldata[entrynum][DEAD]): # If we have found a cell which is still alive livingcells = livingcells + 1 # Increase counter of living cells totalscale = totalscale + float(celldata[entrynum][SCALE]) # Add to scalefactor sum return totalscale / livingcells # power # Plots power spectra of the network separately for E and I cells, obtained from the neural spike trains def power(filepath): binsize = 5.0 # ms for bin size (5ms = 200Hz sampling rate) Icells = getIcellIDs(filepath) # Make list of I cells # Load each segment of saved spike data for i in range(numsegs()): # For each data segment in the file data = read(i) # Read NEURON data for segment i # Create empty numpy array for number of spikes-per-bin firstspk = data[0,0] lastspk = data[0,len(data[1])-1] timecovered = lastspk-firstspk numbins = int(timecovered / binsize) spikesperbinI = np.zeros( numbins ) spikesperbinE = np.zeros( numbins ) # Make MUA time series vector by counting all population spikes during each 'binsize' ms print "Binning file segment %d into %d bins (window size = %f ms)" % (i, numbins, binsize) for j in range(len(data[1])): spiketime = data[0,j] bin = int(float((spiketime % (numbins * binsize)) / binsize)) # Find bin number into which this spike should go cellID = data[1,j] if cellID in Icells: # Check if this is an I or E cell spikesperbinI[bin] = spikesperbinI[bin] + 1 else: spikesperbinE[bin] = spikesperbinE[bin] + 1 # Find and subtract mean from the MUA vector meanmuaI = float(np.sum(spikesperbinI) / len(spikesperbinI)) spikesperbinI = spikesperbinI - meanmuaI meanmuaE = float(np.sum(spikesperbinE) / len(spikesperbinE)) spikesperbinE = spikesperbinE - meanmuaE # Run MTSpec's spectral analysis on the MUA time series # MTSpec documentation: http://krischer.github.com/mtspec/mtspec.multitaper.mtspec.html#mtspec.multitaper.mtspec delta = 1.0 / float(1000.0/binsize) # 1 / sample rate (Hz) time_bandwidth = 4 nfft = None number_of_tapers = None quadratic = False adaptive = True verbose = False optional_output = False statistics = False rshape = False fcrit = False print "Running mtspec for I cells" #spec, freq = multitaper.mtspec(spikesperbin, delta, time_bandwidth, nfft, quadratic, adaptive, verbose, optional_output, statistics, rshape, fcrit) specI, freqI = mtspec(spikesperbinI, delta, time_bandwidth) print "Running mtspec for E cells" specE, freqE = mtspec(spikesperbinE, delta, time_bandwidth) # Save to file for further examination #numpy.savetxt("%s/power-%d-muaE.txt" % (filepath, i), spikesperbinE) #numpy.savetxt("%s/power-%d-specE.txt" % (filepath, i), specE) #numpy.savetxt("%s/power-%d-freqE.txt" % (filepath, i), freqE) print "Generating plots\n" # Create save directory if not already present if (not os.path.isdir("%s/power" % filepath)): os.mkdir("%s/power" % filepath) # I power # Draw labels pyplot.xlabel("Frequency (Hz)") pyplot.ylabel("Normalised power") if titleson: pyplot.suptitle(filepath) pyplot.plot(freqI,specI) # Set axes limits to appropriate values pyplot.ylim(ymax=0.08, ymin=0) pyplot.xlim(xmax=100, xmin=0) pyplot.savefig("%s/power/I-%d.pdf" % (filepath, i), dpi=300, format="pdf") pyplot.clf() # E power # Draw labels pyplot.xlabel("Frequency (Hz)") pyplot.ylabel("Normalised power") if titleson: pyplot.suptitle(filepath) pyplot.plot(freqE,specE) # Set axes limits to appropriate values pyplot.ylim(ymax=0.08, ymin=0) pyplot.xlim(xmax=100, xmin=0) pyplot.savefig("%s/power/E-%d.pdf" % (filepath, i), dpi=300, format="pdf") pyplot.clf() # Save processed plot data to a file for later analysis np.savez("%s/power/I-%d.npz" % (filepath, i), x=freqI, y=specI) np.savez("%s/power/E-%d.npz" % (filepath, i), x=freqE, y=specE) print "Done" def nte(filepath,segs): # Plot normalised transfer entropy (nTE) # Uses a time-binned MUA vector over two given segments of data, for comparison # Input requires IDs of two data segments # Output as a pair of line plots comparing nTEs per population, for the two segments h.usetable_infot = 0 # Turn off log lookup tables binsize = 10.0 # ms for bin size netscale = round(auxvars.numcells / 470) poplabels = ['E6' 'I6' 'I6L' 'E5B' 'E5R' 'I5' 'I5L' 'E4' 'I4' 'I4L' 'E2' 'I2' 'I2L'] popsizes = np.array([59, 25, 13, 17, 65, 25, 13, 30, 20, 14, 150, 25, 13]) # First element was 60, but cell 0 always seems to be missing popsizes *= netscale # Scale-up if we have > 470 cells # Load each segment of spike data -- each segment will be plotted as a separate line for istring in segs: i = int(istring) # Convert char to int print "Reading segment %d" % i # For each data segment in the list data = read(i) # Read NEURON data for segment i # Create empty numpy array for number of spikes-per-bin firstspk = data[0,0] lastspk = data[0,len(data[1])-1] timecovered = lastspk-firstspk numbins = int(timecovered / binsize) spikesperbin = np.zeros( [auxvars.numcells, numbins] ) # x,y matrix (== cellid, spikes per bin) # Make MUA time series vector by counting all population spikes during each 'binsize' ms print "Binning file segment %d into %d bins (window size = %f ms)" % (i, numbins, binsize) for j in range(len(data[1])): spiketime = data[0,j] cellid = data[1,j] bin = int(float((spiketime % (numbins * binsize)) / binsize)) # Find bin number into which this spike should go spikesperbin[cellid,bin] += 1 # Combine cells from each population together allpopMUAs = [] currentpopcell = 0 for popnum in range(len(popsizes)): finalpopcell = currentpopcell + popsizes[popnum] - 1 popMUA = np.sum(spikesperbin[currentpopcell:finalpopcell, :], axis=0) currentpopcell += popsizes[popnum] allpopMUAs.append(popMUA) # append all popMUAs into one structure print allpopMUAs print "\n" normte = np.zeros(len(popsizes)) # For each population, obtain its MUA against all other populations for popnum in range(len(allpopMUAs)): allotherpops = copy.deepcopy(allpopMUAs) thispop = allotherpops.pop(popnum) # Gives us this population, AND removes it from 'remaining' list thispopvec = h.Vector(thispop) # Convert thispop to a hoc Vector # Now, for each 'other' population, compare it to 'thispop' for otherpop in allotherpops: # Convert otherpop to a hoc Vector otherpopvec = h.Vector(otherpop) normtevec = h.normte(thispopvec, otherpopvec, 30) # Try 30 shuffles #normtevec.printf("%8.4f\n") #print "using %f, popnum %d" % (normtevec.x[2], popnum) normte[popnum] += normtevec.x[2] print normte print "\n" pyplot.plot(range(len(popsizes)), normte, linestyle='-', marker='x') # Draw labels xlocations=range(len(popsizes)) # Allow plot to be spaced equally on x-axis, independent of the value pyplot.xticks(xlocations,poplabels) # Display frequency/weight values over the x tick locations pyplot.xlabel("Population") pyplot.ylabel("nTE") if titleson: pyplot.suptitle(filepath) # Save at 300 DPI as 'filepath/deletionscale.pdf' pyplot.savefig("%s/nte.pdf" % filepath, dpi=300, format="pdf") # Save processed plot data to a file for later analysis #np.savez("%s/deletionscale" % filepath, x=xaxis, y=yaxis) pyplot.clf() # Clear figure for next plot print "Done" ################# COMMAND LINE PROCESSING ################# def processargs(args): # Checks list of arguments, and returns a list of cells from either a comma-separated list or # a hyphenated range (e.g. 0,1,3,7,9 or 0-10) if len(args.split()) > 3: print "\n\nERROR: Too many arguments! Permitted: 'noinhib', 'cell1,cell2,..,celln', or '0-100'" if len(args.split()) == 1 or (len(args.split()) == 2 and 'noinhib' in args): # No arguments (or just 'noinhib') supplied return "" # Return an empty list else: # A cell / cell range / time range argument was supplied... args = args.split() args.pop(0) # Remove command name from the front of the list of args if "-" in args[0]: # We're dealing with a request for a range of cells args = args[0].split("-") # Split dash-separated range cellslist = [] for i in range(int(args[0]), int(args[1])+1): cellslist.append(i) return cellslist else: args = args[0].split(",") # Split comma-separated list return args def makegraph(graphtype): if "noinhib" in graphtype: noinhib = True else: noinhib = False command = graphtype.split()[0] if command == "quit": sys.exit() elif command == "raster": if len(graphtype.split()) == 2: rasterrange = graphtype.split()[1] rasterrange = rasterrange.split('-') rasterfrom = rasterrange[0] rasterto = rasterrange[1] else: rasterfrom = 0 rasterto = 0 raster(filepath, rasterfrom, rasterto) elif command == "info": info(filepath, processargs(graphtype), noinhib) elif command == "scale": scale(filepath, processargs(graphtype), noinhib) elif command == "activity": activity(filepath, processargs(graphtype), noinhib) elif command == "power": power(filepath) elif command == "deletionscale": deletionscale(filepath, noinhib) elif command == "all": scale(filepath, "", noinhib) activity(filepath, "", noinhib) deletionscale(filepath, noinhib) raster(filepath, 0, 0) info(filepath, "", noinhib) power(filepath) elif command == "titles": global titleson titleson = not titleson print "\n\n\nTitles on: %s" % titleson elif command == "detail": global detail detail = int(graphtype.split()[1]) print "\nSetting detail level to %d" % detail elif command == "nte": nte(filepath, processargs(graphtype)) else: print "\nUnknown graph type or option '%s'" % command ############### MAIN METHOD ################ # Imports import sys import readline import numpy as np import string import os import linecache import matplotlib matplotlib.use('agg') # Prevent pyplot from trying to open a graphical front-end on headless clients from matplotlib import pyplot from mtspec import * # Multitaper spectral analysis (MTSpec library from http://pypi.python.org/pypi/mtspec) import copy # Check filename was supplied if len(sys.argv) < 2: print "Usage:\npython plot.py [option arg:] [plot_type args:] [plot_type args:] [...:]" print "For a list of plot types and arguments, run interactively (i.e. python plot.py )" sys.exit() if len(sys.argv) >= 3: plot_as_arg = True args = sys.argv[2:] # Get all plot type args graphtypelist = "" for item in args: graphtypelist += item + " " # Concatenate the args string back together (including a space between each item) else: plot_as_arg = False # Handle NEURON imports print "Loading NEURON/grvec routines... \n\n" from neuron import h # Load grvec / intfzip simulation h('xopen("setup.hoc")') h('xopen("nrnoc.hoc")') h('load_file("init.hoc")') h('load_file("grvec.hoc")') h('load_file("labels.hoc")') h('load_file("infot.hoc")') # Load file global auxvars # Vars from the vars file are held here, to be accessible by all methods # Set index numbers for the aux file data lines (format: CELL,TYPE,ACTIVITY(firing rate),POSTSYNACTIVITY,TARGET,SCALE,DEAD) global CELLID CELLID = 0 global TYPE TYPE = 1 global FIRINGRATE FIRINGRATE = 3 global ACTIVITY ACTIVITY = 2 global TARGET TARGET = 4 global SCALE SCALE = 5 global DEAD DEAD = 6 global titleson titleson = False global detail detail = 1 # Take every nth point for plotting (set == 1 to plot every point, or > 1 to reduce points) filepath = sys.argv[1] print "\nLoading data from %s" % filepath loadspks(filepath) auxvars = loadvars(filepath) print "Done." if plot_as_arg: for graphtype in graphtypelist.split(":"): # Make all the requested graph types, then quit makegraph(graphtype) sys.exit() else: while 1: # Get user's graph preference interactively print "\n\n\n" print filepath print "titles toggle graph titles (use 'off' for printing; default off)" print "detail set level of detail for plots (every 'n'th point; default=1)" print "\n" print "all [noinhib] draw all plots [excluding inhibitory cells]" print "raster [starttime-endtime] basic spike raster [between given ms times; default=ALL]" print "info [cells | noinhib] information contribution [comma-separated list of individual cells or range of cells e.g. 0-20 | without inhibitory cells]" print "scale [cell1s | noinhib] scale factor values [comma-separated list of individual cells or range of cells e.g. 0-20 | without inhibitory cells]" print "activity [cells | noinhib] activity values [comma-separated list of individual cells or range of cells e.g. 0-20 | without inhibitory cells]" print "deletionscale scalefactors of cells at time-of-deletion, plotted against total time" print "power oscillatory power spectra" print "quit exit program" graphtype = raw_input("\ne.g. \"raster 100000-150000\", \"scale 5,23,190\", \"info 0-49\", or \"activity no-inhib\"\n>") makegraph(graphtype)