# plotdeletions.py # Mark Rowan, School of Computer Science, University of Birmingham, UK # May 2013 # For a given directory containing experiments, each containing multiple runs, # obtain the proportion of deleted cells (I or E) and plot an error graph # E.g. for experiment 'ADproswt', containing directories 'proswt0.5', 'proswt2', # each of which contains runs 1-20: find the proportion of dead E/I cells from 'aux' file # and plot against the neurostimulation weight value, obtained from 'vars' file. # Usage: python plotdeletions.py # ---------------------------------------------------------------------- # 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! Ignoring this run.\n" % filename #sys.exit() auxvars = 0 return auxvars # readlastauxdata(filepath) # Reads the last "t = ..." segment from a given aux file and returns it as a list of lists def readlastauxdata(filepath): filename = filepath + "/aux" datalist = [] # Create empty list if os.path.exists(filename): # Search to end of file to find line number (nasty hack, but easy to code) with open(filename) as f: for numlines, l in enumerate(f): pass numlines +=2 # Account for EOF (1 line) and "t = ..." (1 line) #print "Number of lines in aux file = %d" % numlines curlineno = numlines - auxvars.numcells # Number of lines to backup while curlineno < numlines: curline = linecache.getline(filename, curlineno) splitline = curline.split(",") # Read this line from the file datalist.append(splitline) # datalist contains cells 0->numcells curlineno += 1 else: print "ERROR: No aux file found at path %s!\n" % filename sys.exit() #print datalist return datalist # There is no sanity-checking here. Ideally we should at least ensure that each # line we read from the file is actually a data-line with the correct format, # rather than the EOF or "t = ..." markers (in case of a truncated file, or similar) # getIcellIDs(auxdata) # Returns a list of all cell IDs for I-cells (so they can be explicitly plotted / ignored) def getIcellIDs(auxdata): cell_list = [] for i in range(auxvars.numcells): celltype = int(auxdata[i][TYPE]) if (h.strm(h.CTYP.o(celltype).s,"^I")): # Cell is I type cell_list.append(int(auxdata[i][CELLID])) return cell_list # getEcellIDs(auxdata) # Returns a list of all cell IDs for E-cells (so they can be explicitly plotted / ignored) def getEcellIDs(auxdata): cell_list = [] for i in range(auxvars.numcells): celltype = int(auxdata[i][TYPE]) if (h.strm(h.CTYP.o(celltype).s,"^I")) == False: # Cell is E type cell_list.append(int(auxdata[i][CELLID])) return cell_list # ---------------------------------------------------------------------- # 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 # Check filename was supplied if len(sys.argv) < 2: print "Usage:\npython plotdeletions.py " sys.exit() # 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")') # 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 # MAIN LOOP filepath = sys.argv[1] # Describes the top-level experiment directory # Determine whether we are plotting frequency or weight plotfreq = 0 plotwt = 0 if "freq" in filepath: plotfreq = 1 else: if "wt" in filepath: plotwt = 1 # Otherwise, we are plotting stimulation start time # Get list of sub-directories for each experimental parameter value dirlist = [o for o in os.listdir(filepath) if os.path.isdir(os.path.join(filepath,o))] print dirlist if "baseline" in dirlist: dirlist.remove("baseline") else: print "'baseline' dir not found. Please add a symlink to a baseline comparison dir'" sys.exit # Pre-allocate global graph x/y axes yaxisI = np.zeros([2,len(dirlist)]) # Dimension 1 is mean, dimension 2 is std yaxisE = np.zeros([2,len(dirlist)]) # Dimension 1 is mean, dimension 2 is std xaxis = np.zeros([1,len(dirlist)]) # For each experiment's sub-directory (i.e. wt/freq value on the x-axis) for exptdir in dirlist: exptdirpath = "%s/%s" % (filepath, exptdir) rundirlist = [o for o in os.listdir(exptdirpath) if os.path.isdir(os.path.join(exptdirpath,o))] print rundirlist exptnum = dirlist.index(exptdir) # Find indexof(exptdir) # Pre-allocate y axis arrays for *this run* thisyaxisI = np.array([]) thisyaxisE = np.array([]) # For each experimental run (seeds 1-20), to get error bars for rundir in rundirlist: print "\nrundir %s" % rundir rundirpath = "%s/%s" % (exptdirpath, rundir) print "rundirpath %s" % rundirpath runnum = rundirlist.index(rundir) # Find indexof(rundir) print "Loading data for run %s of %s from dir '%s'" % (runnum, len(rundirlist), rundir) auxvars = loadvars(rundirpath) # Get vars (for wt/freq) if auxvars: # Only if the 'vars' file was successfully loaded... # Obtain the wt/freq/time value from 'vars' -> xaxis if plotfreq: xaxis[0,exptnum]=auxvars.prosfreq else: if plotwt: xaxis[0,exptnum]=auxvars.proswt else: xaxis[0,exptnum]=float(auxvars.prosthesisstart)/1000/60/60 # Read auxdata auxdata = readlastauxdata(rundirpath) # Get list of cells listofIcells = getIcellIDs(auxdata) listofEcells = getEcellIDs(auxdata) # For each I cell deadIproportion = 0 for Icell in listofIcells: # Obtain the 'dead' flag value from 'aux' at end of run if int(auxdata[Icell][DEAD]): # If cell is dead, add to proportion of dead cells for this run deadIproportion += 1 # float(1.0/len(listofIcells)) # For each E cell deadEproportion = 0 for Ecell in listofEcells: # Obtain the 'dead' flag value from 'aux' at end of run if int(auxdata[Ecell][DEAD]): # If cell is dead, add to proportion of dead cells for this run deadEproportion += 1 # float(1.0/len(listofEcells)) # Compare to baseline for this particular run # Find baseline level of deletion for this particular run BLrundirpath = "%s/baseline/%s" % (filepath, rundir) print "Loading baseline data from %s" % BLrundirpath # Read auxdata BLauxdata = readlastauxdata(BLrundirpath) # For each I cell BLdeadIproportion = 0 for Icell in listofIcells: # Obtain the 'dead' flag value from 'aux' at end of run if int(BLauxdata[Icell][DEAD]): # If cell is dead, add to proportion of dead cells for this run BLdeadIproportion += 1 # float(1.0/len(listofIcells)) # For each E cell BLdeadEproportion = 0 for Ecell in listofEcells: # Obtain the 'dead' flag value from 'aux' at end of run if int(BLauxdata[Ecell][DEAD]): # If cell is dead, add to proportion of dead cells for this run BLdeadEproportion += 1 # float(1.0/len(listofEcells)) thisyaxisI = np.hstack([thisyaxisI,deadIproportion - BLdeadIproportion]) thisyaxisE = np.hstack([thisyaxisE,deadEproportion - BLdeadEproportion]) print thisyaxisE print thisyaxisI # Stack this run's y axis values onto the global graph axes yaxisI[0,exptnum] = np.mean(thisyaxisI) yaxisI[1,exptnum] = np.std(thisyaxisI) yaxisE[0,exptnum] = np.mean(thisyaxisE) yaxisE[1,exptnum] = np.std(thisyaxisE) print xaxis print yaxisI print yaxisE # Sort data xsorted = xaxis[0][np.argsort(xaxis[0])] yisorted = yaxisI[0][np.argsort(xaxis[0])] yierrsorted = yaxisI[1][np.argsort(xaxis[0])] yesorted = yaxisE[0][np.argsort(xaxis[0])] yeerrsorted = yaxisE[1][np.argsort(xaxis[0])] # Save processed plot data to a file for later analysis #np.savez("%s/avgdeadnum.npz" % filepath, x=xsorted, yi=yisorted, ye=yesorted, yierr=yierrsorted, yeerr=yeerrsorted) xlocations=np.arange(0,len(xsorted)) # Allow plot to be spaced equally on x-axis, independent of the value pyplot.axhline(y=0, linewidth=2, color='0.5', linestyle='dashed') # Draw dashed line at y=0 pyplot.errorbar(xlocations, yisorted, yerr=yierrsorted, ecolor='b', elinewidth=2, linestyle='None', marker='x', markersize=10, markeredgewidth=3, markeredgecolor='b') pyplot.errorbar(xlocations, yesorted, yerr=yeerrsorted, ecolor='r', elinewidth=2, linestyle='None', marker='x', markersize=10, markeredgewidth=3, markeredgecolor='r') pyplot.xticks(xlocations,xsorted) # Display frequency/weight values over the x tick locations x1,x2,y1,y2 = pyplot.axis() pyplot.xlim(xmin = x1-1, xmax = x2+1) # Give a bit of space either side of the data #pyplot.ylabel('Change in proportion of dead cells') pyplot.ylabel('Change in # of dead cells') if plotfreq: pyplot.xlabel('Stimulation frequency (Hz)') else: if plotwt: pyplot.xlabel('Stimulation weight') else: pyplot.xlabel('Stimulation start time (hours)') matplotlib.rcParams.update({'font.size': 16}) pyplot.savefig("%s/avgdeadnum.pdf" % filepath, dpi=300, format="pdf")