From: devicerandom Date: Thu, 24 Apr 2008 13:49:48 +0000 (+0000) Subject: Initial SVN upload X-Git-Tag: 0.8.3 X-Git-Url: http://git.tremily.us/?a=commitdiff_plain;h=e11cbc5fa7c94d76ad41942f7f44758130869e48;p=hooke.git Initial SVN upload --- diff --git a/CHANGELOG b/CHANGELOG new file mode 100755 index 0000000..bc807fa --- /dev/null +++ b/CHANGELOG @@ -0,0 +1,391 @@ +0.8.4 +(2008-x-x) + PLUGINS: + macro.py: + hooke does not crash if it doesn't have permissions to create the folder + +0.8.3 +(2008-04-16) + PLUGINS: + generalvclamp.py: + fixed autopeak header + fixed autopeak slope (now unwanted slope values are discarded) + +0.8.2 +(2008-04-10) + PLUGINS: + flatfilts.py: + convfilt does not crash if a file is not a curve + generalvclamp.py: + autopeak now saves curve data correctly + autopeak now generates a dummy note (so that copylog/notelog is aware you measured the curve) + +0.8.1 +(2008-04-07) + PLUGINS: + generalvclamp.py: + fixed DeprecationWarning in flatten + flatfilts.py + convfilt now working + + +0.8.0: +(2008-04-04) + hooke.py: + sanity check of CLI plugins to avoid function overloading at startup + hooke_cli.py ; libhooke.py: + now playlists keep the index (when you reload the playlist, it starts from the + last observed curve) + updated plot to use _send_plot() + hooke.conf accepts lists as arguments for variables in + txt, export now have consistent argument order (thanks to A.G.Casado for pointing me that) + txt crashes no more if no filename is given (thanks to A.G.Casado for pointing me that) + libhookecurve.py: + added add_set() , remove_set() methods to make life easier for plugin writers + procplots.py: + plotmanip_correct() works with new picoforce.py deflection output (see) + PLUGINS: + fit.py: + updated wlc to use _send_plot() + wlc noauto now keeps the contact point + wlc reclick to click again the contact point + temperature now set in hooke.conf + generalvclamp.py: + implemented slope (thanks to Marco Brucale) + implemented autopeak + flatfilts.py: + convfilt,peaks use flattened curve + macro.py: + (new) added macro plugin (thanks to Alberto Gomez Casado) + DRIVERS: + picoforce.py: + fixed trigger bug! (thanks to Alberto Gomez Casado) + better deflection output (separated extension,retraction) + +0.7.5: +(2008-03-27) + hooke_cli.py: + removed outdated size command + PLUGINS: + generalvclamp.py: + implemented flatten + DRIVERS: + added tutorialdriver.py driver + csvdriver.py: + fixed (forgot close_all() method) + +0.7.4: +(2008-03-19) + added csvdriver driver + hooke_cli.py: + fixed plot manipulators handling (now it's safe to comment a + plot manipulator on hooke.conf) + PLUGINS: + fit.py: + fixed possible crash when clicking two times the same point on wlc + +0.7.3: +(2008-01-10) + hooke_cli.py: + fixed crash on copylog + PLUGINS: + massanalysis.py: + Initial release + tutorial.py: + Tutorial plugin, initial release + +0.7.2.1: +(2007-11-30) + PLUGINS: + flatfilt.py: + fixed crash on Windows + +0.7.2: +(2007-11-29) + hooke.py: + new configuration variable hookedir + hooke_cli.py: + copylog now checks if the destination is a real directory + fixed crashes in set + PLUGINS: + generalvclamp.py: + fixed a crash in forcebase when picking two times the same point + flatfilt.py: + fixed crash due to convfilt.conf impossible to load + initial implementation of the blind window for convfilt + initial data set maps (NOT FINISHED) + +0.7.1: +(2007-11-26) + PLUGINS: + flatfilts.py: + fixed possible crash in convfilt + implemented configuration file convfilt.conf + convfilt defaults are now 5 peaks 5 times more the noise absdev + implemented convconf + implemented setconf + libpeakspot.py: + fixed:now it really uses noise_absdev + +0.7.0: +(2007-11-15) + hooke_cli.py: + implemented _send_plot() helper API function + PLUGINS: + generalvclamp.py: + fixed forcebase to work with subtplot + flatfilts.py: + implemented convfilt! + added libpeakspot.py (helping library for convolution filter) + +0.6.5: +(2007-11-06) + hooke_cli.py, hooke.py: + plateau and contact (unmaintained) deleted and scheduled for re-release in generalvramp + implemented _measure_N_points() + PLUGINS: + generalvclamp.py: + implemented forcebase + fit.py: + wlc now accepts and uses temperature as an argument + wlc has been cleaned and uses new APIs + +0.6.4: +(2007-10-23) + hooke_cli.py, libhooke.py: + implemented support for defining order of plotmanip methods in hooke.conf + hooke_cli.py: + implemented delta + implemented point + attempted fix to bug 0033 (notelog crashing Hooke when using Unicode characters) + PLUGINS: + generalvramp.py: + began to move velocity ramp force spectroscopy-specific things in separate plugin + procplots.py: + added detriggerize; "set detrigger" 0/1 disables/enables it. + DRIVERS: + picoforce.py: + removed detriggerize() from driver + +0.6.3: +(2007-10-02) + hooke_cli.py: + rewritten txt command, now working + DRIVERS: + picoforce.py: + implemented detriggerize() to bypass the Picoforce trigger bug + PLUGINS: + superimpose.py: + implemented plotavgimpose + +0.6.2: +(2007-09-27) + hooke_cli.py: + fixed error handling in notelog + smarter handling of directory names in genlist + unexpected error handling in do_plot() + hooke.py: + implemented GetDisplayedPlot event and handlers + PLUGINS: + fit.py: + fixed (bug 0029) about replotting of wlc on a subtplot curve + multiple fitting displayed (to refine...) + +0.6.1: +(2007-08-06) + libhooke.py , hooke.py: + initial support for workdir configuration variable + libhooke.py: + fixed Driver() etc. semantics for gracefully handling unrecognized plots + hooke_cli.py: + fixed export namehandling + fixed plot error handling + PLUGINS: + flatfilts.py: + fixed memory leak + generalclamp.py: + fixed step command + +0.6.0 "Anko": +(2007-07-25) + hooke.py: + initial plugin support for the gui + wlc fitting now 100% plugin + measure_points replaces measure_couple etc. and provides much better extensibility + hooke_cli.py: + curves are sorted at beginning + PLUGINS: + procplots.py: + fft now allows for user selection of curve segment; select the plot; etc. + fit.py: + added gui section of plugin, now completely independent + fixed bug of wlc output + superimpose.py: + new plugin for superimposition of curve segments (still in development) + generalclamp.py: + all clamp commands now in a single plugin + implemented step + +0.5.4: +(2007-06-15) + procplots.py: + fixed fft crash with Numpy 1.0.1 + hooke.py: + fixed crashes if plot.scatter[] was empty + fixed management of multiple plots (bug #0025) + hooke_cli.py + fixed zpiezo error in measurement + hemingclamp.py, picoforce.py: + implemented close_all() method in drivers to avoid too many open files error + flatfilts.py: + fixed memory leak +0.5.3: +(2007-06-06) + wlc.py, hooke.py: + fixing and cleaning fit code: now the fit is part of a PlotObject and 100% coded in wlc.py + plotting of the wlc.py clicked points also begin to be part of a PlotObject + management of 'scatter' style property of plots + hooke_cli.py + fixed measuring error in defl, zpiezo + flatfilts.py: + slightly optimized has_features() routine + procplots.py: + fixed derivplot for every number of vectors + fixed possible crash of subtplot if applied on a file with != 2 plots + added fft command + libhookecurve.py: + fixed xaxis, yaxis for non-default plots: now defined from PlotObject + PlotObject now defines a styles[] vector +0.5.2: +(2007-05-21) + versioning a bit cleaned + fixed bug in hemingclamp.py preventing filename to appear + fixed wxversion problem for 2.8 + fixed too many open files bug (bug 0024) + added index command +0.5.1: +(2007-05-09) + using wxversion to choose from multiple wx versions + fixed old dependencies remaining +0.5.0 "Ingyo": +(2007-05-03) + general code updating and rewriting due to plugin support/better plot management + hooke.py: + initial plugin architecture for the command line. + initial plugin architecture for file drivers + initial plugin architecture for processing plots + export can now export both top and bottom plot (not together) + hooke_cli.py: + wlc fitting moved to fit.py plugin + flatfilt moved to flatfilts.py plugin + subtplot, derivplot moved to procplots.py plugin + double plot temporarily fixed for previous commands + export can now export both top and bottom plot (not together) + +0.4.1: +(2007-02-13) + hooke_cli.py: + double plot now default for clamp experiments + libhooke.py: + fixed bug that prevented flatfilt to work + (maybe) fixed memory leak in flatfilt + +0.4.0 "Hanzei": +(2007-02-08) + general code updating and rewriting due to double plot/force clamp supports + hooke.py: + initial dummy menu sketch + hooke.py, hooke_cli.py: + first general support in code for double plot: + - derivplot now in separate plot + - implemented show and close commands + - all functions should be double plot-aware + - clicking a point is double plot-aware + libhooke.py, hooke_cli.py: + general code cleanup: vectors_to_plot(), subtract_plot(), find_contact_point() and derivplot routines are now methods of class HookeCurve + hooke_cli.py: + implemented quit (alias of exit) + implemented version + libhooke.py, hooke.py, hooke_cli.py: + initial support for force clamp experiments: + - hemingclamp driver supported + - "experiment" flag describes what kind of experiment is a curve + - time, zpiezo, defl commands implemented + libhemingclamp.py: + inital release. + +0.3.1: + hooke.py: + fixed stupid bug in plateau + fixed bug in derivplot and subtplot not taking into account xaxes/yaxes variables +0.3.0: + from now on, all changelog is stored in CHANGELOG + hooke.py, libhooke.py, hooke_cli.py: + fixed plot and flatfilt crash when processing corrupt files + flatfilt output now more verbose + implemented system (execute an external OS command) + implemented copylog (copies annotated curves to a given directory) (todo 0033) + initial txt implementation (exports the current curve as a text file) (todo 0023) + fixed exit behaviour (bug 0013) + xaxes and yaxes variables now control visualization of plot (todo 0018) + new (better) contact point algorithm + workaround for the picoforce trigger bug +0.2.2 : + hooke.py, hooke_cli.py, libhooke.py: + support for fixed persistent length in WLC +0.2.1 : + hooke.py , libhooke.py: + fixed 'wlc noauto' bug (0012) preventing correct contact point to be used +0.2.0 : + hooke_cli.py: + implemented getlist (alias of genlist) + implemented contact (to plot the contact point) + fixed bug 0001 (Hooke crashes when opening a non-pf file) + fixed bug 0008 (Hooke crashes when generating a playlist with malformed namefiles/nonexistent files) + now the plot is refreshed after a "set" command (todo 0014) + wlc fit can use the (new) automatic contact point detection (old behaviour is preserved with "noauto" option) + hooke.py: + fixed versioning printing + complete refactoring of contact point routines + wlc fit adapted to use the (new) automatic contact point detection + wlc fit code a bit cleaned; parts moved to libhooke.py + libhooke.py: + new contact point algorithm (new algorithm) + wlc fit now uses a fancier domain (from contact point to a bit more than last point); initial chunk preparation section moved from hooke.py + + +OLDER CHANGELOGS: + +hooke.py: +0.1.1 : + From now on, all changelog is stored in hooke.py + hooke_cli.py: + corrected bug 0010 (addtolist bug), alerts when hitting start/end of playlist +2006_09_15_devel=0.1.0: initial WLC fit support. We hit 0.1 milestone :D +2006_08_28_devel: refactoring of plot interaction +2006_06_14_devel: fixed libhooke calls +2006_06_08_devel: initial automatic contact point finding +2006_05_30_devel: configuration file support + +hooke_cli.py: +0.1.1 : from now on, all changelog is in hooke.py +2006_09_15_devel: implemented wlc; 0.1.0 milestone. +2006_08_28_devel: refactoring of plot interaction +2006_07_23_devel: implemented note; implemented flatfilt; implemented notelog; exit now warns if playlist/notes + have not been saved. +2006_07_18_devel: implemented subtplot; bug 0007 ("cd" crashing) fixed +2006_06_16_devel: moved math helper functions in libhooke.py +2006_06_14_devel: fixed "jump" output; fixed "exit" (now it works!); fixed off-by-one bug in deflection-correction +2006_06_08_devel: fixed "loadlist" output; +2006_05_30_devel: initial configuration file support; added "set" command; initial deflection-correction support; added "ls" command as an alias of "dir" +2006_05_23_devel: rewriting of playlist-handling code due to major rewrite of hooke_playlist.py + +libhooke.py +0.1.1 : from now on, all changelog is in hooke.py +2006_09_15_devel : initial WLC support +2006_09_14_devel : initial support for Hemingway velocity clamp files, minor refactorings +2006_07_22_devel : implemented math function has_features +2006_06_16_devel : math functions moved here +2006_06_08_devel : hooke_playlist.py becomes libhooke.py +2006_05_30_devel : support for deflection in HookeCurve +2006_05_29_devel : Initial configuration file support +2006_05_23_devel : Major rewrite. Fixed bug 0002 \ No newline at end of file diff --git a/convfilt.conf b/convfilt.conf new file mode 100644 index 0000000..ef5254c --- /dev/null +++ b/convfilt.conf @@ -0,0 +1,20 @@ + + + + + + + + + + \ No newline at end of file diff --git a/csvdriver.py b/csvdriver.py new file mode 100644 index 0000000..4222255 --- /dev/null +++ b/csvdriver.py @@ -0,0 +1,75 @@ +#!/usr/bin/env python + +''' +csvdriver.py + +Simple driver to read general comma-separated values in Hooke + +Columns are read this way: + +X1 , Y1 , X2 , Y2 , X3 , Y3 ... + +If the number of columns is odd, the last column is ignored. + +(c)Massimo Sandal, 2008 +''' + +import libhookecurve as lhc +import libhooke as lh +import csv + +class csvdriverDriver(lhc.Driver): + + def __init__(self, filename): + + self.filedata = open(filename,'r') + self.data = list(self.filedata) + self.filedata.close() + + self.filetype = 'generic' + self.experiment = '' + + self.filename=filename + + def is_me(self): + myfile=file(self.filename) + headerline=myfile.readlines()[0] + myfile.close() + + #using a custom header makes things much easier... + #(looking for raw CSV data is at strong risk of confusion) + if headerline[:-1]=='Hooke data': + return True + else: + return False + + def close_all(self): + self.filedata.close() + + def default_plots(self): + rrows=csv.reader(self.data) + rows=list(rrows) #transform the csv.reader iterator in a normal list + columns=lh.transposed2(rows[1:]) + + main_plot=lhc.PlotObject() + main_plot.vectors=[] + + for index in range(0,len(columns),2): + main_plot.vectors.append([]) + temp_x=columns[index] + temp_y=columns[index+1] + + #convert to float (the csv gives strings) + temp_x=[float(item) for item in temp_x] + temp_y=[float(item) for item in temp_y] + + main_plot.vectors[-1].append(temp_x) + main_plot.vectors[-1].append(temp_y) + + main_plot.units=['x','y'] + main_plot.title=self.filename + main_plot.destination=0 + + return [main_plot] + + \ No newline at end of file diff --git a/fit.py b/fit.py new file mode 100755 index 0000000..b8f4717 --- /dev/null +++ b/fit.py @@ -0,0 +1,440 @@ +#!/usr/bin/env python + +''' +FIT + +Force spectroscopy curves basic fitting plugin. +Licensed under the GNU GPL version 2 + +Non-standard Dependencies: +procplots.py (plot processing plugin) +''' +from libhooke import WX_GOOD, ClickedPoint +import wxversion +wxversion.select(WX_GOOD) +#from wx import PostEvent +#from wx.lib.newevent import NewEvent +import scipy +import numpy as np +import copy +import Queue + +global measure_wlc +global EVT_MEASURE_WLC + +#measure_wlc, EVT_MEASURE_WLC = NewEvent() + +global events_from_fit +events_from_fit=Queue.Queue() #GUI ---> CLI COMMUNICATION + + +class fitCommands: + + def _plug_init(self): + self.wlccurrent=None + self.wlccontact_point=None + self.wlccontact_index=None + + def wlc_fit(self,clicked_points,xvector,yvector, pl_value, T=293): + ''' + Worm-like chain model fitting. + The function is the simple polynomial worm-like chain as proposed by C.Bustamante, J.F.Marko, E.D.Siggia + and S.Smith (Science. 1994 Sep 9;265(5178):1599-600.) + ''' + + '''clicked_points[0] = contact point (calculated or hand-clicked) + clicked_points[1] and [2] are edges of chunk''' + + #STEP 1: Prepare the vectors to apply the fit. + + #indexes of the selected chunk + first_index=min(clicked_points[1].index, clicked_points[2].index) + last_index=max(clicked_points[1].index, clicked_points[2].index) + + #getting the chunk and reverting it + xchunk,ychunk=xvector[first_index:last_index],yvector[first_index:last_index] + xchunk.reverse() + ychunk.reverse() + #put contact point at zero and flip around the contact point (the fit wants a positive growth for extension and force) + xchunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xchunk] + ychunk_corr_up=[-(y-clicked_points[0].graph_coords[1]) for y in ychunk] + #make them arrays + xchunk_corr_up=scipy.array(xchunk_corr_up) + ychunk_corr_up=scipy.array(ychunk_corr_up) + + + #STEP 2: actually do the fit + + #Find furthest point of chunk and add it a bit; the fit must converge + #from an excess! + xchunk_high=max(xchunk_corr_up) + xchunk_high+=(xchunk_high/10) + + #Here are the linearized start parameters for the WLC. + #[lambd=1/Lo , pii=1/P] + + p0=[(1/xchunk_high),(1/(3.5e-10))] + p0_plfix=[(1/xchunk_high)] + + def residuals(params,y,x,T): + ''' + Calculates the residuals of the fit + ''' + lambd, pii=params + + Kb=(1.38065e-23) + #T=293 + therm=Kb*T + + err = y-( (therm*pii/4) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd)) ) + + return err + + def wlc_eval(x,params,pl_value,T): + ''' + Evaluates the WLC function + ''' + if not pl_value: + lambd, pii = params + else: + lambd = params + + if pl_value: + pii=1/pl_value + + Kb=(1.38065e-23) #boltzmann constant + #T=293 #temperature FIXME:should be user-modifiable! + therm=Kb*T #so we have thermal energy + + return ( (therm*pii/4.0) * (((1-(x*lambd))**-2.0) - 1 + (4.0*x*lambd)) ) + + def residuals_plfix(params, y, x, pii, T): + ''' + Calculates the residuals of the fit, if we have the persistent length from an external source + ''' + lambd=params + + Kb=(1.38065e-23) + therm=Kb*T + + err = y-( (therm*pii/4) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd)) ) + + return err + + #make the fit! and obtain params + if pl_value: + plsq=scipy.optimize.leastsq(residuals_plfix, p0_plfix, args=(ychunk_corr_up,xchunk_corr_up,1/pl_value,T)) + else: + plsq=scipy.optimize.leastsq(residuals, p0, args=(ychunk_corr_up,xchunk_corr_up,T)) + + + #STEP 3: plotting the fit + + #obtain domain to plot the fit - from contact point to last_index plus 20 points + thule_index=last_index+10 + if thule_index > len(xvector): #for rare cases in which we fit something at the END of whole curve. + thule_index = len(xvector) + #reverse etc. the domain + xfit_chunk=xvector[clicked_points[0].index:thule_index] + xfit_chunk.reverse() + xfit_chunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xfit_chunk] + xfit_chunk_corr_up=scipy.array(xfit_chunk_corr_up) + + #the fitted curve: reflip, re-uncorrect + yfit=wlc_eval(xfit_chunk_corr_up, plsq[0],pl_value,T) + yfit_down=[-y for y in yfit] + yfit_corr_down=[y+clicked_points[0].graph_coords[1] for y in yfit_down] + + #get out true fit paramers + fit_out=plsq[0] + try: + fit_out=[(1.0/x) for x in fit_out] + except TypeError: #if we fit only 1 parameter, we have a float and not a list in output. + fit_out=[(1.0/fit_out)] + + return fit_out, yfit_corr_down, xfit_chunk + + + def do_wlc(self,args): + ''' + WLC + (fit plugin) + Fits a worm-like chain entropic rise to a given chunk of the curve. + + First you have to click a contact point. + Then you have to click the two edges of the data you want to fit. + The function is the simple polynomial worm-like chain as proposed by + C.Bustamante, J.F.Marko, E.D.Siggia and S.Smith (Science. 1994 + Sep 9;265(5178):1599-600.) + + Arguments: + pl=[value] : Use a fixed persistent length for the fit. If pl is not given, + the fit will be a 2-variable + fit. DO NOT put spaces between 'pl', '=' and the value. + The value must be in meters. + Scientific notation like 0.35e-9 is fine. + + t=[value] : Use a user-defined temperature. The value must be in + kelvins; by default it is 293 K. + DO NOT put spaces between 't', '=' and the value. + + noauto : allows for clicking the contact point by + hand (otherwise it is automatically estimated) the first time. + If subsequent measurements are made, the same contact point + clicked is used + + reclick : redefines by hand the contact point, if noauto has been used before + but the user is unsatisfied of the previously choosen contact point. + --------- + Syntax: wlc [pl=(value)] [t=value] [noauto] + ''' + pl_value=None + T=self.config['temperature'] + for arg in args.split(): + #look for a persistent length argument. + if 'pl=' in arg: + pl_expression=arg.split('=') + pl_value=float(pl_expression[1]) #actual value + #look for a T argument. FIXME: spaces are not allowed between 'pl' and value + if ('t=' in arg[0:2]) or ('T=' in arg[0:2]): + t_expression=arg.split('=') + T=float(t_expression[1]) + + #use the currently displayed plot for the fit + displayed_plot=self._get_displayed_plot() + + #handle contact point arguments correctly + if 'reclick' in args.split(): + print 'Click contact point' + contact_point=self._measure_N_points(N=1, whatset=1)[0] + contact_point_index=contact_point.index + self.wlccontact_point=contact_point + self.wlccontact_index=contact_point.index + self.wlccurrent=self.current.path + elif 'noauto' in args.split(): + if self.wlccontact_index==None or self.wlccurrent != self.current.path: + print 'Click contact point' + contact_point=self._measure_N_points(N=1, whatset=1)[0] + contact_point_index=contact_point.index + self.wlccontact_point=contact_point + self.wlccontact_index=contact_point.index + self.wlccurrent=self.current.path + else: + contact_point=self.wlccontact_point + contact_point_index=self.wlccontact_index + else: + cindex=self.find_contact_point() + contact_point=ClickedPoint() + contact_point.absolute_coords=displayed_plot.vectors[1][0][cindex], displayed_plot.vectors[1][1][cindex] + contact_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1]) + contact_point.is_marker=True + + print 'Click edges of chunk' + points=self._measure_N_points(N=2, whatset=1) + points=[contact_point]+points + params, yfit, xfit = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T) + try: + params, yfit, xfit = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T) + except: + print 'Fit not possible. Probably wrong interval -did you click two *different* points?' + return + + print 'Contour length: ',params[0]*(1.0e+9),' nm' + if len(params)==2: #if we did choose 2-value fit + print 'Persistent length: ',params[1]*(1.0e+9),' nm' + + #add the clicked points in the final PlotObject + clickvector_x, clickvector_y=[], [] + for item in points: + clickvector_x.append(item.graph_coords[0]) + clickvector_y.append(item.graph_coords[1]) + + #create a custom PlotObject to gracefully plot the fit along the curves + + fitplot=copy.deepcopy(displayed_plot) + fitplot.add_set(xfit,yfit) + fitplot.add_set(clickvector_x,clickvector_y) + + if fitplot.styles==[]: + fitplot.styles=[None,None,None,'scatter'] + else: + fitplot.styles+=[None,'scatter'] + + self._send_plot([fitplot]) + + def find_contact_point(self): + ''' + Finds the contact point on the curve. + + The current algorithm (thanks to Francesco Musiani, francesco.musiani@unibo.it and Massimo Sandal) is: + - take care of the PicoForce trigger bug - exclude retraction portions with too high standard deviation + - fit the second half of the retraction curve to a line + - if the fit is not almost horizontal, take a smaller chunk and repeat + - otherwise, we have something horizontal + - so take the average of horizontal points and use it as a baseline + + Then, start from the rise of the retraction curve and look at the first point below the + baseline. + + FIXME: should be moved, probably to generalvclamp.py + ''' + outplot=self.subtract_curves(1) + xret=outplot.vectors[1][0] + ydiff=outplot.vectors[1][1] + + xext=self.plots[0].vectors[0][0] + yext=self.plots[0].vectors[0][1] + xret2=self.plots[0].vectors[1][0] + yret=self.plots[0].vectors[1][1] + + #taking care of the picoforce trigger bug: we exclude portions of the curve that have too much + #standard deviation. yes, a lot of magic is here. + monster=True + monlength=len(xret)-int(len(xret)/20) + finalength=len(xret) + while monster: + monchunk=scipy.array(ydiff[monlength:finalength]) + if abs(scipy.stats.std(monchunk)) < 2e-10: + monster=False + else: #move away from the monster + monlength-=int(len(xret)/50) + finalength-=int(len(xret)/50) + + + #take half of the thing + endlength=int(len(xret)/2) + + ok=False + while not ok: + xchunk=yext[endlength:monlength] + ychunk=yext[endlength:monlength] + regr=scipy.stats.linregress(xchunk,ychunk)[0:2] + #we stop if we found an almost-horizontal fit or if we're going too short... + #FIXME: 0.1 and 6 here are "magic numbers" (although reasonable) + if (abs(regr[1]) > 0.1) and ( endlength < len(xret)-int(len(xret)/6) ) : + endlength+=10 + else: + ok=True + + ymean=scipy.mean(ychunk) #baseline + + index=0 + point = ymean+1 + + #find the first point below the calculated baseline + while point > ymean: + try: + point=yret[index] + index+=1 + except IndexError: + #The algorithm didn't find anything below the baseline! It should NEVER happen + index=0 + return index + + + + def find_contact_point2(self, debug=False): + ''' + TO BE DEVELOPED IN THE FUTURE + Finds the contact point on the curve. + + FIXME: should be moved, probably to generalvclamp.py + ''' + outplot=self.subtract_curves(1) + xret=outplot.vectors[1][0] + ydiff=outplot.vectors[1][1] + + raw_plot=self.current.curve.default_plots()[0] + + '''xext=self.plots[0].vectors[0][0] + yext=self.plots[0].vectors[0][1] + xret2=self.plots[0].vectors[1][0] + yret=self.plots[0].vectors[1][1] + ''' + xext=raw_plot.vectors[0][0] + yext=raw_plot.vectors[0][1] + xret2=raw_plot.vectors[1][0] + yret=raw_plot.vectors[1][1] + #taking care of the picoforce trigger bug: we exclude portions of the curve that have too much + #standard deviation. yes, a lot of magic is here. + + monlength=len(xext) + #take half of the thing + endlength=int(len(xext)/2) + + ok=False + xchunk=xext[endlength:monlength] + ychunk=yext[endlength:monlength] + regr=scipy.polyfit(xchunk,ychunk,1)[0:2] + ''' + while not ok: + xchunk=yext[endlength:monlength] + ychunk=yext[endlength:monlength] + print len(xchunk) + #regr=scipy.stats.linregress(xchunk,ychunk)[0:2] + regr=scipy.polyfit(xchunk,ychunk,1)[0:2] + #we stop if we found an almost-horizontal fit or if we're going too short... + #FIXME: 0.1 and 6 here are "magic numbers" (although reasonable) + if (abs(regr[1]) > 0.1) and ( endlength < len(xret)-int(len(xret)/6) ) : + endlength+=10 + else: + ok=True + ''' + + ''' + ymean=scipy.mean(ychunk) #baseline + + index=0 + point = ymean+1 + + #find the first point below the calculated baseline + while point > ymean: + try: + point=yret[index] + index+=1 + except IndexError: + #The algorithm didn't find anything below the baseline! It should NEVER happen + index=0 + ''' + #fit the contact line + n_contact=100 + x_contact_chunk=xret2[0:n_contact] + y_contact_chunk=yret[0:n_contact] + + regr_contact=scipy.polyfit(x_contact_chunk,y_contact_chunk,1)[0:2] + x_intercept=(regr_contact[1]-regr[1])/(regr[0]-regr_contact[0]) + y_intercept=(regr_contact[0]*x_intercept + regr_contact[1]) + + #now, exploit a ClickedPoint instance to calculate index... + dummy=ClickedPoint() + dummy.absolute_coords=(x_intercept,y_intercept) + dummy.find_graph_coords(xret2,yret) + + if debug: + return dummy.index, regr, regr_contact + else: + return dummy.index + + + def x_do_contact(self,args): + ''' + DEBUG COMMAND to be activated in the future + ''' + index,regr,regr_contact=self.find_contact_point2(debug=True) + print regr + print regr_contact + raw_plot=self.current.curve.default_plots()[0] + xret=raw_plot.vectors[0][0] + #nc_line=[(item*regr[0])+regr[1] for item in x_nc] + nc_line=scipy.polyval(regr,xret) + c_line=scipy.polyval(regr_contact,xret) + + + contact_plot=self.current.curve.default_plots()[0] + contact_plot.add_set(xret, nc_line) + contact_plot.add_set(xret, c_line) + contact_plot.styles=[None,None,None,None] + #contact_plot.styles.append(None) + contact_plot.destination=1 + self._send_plot([contact_plot]) + \ No newline at end of file diff --git a/flatfilts.py b/flatfilts.py new file mode 100755 index 0000000..27e526f --- /dev/null +++ b/flatfilts.py @@ -0,0 +1,421 @@ +#!/usr/bin/env python + +''' +FLATFILTS + +Force spectroscopy curves filtering of flat curves +Licensed under the GNU LGPL version 2 + +Other plugin dependencies: +procplots.py (plot processing plugin) +''' +from libhooke import WX_GOOD +import wxversion +wxversion.select(WX_GOOD) + +import xml.dom.minidom + +import wx +import scipy +import numpy +from numpy import diff + +#import pickle + +import libpeakspot as lps +import libhookecurve as lhc + + +class flatfiltsCommands: + + def _plug_init(self): + #configurate convfilt variables + convfilt_configurator=ConvfiltConfig() + + #different OSes have different path conventions + if self.config['hookedir'][0]=='/': + slash='/' #a Unix or Unix-like system + else: + slash='\\' #it's a drive letter, we assume it's Windows + + self.convfilt_config=convfilt_configurator.load_config(self.config['hookedir']+slash+'convfilt.conf') + + def do_flatfilt(self,args): + ''' + FLATFILT + (flatfilts.py) + Filters out flat (featureless) curves of the current playlist, + creating a playlist containing only the curves with potential + features. + ------------ + Syntax: + flatfilt [min_npks min_deviation] + + min_npks = minmum number of points over the deviation + (default=4) + + min_deviation = minimum signal/noise ratio + (default=9) + + If called without arguments, it uses default values, that + should work most of the times. + ''' + median_filter=7 + min_npks=4 + min_deviation=9 + + args=args.split(' ') + if len(args) == 2: + min_npks=int(args[0]) + min_deviation=int(args[1]) + else: + pass + + print 'Processing playlist...' + notflat_list=[] + + c=0 + for item in self.current_list: + c+=1 + + try: + notflat=self.has_features(item, median_filter, min_npks, min_deviation) + print 'Curve',item.path, 'is',c,'of',len(self.current_list),': features are ',notflat + except: + notflat=False + print 'Curve',item.path, 'is',c,'of',len(self.current_list),': cannot be filtered. Probably unable to retrieve force data from corrupt file.' + + if notflat: + item.features=notflat + item.curve=None #empty the item object, to further avoid memory leak + notflat_list.append(item) + + if len(notflat_list)==0: + print 'Found nothing interesting. Check your playlist, could be a bug or criteria could be too much stringent' + return + else: + print 'Found ',len(notflat_list),' potentially interesting curves' + print 'Regenerating playlist...' + self.pointer=0 + self.current_list=notflat_list + self.current=self.current_list[self.pointer] + self.do_plot(0) + + def has_features(self,item,median_filter,min_npks,min_deviation): + ''' + decides if a curve is flat enough to be rejected from analysis: it sees if there + are at least min_npks points that are higher than min_deviation times the absolute value + of noise. + + Algorithm original idea by Francesco Musiani, with my tweaks and corrections. + ''' + retvalue=False + + item.identify(self.drivers) + #we assume the first is the plot with the force curve + #do the median to better resolve features from noise + flat_plot=self.plotmanip_median(item.curve.default_plots()[0], item, customvalue=median_filter) + flat_vects=flat_plot.vectors + item.curve.close_all() + #needed to avoid *big* memory leaks! + del item.curve + del item + + #absolute value of derivate + yretdiff=diff(flat_vects[1][1]) + yretdiff=[abs(value) for value in yretdiff] + #average of derivate values + diffmean=numpy.mean(yretdiff) + yretdiff.sort() + yretdiff.reverse() + c_pks=0 + for value in yretdiff: + if value/diffmean > min_deviation: + c_pks+=1 + else: + break + + if c_pks>=min_npks: + retvalue = c_pks + + del flat_plot, flat_vects, yretdiff + + return retvalue + + ################################################################ + #-----CONVFILT------------------------------------------------- + #-----Convolution-based peak recognition and filtering. + #Requires the libpeakspot.py library + + def has_peaks(self, plot, abs_devs): + ''' + Finds peak position in a force curve. + FIXME: should be moved in libpeakspot.py + ''' + xret=plot.vectors[1][0] + yret=plot.vectors[1][1] + #Calculate convolution. + convoluted=lps.conv_dx(yret, self.convfilt_config['convolution']) + + #surely cut everything before the contact point + cut_index=self.find_contact_point() + + #cut even more, before the blind window + start_x=xret[cut_index] + blind_index=0 + for value in xret[cut_index:]: + if abs((value) - (start_x)) > self.convfilt_config['blindwindow']*(10**-9): + break + blind_index+=1 + cut_index+=blind_index + + #do the dirty convolution-peak finding stuff + noise_level=lps.noise_absdev(convoluted[cut_index:], self.convfilt_config['positive'], self.convfilt_config['maxcut'], self.convfilt_config['stable']) + above=lps.abovenoise(convoluted,noise_level,cut_index,abs_devs) + peak_location,peak_size=lps.find_peaks(above) + + #take the maximum + for i in range(len(peak_location)): + peak=peak_location[i] + maxpk=min(yret[peak-10:peak+10]) + index_maxpk=yret[peak-10:peak+10].index(maxpk)+(peak-10) + peak_location[i]=index_maxpk + + return peak_location,peak_size + + + def exec_has_peaks(self,item,abs_devs): + ''' + encapsulates has_peaks for the purpose of correctly treating the curve objects in the convfilt loop, + to avoid memory leaks + ''' + item.identify(self.drivers) + #we assume the first is the plot with the force curve + plot=item.curve.default_plots()[0] + + if 'flatten' in self.config['plotmanips']: + #If flatten is present, use it for better recognition of peaks... + flatten=self._find_plotmanip('flatten') #extract flatten plot manipulator + plot=flatten(plot, item, customvalue=1) + + peak_location,peak_size=self.has_peaks(plot,abs_devs) + #close all open files + item.curve.close_all() + #needed to avoid *big* memory leaks! + del item.curve + del item + return peak_location, peak_size + + #------------------------ + #------commands---------- + #------------------------ + def do_peaks(self,args): + ''' + PEAKS + (flatfilts.py) + Test command for convolution filter / test. + ---- + Syntax: peaks [deviations] + absolute deviation = number of times the convolution signal is above the noise absolute deviation. + Default is 5. + ''' + if len(args)==0: + args=self.convfilt_config['mindeviation'] + + + + try: + abs_devs=float(args) + except: + pass + + defplots=self.current.curve.default_plots()[0] #we need the raw, uncorrected plots + + if 'flatten' in self.config['plotmanips']: + flatten=self._find_plotmanip('flatten') #extract flatten plot manipulator + defplots=flatten(defplots, self.current) + else: + print 'You have the flatten plot manipulator not loaded. Enabling it could give you better results.' + + peak_location,peak_size=self.has_peaks(defplots,abs_devs) + print 'Found '+str(len(peak_location))+' peaks.' + #print peak_location + + #if no peaks, we have nothing to plot. exit. + if len(peak_location)==0: + return + + #otherwise, we plot the peak locations. + xplotted_ret=self.plots[0].vectors[1][0] + yplotted_ret=self.plots[0].vectors[1][1] + xgood=[xplotted_ret[index] for index in peak_location] + ygood=[yplotted_ret[index] for index in peak_location] + + recplot=self._get_displayed_plot() + recplot.vectors.append([xgood,ygood]) + if recplot.styles==[]: + recplot.styles=[None,None,'scatter'] + else: + recplot.styles+=['scatter'] + + self._send_plot([recplot]) + + def do_convfilt(self,args): + ''' + CONVFILT + (flatfilts.py) + Filters out flat (featureless) curves of the current playlist, + creating a playlist containing only the curves with potential + features. + ------------ + Syntax: + convfilt [min_npks min_deviation] + + min_npks = minmum number of peaks + (default=4) + + min_deviation = minimum signal/noise ratio *in the convolution* + (default=4) + + If called without arguments, it uses default values. + ''' + + min_npks=self.convfilt_config['minpeaks'] + min_deviation=self.convfilt_config['mindeviation'] + + args=args.split(' ') + if len(args) == 2: + min_npks=int(args[0]) + min_deviation=int(args[1]) + else: + pass + + print 'Processing playlist...' + print '(Please wait)' + notflat_list=[] + + + c=0 + for item in self.current_list: + c+=1 + + try: + peak_location,peak_size=self.exec_has_peaks(item,min_deviation) + if len(peak_location)>=min_npks: + isok='+' + else: + isok='' + print 'Curve',item.path, 'is',c,'of',len(self.current_list),': found '+str(len(peak_location))+' peaks.'+isok + except: + peak_location,peak_size=[],[] + print 'Curve',item.path, 'is',c,'of',len(self.current_list),': cannot be filtered. Probably unable to retrieve force data from corrupt file.' + + if len(peak_location)>=min_npks: + item.peak_location=peak_location + item.peak_size=peak_size + item.curve=None #empty the item object, to further avoid memory leak + notflat_list.append(item) + + #Warn that no flattening had been done. + if not ('flatten' in self.config['plotmanips']): + print 'Flatten manipulator was not found. Processing was done without flattening.' + print 'Try to enable it in your configuration file for better results.' + + if len(notflat_list)==0: + print 'Found nothing interesting. Check your playlist, could be a bug or criteria could be too much stringent' + return + else: + print 'Found ',len(notflat_list),' potentially interesting curves' + print 'Regenerating playlist...' + self.pointer=0 + self.current_list=notflat_list + self.current=self.current_list[self.pointer] + self.do_plot(0) + + def do_convconf(self,args): + ''' + CONVCONFIG + (flatfilts.py) + Prints the current convfilt configuration variables. + ------ + Syntax: convconfig + ''' + print self.convfilt_config + + def do_setconv(self,args): + ''' + SETCONV + (flatfilts.py) + Sets the convfilt configuration variables + ------ + Syntax: setconv variable value + ''' + args=args.split() + try: + self.convfilt_config[args[0]]=eval(args[1]) + except NameError: + self.convfilt_config[args[0]]=args[1] + + +######################### +#HANDLING OF CONFIGURATION FILE +class ConvfiltConfig: + ''' + Handling of convfilt configuration file + + Mostly based on the simple-yet-useful examples of the Python Library Reference + about xml.dom.minidom + + FIXME: starting to look a mess, should require refactoring + ''' + + def __init__(self): + self.config={} + + + def load_config(self, filename): + myconfig=file(filename) + #the following 3 lines are needed to strip newlines. otherwise, since newlines + #are XML elements too, the parser would read them (and re-save them, multiplying + #newlines...) + #yes, I'm an XML n00b + the_file=myconfig.read() + the_file_lines=the_file.split('\n') + the_file=''.join(the_file_lines) + + self.config_tree=xml.dom.minidom.parseString(the_file) + + def getText(nodelist): + #take the text from a nodelist + #from Python Library Reference 13.7.2 + rc = '' + for node in nodelist: + if node.nodeType == node.TEXT_NODE: + rc += node.data + return rc + + def handleConfig(config): + noiseabsdev_elements=config.getElementsByTagName("noise_absdev") + convfilt_elements=config.getElementsByTagName("convfilt") + handleAbsdev(noiseabsdev_elements) + handleConvfilt(convfilt_elements) + + def handleAbsdev(noiseabsdev_elements): + for element in noiseabsdev_elements: + for attribute in element.attributes.keys(): + self.config[attribute]=element.getAttribute(attribute) + + def handleConvfilt(convfilt_elements): + for element in convfilt_elements: + for attribute in element.attributes.keys(): + self.config[attribute]=element.getAttribute(attribute) + + handleConfig(self.config_tree) + #making items in the dictionary machine-readable + for item in self.config.keys(): + try: + self.config[item]=eval(self.config[item]) + except NameError: #if it's an unreadable string, keep it as a string + pass + + return self.config \ No newline at end of file diff --git a/generalclamp.py b/generalclamp.py new file mode 100644 index 0000000..9390146 --- /dev/null +++ b/generalclamp.py @@ -0,0 +1,101 @@ +#!/usr/bin/env python + +''' +GENERALCLAMP.py + +Plugin regarding general force clamp measurements +''' +from libhooke import WX_GOOD, ClickedPoint +import wxversion +wxversion.select(WX_GOOD) +from wx import PostEvent + +class generalclampCommands: + + + def do_showdefl(self,args): + ''' + SHOWDEFL + Shows the deflection plot for a force clamp curve. + Use 'close' to close the plot. + --- + Syntax: showdefl + ''' + if self.current.curve.experiment != 'clamp': + print 'This command makes no sense for a non-force clamp experiment!' + else: + self.current.vectors_to_plot(self.config['correct'],self.config['medfilt'],yclamp='defl') + plot_graph=self.list_of_events['plot_graph'] + wx.PostEvent(self.frame,plot_graph(current=self.current,xaxes=self.config['xaxes'],yaxes=self.config['yaxes'], destination=1)) + + def do_time(self,args): + ''' + TIME + Measure the time difference (in seconds) between two points + Implemented only for force clamp + ---- + Syntax: time + ''' + if self.current.curve.experiment == 'clamp': + print 'Click two points.' + points=self._measure_N_points(N=2) + time=abs(points[0].graph_coords[0]-points[1].graph_coords[0]) + print str(time)+' s' + else: + print 'This command makes no sense for a non-force clamp experiment.' + + def do_zpiezo(self,args): + ''' + ZPIEZO + Measure the zpiezo difference (in nm) between two points + Implemented only for force clamp + ---- + Syntax: zpiezo + ''' + if self.current.curve.experiment == 'clamp': + print 'Click two points.' + points=self._measure_N_points(N=2) + zpiezo=abs(points[0].graph_coords[1]-points[1].graph_coords[1]) + print str(zpiezo*(10**9))+' nm' + else: + print 'This command makes no sense for a non-force clamp experiment.' + + def do_defl(self,args): + ''' + DEFL + Measure the deflection difference (in nm) between two points + Implemented only for force clamp + NOTE: It makes sense only on the time VS defl plot; it is still not masked for the other plot... + ----- + Syntax: defl + ''' + if self.current.curve.experiment == 'clamp': + print 'Click two points.' + points=self._measure_N_points(N=2) + defl=abs(points[0].graph_coords[1]-points[1].graph_coords[1]) + print str(defl*(10**12))+' pN' + else: + print 'This command makes no sense for a non-force clamp experiment.' + + def do_step(self,args): + ''' + STEP + + Measures the length and time duration of a time-Z step + ----- + Syntax: step + ''' + if self.current.curve.experiment == 'clamp': + print 'Click three points in this fashion:' + print ' (0)-------(1)' + print ' |' + print ' |' + print ' (2)----------' + points=self._measure_N_points(N=3,whatset=0) + dz=abs(points[2].graph_coords[1]-points[1].graph_coords[1])*(10e+8) + dt=abs(points[1].graph_coords[0]-points[0].graph_coords[0]) + print 'dZ: ',dz,' nm' + print 'dT: ',dt,' s' + + else: + print 'This command makes no sense for a non-force clamp experiment.' \ No newline at end of file diff --git a/generalvclamp.py b/generalvclamp.py new file mode 100644 index 0000000..4a8cd66 --- /dev/null +++ b/generalvclamp.py @@ -0,0 +1,498 @@ +#!/usr/bin/env python + +''' +generalvclamp.py + +Plugin regarding general velocity clamp measurements +''' + +from libhooke import WX_GOOD, ClickedPoint +import wxversion +wxversion.select(WX_GOOD) +from wx import PostEvent +import numpy as np +import scipy as sp +import copy +import os.path +import time + +import warnings +warnings.simplefilter('ignore',np.RankWarning) + + +class generalvclampCommands: + + def _plug_init(self): + self.basecurrent=None + self.basepoints=None + self.autofile='' + + def do_distance(self,args): + ''' + DISTANCE + (generalvclamp.py) + Measure the distance (in nm) between two points. + For a standard experiment this is the delta X distance. + For a force clamp experiment this is the delta Y distance (actually becomes + an alias of zpiezo) + ----------------- + Syntax: distance + ''' + if self.current.curve.experiment == 'clamp': + print 'You wanted to use zpiezo perhaps?' + return + else: + dx,unitx,dy,unity=self._delta(set=1) + print str(dx*(10**9))+' nm' + + + def do_force(self,args): + ''' + FORCE + (generalvclamp.py) + Measure the force difference (in pN) between two points + --------------- + Syntax: force + ''' + if self.current.curve.experiment == 'clamp': + print 'This command makes no sense for a force clamp experiment.' + return + dx,unitx,dy,unity=self._delta(set=1) + print str(dy*(10**12))+' pN' + + + def do_forcebase(self,args): + ''' + FORCEBASE + (generalvclamp.py) + Measures the difference in force (in pN) between a point and a baseline + took as the average between two points. + + The baseline is fixed once for a given curve and different force measurements, + unless the user wants it to be recalculated + ------------ + Syntax: forcebase [rebase] + rebase: Forces forcebase to ask again the baseline + max: Instead of asking for a point to measure, asks for two points and use + the maximum peak in between + ''' + rebase=False #if true=we select rebase + maxpoint=False #if true=we measure the maximum peak + + plot=self._get_displayed_plot() + whatset=1 #fixme: for all sets + if 'rebase' in args or (self.basecurrent != self.current.path): + rebase=True + if 'max' in args: + maxpoint=True + + if rebase: + print 'Select baseline' + self.basepoints=self._measure_N_points(N=2, whatset=whatset) + self.basecurrent=self.current.path + + if maxpoint: + print 'Select two points' + points=self._measure_N_points(N=2, whatset=whatset) + boundpoints=[points[0].index, points[1].index] + boundpoints.sort() + try: + y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]]) + except ValueError: + print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?' + else: + print 'Select point to measure' + points=self._measure_N_points(N=1, whatset=whatset) + #whatplot=points[0].dest + y=points[0].graph_coords[1] + + #fixme: code duplication + boundaries=[self.basepoints[0].index, self.basepoints[1].index] + boundaries.sort() + to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average + + avg=np.mean(to_average) + forcebase=abs(y-avg) + print str(forcebase*(10**12))+' pN' + + + def plotmanip_flatten(self, plot, current, customvalue=False): + ''' + Subtracts a polynomial fit to the non-contact part of the curve, as to flatten it. + the best polynomial fit is chosen among polynomials of degree 1 to n, where n is + given by the configuration file or by the customvalue. + + customvalue= int (>0) --> starts the function even if config says no (default=False) + ''' + + #not a smfs curve... + if current.curve.experiment != 'smfs': + return plot + + #only one set is present... + if len(self.plots[0].vectors) != 2: + return plot + + #config is not flatten, and customvalue flag is false too + if (not self.config['flatten']) and (not customvalue): + return plot + + max_exponent=12 + delta_contact=0 + + if customvalue: + max_cycles=customvalue + else: + max_cycles=self.config['flatten'] #Using > 1 usually doesn't help and can give artefacts. However, it could be useful too. + + contact_index=self.find_contact_point() + valn=[[] for item in range(max_exponent)] + yrn=[0.0 for item in range(max_exponent)] + errn=[0.0 for item in range(max_exponent)] + + for i in range(int(max_cycles)): + x_ext=plot.vectors[0][0][contact_index+delta_contact:] + y_ext=plot.vectors[0][1][contact_index+delta_contact:] + x_ret=plot.vectors[1][0][contact_index+delta_contact:] + y_ret=plot.vectors[1][1][contact_index+delta_contact:] + for exponent in range(max_exponent): + try: + valn[exponent]=sp.polyfit(x_ext,y_ext,exponent) + yrn[exponent]=sp.polyval(valn[exponent],x_ret) + errn[exponent]=sp.sqrt(sum((yrn[exponent]-y_ext)**2)/float(len(y_ext))) + except TypeError: + print 'Cannot flatten!' + return plot + + best_exponent=errn.index(min(errn)) + + #extension + ycorr_ext=y_ext-yrn[best_exponent]+y_ext[0] #noncontact part + yjoin_ext=np.array(plot.vectors[0][1][0:contact_index+delta_contact]) #contact part + #retraction + ycorr_ret=y_ret-yrn[best_exponent]+y_ext[0] #noncontact part + yjoin_ret=np.array(plot.vectors[1][1][0:contact_index+delta_contact]) #contact part + + ycorr_ext=np.concatenate((yjoin_ext, ycorr_ext)) + ycorr_ret=np.concatenate((yjoin_ret, ycorr_ret)) + + plot.vectors[0][1]=list(ycorr_ext) + plot.vectors[1][1]=list(ycorr_ret) + + return plot + + #---SLOPE--- + def do_slope(self,args): + ''' + SLOPE + (generalvclamp.py) + Measures the slope of a delimited chunk on the return trace. + The chunk can be delimited either by two manual clicks, or have + a fixed width, given as an argument. + --------------- + Syntax: slope [width] + The facultative [width] parameter specifies how many + points will be considered for the fit. If [width] is + specified, only one click will be required. + (c) Marco Brucale, Massimo Sandal 2008 + ''' + + # Reads the facultative width argument + try: + fitspan=int(args) + except: + fitspan=0 + + # Decides between the two forms of user input, as per (args) + if fitspan == 0: + # Gets the Xs of two clicked points as indexes on the current curve vector + print 'Click twice to delimit chunk' + clickedpoints=[] + points=self._measure_N_points(N=2,whatset=1) + clickedpoints=[points[0].index,points[1].index] + clickedpoints.sort() + else: + print 'Click once on the leftmost point of the chunk (i.e.usually the peak)' + clickedpoints=[] + points=self._measure_N_points(N=1,whatset=1) + clickedpoints=[points[0].index-fitspan,points[0].index] + + # Calls the function linefit_between + parameters=[0,0,[],[]] + parameters=self.linefit_between(clickedpoints[0],clickedpoints[1]) + + # Outputs the relevant slope parameter + print 'Slope:' + print str(parameters[0]) + + # Makes a vector with the fitted parameters and sends it to the GUI + xtoplot=parameters[2] + ytoplot=[] + x=0 + for x in xtoplot: + ytoplot.append((x*parameters[0])+parameters[1]) + + clickvector_x, clickvector_y=[], [] + for item in points: + clickvector_x.append(item.graph_coords[0]) + clickvector_y.append(item.graph_coords[1]) + + lineplot=self._get_displayed_plot(0) #get topmost displayed plot + + lineplot.add_set(xtoplot,ytoplot) + lineplot.add_set(clickvector_x, clickvector_y) + + if lineplot.styles==[]: + lineplot.styles=[None,None,None,'scatter'] + else: + lineplot.styles+=[None,'scatter'] + + self._send_plot([lineplot]) + + def linefit_between(self,index1,index2): + ''' + Creates two vectors (xtofit,ytofit) slicing out from the + current return trace a portion delimited by the two indexes + given as arguments. + Then does a least squares linear fit on that slice. + Finally returns [0]=the slope, [1]=the intercept of the + fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors + used for the fit. + (c) Marco Brucale, Massimo Sandal 2008 + ''' + # Translates the indexes into two vectors containing the x,y data to fit + xtofit=self.plots[0].vectors[1][0][index1:index2] + ytofit=self.plots[0].vectors[1][1][index1:index2] + + # Does the actual linear fitting (simple least squares with numpy.polyfit) + linefit=[] + linefit=np.polyfit(xtofit,ytofit,1) + + return (linefit[0],linefit[1],xtofit,ytofit) + +#==================== +#AUTOMATIC ANALYSES +#==================== + def do_autopeak(self,args): + ''' + AUTOPEAK + (generalvclamp.py) + Automatically performs a number of analyses on the peaks of the given curve. + Currently it automatically: + - fits peaks with WLC function + - measures peak maximum forces with a baseline + - measures slope in proximity of peak maximum + Requires flatten plotmanipulator , fit.py plugin , flatfilts.py plugin with convfilt + + Syntax: + autopeak [rebase] [pl=value] [t=value] [noauto] [reclick] + + rebase : Re-asks baseline interval + + pl=[value] : Use a fixed persistent length for the fit. If pl is not given, + the fit will be a 2-variable + fit. DO NOT put spaces between 'pl', '=' and the value. + The value must be in meters. + Scientific notation like 0.35e-9 is fine. + + t=[value] : Use a user-defined temperature. The value must be in + kelvins; by default it is 293 K. + DO NOT put spaces between 't', '=' and the value. + + noauto : allows for clicking the contact point by + hand (otherwise it is automatically estimated) the first time. + If subsequent measurements are made, the same contact point + clicked the first time is used + + reclick : redefines by hand the contact point, if noauto has been used before + but the user is unsatisfied of the previously choosen contact point. + + When you first issue the command, it will ask for the filename. If you are giving the filename + of an existing file, autopeak will resume it and append measurements to it. If you are giving + a new filename, it will create the file and append to it until you close Hooke. + + Useful variables (to set with SET command): + + temperature= temperature of the system for wlc fit (in K) + auto_fit_points = number of points to fit before the peak maximum, for wlc + auto_slope_span = number of points on which measure the slope, for slope + ''' + + pl_value=None + T=self.config['temperature'] + + fit_points=int(self.config['auto_fit_points']) + slope_span=int(self.config['auto_slope_span']) + + delta_force=10 + rebase=False #if true=we select rebase + + #Pick up plot + displayed_plot=self._get_displayed_plot(0) + + if self.current.curve.experiment != 'smfs' or len(displayed_plot.vectors) < 2: + print 'Cannot work on this curve.' + return + + #Look for custom persistent length / custom temperature + for arg in args.split(): + #look for a persistent length argument. + if 'pl=' in arg: + pl_expression=arg.split('=') + pl_value=float(pl_expression[1]) #actual value + #look for a T argument. FIXME: spaces are not allowed between 'pl' and value + if ('t=' in arg[0:2]) or ('T=' in arg[0:2]): + t_expression=arg.split('=') + T=float(t_expression[1]) + + #Handle contact point arguments + #FIXME: code duplication + if 'reclick' in args.split(): + print 'Click contact point' + contact_point=self._measure_N_points(N=1, whatset=1)[0] + contact_point_index=contact_point.index + self.wlccontact_point=contact_point + self.wlccontact_index=contact_point.index + self.wlccurrent=self.current.path + elif 'noauto' in args.split(): + if self.wlccontact_index==None or self.wlccurrent != self.current.path: + print 'Click contact point' + contact_point=self._measure_N_points(N=1, whatset=1)[0] + contact_point_index=contact_point.index + self.wlccontact_point=contact_point + self.wlccontact_index=contact_point.index + self.wlccurrent=self.current.path + else: + contact_point=self.wlccontact_point + contact_point_index=self.wlccontact_index + else: + #Automatically find contact point + cindex=self.find_contact_point() + contact_point=ClickedPoint() + contact_point.absolute_coords=displayed_plot.vectors[1][0][cindex], displayed_plot.vectors[1][1][cindex] + contact_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1]) + contact_point.is_marker=True + + + #Pick up force baseline + whatset=1 #fixme: for all sets + if 'rebase' in args or (self.basecurrent != self.current.path): + rebase=True + if rebase: + print 'Select baseline' + self.basepoints=self._measure_N_points(N=2, whatset=whatset) + self.basecurrent=self.current.path + boundaries=[self.basepoints[0].index, self.basepoints[1].index] + boundaries.sort() + to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average + avg=np.mean(to_average) + + #Find peaks. + defplot=self.current.curve.default_plots()[0] + flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip + defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks + peak_location,peak_size=self.has_peaks(defplot, self.convfilt_config['mindeviation']) + + #Create a new plot to send + fitplot=copy.deepcopy(displayed_plot) + + #Initialize data vectors + c_lengths=[] + p_lengths=[] + forces=[] + slopes=[] + + #Cycle between peaks and do analysis. + for peak in peak_location: + #Do WLC fits. + #FIXME: clean wlc fitting, to avoid this clickedpoint hell + #-create a clicked point for the peak point + peak_point=ClickedPoint() + peak_point.absolute_coords=displayed_plot.vectors[1][0][peak], displayed_plot.vectors[1][1][peak] + peak_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1]) + #-create a clicked point for the other fit point + other_fit_point=ClickedPoint() + other_fit_point.absolute_coords=displayed_plot.vectors[1][0][peak-fit_points], displayed_plot.vectors[1][1][peak-fit_points] + other_fit_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1]) + #do the fit + points=[contact_point, peak_point, other_fit_point] + params, yfit, xfit = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T) + #save wlc values (nm) + c_lengths.append(params[0]*(1.0e+9)) + if len(params)==2: #if we did choose 2-value fit + p_lengths.append(params[1]*(1.0e+9)) + else: + p_lengths.append(pl_value) + #Add WLC fit lines to plot + fitplot.add_set(xfit,yfit) + if len(fitplot.styles)==0: + fitplot.styles=[] + else: + fitplot.styles.append(None) + + #Measure forces + delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force] + y=min(delta_to_measure) + #save force values (pN) + forces.append(abs(y-avg)*(1.0e+12)) + + #Measure slopes + slope=self.linefit_between(peak-slope_span,peak)[0] + slopes.append(slope) + + #Show wlc fits and peak locations + self._send_plot([fitplot]) + self.do_peaks('') + + #Ask the user what peaks to ignore from analysis. + print 'Peaks to ignore (0,1...n from contact point,return to take all)' + print 'N to discard measurement' + exclude_raw=raw_input('Input:') + if exclude_raw=='N': + print 'Discarded.' + return + if not exclude_raw=='': + exclude=exclude_raw.split(',') + try: + exclude=[int(item) for item in exclude] + for i in exclude: + c_lengths[i]=None + p_lengths[i]=None + forces[i]=None + slopes[i]=None + except: + print 'Bad input, taking all...' + #Clean data vectors from ignored peaks + c_lengths=[item for item in c_lengths if item != None] + p_lengths=[item for item in p_lengths if item != None] + forces=[item for item in forces if item != None] + slopes=[item for item in slopes if item != None] + print 'contour (nm)',c_lengths + print 'p (nm)',p_lengths + print 'forces (pN)',forces + print 'slopes (N/m)',slopes + + #Save file info + if self.autofile=='': + self.autofile=raw_input('Filename? (return to ignore) ') + if self.autofile=='': + print 'Not saved.' + return + + if not os.path.exists(self.autofile): + f=open(self.autofile,'w+') + f.write('Analysis started '+time.asctime()+'\n') + f.write('----------------------------------------\n') + f.write('; Contour length (nm) ; Persistence length (nm) ; Max.Force (pN) ; Slope (N/m) \n') + f.close() + + print 'Saving...' + f=open(self.autofile,'a+') + + f.write(self.current.path+'\n') + for i in range(len(c_lengths)): + f.write(' ; '+str(c_lengths[i])+' ; '+str(p_lengths[i])+' ; '+str(forces[i])+' ; '+str(slopes[i])+'\n') + + f.close() + self.do_note('autopeak') + \ No newline at end of file diff --git a/hemingclamp.py b/hemingclamp.py new file mode 100755 index 0000000..bc2bc1d --- /dev/null +++ b/hemingclamp.py @@ -0,0 +1,123 @@ +#!/usr/bin/env python + +''' +libhemingclamp.py + +Library for interpreting Hemingway force spectroscopy files. + +Copyright (C) 2006 Massimo Sandal (University of Bologna, Italy) + +This program is released under the GNU General Public License version 2. +''' +__version__='2007_02_15_devel' + +__changelog__=''' +2007_02_15: fixed time counter with my counter +2007_02_07: Initial implementation +''' +import string +import libhookecurve as lhc + +def hemingclamp_magic(filepath): + ''' + we define our magic heuristic for HemingClamp files + ''' + myfile=file(filepath) + headerlines=myfile.readlines()[0:3] + if headerlines[0][0:10]=='#Hemingway' and headerlines[1][0:19]=='#Experiment: FClamp': + return True + else: + return False + +class DataChunk(list): + '''Dummy class to provide ext and ret methods to the data list. + In this case ext and self can be equal. + ''' + + def ext(self): + return self + + def ret(self): + return self + +class hemingclampDriver(lhc.Driver): + + def __init__(self, filename): + + self.filedata = open(filename,'r') + self.data = self.filedata.readlines()[6:] + self.filedata.close() + + self.filetype = 'hemingclamp' + self.experiment = 'clamp' + + self.filename=filename + + def __del__(self): + self.filedata.close() + + def is_me(self): + ''' + we define our magic heuristic for HemingClamp files + ''' + myfile=file(self.filename) + headerlines=myfile.readlines()[0:3] + myfile.close() + if headerlines[0][0:10]=='#Hemingway' and headerlines[1][0:19]=='#Experiment: FClamp': + return True + else: + return False + + def _getdata_all(self): + time = [] + zpiezo = [] + defl = [] + + for i in self.data: + temp = string.split(i) + #time.append(float(temp[0])*(1.0e-3)) + zpiezo.append(float(temp[2])*(1.0e-9)) + defl.append(float(temp[3])*(1.0e-9)) + + #we rebuild the time counter assuming 1 point = 1 millisecond + c=0.0 + for z in zpiezo: + time.append(c) + c+=(1.0e-3) + + return time,zpiezo,defl + + def time(self): + return DataChunk(self._getdata_all()[0]) + + def zpiezo(self): + return DataChunk(self._getdata_all()[1]) + + def deflection(self): + return DataChunk(self._getdata_all()[2]) + + def close_all(self): + ''' + Explicitly closes all files + ''' + self.filedata.close() + + def default_plots(self): + main_plot=lhc.PlotObject() + defl_plot=lhc.PlotObject() + + time=self.time() + zpiezo=self.zpiezo() + deflection=self.deflection() + + main_plot.vectors=[[time,zpiezo]] + main_plot.units=['seconds','meters'] + main_plot.destination=0 + main_plot.title=self.filename + + defl_plot.vectors=[[time,deflection]] + defl_plot.units=['seconds','Newtons'] + defl_plot.destination=1 + + return [main_plot, defl_plot] + \ No newline at end of file diff --git a/hooke.conf b/hooke.conf new file mode 100755 index 0000000..2fc002d --- /dev/null +++ b/hooke.conf @@ -0,0 +1,51 @@ + + + + + + + + + your_own_directory + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/hooke.jpg b/hooke.jpg new file mode 100755 index 0000000..fca9e88 Binary files /dev/null and b/hooke.jpg differ diff --git a/hooke.py b/hooke.py new file mode 100755 index 0000000..7e29328 --- /dev/null +++ b/hooke.py @@ -0,0 +1,752 @@ +#!/usr/bin/env python + +''' +HOOKE - A force spectroscopy review & analysis tool + +(C) 2008 Massimo Sandal + +Copyright (C) 2008 Massimo Sandal (University of Bologna, Italy). + +This program is released under the GNU General Public License version 2. +''' + +from libhooke import HOOKE_VERSION +from libhooke import WX_GOOD + +import os + +import wxversion +wxversion.select(WX_GOOD) +import wx +import wxmpl +from wx.lib.newevent import NewEvent + +import matplotlib.numerix as nx +import scipy as sp + +from threading import * +import Queue + +from hooke_cli import HookeCli +from libhooke import * +import libhookecurve as lhc + +#import file versions, just to know with what we're working... +from hooke_cli import __version__ as hookecli_version + +global __version__ +global events_from_gui +global config +global CLI_PLUGINS +global GUI_PLUGINS +global LOADED_PLUGINS +global PLOTMANIP_PLUGINS +global FILE_DRIVERS + +__version__=HOOKE_VERSION[0] +__release_name__=HOOKE_VERSION[1] + +events_from_gui=Queue.Queue() #GUI ---> CLI COMMUNICATION + +print 'Starting Hooke.' +#CONFIGURATION FILE PARSING +config_obj=HookeConfig() +config=config_obj.load_config('hooke.conf') + +#IMPORTING PLUGINS + +CLI_PLUGINS=[] +GUI_PLUGINS=[] +PLOTMANIP_PLUGINS=[] +LOADED_PLUGINS=[] + +plugin_commands_namespaces=[] +plugin_gui_namespaces=[] +for plugin_name in config['plugins']: + try: + plugin=__import__(plugin_name) + try: + eval('CLI_PLUGINS.append(plugin.'+plugin_name+'Commands)') #take Command plugin classes + plugin_commands_namespaces.append(dir(eval('plugin.'+plugin_name+'Commands'))) + except: + pass + try: + eval('GUI_PLUGINS.append(plugin.'+plugin_name+'Gui)') #take Gui plugin classes + plugin_gui_namespaces.append(dir(eval('plugin.'+plugin_name+'Gui'))) + except: + pass + except ImportError: + print 'Cannot find plugin ',plugin_name + else: + LOADED_PLUGINS.append(plugin_name) + print 'Imported plugin ',plugin_name + +#eliminate names common to all namespaces +for i in range(len(plugin_commands_namespaces)): + plugin_commands_namespaces[i]=[item for item in plugin_commands_namespaces[i] if (item != '__doc__' and item != '__module__' and item != '_plug_init')] +#check for conflicts in namespaces between plugins +#FIXME: only in commands now, because I don't have Gui plugins to check +#FIXME: how to check for plugin-defined variables (self.stuff) ?? +plugin_commands_names=[] +whatplugin_defines=[] +plugin_gui_names=[] +for namespace,plugin_name in zip(plugin_commands_namespaces, config['plugins']): + for item in namespace: + if item in plugin_commands_names: + i=plugin_commands_names.index(item) #we exploit the fact index gives the *first* occurrence of a name... + print 'Error. Plugin ',plugin_name,' defines a function already defined by ',whatplugin_defines[i],'!' + print 'This should not happen. Please disable one or both plugins and contact the plugin authors to solve the conflict.' + print 'Hooke cannot continue.' + exit() + else: + plugin_commands_names.append(item) + whatplugin_defines.append(plugin_name) + + +config['loaded_plugins']=LOADED_PLUGINS #FIXME: kludge -this should be global but not in config! +#IMPORTING DRIVERS +#FIXME: code duplication +FILE_DRIVERS=[] +LOADED_DRIVERS=[] +for driver_name in config['drivers']: + try: + driver=__import__(driver_name) + try: + eval('FILE_DRIVERS.append(driver.'+driver_name+'Driver)') + except: + pass + except ImportError: + print 'Cannot find driver ',driver_name + else: + LOADED_DRIVERS.append(driver_name) + print 'Imported driver ',driver_name +config['loaded_drivers']=LOADED_DRIVERS + +#LIST OF CUSTOM WX EVENTS FOR CLI ---> GUI COMMUNICATION +#FIXME: do they need to be here? +list_of_events={} + +plot_graph, EVT_PLOT = NewEvent() +list_of_events['plot_graph']=plot_graph + +plot_contact, EVT_PLOT_CONTACT = NewEvent() +list_of_events['plot_contact']=plot_contact + +measure_points, EVT_MEASURE_POINTS = NewEvent() +list_of_events['measure_points']=measure_points + +export_image, EVT_EXPORT_IMAGE = NewEvent() +list_of_events['export_image']=export_image + +close_plot, EVT_CLOSE_PLOT = NewEvent() +list_of_events['close_plot'] = close_plot + +show_plots, EVT_SHOW_PLOTS = NewEvent() +list_of_events['show_plots'] = show_plots + +get_displayed_plot, EVT_GET_DISPLAYED_PLOT = NewEvent() +list_of_events['get_displayed_plot'] = get_displayed_plot +#------------ + +class CliThread(Thread): + + def __init__(self,frame,list_of_events): + Thread.__init__(self) + + #here we have to put temporary references to pass to the cli object. + self.frame=frame + self.list_of_events=list_of_events + + self.debug=0 #to be used in the future + + def run(self): + print '\n\nThis is Hooke, version',__version__ , __release_name__ + print + print '(c) Massimo Sandal, 2006. Released under the GNU General Public License Version 2' + print 'Hooke is Free software.' + print '----' + print '' + + def make_command_class(*bases): + #FIXME: perhaps redundant + return type(HookeCli)("HookeCliPlugged", bases + (HookeCli,), {}) + cli = make_command_class(*CLI_PLUGINS)(self.frame,self.list_of_events,events_from_gui,config,FILE_DRIVERS) + cli.cmdloop() + +''' +GUI CODE + +FIXME: put it in a separate module in the future? +''' +class MainMenuBar(wx.MenuBar): + ''' + Creates the menu bar + ''' + def __init__(self): + wx.MenuBar.__init__(self) + '''the menu description. the key of the menu is XX&Menu, where XX is a number telling + the order of the menus on the menubar. + &Menu is the Menu text + the corresponding argument is ('&Item', 'itemname'), where &Item is the item text and itemname + the inner reference to use in the self.menu_items dictionary. + + See create_menus() to see how it works + + Note: the mechanism on page 124 of "wxPython in Action" is less awkward, maybe, but I want + binding to be performed later. Perhaps I'm wrong :) + ''' + + self.menu_desc={'00&File':[('&Open playlist','openplaymenu'),('&Exit','exitmenu')], + '01&Edit':[('&Export text...','exporttextmenu'),('&Export image...','exportimagemenu')], + '02&Help':[('&About Hooke','aboutmenu')]} + self.create_menus() + + def create_menus(self): + ''' + Smartish routine to create the menu from the self.menu_desc dictionary + Hope it's a workable solution for the future. + ''' + self.menus=[] #the menu objects to append to the menubar + self.menu_items={} #the single menu items dictionary, to bind to events + + names=self.menu_desc.keys() #we gotta sort, because iterating keys goes in odd order + names.sort() + + for name in names: + self.menus.append(wx.Menu()) + for menu_item in self.menu_desc[name]: + self.menu_items[menu_item[1]]=self.menus[-1].Append(-1, menu_item[0]) + + for menu,name in zip(self.menus,names): + self.Append(menu,name[2:]) + +class MainPanel(wx.Panel): + def __init__(self,parent,id): + + wx.Panel.__init__(self,parent,id) + self.splitter = wx.SplitterWindow(self) + +ID_FRAME=100 +class MainWindow(wx.Frame): + '''we make a frame inheriting wx.Frame and setting up things on the init''' + def __init__(self,parent,id,title): + + #----------------------------- + #WX WIDGETS INITIALIZATION + + wx.Frame.__init__(self,parent,ID_FRAME,title,size=(800,600),style=wx.DEFAULT_FRAME_STYLE|wx.NO_FULL_REPAINT_ON_RESIZE) + + self.mainpanel=MainPanel(self,-1) + self.cpanels=[] + + self.cpanels.append(wx.Panel(self.mainpanel.splitter,-1)) + self.cpanels.append(wx.Panel(self.mainpanel.splitter,-1)) + + self.statusbar=wx.StatusBar(self,-1) + self.SetStatusBar(self.statusbar) + + self.mainmenubar=MainMenuBar() + self.SetMenuBar(self.mainmenubar) + + self.controls=[] + self.figures=[] + self.axes=[] + + #This is our matplotlib plot + self.controls.append(wxmpl.PlotPanel(self.cpanels[0],-1)) + self.controls.append(wxmpl.PlotPanel(self.cpanels[1],-1)) + #These are our figure and axes, so to have easy references + #Also, we initialize + self.figures=[control.get_figure() for control in self.controls] + self.axes=[figure.gca() for figure in self.figures] + + self.cpanels[1].Hide() + self.mainpanel.splitter.Initialize(self.cpanels[0]) + + self.sizer_dance() #place/size the widgets + + self.controls[0].SetSize(self.cpanels[0].GetSize()) + self.controls[1].SetSize(self.cpanels[1].GetSize()) + + #------------------------------------------- + #NON-WX WIDGETS INITIALIZATION + + #Flags. + self.click_plot=0 + + #FIXME: These could become a single flag with different (string?) values + self.on_measure_distance=False + self.on_measure_force=False + + self.plot_fit=False + + #Number of points to be clicked + self.num_of_points = 2 + + #Data. + self.current_x_ext=[[],[]] + self.current_y_ext=[[],[]] + self.current_x_ret=[[],[]] + self.current_y_ret=[[],[]] + + self.current_x_unit=[None,None] + self.current_y_unit=[None,None] + + #Initialize xaxes, yaxes + #FIXME: should come from config + self.current_xaxes=0 + self.current_yaxes=0 + + #Other + + + self.index_buffer=[] + + self.clicked_points=[] + + self.events_from_gui = events_from_gui + + ''' + This dictionary keeps all the flags and the relative functon names that + have to be called when a point is clicked. + That is: + - if point is clicked AND foo_flag=True + - foo() + + Conversely, foo_flag is True if a corresponding event is launched by the CLI. + + self.ClickedPoints() takes care of handling this + ''' + + self.click_flags_functions={'measure_points':[False, 'MeasurePoints']} + + #Custom events from CLI --> GUI functions! + #FIXME: Should use the self.Bind() syntax + EVT_PLOT(self, self.PlotCurve) + EVT_PLOT_CONTACT(self, self.PlotContact) + EVT_GET_DISPLAYED_PLOT(self, self.OnGetDisplayedPlot) + EVT_MEASURE_POINTS(self, self.OnMeasurePoints) + EVT_EXPORT_IMAGE(self,self.ExportImage) + EVT_CLOSE_PLOT(self, self.OnClosePlot) + EVT_SHOW_PLOTS(self, self.OnShowPlots) + + #This event and control decide what happens when I click on the plot 0. + wxmpl.EVT_POINT(self, self.controls[0].GetId(), self.ClickPoint0) + wxmpl.EVT_POINT(self, self.controls[1].GetId(), self.ClickPoint1) + + #RUN PLUGIN-SPECIFIC INITIALIZATION + #make sure we execute _plug_init() for every command line plugin we import + for plugin_name in config['plugins']: + try: + plugin=__import__(plugin_name) + try: + eval('plugin.'+plugin_name+'Gui._plug_init(self)') + pass + except AttributeError: + pass + except ImportError: + pass + + + + #WX-SPECIFIC FUNCTIONS + def sizer_dance(self): + ''' + adjust size and placement of wxpython widgets. + ''' + self.splittersizer = wx.BoxSizer(wx.VERTICAL) + self.splittersizer.Add(self.mainpanel.splitter, 1, wx.EXPAND) + + self.plot1sizer = wx.BoxSizer() + self.plot1sizer.Add(self.controls[0], 1, wx.EXPAND) + + self.plot2sizer = wx.BoxSizer() + self.plot2sizer.Add(self.controls[1], 1, wx.EXPAND) + + self.panelsizer=wx.BoxSizer() + self.panelsizer.Add(self.mainpanel, -1, wx.EXPAND) + + self.cpanels[0].SetSizer(self.plot1sizer) + self.cpanels[1].SetSizer(self.plot2sizer) + + self.mainpanel.SetSizer(self.splittersizer) + self.SetSizer(self.panelsizer) + + def binding_dance(self): + self.Bind(wx.EVT_MENU, self.OnOpenPlayMenu, self.menubar.menu_items['openplaymenu']) + self.Bind(wx.EVT_MENU, self.OnExitMenu, self.menubar.menu_items['exitmenu']) + self.Bind(wx.EVT_MENU, self.OnExportText, self.menubar.menu_items['exporttextmenu']) + self.Bind(wx.EVT_MENU, self.OnExportImage, self.menubar.menu_items['exportimagemenu']) + self.Bind(wx.EVT_MENU, self.OnAboutMenu, self.menubar.menu_items['aboutmenu']) + + # DOUBLE PLOT MANAGEMENT + #---------------------- + def show_both(self): + ''' + Shows both plots. + ''' + self.mainpanel.splitter.SplitHorizontally(self.cpanels[0],self.cpanels[1]) + self.mainpanel.splitter.SetSashGravity(0.5) + self.mainpanel.splitter.SetSashPosition(300) #FIXME: we should get it and restore it + self.mainpanel.splitter.UpdateSize() + + def close_plot(self,plot): + ''' + Closes one plot - only if it's open + ''' + if not self.cpanels[plot].IsShown(): + return + if plot != 0: + self.current_plot_dest = 0 + else: + self.current_plot_dest = 1 + self.cpanels[plot].Hide() + self.mainpanel.splitter.Unsplit(self.cpanels[plot]) + self.mainpanel.splitter.UpdateSize() + + + def OnClosePlot(self,event): + self.close_plot(event.to_close) + + def OnShowPlots(self,event): + self.show_both() + + + #FILE MENU FUNCTIONS + #-------------------- + def OnOpenPlayMenu(self, event): + pass + + def OnExitMenu(self,event): + pass + + def OnExportText(self,event): + pass + + def OnExportImage(self,event): + pass + + def OnAboutMenu(self,event): + pass + + #PLOT INTERACTION + #---------------- + def PlotCurve(self,event): + ''' + plots the current ext,ret curve. + ''' + dest=0 + + #FIXME: BAD kludge following. There should be a well made plot queue mechanism, with replacements etc. + #--- + #If we have only one plot in the event, we already have one in self.plots and this is a secondary plot, + #do not erase self.plots but append the new plot to it. + if len(event.plots) == 1 and event.plots[0].destination != 0 and len(self.plots) == 1: + self.plots.append(event.plots[0]) + #if we already have two plots and a new secondary plot comes, we substitute the previous + if len(event.plots) == 1 and event.plots[0].destination != 0 and len(self.plots) > 1: + self.plots[1] = event.plots[0] + else: + self.plots = event.plots + + #FIXME. Should be in PlotObject, somehow + c=0 + for plot in self.plots: + if self.plots[c].styles==[]: + self.plots[c].styles=[None for item in plot.vectors] + + for plot in self.plots: + ''' + MAIN LOOP FOR ALL PLOTS (now only 2 are allowed but...) + ''' + if 'destination' in dir(plot): + dest=plot.destination + + #if the requested panel is not shown, show it + if not ( self.cpanels[dest].IsShown() ): + self.show_both() + + self.axes[dest].hold(False) + self.current_vectors=plot.vectors + self.current_title=plot.title + self.current_plot_dest=dest #let's try this way to take into account the destination plot... + + c=0 + for vectors_to_plot in self.current_vectors: + if len(vectors_to_plot)==2: #3d plots are to come... + if len(plot.styles) > 0 and plot.styles[c] == 'scatter': + self.axes[dest].scatter(vectors_to_plot[0],vectors_to_plot[1]) + else: + self.axes[dest].plot(vectors_to_plot[0],vectors_to_plot[1]) + + self.axes[dest].hold(True) + c+=1 + else: + pass + + #FIXME: tackles only 2d plots + self.axes[dest].set_xlabel(plot.units[0]) + self.axes[dest].set_ylabel(plot.units[1]) + + #FIXME: set smaller fonts + self.axes[dest].set_title(plot.title) + + if plot.xaxes: + #swap X axis + xlim=self.axes[dest].get_xlim() + self.axes[dest].set_xlim((xlim[1],xlim[0])) + if plot.yaxes: + #swap Y axis + ylim=self.axes[dest].get_ylim() + self.axes[dest].set_ylim((ylim[1],ylim[0])) + + self.controls[dest].draw() + + + def PlotContact(self,event): + ''' + plots the contact point + ''' + self.axes[0].hold(True) + self.current_contact_index=event.contact_index + + #now we fake a clicked point + self.clicked_points.append(ClickedPoint()) + self.clicked_points[-1].absolute_coords=self.current_x_ret[dest][self.current_contact_index], self.current_y_ret[dest][self.current_contact_index] + self.clicked_points[-1].is_marker=True + + self._replot() + self.clicked_points=[] + + + def ClickPoint0(self,event): + self.current_plot_dest=0 + self.ClickPoint(event) + def ClickPoint1(self,event): + self.current_plot_dest=1 + self.ClickPoint(event) + + def ClickPoint(self,event): + ''' + this function decides what to do when we receive a left click on the axes. + We trigger other functions: + - the action chosen by the CLI sends an event + - the event raises a flag : self.click_flags_functions['foo'][0] + - the raised flag wants the function in self.click_flags_functions[1] to be called after a click + ''' + for key, value in self.click_flags_functions.items(): + if value[0]: + eval('self.'+value[1]+'(event)') + + def OnMeasurePoints(self,event): + ''' + trigger flags to measure N points + ''' + self.click_flags_functions['measure_points'][0]=True + if 'num_of_points' in dir(event): + self.num_of_points=event.num_of_points + + + def MeasurePoints(self,event,current_set=1): + dest=self.current_plot_dest + try: + current_set=event.set + except AttributeError: + pass + + #find the current plot matching the clicked destination + plot=self._plot_of_dest() + if len(plot.vectors)-1 < current_set: #what happens if current_set is 1 and we have only 1 vector? + current_set=current_set-len(plot.vectors) + + xvector=plot.vectors[current_set][0] + yvector=plot.vectors[current_set][1] + + self.clicked_points.append(ClickedPoint()) + self.clicked_points[-1].absolute_coords=event.xdata, event.ydata + self.clicked_points[-1].find_graph_coords(xvector,yvector) + self.clicked_points[-1].is_marker=True + self.clicked_points[-1].is_line_edge=True + self.clicked_points[-1].dest=dest + + self._replot() + + if len(self.clicked_points)==self.num_of_points: + self.events_from_gui.put(self.clicked_points) + #restore to default state: + self.clicked_points=[] + self.click_flags_functions['measure_points'][0]=False + + + def OnGetDisplayedPlot(self,event): + if 'dest' in dir(event): + self.GetDisplayedPlot(event.dest) + else: + self.GetDisplayedPlot(self.current_plot_dest) + + def GetDisplayedPlot(self,dest): + ''' + returns to the CLI the currently displayed plot for the given destination + ''' + displayed_plot=self._plot_of_dest(dest) + events_from_gui.put(displayed_plot) + + def ExportImage(self,event): + ''' + exports an image as a file. + Current supported file formats: png, eps + (matplotlib docs say that jpeg should be supported too, but with .jpg it doesn't work for me!) + ''' + #dest=self.current_plot_dest + dest=event.dest + filename=event.name + self.figures[dest].savefig(filename) + + def _find_nearest_point(self, mypoint, dataset=1): + ''' + Given a clicked point on the plot, finds the nearest point in the dataset (in X) that + corresponds to the clicked point. + ''' + dest=self.current_plot_dest + + xvector=plot.vectors[dataset][0] + yvector=plot.vectors[dataset][1] + + #Ye Olde sorting algorithm... + #FIXME: is there a better solution? + index=0 + best_index=0 + best_diff=10^9 #hope we never go over this magic number :( + for point in xvector: + diff=abs(point-mypoint) + if diff 0 and plot.styles[c]=='scatter': + self.axes[dest].scatter(plotset[0], plotset[1]) + else: + self.axes[dest].plot(plotset[0], plotset[1]) + c+=1 + #plot points we have clicked + for item in self.clicked_points: + if item.is_marker: + if item.graph_coords==(None,None): #if we have no graph coords, we display absolute coords + self.axes[dest].scatter([item.absolute_coords[0]],[item.absolute_coords[1]]) + else: + self.axes[dest].scatter([item.graph_coords[0]],[item.graph_coords[1]]) + + if self.plot_fit: + print 'DEBUGGING WARNING: use of self.plot_fit is deprecated!' + self.axes[dest].plot(self.plot_fit[0],self.plot_fit[1]) + + self.axes[dest].hold(True) + #set old axes again + self.axes[dest].set_xlim(xlim) + self.axes[dest].set_ylim(ylim) + #set title and names again... + self.axes[dest].set_title(self.current_title) + self.axes[dest].set_xlabel(plot.units[0]) + self.axes[dest].set_ylabel(plot.units[1]) + #and redraw! + self.controls[dest].draw() + + +class MySplashScreen(wx.SplashScreen): + """ + Create a splash screen widget. + That's just a fancy addition... every serious application has a splash screen! + """ + def __init__(self, frame): + # This is a recipe to a the screen. + # Modify the following variables as necessary. + #aBitmap = wx.Image(name = "wxPyWiki.jpg").ConvertToBitmap() + aBitmap=wx.Image(name='hooke.jpg').ConvertToBitmap() + splashStyle = wx.SPLASH_CENTRE_ON_SCREEN | wx.SPLASH_TIMEOUT + splashDuration = 2000 # milliseconds + splashCallback = None + # Call the constructor with the above arguments in exactly the + # following order. + wx.SplashScreen.__init__(self, aBitmap, splashStyle, + splashDuration, None, -1) + wx.EVT_CLOSE(self, self.OnExit) + self.frame=frame + wx.Yield() + + def OnExit(self, evt): + self.Hide() + + self.frame.Show() + # The program will freeze without this line. + evt.Skip() # Make sure the default handler runs too... + + +#------------------------------------------------------------------------------ + +def main(): + + #save the directory where Hooke is located + config['hookedir']=os.getcwd() + + #now change to the working directory. + try: + os.chdir(config['workdir']) + except OSError: + print "Warning: Invalid work directory." + + app=wx.PySimpleApp() + + def make_gui_class(*bases): + return type(MainWindow)("MainWindowPlugged", bases + (MainWindow,), {}) + + main_frame = make_gui_class(*GUI_PLUGINS)(None, -1, ('Hooke '+__version__)) + + #FIXME. The frame.Show() is called by the splashscreen here! Ugly as hell. + + mysplash=MySplashScreen(main_frame) + mysplash.Show() + + my_cmdline=CliThread(main_frame, list_of_events) + my_cmdline.start() + + + app.MainLoop() + +main() diff --git a/hooke_cli.py b/hooke_cli.py new file mode 100755 index 0000000..8ed2ca0 --- /dev/null +++ b/hooke_cli.py @@ -0,0 +1,871 @@ +#!/usr/bin/env python + +''' +hooke_cli.py + +Command line module of Hooke. + +Copyright (C) 2006 Massimo Sandal (University of Bologna, Italy). + +This program is released under the GNU General Public License version 2. +''' + + +from libhooke import * #FIXME +import libhookecurve as lhc + +from libhooke import WX_GOOD +from libhooke import HOOKE_VERSION + +import wxversion +wxversion.select(WX_GOOD) +import wx + +from wx.lib.newevent import NewEvent +from matplotlib.numerix import * #FIXME + +import xml.dom.minidom +import sys, os, os.path, glob, shutil +import Queue +import cmd +import time + +global __version__ +global __codename__ +global __releasedate__ +__version__ = HOOKE_VERSION[0] +__codename__ = HOOKE_VERSION[1] +__releasedate__ = HOOKE_VERSION[2] + +from matplotlib import __version__ as mpl_version +from wx import __version__ as wx_version +from wxmpl import __version__ as wxmpl_version +from scipy import __version__ as scipy_version +from numpy import __version__ as numpy_version +from sys import version as python_version +import platform + + +class HookeCli(cmd.Cmd): + + def __init__(self,frame,list_of_events,events_from_gui,config,drivers): + cmd.Cmd.__init__(self) + + self.prompt = 'hooke: ' + + + self.current_list=[] #the playlist we're using + + self.current=None #the current curve under analysis. + self.plots=None + ''' + The actual hierarchy of the "current curve" is a bit complex: + + self.current = the lhc.HookeCurve container object of the current curve + self.current.curve = the current "real" curve object as defined in the filetype driver class + self.current.curve.default_plots() = the default plots of the filetype driver. + + The plot objects obtained by mean of self.current.curve.default_plots() + then undergoes modifications by the plotmanip + modifier functions. The modified plot is saved in self.plots and used if needed by other functions. + ''' + + + self.pointer=0 #a pointer to navigate the current list + + #Things that come from outside + self.frame=frame #the wx frame we refer to + self.list_of_events=list_of_events #a list of wx events we use to interact with the GUI + self.events_from_gui=events_from_gui #the Queue object we use to have messages from the GUI + self.config=config #the configuration dictionary + self.drivers=drivers #the file format drivers + + #get plot manipulation functions + plotmanip_functions=[] + for object_name in dir(self): + if object_name[0:9]=='plotmanip': + plotmanip_functions.append(getattr(self,object_name)) + #put plotmanips in order + self.plotmanip=[None for item in self.config['plotmanips']] + for item in plotmanip_functions: + namefunction=item.__name__[10:] + if namefunction in self.config['plotmanips']: + nameindex=self.config['plotmanips'].index(namefunction) #index of function in plotmanips config + self.plotmanip[nameindex] = item + else: + pass + + + self.playlist_saved=0 + self.playlist_name='' + self.notes_saved=1 + + #Data that must be saved in the playlist, related to the whole playlist (not individual curves) + self.playlist_generics={} + + #make sure we execute _plug_init() for every command line plugin we import + for plugin_name in self.config['plugins']: + try: + plugin=__import__(plugin_name) + try: + eval('plugin.'+plugin_name+'Commands._plug_init(self)') + except AttributeError: + pass + except ImportError: + pass + + +#HELPER FUNCTIONS +#Everything sending an event should be here + def _measure_N_points(self, N, whatset=1): + ''' + general helper function for N-points measures + ''' + measure_points=self.list_of_events['measure_points'] + wx.PostEvent(self.frame, measure_points(num_of_points=N, set=whatset)) + while 1: + try: + points=self.frame.events_from_gui.get() + break + except Empty: + pass + return points + + def _get_displayed_plot(self,dest=0): + ''' + returns the currently displayed plot. + ''' + wx.PostEvent(self.frame, self.list_of_events['get_displayed_plot'](dest=dest)) + while 1: + try: + displayed_plot=self.events_from_gui.get() + except Empty: + pass + if displayed_plot: + break + return displayed_plot + + def _send_plot(self,plots): + ''' + sends a plot to the GUI + ''' + wx.PostEvent(self.frame, self.list_of_events['plot_graph'](plots=plots)) + return + + def _find_plotmanip(self, name): + ''' + returns a plot manipulator function from its name + ''' + return self.plotmanip[self.config['plotmanips'].index(name)] + +#HERE COMMANDS BEGIN + + def help_set(self): + print ''' +SET +Sets a local configuration variable +------------- +Syntax: set [variable] [value] + ''' + def do_set(self,args): + #FIXME: some variables in self.config should be hidden or intelligently configurated... + args=args.split() + if len(args)==0: + print 'You must specify a variable and a value' + print 'Available variables:' + print self.config.keys() + return + if args[0] not in self.config.keys(): + print 'This is not an internal Hooke variable!' + return + if len(args)==1: + #FIXME:we should reload the config file and reset the config value + print self.config[args[0]] + return + key=args[0] + try: #try to have a numeric value + value=float(args[1]) + except ValueError: #if it cannot be converted to float, it's None, or a string... + if value.lower()=='none': + value=None + else: + value=args[1] + + self.config[key]=value + self.do_plot(0) + +#PLAYLIST MANAGEMENT AND NAVIGATION +#------------------------------------ + + def help_loadlist(self): + print ''' +LOADLIST +Loads a file playlist +----------- +Syntax: loadlist [playlist file] + ''' + def do_loadlist(self, args): + #checking for args: if nothing is given as input, we warn and exit. + while len(args)==0: + args=raw_input('File to load?') + + arglist=args.split() + play_to_load=arglist[0] + + #We assume a Hooke playlist has the extension .hkp + if play_to_load[-4:] != '.hkp': + play_to_load+='.hkp' + + try: + playxml=PlaylistXML() + self.current_list, self.playlist_generics=playxml.load(play_to_load) + self.current_playxml=playxml + except IOError: + print 'File not found.' + return + + print 'Loaded %s curves' %len(self.current_list) + + if 'pointer' in self.playlist_generics.keys(): + self.pointer=int(self.playlist_generics['pointer']) + else: + #if no pointer is found, set the current curve as the first curve of the loaded playlist + self.pointer=0 + print 'Starting at curve ',self.pointer + + self.current=self.current_list[self.pointer] + + #resets saved/notes saved state + self.playlist_saved=0 + self.playlist_name='' + self.notes_saved=0 + + self.do_plot(0) + + + def help_genlist(self): + print ''' +GENLIST +Generates a file playlist. +Note it doesn't *save* it: see savelist for this. + +If [input files] is a directory, it will use all files in the directory for playlist. +So: +genlist dir +genlist dir/ +genlist dir/*.* + +are all equivalent syntax. +------------ +Syntax: genlist [input files] + +''' + def do_genlist(self,args): + #args list is: input path, output name + if len(args)==0: + args=raw_input('Input files?') + + arglist=args.split() + list_path=arglist[0] + + #if it's a directory, is like /directory/*.* + #FIXME: probably a bit kludgy. + if os.path.isdir(list_path): + if platform.system == 'Windows': + SLASH="\\" + else: + SLASH="/" + if list_path[-1] == SLASH: + list_path=list_path+'*.*' + else: + list_path=list_path+SLASH+'*.*' + + #expanding correctly the input list with the glob module :) + list_files=glob.glob(list_path) + list_files.sort() + + self.current_list=[] + for item in list_files: + try: + self.current_list.append(lhc.HookeCurve(os.path.abspath(item))) + except: + pass + + self.pointer=0 + if len(self.current_list)>0: + self.current=self.current_list[self.pointer] + else: + print 'Empty list!' + return + + #resets saved/notes saved state + self.playlist_saved=0 + self.playlist_name='' + self.notes_saved=0 + + self.do_plot(0) + + + def do_savelist(self,args): + ''' + SAVELIST + Saves the current file playlist on disk. + ------------ + Syntax: savelist [filename] + ''' + while len(args)==0: + args=raw_input('Input files?') + + output_filename=args + + self.playlist_generics['pointer']=self.pointer + + #autocomplete filename if not specified + if output_filename[-4:] != '.hkp': + output_filename+='.hkp' + + playxml=PlaylistXML() + playxml.export(self.current_list, self.playlist_generics) + playxml.save(output_filename) + + #remembers we have saved playlist + self.playlist_saved=1 + + def help_addtolist(self): + print ''' +ADDTOLIST +Adds a file to the current playlist +-------------- +Syntax: addtolist [filename] +''' + def do_addtolist(self,args): + #args list is: input path + if len(args)==0: + print 'You must give the input filename you want to add' + self.help_addtolist() + return + + filenames=glob.glob(args) + + for filename in filenames: + self.current_list.append(lhc.HookeCurve(os.path.abspath(filename))) + #we need to save playlist + self.playlist_saved=0 + + def help_printlist(self): + print ''' +PRINTLIST +Prints the list of curves in the current playlist +------------- +Syntax: printlist +''' + def do_printlist(self,args): + for item in self.current_list: + print item.path + + + def help_jump(self): + print ''' +JUMP +Jumps to a given curve. +------ +Syntax: jump {$curve} + +If the curve is not in the current playlist, it politely asks if we want to add it. + ''' + def do_jump(self,filename): + ''' + jumps to the curve with the given filename. + if the filename is not in the playlist, it asks if we must add it or not. + ''' + + if filename=='': + filename=raw_input('Jump to?') + + filepath=os.path.abspath(filename) + print filepath + + c=0 + item_not_found=1 + while item_not_found: + try: + + if self.current_list[c].path == filepath: + self.pointer=c + self.current=self.current_list[self.pointer] + item_not_found=0 + self.do_plot(0) + else: + c+=1 + except IndexError: + #We've found the end of the list. + answer=raw_input('Curve not found in playlist. Add it to list?') + if answer.lower()[0]=='y': + try: + self.do_addtolist(filepath) + except: + print 'Curve file not found.' + return + self.current=self.current_list[-1] + self.pointer=(len(current_list)-1) + self.do_plot(0) + + item_not_found=0 + + + def do_index(self,args): + ''' + INDEX + Prints the index of the current curve in the list + ----- + Syntax: index + ''' + print self.pointer+1, 'of', len(self.current_list) + + + def help_next(self): + print ''' +NEXT +Go the next curve in the playlist. +If we are at the last curve, we come back to the first. +----- +Syntax: next, n + ''' + def do_next(self,args): + self.current.curve.close_all() + if self.pointer == (len(self.current_list)-1): + self.pointer=0 + print 'Playlist finished; back to first curve.' + else: + self.pointer+=1 + + self.current=self.current_list[self.pointer] + self.do_plot(0) + + + def help_n(self): + self.help_next() + def do_n(self,args): + self.do_next(args) + + def help_previous(self,args): + print ''' +PREVIOUS +Go to the previous curve in the playlist. +If we are at the first curve, we jump to the last. +------- +Syntax: previous, p + ''' + def do_previous(self,args): + self.current.curve.close_all() + if self.pointer == 0: + self.pointer=(len(self.current_list)-1) + print 'Start of playlist; jump to last curve.' + else: + self.pointer-=1 + + self.current=self.current_list[self.pointer] + self.do_plot(args) + + + def help_p(self): + self.help_previous() + def do_p(self,args): + self.do_previous(args) + + +#PLOT INTERACTION COMMANDS +#------------------------------- + def help_plot(self): + print ''' +PLOT +Plots the current force curve +------- +Syntax: plot + ''' + def do_plot(self,args): + + self.current.identify(self.drivers) + self.plots=self.current.curve.default_plots() + try: + self.plots=self.current.curve.default_plots() + except Exception, e: + print 'Unexpected error occurred in do_plot().' + print e + return + + #apply the plotmanip functions eventually present + nplots=len(self.plots) + c=0 + while c 1: + dest=int(args[1]) + + export_image=self.list_of_events['export_image'] + wx.PostEvent(self.frame, export_image(name=name, dest=dest)) + + + def help_txt(self): + print ''' +TXT +Saves the current curve as a text file +Columns are, in order: +X1 , Y1 , X2 , Y2 , X3 , Y3 ... + +------------- +Syntax: txt [filename] {plot to export} + ''' + def do_txt(self,args): + + def transposed2(lists, defval=0): + ''' + transposes a list of lists, i.e. from [[a,b,c],[x,y,z]] to [[a,x],[b,y],[c,z]] without losing + elements + (by Zoran Isailovski on the Python Cookbook online) + ''' + if not lists: return [] + return map(lambda *row: [elem or defval for elem in row], *lists) + + whichplot=0 + args=args.split() + if len(args)==0: + filename=raw_input('Filename?') + else: + filename=args[0] + try: + whichplot=int(args[1]) + except: + pass + + columns=[] + for dataset in self.plots[whichplot].vectors: + for i in range(0,len(dataset)): + columns.append([]) + for value in dataset[i]: + columns[-1].append(str(value)) + + rows=transposed2(columns, 'nan') + rows=[' , '.join(item) for item in rows] + text='\n'.join(rows) + + txtfile=open(filename,'w+') + txtfile.write(text) + txtfile.close() + + + #LOGGING, REPORTING, NOTETAKING + def help_note(self): + print ''' +NOTE +Writes or displays a note about the current curve. +If [anything] is empty, it displays the note, otherwise it adds a note. +The note is then saved in the playlist if you issue a savelist command +--------------- +Syntax: note [anything] +''' + def do_note(self,args): + if args=='': + print self.current_list[self.pointer].notes + else: + #bypass UnicodeDecodeError troubles + try: + args=args.decode('ascii') + except: + args=args.decode('ascii','ignore') + if len(args)==0: + args='?' + + self.current_list[self.pointer].notes=args + self.notes_saved=0 + + def help_notelog(self): + print ''' +NOTELOG +Writes a log of the notes taken during the session for the current +playlist +-------------- +Syntax notelog [filename] +''' + def do_notelog(self,args): + + if len(args)==0: + args=raw_input('Notelog filename?') + + note_lines='Notes taken at '+time.asctime()+'\n' + for item in self.current_list: + if len(item.notes)>0: + #FIXME: log should be justified + #FIXME: file path should be truncated... + note_string=(item.path+' | '+item.notes+'\n') + note_lines+=note_string + + try: + f=open(args,'a+') + f.write(note_lines) + f.close() + except IOError, (ErrorNumber, ErrorMessage): + print 'Error: notes cannot be saved. Catched exception:' + print ErrorMessage + + self.notes_saved=1 + + def help_copylog(self): + print ''' +COPYLOG +Moves the annotated curves to another directory +----------- +Syntax copylog [directory] + ''' + def do_copylog(self,args): + + if len(args)==0: + args=raw_input('Destination directory?') + + mydir=os.path.abspath(args) + if not os.path.isdir(mydir): + print 'Destination is not a directory.' + return + + for item in self.current_list: + if len(item.notes)>0: + try: + shutil.copy(item.path, mydir) + except OSError: + print 'OSError. Cannot copy file. Perhaps you gave me a wrong directory?' + + +#OS INTERACTION COMMANDS +#----------------- + def help_dir(self): + print ''' +DIR, LS +Lists the files in the directory +--------- +Syntax: dir [path] + ls [path] + ''' + def do_dir(self,args): + + if len(args)==0: + args='*' + print glob.glob(args) + + def help_ls(self): + self.help_dir(self) + def do_ls(self,args): + self.do_dir(args) + + def help_pwd(self): + print ''' +PWD +Gives the current working directory. +------------ +Syntax: pwd + ''' + def do_pwd(self,args): + print os.getcwd() + + def help_cd(self): + print ''' +CD +Changes the current working directory +----- +Syntax: cd + ''' + def do_cd(self,args): + mypath=os.path.abspath(args) + try: + os.chdir(mypath) + except OSError: + print 'I cannot access that directory.' + + + def help_system(self): + print ''' +SYSTEM +Executes a system command line and reports the output +----- +Syntax system [command line] + ''' + pass + def do_system(self,args): + waste=os.system(args) + + def do_debug(self,args): + ''' + this is a dummy command where I put debugging things + ''' + print self.config['plotmanips'] + pass + + def help_current(self): + print ''' +CURRENT +Prints the current curve path. +------ +Syntax: current + ''' + def do_current(self,args): + print self.current.path + + def do_version(self,args): + ''' + VERSION + ------ + Prints the current version and codename, plus library version. Useful for debugging. + ''' + print 'Hooke '+__version__+' ('+__codename__+')' + print 'Released on: '+__releasedate__ + print '---' + print 'Python version: '+python_version + print 'WxPython version: '+wx_version + print 'wxMPL version: '+wxmpl_version + print 'Matplotlib version: '+mpl_version + print 'SciPy version: '+scipy_version + print 'NumPy version: '+numpy_version + print '---' + print 'Platform: '+str(platform.uname()) + print '---' + print 'Loaded plugins:',self.config['loaded_plugins'] + + def help_exit(self): + print ''' +EXIT, QUIT +Exits the program cleanly. +------ +Syntax: exit +Syntax: quit +''' + def do_exit(self,args): + we_exit='N' + + if (not self.playlist_saved) or (not self.notes_saved): + we_exit=raw_input('You did not save your playlist and/or notes. Exit?') + else: + we_exit=raw_input('Exit?') + + if we_exit[0].upper()=='Y': + wx.CallAfter(self.frame.Close) + sys.exit(0) + else: + return + + def help_quit(self): + self.help_exit() + def do_quit(self,args): + self.do_exit(args) + + + +if __name__ == '__main__': + mycli=HookeCli(0) + mycli.cmdloop() \ No newline at end of file diff --git a/libhooke.py b/libhooke.py new file mode 100755 index 0000000..65b9636 --- /dev/null +++ b/libhooke.py @@ -0,0 +1,317 @@ +#!/usr/bin/env python + +''' +libhooke.py + +General library of internal objects and utilities for Hooke. + +Copyright (C) 2006 Massimo Sandal (University of Bologna, Italy). +With algorithms contributed by Francesco Musiani (University of Bologna, Italy) + +This program is released under the GNU General Public License version 2. +''' + + + +import libhookecurve as lhc + +import scipy +import scipy.signal +import scipy.optimize +import scipy.stats +import numpy +import xml.dom.minidom +import os +import string +import csv + +HOOKE_VERSION=['0.8.3', 'Seinei', '2008-04-16'] +WX_GOOD=['2.6','2.8'] + +class PlaylistXML: + ''' + This module allows for import/export of an XML playlist into/out of a list of HookeCurve objects + ''' + + def __init__(self): + + self.playlist=None #the DOM object representing the playlist data structure + self.playpath=None #the path of the playlist XML file + self.plaything=None + self.hidden_attributes=['curve'] #This list contains hidden attributes that we don't want to go into the playlist. + + def export(self, list_of_hooke_curves, generics): + ''' + Creates an initial playlist from a list of files. + A playlist is an XML document with the following syntaxis: + + + + + ''' + + #create the output playlist, a simple XML document + impl=xml.dom.minidom.getDOMImplementation() + #create the document DOM object and the root element + newdoc=impl.createDocument(None, "playlist",None) + top_element=newdoc.documentElement + + #save generics variables + playlist_generics=newdoc.createElement("generics") + top_element.appendChild(playlist_generics) + for key in generics.keys(): + newdoc.createAttribute(key) + playlist_generics.setAttribute(key,str(generics[key])) + + #save curves and their attributes + for item in list_of_hooke_curves: + #playlist_element=newdoc.createElement("curve") + playlist_element=newdoc.createElement("element") + top_element.appendChild(playlist_element) + for key in item.__dict__: + if not (key in self.hidden_attributes): + newdoc.createAttribute(key) + playlist_element.setAttribute(key,str(item.__dict__[key])) + + self.playlist=newdoc + + def load(self,filename): + ''' + loads a playlist file + ''' + myplay=file(filename) + self.playpath=filename + + #the following 3 lines are needed to strip newlines. otherwise, since newlines + #are XML elements too (why?), the parser would read them (and re-save them, multiplying + #newlines...) + #yes, I'm an XML n00b + the_file=myplay.read() + the_file_lines=the_file.split('\n') + the_file=''.join(the_file_lines) + + self.playlist=xml.dom.minidom.parseString(the_file) + + #inner parsing functions + def handlePlaylist(playlist): + list_of_files=playlist.getElementsByTagName("element") + generics=playlist.getElementsByTagName("generics") + return handleFiles(list_of_files), handleGenerics(generics) + + def handleGenerics(generics): + generics_dict={} + if len(generics)==0: + return generics_dict + + for attribute in generics[0].attributes.keys(): + generics_dict[attribute]=generics[0].getAttribute(attribute) + return generics_dict + + def handleFiles(list_of_files): + new_playlist=[] + for myfile in list_of_files: + #rebuild a data structure from the xml attributes + the_curve=lhc.HookeCurve(myfile.getAttribute('path')) + for attribute in myfile.attributes.keys(): #extract attributes for the single curve + the_curve.__dict__[attribute]=myfile.getAttribute(attribute) + new_playlist.append(the_curve) + + return new_playlist #this is the true thing returned at the end of this function...(FIXME: clarity) + + return handlePlaylist(self.playlist) + + + def save(self,output_filename): + ''' + saves the playlist in a XML file. + ''' + outfile=file(output_filename,'w') + self.playlist.writexml(outfile,indent='\n') + outfile.close() + + +class HookeConfig: + ''' + Handling of Hooke configuration file + + Mostly based on the simple-yet-useful examples of the Python Library Reference + about xml.dom.minidom + + FIXME: starting to look a mess, should require refactoring + ''' + + def __init__(self): + self.config={} + self.config['plugins']=[] + self.config['drivers']=[] + self.config['plotmanips']=[] + + def load_config(self, filename): + myconfig=file(filename) + + #the following 3 lines are needed to strip newlines. otherwise, since newlines + #are XML elements too, the parser would read them (and re-save them, multiplying + #newlines...) + #yes, I'm an XML n00b + the_file=myconfig.read() + the_file_lines=the_file.split('\n') + the_file=''.join(the_file_lines) + + self.config_tree=xml.dom.minidom.parseString(the_file) + + def getText(nodelist): + #take the text from a nodelist + #from Python Library Reference 13.7.2 + rc = '' + for node in nodelist: + if node.nodeType == node.TEXT_NODE: + rc += node.data + return rc + + def handleConfig(config): + display_elements=config.getElementsByTagName("display") + plugins_elements=config.getElementsByTagName("plugins") + drivers_elements=config.getElementsByTagName("drivers") + workdir_elements=config.getElementsByTagName("workdir") + plotmanip_elements=config.getElementsByTagName("plotmanips") + handleDisplay(display_elements) + handlePlugins(plugins_elements) + handleDrivers(drivers_elements) + handleWorkdir(workdir_elements) + handlePlotmanip(plotmanip_elements) + + def handleDisplay(display_elements): + for element in display_elements: + for attribute in element.attributes.keys(): + self.config[attribute]=element.getAttribute(attribute) + + def handlePlugins(plugins): + for plugin in plugins[0].childNodes: + try: + self.config['plugins'].append(str(plugin.tagName)) + except: #if we allow fancy formatting of xml, there is a text node, so tagName fails for it... + pass + #FIXME: code duplication + def handleDrivers(drivers): + for driver in drivers[0].childNodes: + try: + self.config['drivers'].append(str(driver.tagName)) + except: #if we allow fancy formatting of xml, there is a text node, so tagName fails for it... + pass + + def handlePlotmanip(plotmanips): + for plotmanip in plotmanips[0].childNodes: + try: + self.config['plotmanips'].append(str(plotmanip.tagName)) + except: #if we allow fancy formatting of xml, there is a text node, so tagName fails for it... + pass + + def handleWorkdir(workdir): + wdir=getText(workdir[0].childNodes) + self.config['workdir']=wdir.strip() + + handleConfig(self.config_tree) + #making items in the dictionary more machine-readable + for item in self.config.keys(): + try: + self.config[item]=float(self.config[item]) + except TypeError: #we are dealing with a list, probably. keep it this way. + try: + self.config[item]=eval(self.config[item]) + except: #not a list, not a tuple, probably a string? + pass + except ValueError: #if we can't get it to a number, it must be None or a string + if string.lower(self.config[item])=='none': + self.config[item]=None + else: + pass + + return self.config + + + def save_config(self, config_filename): + print 'Not Implemented.' + pass + + +class ClickedPoint: + ''' + this class defines what a clicked point on the curve plot is + ''' + def __init__(self): + + self.is_marker=None #boolean ; decides if it is a marker + self.is_line_edge=None #boolean ; decides if it is the edge of a line (unused) + self.absolute_coords=(None,None) #(float,float) ; the absolute coordinates of the clicked point on the graph + self.graph_coords=(None,None) #(float,float) ; the coordinates of the plot that are nearest in X to the clicked point + self.index=None #integer ; the index of the clicked point with respect to the vector selected + self.dest=None #0 or 1 ; 0=top plot 1=bottom plot + + + def find_graph_coords(self, xvector, yvector): + ''' + Given a clicked point on the plot, finds the nearest point in the dataset (in X) that + corresponds to the clicked point. + ''' + + #Ye Olde sorting algorithm... + index=0 + best_index=0 + best_diff=10^9 #FIXME:hope we never go over this magic number, doing best_diff=max(xvector)-min(xvector) can be better... + for point in xvector: + diff=abs(point-self.absolute_coords[0]) + if diff1: + for index in range(len(peaks_location)-1): + if peaks_location[index+1]-peaks_location[index] < seedouble: + temp_location=peaks_location[:index]+peaks_location[index+1:] + if temp_location != []: + peaks_location=temp_location + + return peaks_location,peaks_size \ No newline at end of file diff --git a/macro.py b/macro.py new file mode 100644 index 0000000..acb25ce --- /dev/null +++ b/macro.py @@ -0,0 +1,231 @@ +#!/usr/bin/env python + +''' +COMMAND MACRO PLUGIN FOR HOOKE + +Records, saves and executes batches of commands +(c)Alberto Gomez-Casado 2008 +''' + +import libhookecurve as lhc +import os.path + +class macroCommands: + + currentmacro=[] + pause=0 + auxprompt=[] + macrodir=[] + + + def _plug_init(self): + self.currentmacro=[] + self.auxprompt=self.prompt + self.macrodir=os.curdir + if not os.path.exists(os.path.join(self.macrodir,'macros')): + try: + os.mkdir('macros') + self.macrodir=os.path.join(self.macrodir,'macros') + except: + print 'Warning: cannot create macros folder.' + print 'Probably you do not have permissions in your Hooke folder, use macro at your own risk.' + + def collect(self): + + print 'Enter STOP / PAUSE to go back to normal mode\nUNDO to remove last command' + line=[] + while not(line=='STOP' or line=='PAUSE'): + line=raw_input('hooke (macroREC): ') + if line=='PAUSE': + self.pause=1 + self.prompt='hooke (macroPAUSE): ' + break + if line=='STOP': + self.prompt=self.auxprompt + self.do_recordmacro('stop') + break + if line=='UNDO': + self.currentmacro.pop() + continue + param=line.split() + + #FIXME check if accessing param[2] when it doesnt exist breaks something + if param[0] =='export': + exportline=param[0]+' __curve__ ' + if len(param)==3: + exportline=exportline+param[2] + self.currentmacro.append(exportline) + self.onecmd(line) + continue + + if param[0] =='txt': + exportline=param[0] + if len(param)==3: + exportline=exportline+' '+param[2] + exportline=exportline+'__curve__' + self.currentmacro.append(exportline) + self.onecmd(line) + continue + + self.onecmd(line) + + self.currentmacro.append(line) + + + def do_recordmacro(self, args): + '''RECORDMACRO + Stores input commands to create script files + ------- + Syntax: recordmacro [start / stop] + If a macro is currently paused start resumes recording + ''' + + + if len(args)==0: + args='start' + + if args=='stop': + self.pause=0 + self.prompt=self.auxprompt + if len(self.currentmacro) != 0: + answer=raw_input('Do you want to save this macro? ') + if answer[0].lower() == 'y': + self.do_savemacro('') + else: + print 'Macro discarded' + self.currentmacro=[] + else: + print 'Macro was empty' + + if args=='start': + + if self.pause==1: + self.pause=0 + self.collect() + else: + if len(self.currentmacro) != 0: + answer=raw_input('Another macro is already beign recorded\nDo you want to save it?') + if answer[0].lower() == 'y': + self.do_savemacro('') + else: + print 'Old macro discarded, you can start recording the new one' + + self.currentmacro=[] + self.collect() + + + def do_savemacro(self, macroname): + + '''SAVEMACRO + Saves previously recorded macro into a script file for future use + ------- + Syntax: savemacro [macroname] + If no macroname is supplied one will be interactively asked + ''' + + saved_ok=0 + if self.currentmacro==None: + print 'No macro is being recorded!' + return 0 + if len(macroname)==0: + macroname=raw_input('Enter new macro name: ') + if len(macroname) == 0: + print 'Invalid name' + + macroname=os.path.join(self.macrodir,macroname+'.hkm') + if os.path.exists(macroname): + overwrite=raw_input('That name is in use, overwrite?') + if overwrite[0].lower()!='y': + print 'Cancelled save' + return 0 + txtfile=open(macroname,'w+') + self.currentmacro='\n'.join(self.currentmacro) + txtfile.write(self.currentmacro) + txtfile.close() + print 'Saved on '+macroname + self.currentmacro=[] + + def do_execmacro (self, args): + + '''EXECMACRO + Loads a macro and executes it over current curve / playlist + ----- + Syntax: execmacro macroname [playlist] [v] + + macroname.hkm should be present at [hooke]/macros directory + By default the macro will be executed over current curve + passing 'playlist' word as second argument executes macroname + over all curves + By default curve(s) will be processed silently, passing 'v' + as second/third argument will print each command that is + executed + + Note that macros applied to playlists should end by export + commands so the processed curves are not lost + ''' + verbose=0 + cycle=0 + curve=None + + if len(self.currentmacro) != 0: + print 'Warning!: you are calling a macro while recording other' + if len(args) == 0: + print 'You must provide a macro name' + return 0 + args=args.split() + + print 'args ' + ' '.join(args) + + if len(args)>1: + if args[1] == 'playlist': + cycle=1 + print 'Remember! macros applied over playlists should include export orders' + if len(args)>2 and args[2] == 'v': + verbose=1 + else: + if args[1] == 'v': + verbose=1 + print cycle + print verbose + + macropath=os.path.join(self.macrodir,args[0]+'.hkm') + if not os.path.exists(macropath): + print 'Could not find a macro with that name' + return 0 + txtfile=open(macropath) + if cycle ==1: + #print self.current_list + for item in self.current_list: + self.current=item + self.do_plot(0) + + for command in txtfile: + + if verbose==1: + print 'Executing command '+command + testcmd=command.split() + w=0 + for word in testcmd: + if word=='__curve__': + testcmd[w]=os.path.splitext(os.path.basename(item.path))[0] + w=w+1 + self.onecmd(' '.join(testcmd)) + self.current.curve.close_all() + txtfile.seek(0) + else: + for command in txtfile: + testcmd=command.split() + w=0 + if verbose==1: + print 'Executing command '+command + for word in testcmd: + if word=='__curve__': + testcmd[w]=os.path.splitext(os.path.basename(self.current.path))[0] + w=w+1 + self.onecmd(' '.join(testcmd)) + + + + + + diff --git a/massanalysis.py b/massanalysis.py new file mode 100644 index 0000000..e3d3830 --- /dev/null +++ b/massanalysis.py @@ -0,0 +1,143 @@ +#!/usr/bin/env python + +''' +massanalysis.py + +Global analysis of force curves with various parameters + +Requires: +libpeakspot.py +flatfilts.py +''' + + +import libpeakspot as lps +import libhookecurve as lhc +import libhooke as lh +import numpy as np + +import csv + +class massanalysisCommands: + + def _plug_init(self): + self.mass_variables={} + self.interesting_variables=['curve','firstpeak_distance','lastpeak_distance','Npeaks','median_distance','mean_distance'] + self._clean_data() + + def _clean_data(self): + for variable in self.interesting_variables: + self.mass_variables[variable]=[] + + def peak_position_from_contact(self, item, locations): + ''' + calculates X distance of a peak from the contact point + ''' + item.identify(self.drivers) + + real_positions=[] + cut_index=self.find_contact_point() + + #we assume the first is the plot with the force curve + plot=item.curve.default_plots()[0] + xret=plot.vectors[1][0] + + start_x=xret[cut_index] + + real_positions=[abs((xret[index])-(start_x)) for index in locations] + #close all open files + item.curve.close_all() + #needed to avoid *big* memory leaks! + del item.curve + del item + return real_positions + + def do_maplist(self,args): + ''' + MAPLIST + (flatfilts.py) + ---- + pass + ''' + self._clean_data() #if we recall it, clean previous data! + min_deviation=self.convfilt_config['mindeviation'] + + + c=0 + for item in self.current_list: + try: + peak_location,peak_size=self.exec_has_peaks(item, min_deviation) + real_positions=self.peak_position_from_contact(item, peak_location) + + self.mass_variables['Npeaks'].append(len(peak_location)) + + if len(peak_location) > 1: + self.mass_variables['firstpeak_distance'].append(min(real_positions)) + self.mass_variables['lastpeak_distance'].append(max(real_positions)) + + distancepeaks=[] + for index in range(len(real_positions)-1): + distancepeaks.append(real_positions[index+1]-real_positions[index]) + else: + self.mass_variables['firstpeak_distance'].append(0) + self.mass_variables['lastpeak_distance'].append(0) + + if len(peak_location) > 2: + self.mass_variables['median_distance'].append(np.median(distancepeaks)) + self.mass_variables['mean_distance'].append(np.mean(distancepeaks)) + else: + self.mass_variables['median_distance'].append(0) + self.mass_variables['mean_distance'].append(0) + + print 'curve',c + except SyntaxError: + print 'curve',c,'not mapped' + pass + + c+=1 + + def do_plotmap(self,args): + ''' + ''' + args=args.split() + if len(args)>1: + x=self.mass_variables[args[0]] + y=self.mass_variables[args[1]] + else: + print 'Give me two arguments between those:' + print self.interesting_variables + return + + scattermap=lhc.PlotObject() + scattermap.vectors=[[]] + scattermap.vectors[0].append(x) + scattermap.vectors[0].append(y) + + scattermap.units=[args[0],args[1]] + scattermap.styles=['scatter'] + scattermap.destination=1 + + self._send_plot([scattermap]) + + def do_savemaps(self,args): + ''' + args=filename + ''' + + ''' + def csv_write_cols(data, f): + + #from Bruno Desthuillers on comp.lang.python + + writer = csv.writer(f) + keys = data.keys() + writer.writerow(dict(zip(keys,keys))) + for row in zip(*data.values()): + writer.writerow(dict(zip(keys, row))) + ''' + + f=open(args,'wb') + lh.csv_write_dictionary(f,self.mass_variables) + f.close() + + \ No newline at end of file diff --git a/picoforce.py b/picoforce.py new file mode 100755 index 0000000..7c0e926 --- /dev/null +++ b/picoforce.py @@ -0,0 +1,546 @@ +#!/usr/bin/env python + +''' +libpicoforce.py + +Library for interpreting Picoforce force spectroscopy files. + +Copyright (C) 2006 Massimo Sandal (University of Bologna, Italy). + +This program is released under the GNU General Public License version 2. +''' + +import re, struct +from scipy import arange + +import libhookecurve as lhc + +__version__='0.0.0.20080404' + + +class DataChunk(list): + #Dummy class to provide ext and ret methods to the data list. + + def ext(self): + halflen=(len(self)/2) + return self[0:halflen] + + def ret(self): + halflen=(len(self)/2) + return self[halflen:] + +class picoforceDriver(lhc.Driver): + + #Construction and other special methods + + def __init__(self,filename): + ''' + constructor method + ''' + + self.textfile=file(filename) + self.binfile=file(filename,'rb') + + #The 0,1,2 data chunks are: + #0: D (vs T) + #1: Z (vs T) + #2: D (vs Z) + + + self.filepath=filename + self.debug=False + + self.filetype='picoforce' + self.experiment='smfs' + + + #Hidden methods. These are meant to be used only by API functions. If needed, however, + #they can be called just like API methods. + + def _get_samples_line(self): + ''' + Gets the samples per line parameters in the file, to understand trigger behaviour. + ''' + self.textfile.seek(0) + + samps_expr=re.compile(".*Samps") + + samps_values=[] + for line in self.textfile.readlines(): + if samps_expr.match(line): + try: + samps=int(line.split()[2]) #the third word splitted is the offset (in bytes) + samps_values.append(samps) + except: + pass + + #We raise a flag for the fact we meet an offset, otherwise we would take spurious data length arguments. + + return int(samps_values[0]) + + def _get_chunk_coordinates(self): + ''' + This method gets the coordinates (offset and length) of a data chunk in our + Picoforce file. + + It returns a list containing two tuples: + the first element of each tuple is the data_offset, the second is the corresponding + data size. + + In near future probably each chunk will get its own data structure, with + offset, size, type, etc. + ''' + self.textfile.seek(0) + + offset_expr=re.compile(".*Data offset") + length_expr=re.compile(".*Data length") + + data_offsets=[] + data_sizes=[] + flag_offset=0 + + for line in self.textfile.readlines(): + + if offset_expr.match(line): + offset=int(line.split()[2]) #the third word splitted is the offset (in bytes) + data_offsets.append(offset) + #We raise a flag for the fact we meet an offset, otherwise we would take spurious data length arguments. + flag_offset=1 + + #same for the data length + if length_expr.match(line) and flag_offset: + size=int(line.split()[2]) + data_sizes.append(size) + #Put down the offset flag until the next offset is met. + flag_offset=0 + + return zip(data_offsets,data_sizes) + + def _get_data_chunk(self,whichchunk): + ''' + reads a data chunk and converts it in 16bit signed int. + ''' + offset,size=self._get_chunk_coordinates()[whichchunk] + + + self.binfile.seek(offset) + raw_chunk=self.binfile.read(size) + + my_chunk=[] + for data_position in range(0,len(raw_chunk),2): + data_unit_bytes=raw_chunk[data_position:data_position+2] + #The unpack function converts 2-bytes in a signed int ('h'). + #we use output[0] because unpack returns a 1-value tuple, and we want the number only + data_unit=struct.unpack('h',data_unit_bytes)[0] + my_chunk.append(data_unit) + + return DataChunk(my_chunk) + + def _get_Zscan_info(self,index): + ''' + gets the Z scan informations needed to interpret the data chunk. + These info come from the general section, BEFORE individual chunk headers. + + By itself, the function will parse for three parameters. + (index) that tells the function what to return when called by + exposed API methods. + index=0 : returns Zscan_V_LSB + index=1 : returns Zscan_V_start + index=2 : returns Zscan_V_size + ''' + self.textfile.seek(0) + + ciaoforcelist_expr=re.compile(".*Ciao force") + zscanstart_expr=re.compile(".*@Z scan start") + zscansize_expr=re.compile(".*@Z scan size") + + ciaoforce_flag=0 + theline=0 + for line in self.textfile.readlines(): + if ciaoforcelist_expr.match(line): + ciaoforce_flag=1 #raise a flag: zscanstart and zscansize params to read are later + + if ciaoforce_flag and zscanstart_expr.match(line): + raw_Zscanstart_line=line.split() + + if ciaoforce_flag and zscansize_expr.match(line): + raw_Zscansize_line=line.split() + + Zscanstart_line=[] + Zscansize_line=[] + for itemscanstart,itemscansize in zip(raw_Zscanstart_line,raw_Zscansize_line): + Zscanstart_line.append(itemscanstart.strip('[]()')) + Zscansize_line.append(itemscansize.strip('[]()')) + + Zscan_V_LSB=float(Zscanstart_line[6]) + Zscan_V_start=float(Zscanstart_line[8]) + Zscan_V_size=float(Zscansize_line[8]) + + return (Zscan_V_LSB,Zscan_V_start,Zscan_V_size)[index] + + def _get_Z_magnify_scale(self,whichchunk): + ''' + gets Z scale and Z magnify + Here we get Z scale/magnify from the 'whichchunk' only. + whichchunk=1,2,3 + TODO: make it coherent with data_chunks syntaxis (0,1,2) + + In future, should we divide the *file* itself into chunk descriptions and gain + true chunk data structures? + ''' + self.textfile.seek(0) + + z_scale_expr=re.compile(".*@4:Z scale") + z_magnify_expr=re.compile(".*@Z magnify") + + ramp_size_expr=re.compile(".*@4:Ramp size") + ramp_offset_expr=re.compile(".*@4:Ramp offset") + + occurrences=0 + found_right=0 + + + for line in self.textfile.readlines(): + if z_magnify_expr.match(line): + occurrences+=1 + if occurrences==whichchunk: + found_right=1 + raw_z_magnify_expression=line.split() + else: + found_right=0 + + if found_right and z_scale_expr.match(line): + raw_z_scale_expression=line.split() + if found_right and ramp_size_expr.match(line): + raw_ramp_size_expression=line.split() + if found_right and ramp_offset_expr.match(line): + raw_ramp_offset_expression=line.split() + + return float(raw_z_magnify_expression[5]),float(raw_z_scale_expression[7]), float(raw_ramp_size_expression[7]), float(raw_ramp_offset_expression[7]), float(raw_z_scale_expression[5][1:]) + + + #Exposed APIs. + #These are the methods that are meant to be called from external apps. + + def LSB_to_volt(self,chunknum,voltrange=20): + ''' + Converts the LSB data of a given chunk (chunknum=0,1,2) in volts. + First step to get the deflection and the force. + + SYNTAXIS: + item.LSB_to_volt(chunknum, [voltrange]) + + The voltrange is by default set to 20 V. + ''' + return DataChunk([((float(lsb)/65535)*voltrange) for lsb in self.data_chunks[chunknum]]) + + def LSB_to_deflection(self,chunknum,deflsensitivity=None,voltrange=20): + ''' + Converts the LSB data in deflection (meters). + + SYNTAXIS: + item.LSB_to_deflection(chunknum, [deflection sensitivity], [voltrange]) + + chunknum is the chunk you want to parse (0,1,2) + + The deflection sensitivity by default is the one parsed from the file. + The voltrange is by default set to 20 V. + ''' + if deflsensitivity is None: + deflsensitivity=self.get_deflection_sensitivity() + + lsbvolt=self.LSB_to_volt(chunknum) + return DataChunk([volt*deflsensitivity for volt in lsbvolt]) + + def deflection(self): + ''' + Get the actual force curve deflection. + ''' + deflchunk= self.LSB_to_deflection(2) + return deflchunk.ext(),deflchunk.ret() + + def LSB_to_force(self,chunknum=2,Kspring=None,voltrange=20): + ''' + Converts the LSB data (of deflection) in force (newtons). + + SYNTAXIS: + item.LSB_to_force([chunknum], [spring constant], [voltrange]) + + chunknum is the chunk you want to parse (0,1,2). The chunk used is by default 2. + The spring constant by default is the one parsed from the file. + The voltrange is by default set to 20 V. + ''' + if Kspring is None: + Kspring=self.get_spring_constant() + + lsbdefl=self.LSB_to_deflection(chunknum) + return DataChunk([(meter*Kspring) for meter in lsbdefl]) + + def get_Zscan_V_start(self): + return self._get_Zscan_info(1) + + def get_Zscan_V_size(self): + return self._get_Zscan_info(2) + + def get_Z_scan_sensitivity(self): + ''' + gets Z sensitivity + ''' + self.textfile.seek(0) + + z_sensitivity_expr=re.compile(".*@Sens. Zsens") + + for line in self.textfile.readlines(): + if z_sensitivity_expr.match(line): + z_sensitivity=float(line.split()[3]) + #return it in SI units (that is: m/V, not nm/V) + return z_sensitivity*(10**(-9)) + + def get_Z_magnify(self,whichchunk): + ''' + Gets the Z magnify factor. Normally it is 1, unknown exact use as of 2006-01-13 + ''' + return self._get_Z_magnify_scale(whichchunk)[0] + + def get_Z_scale(self,whichchunk): + ''' + Gets the Z scale. + ''' + return self._get_Z_magnify_scale(whichchunk)[1] + + def get_ramp_size(self,whichchunk): + ''' + Gets the -user defined- ramp size + ''' + return self._get_Z_magnify_scale(whichchunk)[2] + + def get_ramp_offset(self,whichchunk): + ''' + Gets the ramp offset + ''' + return self._get_Z_magnify_scale(whichchunk)[3] + + def get_Z_scale_LSB(self,whichchunk): + ''' + Gets the LSB-to-volt conversion factor of the Z data. + (so called hard-scale in the Nanoscope documentation) + + ''' + return self._get_Z_magnify_scale(whichchunk)[4] + + def get_deflection_sensitivity(self): + ''' + gets deflection sensitivity + ''' + self.textfile.seek(0) + + def_sensitivity_expr=re.compile(".*@Sens. DeflSens") + + for line in self.textfile.readlines(): + if def_sensitivity_expr.match(line): + def_sensitivity=float(line.split()[3]) + break + #return it in SI units (that is: m/V, not nm/V) + return def_sensitivity*(10**(-9)) + + def get_spring_constant(self): + ''' + gets spring constant. + We actually find *three* spring constant values, one for each data chunk (F/t, Z/t, F/z). + They are normally all equal, but we retain all three for future... + ''' + self.textfile.seek(0) + + springconstant_expr=re.compile(".*Spring Constant") + + constants=[] + + for line in self.textfile.readlines(): + if springconstant_expr.match(line): + constants.append(float(line.split()[2])) + + return constants[0] + + def get_Zsensorsens(self): + ''' + gets Zsensorsens for Z data. + + This is the sensitivity needed to convert the LSB data in nanometers for the Z-vs-T data chunk. + ''' + self.textfile.seek(0) + + zsensorsens_expr=re.compile(".*Sens. ZSensorSens") + + for line in self.textfile.readlines(): + if zsensorsens_expr.match(line): + zsensorsens_raw_expression=line.split() + #we must take only first occurrence, so we exit from the cycle immediately + break + + return (float(zsensorsens_raw_expression[3]))*(10**(-9)) + + def Z_data(self): + ''' + returns converted ext and ret Z curves. + They're on the second chunk (Z vs t). + ''' + #Zmagnify_zt=self.get_Z_magnify(2) + #Zscale_zt=self.get_Z_scale(2) + Zlsb_zt=self.get_Z_scale_LSB(2) + #rampsize_zt=self.get_ramp_size(2) + #rampoffset_zt=self.get_ramp_offset(2) + zsensorsens=self.get_Zsensorsens() + + ''' + The magic formula that converts the Z data is: + + meters = LSB * V_lsb_conversion_factor * ZSensorSens + ''' + + #z_curves=[item*Zlsb_zt*zsensorsens for item in self.data_chunks[1].pair['ext']],[item*Zlsb_zt*zsensorsens for item in self.data_chunks[1].pair['ret']] + z_curves=[item*Zlsb_zt*zsensorsens for item in self.data_chunks[1].ext()],[item*Zlsb_zt*zsensorsens for item in self.data_chunks[1].ret()] + z_curves=[DataChunk(item) for item in z_curves] + return z_curves + + def Z_extremes(self): + ''' + returns the extremes of the Z values + ''' + zcurves=self.Z_data() + z_extremes={} + z_extremes['ext']=zcurves[0][0],zcurves[0][-1] + z_extremes['ret']=zcurves[1][0],zcurves[1][-1] + + return z_extremes + + def Z_step(self): + ''' + returns the calculated step between the Z values + ''' + zrange={} + zpoints={} + + z_extremes=self.Z_extremes() + + zrange['ext']=abs(z_extremes['ext'][0]-z_extremes['ext'][1]) + zrange['ret']=abs(z_extremes['ret'][0]-z_extremes['ret'][1]) + + #We must take 1 from the calculated zpoints, or when I use the arange function gives me a point more + #with the step. That is, if I have 1000 points, and I use arange(start,stop,step), I have 1001 points... + #For cleanness, solution should really be when using arange, but oh well... + zpoints['ext']=len(self.Z_data()[0])-1 + zpoints['ret']=len(self.Z_data()[1])-1 + #this syntax must become coherent!! + return (zrange['ext']/zpoints['ext']),(zrange['ret']/zpoints['ret']) + + def Z_domains(self): + ''' + returns the Z domains on which to plot the force data. + + The Z domains are returned as a single long DataChunk() extended list. The extension and retraction part + can be extracted using ext() and ret() methods. + ''' + x1step=self.Z_step()[0] + x2step=self.Z_step()[1] + + try: + xext=arange(self.Z_extremes()['ext'][0],self.Z_extremes()['ext'][1],-x1step) + xret=arange(self.Z_extremes()['ret'][0],self.Z_extremes()['ret'][1],-x2step) + except: + xext=arange(0,1) + xret=arange(0,1) + print 'picoforce.py: Warning. xext, xret domains cannot be extracted.' + + if not (len(xext)==len(xret)): + xext=xret + if self.debug: + #print warning + print "picoforce.py: Warning. Extension and retraction domains have different sizes." + print "length extension: ", len(xext) + print "length retraction: ", len(xret) + print "You cannot trust the resulting curve." + print "Until a solution is found, I substitute the ext domain with the ret domain. Sorry." + + + return DataChunk(xext.tolist()+xret.tolist()) + + def Z_scan_size(self): + return self.get_Zscan_V_size()*self.get_Z_scan_sensitivity() + + def Z_start(self): + return self.get_Zscan_V_start()*self.get_Z_scan_sensitivity() + + def ramp_size(self,whichchunk): + ''' + to be implemented if needed + ''' + raise "Not implemented yet." + + + def ramp_offset(self,whichchunk): + ''' + to be implemented if needed + ''' + raise "Not implemented yet." + + def detriggerize(self, forcext): + ''' + Cuts away the trigger-induced s**t on the extension curve. + DEPRECATED + cutindex=2 + startvalue=forcext[0] + + for index in range(len(forcext)-1,2,-2): + if forcext[index]>startvalue: + cutindex=index + else: + break + + return cutindex + ''' + return 0 + + def is_me(self): + ''' + self-identification of file type magic + ''' + curve_file=file(self.filepath) + header=curve_file.read(30) + curve_file.close() + + if header[2:17] == 'Force file list': #header of a picoforce file + self.data_chunks=[self._get_data_chunk(num) for num in [0,1,2]] + return True + else: + return False + + def close_all(self): + ''' + Explicitly closes all files + ''' + self.textfile.close() + self.binfile.close() + + def default_plots(self): + ''' + creates the default PlotObject + ''' + + + force=self.LSB_to_force() + zdomain=self.Z_domains() + + samples=self._get_samples_line() + #cutindex=0 + #cutindex=self.detriggerize(force.ext()) + + main_plot=lhc.PlotObject() + + main_plot.vectors=[[zdomain.ext()[0:samples], force.ext()[0:samples]],[zdomain.ret()[0:samples], force.ret()[0:samples]]] + main_plot.normalize_vectors() + main_plot.units=['meters','newton'] + main_plot.destination=0 + main_plot.title=self.filepath + + + return [main_plot] \ No newline at end of file diff --git a/procplots.py b/procplots.py new file mode 100755 index 0000000..f1cc0a3 --- /dev/null +++ b/procplots.py @@ -0,0 +1,251 @@ +#!/usr/bin/env python + +''' +PROCPLOTS +Processed plots plugin for force curves. + +Licensed under the GNU GPL version 2 +''' +from libhooke import WX_GOOD +import wxversion +wxversion.select(WX_GOOD) + +import wx +import libhookecurve as lhc +import numpy as np +import scipy as sp +import scipy.signal +import copy + +class procplotsCommands: + + def _plug_init(self): + pass + + def do_derivplot(self,args): + ''' + DERIVPLOT + (procplots.py plugin) + Plots the derivate (actually, the discrete differentiation) of the current force curve + --------- + Syntax: derivplot + ''' + dplot=self.derivplot_curves() + plot_graph=self.list_of_events['plot_graph'] + wx.PostEvent(self.frame,plot_graph(plots=[dplot])) + + def derivplot_curves(self): + ''' + do_derivplot helper function + + create derivate plot curves for force curves. + ''' + dplot=lhc.PlotObject() + dplot.vectors=[] + + for vector in self.plots[0].vectors: + dplot.vectors.append([]) + dplot.vectors[-1].append(vector[0][:-1]) + dplot.vectors[-1].append(np.diff(vector[1])) + + dplot.destination=1 + dplot.units=self.plots[0].units + + return dplot + + def do_subtplot(self,args): + ''' + SUBTPLOT + (procplots.py plugin) + Plots the difference between ret and ext current curve + ------- + Syntax: subtplot + ''' + #FIXME: sub_filter and sub_order must be args + + if len(self.plots[0].vectors) != 2: + print 'This command only works on a curve with two different plots.' + pass + + outplot=self.subtract_curves(sub_order=1) + + plot_graph=self.list_of_events['plot_graph'] + wx.PostEvent(self.frame,plot_graph(plots=[outplot])) + + def subtract_curves(self, sub_order): + ''' + subtracts the extension from the retraction + ''' + xext=self.plots[0].vectors[0][0] + yext=self.plots[0].vectors[0][1] + xret=self.plots[0].vectors[1][0] + yret=self.plots[0].vectors[1][1] + + #we want the same number of points + maxpoints_tot=min(len(xext),len(xret)) + xext=xext[0:maxpoints_tot] + yext=yext[0:maxpoints_tot] + xret=xret[0:maxpoints_tot] + yret=yret[0:maxpoints_tot] + + if sub_order: + ydiff=[yretval-yextval for yretval,yextval in zip(yret,yext)] + else: #reverse subtraction (not sure it's useful, but...) + ydiff=[yextval-yretval for yextval,yretval in zip(yext,yret)] + + outplot=copy.deepcopy(self.plots[0]) + outplot.vectors[0][0], outplot.vectors[1][0] = xext,xret #FIXME: if I use xret, it is not correct! + outplot.vectors[1][1]=ydiff + outplot.vectors[0][1]=[0 for item in outplot.vectors[1][0]] + + return outplot + + +#-----PLOT MANIPULATORS + def plotmanip_median(self, plot, current, customvalue=None): + ''' + does the median of the y values of a plot + ''' + if customvalue: + median_filter=customvalue + else: + median_filter=self.config['medfilt'] + + if median_filter==0: + return plot + + if float(median_filter)/2 == int(median_filter)/2: + median_filter+=1 + + nplots=len(plot.vectors) + c=0 + while cstartvalue: + cutindex=index + else: + break + + plot.vectors[0][0]=plot.vectors[0][0][:cutindex] + plot.vectors[0][1]=plot.vectors[0][1][:cutindex] + + return plot + ''' + + + +#FFT--------------------------- + def fft_plot(self, vector): + ''' + calculates the fast Fourier transform for the selected vector in the plot + ''' + fftplot=lhc.PlotObject() + fftplot.vectors=[[]] + + fftlen=len(vector)/2 #need just 1/2 of length + fftplot.vectors[-1].append(np.arange(1,fftlen).tolist()) + + try: + fftplot.vectors[-1].append(abs(np.fft(vector)[1:fftlen]).tolist()) + except TypeError: #we take care of newer NumPy (1.0.x) + fftplot.vectors[-1].append(abs(np.fft.fft(vector)[1:fftlen]).tolist()) + + + fftplot.destination=1 + + + return fftplot + + + def do_fft(self,args): + ''' + FFT + (procplots.py plugin) + Plots the fast Fourier transform of the selected plot + --- + Syntax: fft [top,bottom] [select] [0,1...] + + By default, fft performs the Fourier transform on all the 0-th data set on the + top plot. + + [top,bottom]: which plot is the data set to fft (default: top) + [select]: you pick up two points on the plot and fft only the segment between + [0,1,...]: which data set on the selected plot you want to fft (default: 0) + ''' + + #args parsing + #whatplot = plot to fft + #whatset = set to fft in the plot + select=('select' in args) + if 'top' in args: + whatplot=0 + elif 'bottom' in args: + whatplot=1 + else: + whatplot=0 + whatset=0 + for arg in args: + try: + whatset=int(arg) + except ValueError: + pass + + if select: + points=self._measure_N_points(N=2, whatset=whatset) + boundaries=[points[0].index, points[1].index] + boundaries.sort() + y_to_fft=self.plots[whatplot].vectors[whatset][1][boundaries[0]:boundaries[1]] #y + else: + y_to_fft=self.plots[whatplot].vectors[whatset][1] #y + + fftplot=self.fft_plot(y_to_fft) + fftplot.units=['frequency', 'power'] + plot_graph=self.list_of_events['plot_graph'] + wx.PostEvent(self.frame,plot_graph(plots=[fftplot])) + \ No newline at end of file diff --git a/superimpose.py b/superimpose.py new file mode 100644 index 0000000..fdd0a33 --- /dev/null +++ b/superimpose.py @@ -0,0 +1,98 @@ +#!/usr/bin/env python + +from libhooke import WX_GOOD +import wxversion +wxversion.select(WX_GOOD) +from wx import PostEvent + +import libhookecurve as lhc +from numpy import arange, mean + +class superimposeCommands: + + def _plug_init(self): + self.imposed=[] + + def do_selimpose(self,args): + ''' + SELIMPOSE (superimpose.py plugin) + Hand-selects the curve portion to superimpose + ''' + #fixme: set to superimpose should be in args + + if args=='clear': + self.imposed=[] + return + + current_set=1 + + points=self._measure_two_points() + boundaries=[points[0].index, points[1].index] + boundaries.sort() + + theplot=self.plots[0] + #append the selected section + self.imposed.append([]) + self.imposed[-1].append(theplot.vectors[1][0][boundaries[0]:boundaries[1]]) #x + self.imposed[-1].append(theplot.vectors[1][1][boundaries[0]:boundaries[1]]) #y + + #align X first point + self.imposed[-1][0] = [item-self.imposed[-1][0][0] for item in self.imposed[-1][0]] + #align Y first point + self.imposed[-1][1] = [item-self.imposed[-1][1][0] for item in self.imposed[-1][1]] + + def do_plotimpose(self,args): + ''' + PLOTIMPOSE (sumperimpose.py plugin) + plots superimposed curves + ''' + imposed_object=lhc.PlotObject() + imposed_object.vectors=self.imposed + print 'Plotting',len(imposed_object.vectors),'imposed curves' + + imposed_object.normalize_vectors() + + imposed_object.units=self.plots[0].units + imposed_object.title='Imposed curves' + imposed_object.destination=1 + + plot_graph=self.list_of_events['plot_graph'] + PostEvent(self.frame,plot_graph(plots=[imposed_object])) + + def do_plotavgimpose(self,args): + ''' + PLOTAVGIMPOSE (superimpose.py plugin) + Plots the average of superimposed curves using a running window + ''' + step=(-5*(10**-10)) + #find extension of each superimposed curve + min_x=[] + for curve in self.imposed: + min_x.append(min(curve[0])) + + #find minimum extension + min_ext_limit=max(min_x) + + x_avg=arange(step,min_ext_limit,step) + y_avg=[] + for value in x_avg: + to_avg=[] + for curve in self.imposed: + for xvalue, yvalue in zip(curve[0],curve[1]): + if xvalue >= (value+step) and xvalue <= (value-step): + to_avg.append(yvalue) + y_avg.append(mean(to_avg)) + + print 'len x len y' + print len(x_avg), len(y_avg) + print y_avg + + avg_object=lhc.PlotObject() + avg_object.vectors=[[x_avg, y_avg]] + avg_object.normalize_vectors() + avg_object.units=self.plots[0].units + avg_object.title="Average curve" + avg_object.destination=1 + + plot_graph=self.list_of_events['plot_graph'] + PostEvent(self.frame,plot_graph(plots=[avg_object])) \ No newline at end of file diff --git a/tutorial.py b/tutorial.py new file mode 100755 index 0000000..cb24d68 --- /dev/null +++ b/tutorial.py @@ -0,0 +1,544 @@ +#!/usr/bin/env python + +''' +TUTORIAL PLUGIN FOR HOOKE + +This plugin contains example commands to teach how to write an Hooke plugin, including description of main Hooke +internals. +(c)Massimo Sandal 2007 +''' + +import libhookecurve as lhc + +import numpy as np + +''' +SYNTAX OF DATA TYPE DECLARATION: + type = type of object + [ type ] = list containing objects of type + {typekey:typearg} = dictionary with keys of type typekey and args of type typearg + ( type ) = tuple containing objects of type +''' + + +class tutorialCommands: + ''' + Here we define the class containing all the Hooke commands we want to define + in the plugin. + + Notice the class name!! + The syntax is filenameCommands. That is, if your plugin is pluggy.py, your class + name is pluggyCommands. + + Otherwise, the class will be ignored by Hooke. + ''' + + def _plug_init(self): + ''' + This is the plugin initialization. + When Hooke starts and the plugin is loaded, this function is executed. + If there is something you need to do when Hooke starts, code it in this function. + ''' + print 'I am the Tutorial plugin initialization!' + + #Here we initialize a local configuration variable; see plotmanip_absvalue() for explanation. + self.config['tutorial_absvalue']=0 + pass + + def do_nothing(self,args): + ''' + This is a boring but working example of an actual Hooke command. + A Hooke command is a function of the xxxxCommands class, which is ALWAYS defined + this way: + + def do_nameofcommand(self,args) + + *do_ is needed to make Hooke understand this function is a command + *nameofcommand is how the command will be called in the Hooke command line. + *self is, well, self + *args is ALWAYS needed (otherwise Hooke will crash executing the command). We will see + later what args is. + + Note that if you now start Hooke with this plugin activated and you type in the Hooke command + line "help nothing" you will see this very text as output. So the help of a command is a + string comment below the function definition, like this one. + + Commands usually return None. + ''' + print 'I am a Hooke command. I do nothing.' + + def do_printargs(self,args): + ''' + This command prints the args you give to it. + args is always a string, that contains everything you write after the command. + So if you issue "mycommand blah blah 12345" args is "blah blah 12345". + + Again, args is needed in the definition even if your command does not use it. + ''' + print 'You gave me those args: '+args + + def help_tutorial(self): + ''' + This is a help function. + If you want a help function for something that is not a command, you can write a help + function like this. Calling "help tutorial" will execute this function. + ''' + print 'You called help_tutorial()' + + def do_environment(self,args): + ''' + This plugin contains a panoramic of the Hooke command line environment variables, + and prints their current value. + ''' + + '''self.current_list + TYPE: [ libhookecurve.HookeCurve ], len=variable + contains the actual playlist of Hooke curve objects. + Each HookeCurve object represents a reference to a data file. + We will see later in detail how do they work. + ''' + print 'current_list length:',len(self.current_list) + print 'current_list 0th:',self.current_list[0] + + '''self.pointer + TYPE: int + contains the index of + the current curve in the playlist + ''' + print 'pointer: ',self.pointer + + '''self.current + TYPE: libhookecurve.HookeCurve + contains the current curve displayed. + We will see later how it works. + ''' + print 'current:',self.current + + '''self.plots + TYPE: [ libhookecurve.PlotObject ], len=1,2 + contains the current default plots. + Each PlotObject contains all info needed to display + the plot: apart from the data vectors, the title, destination + etc. + Usually self.plots[0] is the default topmost plot, self.plots[1] is the + accessory bottom plot. + ''' + print 'plots:',self.plots + + '''self.config + TYPE: { string:anything } + contains the current Hooke configuration variables, in form of a dictionary. + ''' + print 'config:',self.config + + '''self.plotmanip + TYPE: [ function ] + Contains the ordered plot manipulation functions. + These functions are called to modify the default plot by default before it is plotted. + self.plots contains the plot passed through the plot manipulators. + We will see it better later. + *YOU SHOULD NEVER MODIFY THAT* + ''' + print 'plotmanip: ',self.plotmanip + + '''self.drivers + TYPE: [ class ] + Contains the plot reading drivers. + *YOU SHOULD NEVER MODIFY THAT* + ''' + print 'drivers: ',self.drivers + + '''self.frame + TYPE: wx.Frame + Contains the wx Frame of the GUI. + ***NEVER, EVER TOUCH THAT.*** + ''' + print 'frame: ',self.frame + + '''self.list_of_events + TYPE: { string:wx.Event } + Contains the wx.Events to communicate with the GUI. + Usually not necessary to use it, unless you want + to create a GUI plugin. + ''' + print 'list of events:',self.list_of_events + + '''self.events_from_gui + TYPE: Queue.Queue + Contains the Queue where data from the GUI is put. + Usually not necessary to use it, unless you want + to create a GUI plugin. + ''' + print 'events from gui:',self.events_from_gui + + '''self.playlist_saved + TYPE: Int (0/1) ; Boolean + Flag that tells if the playlist has been saved or not. + ''' + print 'playlist saved:',self.playlist_saved + + '''self.playlist_name + TYPE: string + Name of current playlist + ''' + print 'playlist name:',self.playlist_name + + '''self.notes_saved + TYPE: Int (0/1) ; Boolean + Flag that tells if the playlist has been saved or not. + ''' + print 'notes saved:',self.notes_saved + + + def do_myfirstplot(self,args): + ''' + In this function, we see how to create a PlotObject and send it to the screen. + ***Read the code of PlotObject in libhookecurve.py before!***. + ''' + + #We generate some interesting data to plot for this example. + xdata1=np.arange(-5,5,0.1) + xdata2=np.arange(-5,5,0.1) + ydata1=[item**2 for item in xdata1] + ydata2=[item**3 for item in xdata2] + + #Create the object. + #The PlotObject class lives in the libhookecurve library. + myplot=lhc.PlotObject() + myplot.vectors=[[],[]] #Decide we will have two data sets in this plot + ''' + The *data* of the plot live in the .vectors list. + + plot.vectors is a multidimensional array: + plot.vectors[0]=set1 + plot.vectors[1]=set2 + plot.vectors[2]=sett3 + etc. + + 2 curves in a x,y plot are: + [[[x1],[y1]],[[x2],[y2]]] + for example: + x1 y1 x2 y2 + [[[1,2,3,4],[10,20,30,40]],[[3,6,9,12],[30,60,90,120]]] + x1 = self.vectors[0][0] + y1 = self.vectors[0][1] + x2 = self.vectors[1][0] + y2 = self.vectors[1][1] + ''' + #Pour 0-th dataset into myplot: + myplot.vectors[0].append(xdata1) #...x + myplot.vectors[0].append(ydata1) #...then y + + #Pour 1-st dataset into myplot: + myplot.vectors[1].append(xdata2) #...x + myplot.vectors[1].append(ydata2) #...then y + + #Add units to x and y axes + #units=[string, string] + myplot.units=['x axis','y axis'] + + #Where do we want the plot? 0=top, 1=bottom + myplot.destination=1 + + '''Send it to the GUI. + Note that you *have* to encapsulate it into a list, so you + have to send [myplot], not simply myplot. + + You can also send more two plots at once + self.send_plot([plot1,plot2]) + ''' + self._send_plot([myplot]) + + + def do_myfirstscatter(self,args): + ''' + How to draw a scatter plot. + ''' + #We generate some interesting data to plot for this example. + xdata1=np.arange(-5,5,1) + xdata2=np.arange(-5,5,1) + ydata1=[item**2 for item in xdata1] + ydata2=[item**3 for item in xdata2] + + myplot=lhc.PlotObject() + myplot.vectors=[[],[]] + #Pour 0-th dataset into myplot: + myplot.vectors[0].append(xdata1) #...x + myplot.vectors[0].append(ydata1) #...then y + + #Pour 1-st dataset into myplot: + myplot.vectors[1].append(xdata2) #...x + myplot.vectors[1].append(ydata2) #...then y + + #Add units to x and y axes + myplot.units=['x axis','y axis'] + + #Where do we want the plot? 0=top, 1=bottom + myplot.destination=1 + + '''None=standard line plot + 'scatter'=scatter plot + By default, the styles attribute is an empty list. If you + want to define a scatter plot, you must define all other + plots as None or 'scatter', depending on what you want. + + Here we define the second set to be plotted as scatter, + and the first to be plotted as line. + ''' + myplot.styles=[None,'scatter'] + + self._send_plot([myplot]) + + + def do_clickaround(self,args): + ''' + Here we click two points on the curve and take some parameters from the points + we have clicked. + ''' + + ''' + points = self._measure_N_points(N=Int, whatset=Int) + *N = number of points to measure(1...n) + *whatset = data set to measure (0,1...n) + *points = a list of ClickedPoint objects, one for each point requested + ''' + points=self._measure_N_points(N=2,whatset=1) + print 'You clicked the following points.' + + ''' + These are the absolute coordinates of the + point clicked. + [float, float] = x,y + ''' + print 'Absolute coordinates:' + print points[0].absolute_coords + print points[1].absolute_coords + print + + ''' + These are the coordinates of the points + clicked, remapped on the graph. + Hooke looks at the graph point which X + coordinate is next to the X coordinate of + the point measured, and uses that point + as the actual clicked point. + [float, float] = x,y + ''' + print 'Coordinates on the graph:' + print points[0].graph_coords + print points[1].graph_coords + print + + ''' + These are the indexes of the clicked points + on the dataset vector. + ''' + print 'Index of points on the graph:' + print points[0].index + print points[1].index + + + def help_thedifferentplots(self): + ''' + The *three* different default plots you should be familiar with + in Hooke. + + Each plot contains of course the respective data in their + vectors attribute, so here you learn also which data access for + each situation. + ''' + print ''' + 1. THE RAW, CURRENT PLOTS + + self.current + --- + Contains the current libhookecurve.HookeCurve container object. + A HookeCurve object defines only two default attributes: + + * self.current.path = string + The path of the current displayed curve + + * self.current.curve = libhookecurve.Driver + The curve object. This is not only generated by the driver, + this IS a driver instance in itself. + This means that in self.current.curve you can access the + specific driver APIs, if you know them. + + And defines only one method: + * self.current.identify() + Fills in the self.current.curve object. + See in the cycling tutorial. + + ***** + The REAL curve data actually lives in: + --- + * self.current.curve.default_plots() = [ libhooke.PlotObject ] + Contains the raw PlotObject-s, as "spitted out" by the driver, without any + intervention. + This is as close to the raw data as Hooke gets. + + One or two plots can be spit out; they are always enclosed in a list. + ***** + + Methods of self.current.curve are: + --- + + * self.current.curve.is_me() + (Used by identify() only.) + + * self.current.curve.close_all() + Closes all driver open files; see the cycling tutorial. + ''' + + print ''' + 2. THE PROCESSED, DEFAULT PLOT + + The plot that is spitted out by the driver is *not* the usual default plot + that is displayed by calling "plot" at the Hooke prompt. + + This is because the raw, driver-generated plot is usually *processed* by so called + *plot processing* functions. We will see in the tutorial how to define + them. + + For example, in force spectroscopy force curves, raw data are automatically corrected + for deflection. Other data can be, say, filtered by default. + + The default plots are accessible in + self.plots = [ libhooke.PlotObject ] + + self.plots[0] is usually the topmost plot + self.plots[1] is usually the bottom plot (if present) + ''' + + print ''' + 3. THE PLOT DISPLAYED RIGHT NOW. + + Sometimes the plots you are displaying *right now* is different from the previous + two. You may have a fit trace, you may have issued some command that spits out + a custom plot and you want to rework that, whatever. + + You can obtain in any moment the plot currently displayed by Hooke by issuing + + PlotObject = self._get_displayed_plot(dest) + * dest = Int (0/1) + dest=0 : top plot + dest=1 : bottom plot + ''' + + + def do_cycling(self,args): + ''' + Here we cycle through our playlist and print some info on the curves we find. + Cycling through the playlist needs a bit of care to avoid memory leaks and dangling + open files... + + Look at the source code for more information. + ''' + + def things_when_cycling(item): + ''' + We encapsulate here everything has to open the actual curve file. + By doing it all here, we avoid to do acrobacies when deleting objects etc. + in the main loop: we do the dirty stuff here. + ''' + + ''' + identify() + + This method looks for the correct driver in self.drivers to use; + and puts the curve content in the .curve attribute. + Basically, until identify() is called, the HookeCurve object + is just an empty shell. When identify() is called (usually by + the Hooke plot routine), the HookeCurve object is "filled" with + the actual curve. + ''' + + item.identify(self.drivers) + + ''' + After the identify(), item.curve contains the curve, and item.curve.default_plots() behaves exactly like + self.current.curve.default_plots() -but for the given item. + ''' + itplot=item.curve.default_plots() + + print 'length of X1 vector:',len(itplot[0].vectors[0][0]) #just to show something + + ''' + The following three lines are a magic spell you HAVE to do + before closing the function. + (Otherwise you will be plagued by unpredicatable, system-dependent bugs.) + ''' + item.curve.close_all() #Avoid open files dangling + del item.curve #Avoid memory leaks + del item #Just be paranoid. Don't ask. + + return + + + c=0 + for item in self.current_list: + print 'Looking at curve ',c,'of',len(self.current_list) + things_when_cycling(item) + c+=1 + + return + + + + def plotmanip_absvalue(self, plot, current, customvalue=None): + ''' + This function defines a PLOT MANIPULATOR. + A plot manipulator is a function that takes a plot in input, does something to the plot + and returns the modified plot in output. + The function, once plugged, gets automatically called everytime self.plots is updated + + For example, in force spectroscopy force curves, raw data are automatically corrected + for deflection. Other data can be, say, filtered by default. + + To create and activate a plot manipulator you have to: + * Write a function (like this) which name starts with "plotmanip_" (just like commands + start with "do_") + * The function must support four arguments: + self : (as usual) + plot : a plot object + current : (usually not used, deprecated) + customvalue=None : a variable containing custom value(s) you need for your plot manipulators. + * The function must return a plot object. + * Add an entry in hooke.conf: if your function is "plotmanip_something" you will have + to add in the plotmanips section: example + + + + + + + + + Important: Plot manipulators are *in pipe*: each plot manipulator output becomes the input of the next one. + The order in hooke.conf *is the order* in which plot manipulators are connected, so in the example above + we have: + self.current.curve.default_plots() --> detriggerize --> correct --> median --> something --> self.plots + ''' + + ''' + Here we see what is in a configuration variable to enable/disable the plot manipulator as user will using + the Hooke "set" command. + Typing "set tutorial_absvalue 0" disables the plot manipulator; typing "set tutorial_absvalue 1" will enable it. + ''' + if not self.config['tutorial_absvalue']: + return plot + + #We do something to the plot, for demonstration's sake + #If we needed variables, we would have used customvalue. + plot.vectors[0][1]=[abs(i) for i in plot.vectors[0][1]] + plot.vectors[1][1]=[abs(i) for i in plot.vectors[1][1]] + + #Return the plot object. + return plot + + +#TODO IN TUTORIAL: +#how to add lines to an existing plot!! +#peaks +#configuration files +#gui plugins diff --git a/tutorialdriver.py b/tutorialdriver.py new file mode 100644 index 0000000..b16512a --- /dev/null +++ b/tutorialdriver.py @@ -0,0 +1,208 @@ +#!/usr/bin/env python + +''' +tutorialdriver.py + +TUTORIAL DRIVER FOR HOOKE + +Example driver to teach how to write a driver for data types. +(c)Massimo Sandal 2008 +''' + +''' +Here we define a (fake) file format that is read by this driver. The file format is as following: + +TUTORIAL_FILE +PLOT1 +X1 +n1 +n2 +... +nN +Y1 +n1 +n2 +... +nN +X2 +n1 +n2 +.. +nN +Y2 +n1 +n2 +.. +nN +PLOT2 +X1 +... +Y1 +... +X2 +... +Y2 +... +END +that is, two plots with two datasets each. +''' + +import libhookecurve as lhc #We need to import this library to define some essential data types + +class tutorialdriverDriver(lhc.Driver): + ''' + This is a *generic* driver, not a specific force spectroscopy driver. + See the written documentation to see what a force spectroscopy driver must be defined to take into account Hooke facilities. + + Our driver is a class with the name convention nameofthedriverDriver, where "nameofthedriver" is the filename. + The driver must inherit from the parent class lhc.Driver, so the syntax is + class nameofthedriverDriver(lhc.Driver) + ''' + def __init__(self, filename): + ''' + THIS METHOD MUST BE DEFINED. + The __init__ method MUST call the filename, so that it can open the file. + ''' + self.filename=filename #self.filename can always be useful, and should be defined + self.filedata = open(filename,'r') #We open the file + ''' + In this case, we have a data format that is just a list of ASCII values, so we can just divide that in rows, and generate a list + with each item being a row. + Of course if your data files are binary, or follow a different approach, do whatever you need. :) + ''' + self.data = list(self.filedata) + self.filedata.close() #remember to close the file + + '''These are two strings that can be used by Hooke commands/plugins to understand what they are looking at. They have no other + meaning. They have to be somehow defined however - commands often look for those variables. + + self.filetype should contain the name of the exact filetype defined by the driver (so that filetype-specific commands can know + if they're dealing with the correct filetype) + self.experiment should contain instead the type of data involved (for example, various drivers can be used for force-clamp experiments, + but hooke commands could like to know if we're looking at force clamp data, regardless of their origin, and not other + kinds of data) + + Of course, all other variables you like can be defined in the class. + ''' + self.filetype = 'tutorial' + self.experiment = 'generic' + + def is_me(self): + ''' + THIS METHOD MUST BE DEFINED. + RETURNS: Boolean (True or False) + This method must be an heuristic that looks at the file content and decides if the file can be opened by the driver itself. + It returns True if the file opened can be interpreted by the current driver, False otherwise. + Defining this method allows Hooke to understand what kind of files we're looking at automatically. + + We have to re-open/re-close the file here. + ''' + + myfile=open(self.filename, 'r') + headerline=myfile.readlines()[0] #we take the first line + myfile.close() + + ''' + Here, our "magic fingerprint" is the TUTORIAL_FILE header. Of course, depending on the data file, you can have interesting + headers, or patterns, etc. that you can use to guess the data format. What matters is successful recognizing, and returning + a boolean (True/False). + ''' + if headerline[:-1]=='TUTORIAL_FILE': #[:-1], otherwise the return character is included in the line + return True + else: + return False + + def _generate_vectors(self): + ''' + Here we parse the data and generate the raw vectors. This method has only to do with the peculiar file format here, so it's of + no big interest (I just wrote it to present a functional driver). + + Only thing to remember, it can be nice to call methods that are used only "internally" by the driver (or by plugins) with a + "_" prefix, so to have a visual remark. But it's just an advice. + ''' + vectors={'PLOT1':[[],[],[],[]] , 'PLOT2':[[],[],[],[]]} + positions={'X1':0,'Y1':1,'X2':2,'Y2':3} + whatplot='' + pos=0 + for item in self.data: + try: + num=float(item) + vectors[whatplot][pos].append(num) + except ValueError: + if item[:-1]=='PLOT1': + whatplot=item[:-1] + elif item[:-1]=='PLOT2': + whatplot=item[:-1] + elif item[0]=='X' or item[0]=='Y': + pos=positions[item[:-1]] + else: + pass + + return vectors + + def close_all(self): + ''' + THIS METHOD MUST BE DEFINED. + This method is a precaution method that is invoked when cycling to avoid eventually dangling open files. + In this method, all file objects defined in the driver must be closed. + ''' + self.filename.close() + + + def default_plots(self): + ''' + THIS METHOD MUST BE DEFINED. + RETURNS: [ lhc.PlotObject ] or [ lhc.PlotObject, lhc.PlotObject] + + This is the method that returns the plots to Hooke. + It must return a list with *one* or *two* PlotObjects. + + See the libhookecurve.py source code to see how PlotObjects are defined and work in detail. + ''' + gen_vectors=self._generate_vectors() + + #Here we create the main plot PlotObject and initialize its vectors. + main_plot=lhc.PlotObject() + main_plot.vectors=[] + #The same for the other plot. + other_plot=lhc.PlotObject() + other_plot.vectors=[] + + ''' + Now we fill the plot vectors with our data. + set 1 set 2 + The "correct" shape of the vector is [ [[x1,x2,x3...],[y1,y2,y3...]] , [[x1,x2,x3...],[y1,y2,y3...]] ], so we have to put stuff in this way into it. + + The add_set() method takes care of this , just use plot.add_set(x,y). + ''' + main_plot.add_set(gen_vectors['PLOT1'][0],gen_vectors['PLOT1'][1]) + main_plot.add_set(gen_vectors['PLOT1'][2],gen_vectors['PLOT1'][3]) + + other_plot.add_set(gen_vectors['PLOT2'][0],gen_vectors['PLOT2'][1]) + other_plot.add_set(gen_vectors['PLOT2'][2],gen_vectors['PLOT2'][3]) + + ''' + normalize_vectors() trims the vectors, so that if two x/y couples are of different lengths, the latest + points are trimmed (otherwise we have a python error). Always a good idea to run it, to avoid crashes. + ''' + main_plot.normalize_vectors() + other_plot.normalize_vectors() + + ''' + Here we define: + - units: [string, string], define the measure units of X and Y axes + - destination: 0/1 , defines where to plot the plot (0=top, 1=bottom), default=0 + - title: string , the plot title. + + for each plot. + Again, see libhookecurve.py comments for details. + ''' + main_plot.units=['unit of x','unit of y'] + main_plot.destination=0 + main_plot.title=self.filename+' main' + + other_plot.units=['unit of x','unit of y'] + other_plot.destination=1 + other_plot.title=self.filename+' other' + + return [main_plot, other_plot] \ No newline at end of file