1 # Copyright (C) 2010 Alberto Gomez-Casado
2 # W. Trevor King <wking@drexel.edu>
4 # This file is part of Hooke.
6 # Hooke is free software: you can redistribute it and/or
7 # modify it under the terms of the GNU Lesser General Public
8 # License as published by the Free Software Foundation, either
9 # version 3 of the License, or (at your option) any later version.
11 # Hooke is distributed in the hope that it will be useful,
12 # but WITHOUT ANY WARRANTY; without even the implied warranty of
13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 # GNU Lesser General Public License for more details.
16 # You should have received a copy of the GNU Lesser General Public
17 # License along with Hooke. If not, see
18 # <http://www.gnu.org/licenses/>.
20 #FIXME clean this, probably some dependencies are not needed
22 from ..libhooke import WX_GOOD, ClickedPoint
24 wxversion.select(WX_GOOD)
25 from wx import PostEvent
33 warnings.simplefilter('ignore',np.RankWarning)
37 global EVT_MEASURE_WLC
39 global events_from_fit
40 events_from_fit=Queue.Queue() #GUI ---> CLI COMMUNICATION
43 class multifitCommands:
45 def do_multifit(self,args):
49 Presents curves for manual analysis in a comfortable mouse-only fashion.
50 Obtains contour length, persistance length, rupture force and
52 WLC is shown in red, eFJC in blue.
55 multifit [pl=value] [kl=value] [t=value] [slopew=value] [basew=value]
58 pl=[value] and kl=[value]: Use a fixed persistent length (WLC) or Kuhn
59 length (eFJC) for the fit. If pl is not given, the fit will be
61 DO NOT put spaces between 'pl', '=' and the value.
62 eFJC fit works better with a fixed kl
63 The value must be in nanometers.
65 t=[value] : Use a user-defined temperature. The value must be in
66 kelvins; by default it is 293 K.
67 DO NOT put spaces between 't', '=' and the value.
69 slopew and basew : width in points for slope fitting (points to the
70 right of clicked rupture) and base level fitting (points to
71 the left of clicked top of rupture).
72 If slopew is not provided, the routine will ask for a 5th
73 point to fit the slope.
74 DO NOT put spaces between 'slopew' or 'basew', '=' value.
76 slope2p : calculates the slope from the two clicked points,
77 instead of fitting data between them
79 justone : performs the fits over current curve instead of iterating
81 See fit command help for more information on the options and fit
83 NOTE: centerzero plot modifier should be activated (set centerzero 1).
86 #NOTE duplicates a lot of code from do_fit in fit.py, a call to it could be used directly
87 #but it is easier to control the program flow bypassing it
91 T=self.config['temperature']
97 #FIXME spaces are not allowed between pl or t and value
98 for arg in args.split():
99 #look for a persistent length argument.
101 pl_expression=arg.split('=')
102 pl_value=float(pl_expression[1]) #actual value
104 kl_expression=arg.split('=')
105 kl_value=float(kl_expression[1]) #actual value
106 #look for a T argument.
107 if ('t=' in arg) or ('T=' in arg):
108 t_expression=arg.split('=')
109 T=float(t_expression[1])
110 #look for a basew argument.
111 if ('basew=' in arg):
112 basew_expression=arg.split('=')
113 basew=int(basew_expression[1])
114 #look for a slopew argument.
115 if ('slopew=' in arg):
116 slopew_expression=arg.split('=')
117 slopew=int(slopew_expression[1])
118 if('justone' in arg):
120 if('slope2p' in arg):
128 print 'What curve no. you would like to start? (enter for ignoring)'
131 if skip.isdigit()==False:
135 print 'Skipping '+str(skip)+ ' curves'
138 #begin presenting curves for analysis
139 while curveindex <len(self.current_list):
142 curve=self.current_list[curveindex]
147 #give periodically the opportunity to stop the analysis
148 if counter%20==0 and counter>0:
149 print '\nYou checked '+str(counter)+' curves. Do you want to continue?'
163 print '\nCurve '+str(curveindex+1)+' of '+str(len(self.current_list))
164 print 'Click contact point or left end (red area)of the curve to skip'
165 #hightlights an area dedicated to reject current curve
166 sizeret=len(self.plots[0].vectors[1][0])
167 noarea=self.plots[0].vectors[1][0][-sizeret/6:-1]
168 noareaplot=copy.deepcopy(self._get_displayed_plot())
169 #noise dependent slight shift to improve visibility
170 noareaplot.add_set(noarea,np.ones_like(noarea)*5.0*np.std(self.plots[0].vectors[1][1][-sizeret/6:-1]))
171 noareaplot.styles+=['scatter']
172 noareaplot.colors+=["red"]
173 self._send_plot([noareaplot])
174 contact_point=self._measure_N_points(N=1, whatset=1)[0]
175 contact_point_index=contact_point.index
176 noareaplot.remove_set(-1)
177 #remove set leaves some styles disturbing the next graph, fixed by doing do_plot
179 if contact_point_index>5.0/6.0*sizeret:
185 self.wlccontact_point=contact_point
186 self.wlccontact_index=contact_point.index
187 self.wlccurrent=self.current.path
189 print 'Now click two points for the fitting area (one should be the rupture point)'
190 wlcpoints=self._measure_N_points(N=2,whatset=1)
191 print 'And one point of the top of the jump'
192 toppoint=self._measure_N_points(N=1,whatset=1)
193 slopepoint=ClickedPoint()
195 print 'Click a point to calculate slope'
196 slopepoint=self._measure_N_points(N=1,whatset=1)
198 fitpoints=[contact_point]+wlcpoints
199 #use the currently displayed plot for the fit
200 displayed_plot=self._get_displayed_plot()
202 #use both fit functions
204 wlcparams, wlcyfit, wlcxfit, wlcfit_errors,wlc_qstd = self.wlc_fit(fitpoints, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T, return_errors=True )
207 print 'WLC fit not possible'
211 fjcparams, fjcyfit, fjcxfit, fjcfit_errors,fjc_qstd = self.efjc_fit(fitpoints, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],kl_value,T, return_errors=True )
214 print 'eFJC fit not possible'
217 #Measure rupture force
218 ruptpoint=ClickedPoint()
219 if wlcpoints[0].graph_coords[1]<wlcpoints[1].graph_coords[1]:
220 ruptpoint=wlcpoints[0]
222 ruptpoint=wlcpoints[1]
223 tpindex=toppoint[0].index
224 toplevel=np.average(displayed_plot.vectors[1][1][tpindex:tpindex+basew])
225 force=toplevel-ruptpoint.graph_coords[1]
227 #Measure the slope - loading rate
231 xint=ruptpoint.graph_coords[0]-slopepoint[0].graph_coords[0]
232 yint=ruptpoint.graph_coords[1]-slopepoint[0].graph_coords[1]
234 slopeplot=self._get_displayed_plot(0) #get topmost displayed plot
235 slopeplot.add_set([ruptpoint.graph_coords[0],slopepoint[0].graph_coords[0]],[ruptpoint.graph_coords[1],slopepoint[0].graph_coords[1]])
236 if slopeplot.styles==[]:
237 slopeplot.styles=[None,None,'scatter']
238 slopeplot.colors=[None,None,'red']
240 slopeplot.styles+=['scatter']
241 slopeplot.colors+=['red']
243 slope=self._slope([ruptpoint]+slopepoint,slopew)
245 slope=self._slope([ruptpoint],slopew)
247 #plot results (_slope already did)
249 #now we have the fit, we can plot it
250 #add the clicked points in the final PlotObject
251 clickvector_x, clickvector_y=[], []
252 for item in wlcpoints:
253 clickvector_x.append(item.graph_coords[0])
254 clickvector_y.append(item.graph_coords[1])
256 #create a custom PlotObject to gracefully plot the fit along the curves
258 #first fix those irritating zoom-destroying fits
259 lowestpoint=min(displayed_plot.vectors[1][1])
260 for i in np.arange(0,len(wlcyfit)):
261 if wlcyfit[i] < 1.2*lowestpoint:
263 for i in np.arange(0,len(fjcyfit)):
264 if fjcyfit[i] < 1.2*lowestpoint:
267 fitplot=copy.deepcopy(displayed_plot)
268 fitplot.add_set(wlcxfit,wlcyfit)
269 fitplot.add_set(fjcxfit,fjcyfit)
270 fitplot.add_set(clickvector_x,clickvector_y)
272 fitplot.styles+=[None,None,'scatter']
273 fitplot.colors+=["red","blue",None]
275 self._send_plot([fitplot])
278 #present results of measurement
279 if len(wlcparams)==1:
280 wlcparams.append(pl_value*1e-9)
281 if len(fjcparams)==1:
282 fjcparams.append(kl_value*1e-9)
285 fjcfit_nm=[i*(10**9) for i in fjcfit_errors]
286 if len(fjcfit_nm)==1:
291 wlcfit_nm=[i*(10**9) for i in wlcfit_errors]
292 if len(wlcfit_nm)==1:
297 wlc_fitq=wlc_qstd/np.std(displayed_plot.vectors[1][1][-20:-1])
298 fjc_fitq=fjc_qstd/np.std(displayed_plot.vectors[1][1][-20:-1])
301 print 'WLC contour : '+str(1e9*wlcparams[0])+u' \u00b1 '+str(wlcfit_nm[0])+' nm'
302 print 'Per. length : '+str(1e9*wlcparams[1])+u' \u00b1 '+str(wlcfit_nm[1])+' nm'
303 print 'Quality :'+str(wlc_fitq)
305 print 'eFJC contour : '+str(1e9*fjcparams[0])+u' \u00b1 '+str(fjcfit_nm[0])+' nm'
306 print 'Kuhn length (e): '+str(1e9*fjcparams[1])+u' \u00b1 '+str(fjcfit_nm[1])+' nm'
307 print 'Quality :'+str(fjc_fitq)
309 print 'Force : '+str(1e12*force)+' pN'
310 print 'Slope : '+str(slope)+' N/m'
313 #FIXME all drivers should provide retract velocity, in SI or homogeneous units
314 lrate=slope*self.current.curve.retract_velocity
315 print 'Loading rate (UNTESTED):'+str(lrate)
317 print 'Loading rate : n/a'
323 #save accept/discard/repeat
324 print '\n\nDo you want to save these?'
328 if self.autofile=='':
329 self.autofile=raw_input('Results filename? (press enter for a random generated name)')
330 if self.autofile=='':
331 [osf,name]=tempfile.mkstemp(prefix='analysis-',suffix='.txt',text=True,dir=self.config['hookedir'])
332 print 'saving in '+name
335 f=open(self.autofile,'a+')
336 f.write('Analysis started '+time.asctime()+'\n')
337 f.write('----------------------------------------\n')
338 f.write(' File ; WLC Cont. length (nm) ; Sigma WLC cl; Per. Length (nm) ; Sigma pl ; WLC-Q ; eFJC Cont. length (nm) ; Sigma eFJC cl ; Kuhn length (nm); Sigma kl ; FJC-Q ;Force (pN) ; Slope (N/m) ; (untested) Loading rate (pN/s)\n')
341 if not os.path.exists(self.autofile):
342 f=open(self.autofile,'a+')
343 f.write('Analysis started '+time.asctime()+'\n')
344 f.write('----------------------------------------\n')
345 f.write(' File ; WLC Cont. length (nm) ; Sigma WLC cl; Per. Length (nm) ; Sigma pl; WLC-Q ;eFJC Cont. length (nm) ; Sigma eFJC cl ; Kuhn length (nm); Sigma kl ; FJC-Q ;Force (pN) ; Slope (N/m) ; (untested) Loading rate (pN/s)\n')
350 f=open(self.autofile,'a+')
351 f.write(self.current.path)
352 f.write(' ; '+str(1e9*wlcparams[0])+' ; '+str(wlcfit_nm[0])+' ; '+str(1e9*wlcparams[1])+' ; '+str(wlcfit_nm[1])+' ; '+ str(wlc_fitq)+' ; '+str(1e9*fjcparams[0])+' ; '+str(fjcfit_nm[0])+' ; '+str(1e9*fjcparams[1])+' ; '+str(fjcfit_nm[1])+' ; '+str(fjc_fitq)+' ; '+str(1e12*force)+' ; '+ str(slope)+' ; '+str(lrate)+'\n')
355 print '\nWould you like to try again on this same curve?'
361 print 'You saved '+str(savecounter)+' out of '+str(len(self.current_list))+' total curves.'
366 #shows something in the graph to indicate we expect YNclick
367 #FIXME is there an easy way to write in the graph?
368 askplot=self._get_displayed_plot()
369 center=np.min(askplot.vectors[1][0])+15e-9
370 scale=np.std(askplot.vectors[1][1][-50:-1])
371 hscale=np.mean(np.diff(askplot.vectors[1][0][-10:-1]))
372 xhline=[center-20*hscale,center-10*hscale,center+10*hscale,center+20*hscale]
373 ytopline=[15*scale,20*scale,20*scale,15*scale]
374 ybotline=[-15*scale,-20*scale,-20*scale,-15*scale]
375 yvline=np.arange(-25*scale,26*scale,scale/2)
376 xvline=np.ones_like(yvline)*center
377 xshow=np.concatenate((xvline,xhline,xhline))
378 yshow=np.concatenate((yvline,ytopline,ybotline))
379 askplot.add_set(xshow,yshow)
380 askplot.styles+=['scatter']
381 askplot.colors+=["black"]
382 self._send_plot([askplot])
383 print 'Click any point over y=0 for Yes, under it for No'
384 point=self._measure_N_points(N=1,whatset=1)[0]
385 value=point.absolute_coords[1]
386 askplot.remove_set(-1)
387 self._send_plot([askplot])