Find peak height for sequentially numbered potentiostat data files
—————————–
#!/usr/bin/env python """This program finds peak height values (i.e., peak currents) from .txt files containing squarewave voltammogram data (generated by CH Instruments Potentiostats). """ __author__ = "Andrew J. Bonham" __copyright__ = "Copyright 2010" __credits__ = ["Andrew J. Bonham", "Ryan White"] __version__ = 1.06 __maintainer__ = "Andrew J. Bonham" __email__ = "bonham@gmail.com" __status__ = "Production" # Setup: import basic modules we need import csv import os import time import platform import glob import pickle #Setup: import Tkinter for GUI import Tkinter try: #use ttk for native widgets if possible import ttk except: import Tkinter as ttk import tkFileDialog #Setup: import modules for curve fitting and graphing import numpy from scipy.optimize import leastsq import pylab #Define our classes class PeakFinderApp(Tkinter.Tk): '''This is the gui for Peak Finder. This application displays a user interface to select a directory and to specify file formats and output filename, then reports on the program function with a simple progressbar.''' def __init__(self): '''set up the peak finder app GUI object.''' self.directory_manager() #set the default directory to user's home directory Tkinter.Tk.__init__(self) #initialize the app as a GUI object self.window_title = "Peak Finder " + str(__version__) logic = PeakLogic(self) #invoke PeakLogic object to do the actual work self.title(self.window_title) self.geometry('+100+100') #Create the frame for the main window mainframe = ttk.Frame(self) #mainframe.grid_propagate(0) mainframe.grid(column=0, row=0, sticky=(Tkinter.N, Tkinter.W, Tkinter.E, Tkinter.S), padx=5, pady=5) mainframe.columnconfigure(0, weight=1) mainframe.rowconfigure(0, weight=1) self.resizable(Tkinter.FALSE,Tkinter.FALSE) #define our variables in GUI-form and give sane defaults self.filename_ = Tkinter.StringVar(value='output1') self.format_ = Tkinter.StringVar(value='run_') self.output = Tkinter.StringVar() self.longpath_ = Tkinter.StringVar(value=self.mydocs) self.dir_selected = Tkinter.IntVar(value=0) self.progress_ = Tkinter.IntVar(value=0) self.init_potential_ = Tkinter.DoubleVar(value=-0.18) self.final_potential_ = Tkinter.DoubleVar(value=-0.40) self.peak_center_ = Tkinter.DoubleVar(value=-0.322) self.final_edge_ = Tkinter.DoubleVar(value=-0.5) self.init_edge_ = Tkinter.DoubleVar(value=0) #display the entry boxes ttk.Label(mainframe, text=self.window_title, font='helvetica 12 bold').grid(column=1, row=1, sticky=Tkinter.W, columnspan=2) self.filename_entry = ttk.Entry(mainframe, width=12, textvariable=self.filename_).grid(column=2, row=2, sticky=(Tkinter.W, Tkinter.E)) ttk.Label(mainframe, text="Output Filename:").grid(column=1, row=2, sticky=Tkinter.W) self.format_entry = ttk.Entry(mainframe, width=12, textvariable=self.format_).grid(column=2, row=3, sticky=(Tkinter.W, Tkinter.E)) ttk.Label(mainframe, text="Data Filename Format:").grid(column=1, row=3, sticky=Tkinter.W) self.init_potential = ttk.Entry(mainframe, width=12, textvariable=self.init_potential_).grid(column=2, row=9, sticky=(Tkinter.W, Tkinter.E)) ttk.Label(mainframe, text="Initial Potential:").grid(column=1, row=9, sticky=Tkinter.W) self.final_potential = ttk.Entry(mainframe, width=12, textvariable=self.final_potential_).grid(column=2, row=10, sticky=(Tkinter.W, Tkinter.E)) ttk.Label(mainframe, text="Final Potential:").grid(column=1, row=10, sticky=Tkinter.W) self.peak_center = ttk.Entry(mainframe, width=12, textvariable=self.peak_center_).grid(column=2, row=11, sticky=(Tkinter.W, Tkinter.E)) ttk.Label(mainframe, text="Peak Center:").grid(column=1, row=11, sticky=Tkinter.W) ttk.Label(mainframe, text="Peak Location Guess",font='helvetica 10 bold').grid(column=1, row=7, sticky=Tkinter.W) ttk.Label(mainframe, text="(only change if necessary)").grid(column=1, row=8, sticky=Tkinter.W) #Display the buttons ttk.Button(mainframe, text="Find Peaks", command=logic.peakfinder).grid(column=1, row=14, sticky=Tkinter.W) ttk.Button(mainframe, text="Test fit", command=logic.testFit).grid(column=2, row=14, sticky=Tkinter.E) ttk.Button(mainframe, text="Data Dir.", command=self.dir_select).grid(column=1, row=6, sticky=Tkinter.W) #Set up a space to show the output ttk.Label(mainframe, textvariable=self.output).grid(column=1, row=15, sticky=(Tkinter.W, Tkinter.E), columnspan=5) #Set up Outlier correction ttk.Label(mainframe, text="Regions to include (from -> to):").grid(column=1, row=12, sticky=Tkinter.W) self.final_edge = ttk.Entry(mainframe, width=12, textvariable=self.final_edge_).grid(column=1, row=13, sticky=(Tkinter.E)) self.init_edge = ttk.Entry(mainframe, width=12, textvariable=self.init_edge_).grid(column=2, row=13, sticky=(Tkinter.W, Tkinter.E)) #Pad the widgets for prettiness for child in mainframe.winfo_children(): child.grid_configure(padx=5, pady=5) #Display a progressbar self.bar = progressBar(mainframe, width = 100, height = 10) self.bar.update(0) #ttk.Label(mainframe, text="").grid(column=1, row=16, sticky=(Tkinter.W, Tkinter.E)) #add a system menu self.option_add('*tearOff', Tkinter.FALSE) self.menubar = Tkinter.Menu(self) self['menu'] = self.menubar ##self.sysmenubar = Tkinter.Menu(self) ##self.sysmenu = Tkinter.Menu(self.sysmenubar, name='system') ##self.sysmenubar.add_cascade(menu=self.sysmenu, label='sysmenu') ##self.sysmenubar.add_command(label='Test', command=self.destroy) self.menu_file = Tkinter.Menu(self.menubar) self.menubar.add_cascade(menu=self.menu_file, label='File') self.menu_file.add_command(label="Load Params", command=self.loadParams) self.menu_file.add_command(label="Save Params", command=self.saveParams) self.menu_file.add_command(label='Exit', command=self.destroy) self.menu_help = Tkinter.Menu(self.menubar, name='help') self.menubar.add_cascade(menu=self.menu_help, label='Help') self.menu_help.add_command(label='How to Use', command=self.helpPopUp) self.menu_help.add_command(label='About Peak Finder...', command=self.aboutPopUp) self.status = StatusBar(self,0,18) self.status.set("Ready") #Run the program event loop self.mainloop() def directory_manager(self): '''PeakFinderApp.directory_manager() sets the initial directory to the users home directory.''' if platform.system() == 'Windows': #import windows file management from win32com.shell import shell, shellcon self.mydocs = shell.SHGetFolderPath(0, shellcon.CSIDL_PERSONAL, 0, 0) else: #import mac/linux file management self.mydocs = os.getenv("HOME") os.chdir(self.mydocs) def aboutPopUp(self): '''Display a pop-up menu about the program.''' self.aboutDialog = Tkinter.Toplevel(self) self.aboutDialog.resizable(Tkinter.FALSE,Tkinter.FALSE) self.aboutDialog.geometry('+400+100') #Create the frame for the about window aboutframe = ttk.Frame(self.aboutDialog, width='200', height='200') aboutframe.grid_propagate(0) aboutframe.grid(column=0, row=0, sticky=(Tkinter.N, Tkinter.W, Tkinter.E, Tkinter.S), padx=5, pady=5) aboutframe.columnconfigure(0, weight=1) aboutframe.rowconfigure(0, weight=1) ttk.Label(aboutframe, text=self.window_title, font='helvetica 12 bold', anchor="center").grid(column=0, row=0, sticky=(Tkinter.N,Tkinter.W, Tkinter.E)) ttk.Label(aboutframe, text='Voltammogram data analysis\nsoftware for CH Instruments data.\n\nWritten by\n'+str(__author__)+'\nhttp://www.andrewjbonham.com\n'+str(__copyright__)+'\n\n\n',anchor="center", justify="center").grid(column=0, row=1,sticky=Tkinter.N) ttk.Button(aboutframe, text="Close", command=self.aboutDialog.destroy).grid(column=0, row=4, sticky=(Tkinter.S)) self.aboutDialog.mainloop() def helpPopUp(self): '''Display a pop-up menu explaining how to use the program.''' self.helpDialog = Tkinter.Toplevel(self) self.helpDialog.resizable(Tkinter.FALSE,Tkinter.FALSE) self.helpDialog.geometry('+400+100') #Create the frame for the about window helpframe = ttk.Frame(self.helpDialog, width='200', height='240') ##helpframe.grid_propagate(0) helpframe.grid(column=0, row=0, sticky=(Tkinter.N, Tkinter.W, Tkinter.E, Tkinter.S), padx=5, pady=5) helpframe.columnconfigure(0, weight=1) helpframe.rowconfigure(0, weight=1) ttk.Label(helpframe, text=self.window_title+' Help', font='helvetica 12 bold', anchor="center").grid(column=0, row=0, sticky=(Tkinter.N,Tkinter.W, Tkinter.E)) helptext = Tkinter.Text(helpframe, width=40,height=9) helpmessage = "Peak Finder is used to find the peak current for a methylene blue reduction peak for CH instruments voltammogram data files." +\ "The data files must be sequentially numbered in the specified directory and be in text format.\n\n" + \ "The initial and final potential should be adjusted to lie outside the gaussian peak." helptext.insert('1.0',helpmessage) helptext.config(state='disabled', bg='SystemButtonFace', borderwidth=0, font='helvetica 10', wrap='word') helptext.grid(column=0, row=1,sticky=Tkinter.N, padx=10) ttk.Button(helpframe, text="Close", command=self.helpDialog.destroy).grid(column=0, row=4, sticky=(Tkinter.S)) self.helpDialog.mainloop() def dir_select(self,*args): '''PeakFinderApp.dir_select() allows the user to select a directory where datafiles are located.''' try: longpath = tkFileDialog.askdirectory() self.longpath_.set(longpath) self.dir_selected.set(1) return longpath except: pass def saveParams(self): '''Save the analysis parameters.''' if int(self.dir_selected.get()) == 1: storage = [] storage.append(self.filename_.get()) storage.append(self.format_.get()) storage.append(self.longpath_.get()) storage.append(self.init_potential_.get()) storage.append(self.final_potential_.get()) self.file_opt = options = {} options['defaultextension'] = '' # couldn't figure out how this works options['filetypes'] = [('data files', '.dat')] options['initialdir'] = self.longpath_.get() options['initialfile'] = 'myfile.dat' options['parent'] = self options['title'] = 'Save Parameters' storagefile = tkFileDialog.asksaveasfilename(**self.file_opt) if storagefile: f = open(storagefile, 'wb') pickle.dump(storage, f) f.close() else: #Display a directory selection warning message self.warnDialog = Tkinter.Toplevel(self) self.warnDialog.title("Warning") self.warnDialog.resizable(Tkinter.FALSE,Tkinter.FALSE) self.warnDialog.geometry('+120+120') #Create the frame for the about window warnframe = ttk.Frame(self.warnDialog, width='200', height='100') warnframe.grid_propagate(0) warnframe.grid(column=0, row=0, sticky=(Tkinter.N, Tkinter.W, Tkinter.E, Tkinter.S), padx=5, pady=5) warnframe.columnconfigure(0, weight=1) warnframe.rowconfigure(0, weight=1) ttk.Label(warnframe, text="Select Directory", font='helvetica 12 bold', anchor="center").grid(column=0, row=0, sticky=(Tkinter.N,Tkinter.W, Tkinter.E)) ttk.Label(warnframe, text="An analysis directory must be included\nin saved parameters.\n\n",anchor="center", justify="center").grid(column=0, row=1,sticky=Tkinter.N) ttk.Button(warnframe, text="Close", command=self.warnDialog.destroy).grid(column=0, row=4, sticky=(Tkinter.S)) self.warnDialog.mainloop() def loadParams(self): '''Retrieve a file of analysis parameters.''' self.file_opt = options = {} options['defaultextension'] = '' # couldn't figure out how this works options['filetypes'] = [('data files', '.dat')] options['initialdir'] = self.longpath_.get() options['initialfile'] = 'myfile.dat' options['parent'] = self options['title'] = 'Save Parameters' storagefile = tkFileDialog.askopenfilename(**self.file_opt) if storagefile: f = open(storagefile, 'rb') storage = pickle.load(f) f.close() self.filename_.set(storage[0]) self.format_.set(storage[1]) self.longpath_.set(storage[2]) self.init_potential_.set(storage[3]) self.final_potential_.set(storage[4]) self.dir_selected.set(1) class PeakLogic(): '''This is the internal logic of Peak Finder. Peaklogic looks in a specified directory for files in a given naming format, then fits the square wave voltammogram data inside to a non-linear polynomial plus gaussian function, subtracts the polynomial, and reports the peak current. Ultimately, it prints a pretty graph of the data, which can show multiple runs.''' def __init__(self,app): '''Initialize PeakLogic, passing the calling gui object to this class.''' self.app = app def peakfinder(self,*args): '''PeakLogic.peakfinder() is the core function to find and report peak current values. PeakLogic.peakfinder() looks in the user-selected directory for files in the format specified, then sends paths, filenames, and number of trials to thePeakLogic.loopForPeaks function.''' if int(self.app.dir_selected.get()) == 1: # Make sure the user has selected a directory try: #grab the text variables that we need from the gui filename = str(self.app.filename_.get()) longpath = str(self.app.longpath_.get()) format = str(self.app.format_.get()) self.app.status.set("Running...") #set the status bar os.chdir(longpath) #open our self.output file and set it up self.g = open(str(filename) + ".csv", 'w') self.g.write(",".join(('scan number','timestamp','time (s)','peak current')) + '\n') trials = len(glob.glob("%s/%s*.txt" % (longpath, format))) #detect number of runs self.app.bar.set_maxval(int(trials)+10) #set the progressbar size based on trial number timelist, iplist = self.loopForPeaks(longpath,format,trials) #run the peakfinder math self.g.close() #close the output file after it is filed by loopForPeaks self.app.output.set("Wrote output to " + filename + ".csv") #Show the user what was found self.grapher(timelist,iplist,'Seconds','Current',filename) #display a graph of the data return iplist except ValueError: pass def loopForPeaks(self, longpath, format, trials): '''PeakLogic.loopForPeaks() will loop through each file, collecting data and sending it to the peakMath function, as well as extracting time information.''' #clear some lists to hold our data full_x_lists, full_y_lists, timelist, realtimelist = [], [], [], [] startT = 0 try: pylab.close(2) #close test fitting graph if open except: pass for x in range (1, int(trials)+1): #loop through each file try: file = os.path.join(longpath,format + str(x) + ".txt") self.f = csv.reader(open(file),delimiter=',') listfile = list(self.f) x_list, y_list, rx_list = [], [], [] #clear some lists to hold data rt = listfile[0][1].split(' ') #grab the time the trial was run at convT = str(rt[4]).split(':') #split the time into individual parts realtimelist.append(str(rt[4])) tupT = (2010,5,25,int(convT[0]),int(convT[1]),int(convT[2]),1,145,0) #turn the time into full date format pointT = time.mktime(tupT) #convert the full time format into Unix time if x == 1: #if it's the first data point, set the initial time to zero startT = pointT pointTcorr = pointT - startT # subtract initial time to get time since the start of the trial x_list.append(pointTcorr) datalist = listfile[35:] #remove the header rows from the file, leaving just the data for row in datalist: if row != []: #skip empty lines in the file if row[0]: rx_list.append(row[0]) if row[1]: y_list.append(row[1]) full_x_lists.append(rx_list) full_y_lists.append(y_list) timelist.append(pointTcorr) except: pass #print len(full_x_lists) #print len(full_y_lists) iplist = self.peakMath(full_x_lists,full_y_lists) #feed the extracted x and y to PeakLogic.peakMath() #print "received ip list" for i,v,y in zip(iplist,timelist,realtimelist): #write the output csv file self.g.write(str(x) + ',' + str(y) + ','+ str(v) + ',' + str(i) + '\n') return timelist, iplist #return time and peak current for graphing def peakMath(self,listsx,listsy): '''PeakLogic.peakMath() feeds each x-y file to the fittingMath() function, and returns a list of peak currents (ip).''' iplist = [] self.count = 1 for xfile,yfile in zip(listsx,listsy): ip = max(0, self.fittingMath(xfile,yfile,1)) #for weird peak currents, set value to zero iplist.append(ip) self.app.bar.update(self.count) #update the progress bar self.count += 1 return iplist def fittingMath(self,xfile,yfile,flag=1): '''PeakLogic.fittingMath() fits the data to a polynomial and a gaussian, then subtracts the polynomial to find peak current..''' try: center = self.app.peak_center_.get() v0 = [1.5e-6,-1e-7,3e-6,5e-7,center,30,1e-9] #give reasonable starting values for non-linear regression x = numpy.array (xfile, dtype=numpy.float64) y = numpy.array (yfile, dtype=numpy.float64) #fp is full portion with the polynomial plus gaussian fp = lambda v, x: v[0] + (x*v[1])+(v[2]*(x**2))+(v[6]*(x**3)) + v[3]*numpy.exp(-(((x-v[4])**2)/2*v[5]**2)) #pp is the polynomial portion pp = lambda v, x: v[0] + (x*v[1])+(v[2]*(x**2))+(v[6]*(x**3)) #gp is the gaussian portion gp = lambda v, x: v[3]*numpy.exp(-(((x-v[4])**2)/2*v[5]**2)) ep = lambda v, x, y: (pp(v,x)-y) #ep is the error for a polynomial fit eg = lambda v, x, y: (gp(v,x)-y) #eg is the error for a gaussian fit e = lambda v, x, y: (fp(v,x)-y) #e is the error of the full fit from the real data #First, we will fit just the polynomial portion, ignoring the middle values passingx, passingy = self.truncEdges(xfile,yfile) # cut out outliers outx,outy = self.truncList(passingx,passingy) # cut out the middle values and return the edges v, success = leastsq(ep, v0, args=(outx,outy)) #fitting the polynomial y_subbed = [] #clear a list for putting revised y values for a,b in zip(xfile,yfile): #remove polynomial contribution from y-data a = numpy.array (a, dtype=numpy.float64) b = numpy.array (b, dtype=numpy.float64) temp = b - pp(v,a) y_subbed.append(temp) suby = numpy.array (y_subbed, dtype=numpy.float64) #put the result back into an array #Now, we can fit the gaussian portion v, success = leastsq(eg, v, args=(x,suby)) # fit the gaussin to the polynomial-subtracted data ip = max(fp(v,x)-pp(v,x)) #find the max value of the gaussian peak (minus the polynomial) if flag == 1: return ip if flag == 0: return x,y,fp(v,x),pp(v,x) except: return 0 def truncList(self,listx,listy): '''PeakLogic.truncList() removes the central portions of an x-y data list (where we suspect the gaussian is found) and returns the outer edges.''' newx, newy = [], [] #clear lists for spot,h in zip(listx,listy): spot = float(spot) #convert from string h = float(h) #convert from string low = float(self.app.final_potential_.get()) high = float(self.app.init_potential_.get()) if spot < low: #add low values newx.append(spot) newy.append(h) if spot > high: #add high values newx.append(spot) newy.append(h) #convert results back to an array px = numpy.array (newx, dtype=numpy.float64) py = numpy.array (newy, dtype=numpy.float64) return px, py #return partial x and partial y def truncEdges(self,listx,listy): '''PeakLogic.truncEdges() removes outlier regions of known bad signal from an x-y data list and returns the inner edges.''' newx, newy = [], [] #clear lists for spot,h in zip(listx,listy): spot = float(spot) #convert from string h = float(h) #convert from string low = float(self.app.final_edge_.get()) high = float(self.app.init_edge_.get()) if spot > low: #add low values #print "Ping" if spot < high: #print "Pong" newx.append(spot) newy.append(h) #else: #print "Inner Pass" #pass #else: #print "Outer Pass" #pass #convert results back to an array px = numpy.array (newx, dtype=numpy.float64) py = numpy.array (newy, dtype=numpy.float64) return px, py #return partial x and partial y def testFit(self): '''PeakLogic.testFit() performs a non-linear regression fit for the first data point and displays it for the user.''' if int(self.app.dir_selected.get()) == 1: # Make sure the user has selected a directory try: longpath = str(self.app.longpath_.get()) format = str(self.app.format_.get()) self.app.status.set("Running...") file = os.path.join(longpath,format + "1.txt") self.testfile = csv.reader(open(file),delimiter=',') #open the first data file listfile = [] for row in self.testfile: listfile.append(row) #turn the csv.reader object into a list datalist = listfile[35:] #remove the header rows from the file, leaving just the data x_list, y_list = [], [] #clear lists for row in datalist: if row != []: #skip empty lines if row[0]: x_list.append(row[0]) if row[1]: y_list.append(row[1]) x,y,y_fp,y_pp = self.fittingMath(x_list,y_list,flag=0) #send x and y data to fittingMath, request a verbose return self.testGrapher(x,y,y_fp,y_pp) #send fitted data to be graphed except: pass def grapher(self,x,y,labelX,labelY,fileTitle): '''PeakLogic.grapher() displays a nifty graph of the data using matplotlib.''' try: pylab.figure(1) line = pylab.plot(x, y, 'bo', x, y, linewidth=1.0,figure=1) pylab.xlabel(str(labelX)) pylab.ylabel(str(labelY)) pylab.title('Peak Current for ' + str(fileTitle)) pylab.grid(True) pylab.get_current_fig_manager().window.wm_geometry("+400+100") try: pylab.savefig(self.app.filename_.get()+".png") except: pass pylab.show() self.count = self.count + 9 self.app.bar.update(self.count) self.app.status.set("Finished") except: pass def testGrapher(self,x,y,y_fp,y_pp): '''PeakLogic.testGrapher() displays a graph of the test fitting.''' try: try: pylab.close(2) #close previous test fitting if open except: pass pylab.figure(2) pylab.plot(x,y,'ro',label='data') pylab.plot(x,y_pp,label='poly') pylab.plot(x,y_fp,label='full') pylab.legend() pylab.xlabel("Potential") pylab.ylabel("Current") pylab.title("Test regression of first datafile") pylab.get_current_fig_manager().window.wm_geometry("+400+100") pylab.show() self.app.status.set("Finished") except: pass class progressBar(): '''Create a Tkinter Progress bar widget.''' def __init__(self, root, width = 100, height = 10, maxval = 100): '''Initialize progressBar to make the Tkinter progressbar. This uses the canvas widget, since the native Tkinter progressbar doesnt work well.''' self.root = root #attach to the calling object gui self.maxval = float(maxval) self.canvas = Tkinter.Canvas(self.root, width = width, height = height, highlightt = 0, relief = 'ridge', borderwidth = 2) self.canvas.create_rectangle(0, 0, 0, 0, fill = 'blue') self.label = ttk.Label(self.root, text = 'Progress:').grid(column=1, row=16, sticky=Tkinter.W) self.canvas.grid(column=1,row=17, sticky=(Tkinter.W, Tkinter.E), columnspan=3) def set_maxval(self, maxval): '''progressBar.set_maxval() sets the max value of the progressbar.''' self.maxval = float(maxval) def update(self, value = 0): '''progressBar.update() updates the progressbar to a specified value.''' value = max(value, 0) value = min(value, self.maxval) self.canvas.delete(Tkinter.ALL) self.canvas.create_rectangle(0, 0, self.canvas.winfo_width() * value / self.maxval, self.canvas.winfo_reqheight(), fill='blue') self.root.update() class StatusBar(): """ A statusbar (not used) """ def __init__(self, master,columnset,rowset): #self.config(bg="green") self.label = ttk.Label(master, borderwidth=1, relief='ridge', background='#EEEEEE', anchor=Tkinter.W) self.label.grid(column=columnset,row=rowset,sticky=(Tkinter.W,Tkinter.E,Tkinter.S)) def set(self, input): """ Sets the statusbar label @param input: is the new text """ self.label.config(text=input) self.label.update_idletasks() def clear(self,stuff): """ Clears the status bar """ self.label.config(text="") self.label.update_idletasks() #### Party #### if __name__ == '__main__': app = PeakFinderApp()