3 #FIXME clean this, probably some dependencies are not needed
5 from ..libhooke import WX_GOOD, ClickedPoint
7 wxversion.select(WX_GOOD)
8 from wx import PostEvent
16 warnings.simplefilter('ignore',np.RankWarning)
20 global EVT_MEASURE_WLC
22 global events_from_fit
23 events_from_fit=Queue.Queue() #GUI ---> CLI COMMUNICATION
26 class multifitCommands:
28 def do_multifit(self,args):
32 Presents curves for manual analysis in a comfortable mouse-only fashion.
33 Obtains contour length, persistance length, rupture force and
35 WLC is shown in red, eFJC in blue.
38 multifit [pl=value] [kl=value] [t=value] [slopew=value] [basew=value]
41 pl=[value] and kl=[value]: Use a fixed persistent length (WLC) or Kuhn
42 length (eFJC) for the fit. If pl is not given, the fit will be
44 DO NOT put spaces between 'pl', '=' and the value.
45 eFJC fit works better with a fixed kl
46 The value must be in nanometers.
48 t=[value] : Use a user-defined temperature. The value must be in
49 kelvins; by default it is 293 K.
50 DO NOT put spaces between 't', '=' and the value.
52 slopew and basew : width in points for slope fitting (points to the
53 right of clicked rupture) and base level fitting (points to
54 the left of clicked top of rupture).
55 If slopew is not provided, the routine will ask for a 5th
56 point to fit the slope.
57 DO NOT put spaces between 'slopew' or 'basew', '=' value.
59 slope2p : calculates the slope from the two clicked points,
60 instead of fitting data between them
62 justone : performs the fits over current curve instead of iterating
64 See fit command help for more information on the options and fit
66 NOTE: centerzero plot modifier should be activated (set centerzero 1).
69 #NOTE duplicates a lot of code from do_fit in fit.py, a call to it could be used directly
70 #but it is easier to control the program flow bypassing it
74 T=self.config['temperature']
80 #FIXME spaces are not allowed between pl or t and value
81 for arg in args.split():
82 #look for a persistent length argument.
84 pl_expression=arg.split('=')
85 pl_value=float(pl_expression[1]) #actual value
87 kl_expression=arg.split('=')
88 kl_value=float(kl_expression[1]) #actual value
89 #look for a T argument.
90 if ('t=' in arg) or ('T=' in arg):
91 t_expression=arg.split('=')
92 T=float(t_expression[1])
93 #look for a basew argument.
95 basew_expression=arg.split('=')
96 basew=int(basew_expression[1])
97 #look for a slopew argument.
98 if ('slopew=' in arg):
99 slopew_expression=arg.split('=')
100 slopew=int(slopew_expression[1])
101 if('justone' in arg):
103 if('slope2p' in arg):
111 print 'What curve no. you would like to start? (enter for ignoring)'
114 if skip.isdigit()==False:
118 print 'Skipping '+str(skip)+ ' curves'
121 #begin presenting curves for analysis
122 while curveindex <len(self.current_list):
125 curve=self.current_list[curveindex]
130 #give periodically the opportunity to stop the analysis
131 if counter%20==0 and counter>0:
132 print '\nYou checked '+str(counter)+' curves. Do you want to continue?'
146 print '\nCurve '+str(curveindex+1)+' of '+str(len(self.current_list))
147 print 'Click contact point or left end (red area)of the curve to skip'
148 #hightlights an area dedicated to reject current curve
149 sizeret=len(self.plots[0].vectors[1][0])
150 noarea=self.plots[0].vectors[1][0][-sizeret/6:-1]
151 noareaplot=copy.deepcopy(self._get_displayed_plot())
152 #noise dependent slight shift to improve visibility
153 noareaplot.add_set(noarea,np.ones_like(noarea)*5.0*np.std(self.plots[0].vectors[1][1][-sizeret/6:-1]))
154 noareaplot.styles+=['scatter']
155 noareaplot.colors+=["red"]
156 self._send_plot([noareaplot])
157 contact_point=self._measure_N_points(N=1, whatset=1)[0]
158 contact_point_index=contact_point.index
159 noareaplot.remove_set(-1)
160 #remove set leaves some styles disturbing the next graph, fixed by doing do_plot
162 if contact_point_index>5.0/6.0*sizeret:
168 self.wlccontact_point=contact_point
169 self.wlccontact_index=contact_point.index
170 self.wlccurrent=self.current.path
172 print 'Now click two points for the fitting area (one should be the rupture point)'
173 wlcpoints=self._measure_N_points(N=2,whatset=1)
174 print 'And one point of the top of the jump'
175 toppoint=self._measure_N_points(N=1,whatset=1)
176 slopepoint=ClickedPoint()
178 print 'Click a point to calculate slope'
179 slopepoint=self._measure_N_points(N=1,whatset=1)
181 fitpoints=[contact_point]+wlcpoints
182 #use the currently displayed plot for the fit
183 displayed_plot=self._get_displayed_plot()
185 #use both fit functions
187 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 )
190 print 'WLC fit not possible'
194 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 )
197 print 'eFJC fit not possible'
200 #Measure rupture force
201 ruptpoint=ClickedPoint()
202 if wlcpoints[0].graph_coords[1]<wlcpoints[1].graph_coords[1]:
203 ruptpoint=wlcpoints[0]
205 ruptpoint=wlcpoints[1]
206 tpindex=toppoint[0].index
207 toplevel=np.average(displayed_plot.vectors[1][1][tpindex:tpindex+basew])
208 force=toplevel-ruptpoint.graph_coords[1]
210 #Measure the slope - loading rate
214 xint=ruptpoint.graph_coords[0]-slopepoint[0].graph_coords[0]
215 yint=ruptpoint.graph_coords[1]-slopepoint[0].graph_coords[1]
217 slopeplot=self._get_displayed_plot(0) #get topmost displayed plot
218 slopeplot.add_set([ruptpoint.graph_coords[0],slopepoint[0].graph_coords[0]],[ruptpoint.graph_coords[1],slopepoint[0].graph_coords[1]])
219 if slopeplot.styles==[]:
220 slopeplot.styles=[None,None,'scatter']
221 slopeplot.colors=[None,None,'red']
223 slopeplot.styles+=['scatter']
224 slopeplot.colors+=['red']
226 slope=self._slope([ruptpoint]+slopepoint,slopew)
228 slope=self._slope([ruptpoint],slopew)
230 #plot results (_slope already did)
232 #now we have the fit, we can plot it
233 #add the clicked points in the final PlotObject
234 clickvector_x, clickvector_y=[], []
235 for item in wlcpoints:
236 clickvector_x.append(item.graph_coords[0])
237 clickvector_y.append(item.graph_coords[1])
239 #create a custom PlotObject to gracefully plot the fit along the curves
241 #first fix those irritating zoom-destroying fits
242 lowestpoint=min(displayed_plot.vectors[1][1])
243 for i in np.arange(0,len(wlcyfit)):
244 if wlcyfit[i] < 1.2*lowestpoint:
246 for i in np.arange(0,len(fjcyfit)):
247 if fjcyfit[i] < 1.2*lowestpoint:
250 fitplot=copy.deepcopy(displayed_plot)
251 fitplot.add_set(wlcxfit,wlcyfit)
252 fitplot.add_set(fjcxfit,fjcyfit)
253 fitplot.add_set(clickvector_x,clickvector_y)
255 fitplot.styles+=[None,None,'scatter']
256 fitplot.colors+=["red","blue",None]
258 self._send_plot([fitplot])
261 #present results of measurement
262 if len(wlcparams)==1:
263 wlcparams.append(pl_value*1e-9)
264 if len(fjcparams)==1:
265 fjcparams.append(kl_value*1e-9)
268 fjcfit_nm=[i*(10**9) for i in fjcfit_errors]
269 if len(fjcfit_nm)==1:
274 wlcfit_nm=[i*(10**9) for i in wlcfit_errors]
275 if len(wlcfit_nm)==1:
280 wlc_fitq=wlc_qstd/np.std(displayed_plot.vectors[1][1][-20:-1])
281 fjc_fitq=fjc_qstd/np.std(displayed_plot.vectors[1][1][-20:-1])
284 print 'WLC contour : '+str(1e9*wlcparams[0])+u' \u00b1 '+str(wlcfit_nm[0])+' nm'
285 print 'Per. length : '+str(1e9*wlcparams[1])+u' \u00b1 '+str(wlcfit_nm[1])+' nm'
286 print 'Quality :'+str(wlc_fitq)
288 print 'eFJC contour : '+str(1e9*fjcparams[0])+u' \u00b1 '+str(fjcfit_nm[0])+' nm'
289 print 'Kuhn length (e): '+str(1e9*fjcparams[1])+u' \u00b1 '+str(fjcfit_nm[1])+' nm'
290 print 'Quality :'+str(fjc_fitq)
292 print 'Force : '+str(1e12*force)+' pN'
293 print 'Slope : '+str(slope)+' N/m'
296 #FIXME all drivers should provide retract velocity, in SI or homogeneous units
297 lrate=slope*self.current.curve.retract_velocity
298 print 'Loading rate (UNTESTED):'+str(lrate)
300 print 'Loading rate : n/a'
306 #save accept/discard/repeat
307 print '\n\nDo you want to save these?'
311 if self.autofile=='':
312 self.autofile=raw_input('Results filename? (press enter for a random generated name)')
313 if self.autofile=='':
314 [osf,name]=tempfile.mkstemp(prefix='analysis-',suffix='.txt',text=True,dir=self.config['hookedir'])
315 print 'saving in '+name
318 f=open(self.autofile,'a+')
319 f.write('Analysis started '+time.asctime()+'\n')
320 f.write('----------------------------------------\n')
321 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')
324 if not os.path.exists(self.autofile):
325 f=open(self.autofile,'a+')
326 f.write('Analysis started '+time.asctime()+'\n')
327 f.write('----------------------------------------\n')
328 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')
333 f=open(self.autofile,'a+')
334 f.write(self.current.path)
335 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')
338 print '\nWould you like to try again on this same curve?'
344 print 'You saved '+str(savecounter)+' out of '+str(len(self.current_list))+' total curves.'
349 #shows something in the graph to indicate we expect YNclick
350 #FIXME is there an easy way to write in the graph?
351 askplot=self._get_displayed_plot()
352 center=np.min(askplot.vectors[1][0])+15e-9
353 scale=np.std(askplot.vectors[1][1][-50:-1])
354 hscale=np.mean(np.diff(askplot.vectors[1][0][-10:-1]))
355 xhline=[center-20*hscale,center-10*hscale,center+10*hscale,center+20*hscale]
356 ytopline=[15*scale,20*scale,20*scale,15*scale]
357 ybotline=[-15*scale,-20*scale,-20*scale,-15*scale]
358 yvline=np.arange(-25*scale,26*scale,scale/2)
359 xvline=np.ones_like(yvline)*center
360 xshow=np.concatenate((xvline,xhline,xhline))
361 yshow=np.concatenate((yvline,ytopline,ybotline))
362 askplot.add_set(xshow,yshow)
363 askplot.styles+=['scatter']
364 askplot.colors+=["black"]
365 self._send_plot([askplot])
366 print 'Click any point over y=0 for Yes, under it for No'
367 point=self._measure_N_points(N=1,whatset=1)[0]
368 value=point.absolute_coords[1]
369 askplot.remove_set(-1)
370 self._send_plot([askplot])