3 """Force spectroscopy curves basic fitting plugin.
5 Non-standard Dependencies:
6 procplots.py (plot processing plugin)
9 from ..libhooke import WX_GOOD, ClickedPoint
12 wxversion.select(WX_GOOD)
13 #from wx import PostEvent
14 #from wx.lib.newevent import NewEvent
23 global EVT_MEASURE_WLC
25 #measure_wlc, EVT_MEASURE_WLC = NewEvent()
27 global events_from_fit
28 events_from_fit=Queue.Queue() #GUI ---> CLI COMMUNICATION
31 class fitCommands(object):
35 self.wlccontact_point=None
36 self.wlccontact_index=None
39 '''Calculates the average distance from data to fit, scaled by the standard deviation
40 of the free cantilever area (thermal noise)
45 def wlc_fit(self,clicked_points,xvector,yvector, pl_value, T=293, return_errors=False):
47 Worm-like chain model fitting.
48 The function is the simple polynomial worm-like chain as proposed by C.Bustamante, J.F.Marko, E.D.Siggia
49 and S.Smith (Science. 1994 Sep 9;265(5178):1599-600.)
53 clicked_points[0] = contact point (calculated or hand-clicked)
54 clicked_points[1] and [2] are edges of chunk
57 #STEP 1: Prepare the vectors to apply the fit.
58 if pl_value is not None:
59 pl_value=pl_value/(10**9)
61 #indexes of the selected chunk
62 first_index=min(clicked_points[1].index, clicked_points[2].index)
63 last_index=max(clicked_points[1].index, clicked_points[2].index)
65 #getting the chunk and reverting it
66 xchunk,ychunk=xvector[first_index:last_index],yvector[first_index:last_index]
69 #put contact point at zero and flip around the contact point (the fit wants a positive growth for extension and force)
70 xchunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xchunk]
71 ychunk_corr_up=[-(y-clicked_points[0].graph_coords[1]) for y in ychunk]
74 xchunk_corr_up=scipy.array(xchunk_corr_up)
75 ychunk_corr_up=scipy.array(ychunk_corr_up)
78 #STEP 2: actually do the fit
80 #Find furthest point of chunk and add it a bit; the fit must converge
82 xchunk_high=max(xchunk_corr_up)
83 xchunk_high+=(xchunk_high/10)
85 #Here are the linearized start parameters for the WLC.
86 #[lambd=1/Lo , pii=1/P]
88 p0=[(1/xchunk_high),(1/(3.5e-10))]
89 p0_plfix=[(1/xchunk_high)]
92 fixme: remove these comments after testing
94 def dist(px,py,linex,liney):
95 distancesx=scipy.array([(px-x)**2 for x in linex])
96 minindex=np.argmin(distancesx)
97 print px, linex[0], linex[-1]
98 return (py-liney[minindex])**2
100 def f_wlc(params,x,T=T):
102 wlc function for ODR fitting
107 y=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
110 def f_wlc_plfix(params,x,pl_value=pl_value,T=T):
112 wlc function for ODR fitting
118 y=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
122 realdata=scipy.odr.RealData(xchunk_corr_up,ychunk_corr_up)
124 model=scipy.odr.Model(f_wlc_plfix)
125 o = scipy.odr.ODR(realdata, model, p0_plfix)
127 model=scipy.odr.Model(f_wlc)
128 o = scipy.odr.ODR(realdata, model, p0)
130 o.set_job(fit_type=2)
132 fit_out=[(1/i) for i in out.beta]
134 #Calculate fit errors from output standard deviations.
135 #We must propagate the error because we fit the *inverse* parameters!
136 #The error = (error of the inverse)*(value**2)
138 for sd,value in zip(out.sd_beta, fit_out):
139 err_real=sd*(value**2)
140 fit_errors.append(err_real)
142 def wlc_eval(x,params,pl_value,T):
144 Evaluates the WLC function
154 Kb=(1.38065e-23) #boltzmann constant
155 therm=Kb*T #so we have thermal energy
157 return ( (therm*pii/4.0) * (((1-(x*lambd))**-2.0) - 1 + (4.0*x*lambd)) )
160 #STEP 3: plotting the fit
162 #obtain domain to plot the fit - from contact point to last_index plus 20 points
163 thule_index=last_index+10
164 if thule_index > len(xvector): #for rare cases in which we fit something at the END of whole curve.
165 thule_index = len(xvector)
166 #reverse etc. the domain
167 xfit_chunk=xvector[clicked_points[0].index:thule_index]
169 xfit_chunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xfit_chunk]
170 xfit_chunk_corr_up=scipy.array(xfit_chunk_corr_up)
172 #the fitted curve: reflip, re-uncorrect
173 yfit=wlc_eval(xfit_chunk_corr_up, out.beta, pl_value,T)
174 yfit_down=[-y for y in yfit]
175 yfit_corr_down=[y+clicked_points[0].graph_coords[1] for y in yfit_down]
178 #calculate fit quality
180 yqeval=wlc_eval(xchunk_corr_up,out.beta,pl_value,T)
181 #we need to cut the extra from thuleindex
182 for qindex in np.arange(0,len(yqeval)):
183 qsum+=(yqeval[qindex]-ychunk_corr_up[qindex])**2
184 qstd=np.sqrt(qsum/len(ychunk_corr_up))
188 return fit_out, yfit_corr_down, xfit_chunk, fit_errors, qstd
190 return fit_out, yfit_corr_down, xfit_chunk, None, qstd
192 def fjc_fit(self,clicked_points,xvector,yvector, pl_value, T=293, return_errors=False):
194 Freely-jointed chain function
195 ref: C.Ray and B.B. Akhremitchev; http://www.chem.duke.edu/~boris/research/force_spectroscopy/fit_efjc.pdf
198 clicked_points[0] is the contact point (calculated or hand-clicked)
199 clicked_points[1] and [2] are edges of chunk
202 #STEP 1: Prepare the vectors to apply the fit.
203 if pl_value is not None:
204 pl_value=pl_value/(10**9)
206 #indexes of the selected chunk
207 first_index=min(clicked_points[1].index, clicked_points[2].index)
208 last_index=max(clicked_points[1].index, clicked_points[2].index)
210 #getting the chunk and reverting it
211 xchunk,ychunk=xvector[first_index:last_index],yvector[first_index:last_index]
214 #put contact point at zero and flip around the contact point (the fit wants a positive growth for extension and force)
215 xchunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xchunk]
216 ychunk_corr_up=[-(y-clicked_points[0].graph_coords[1]) for y in ychunk]
220 xchunk_corr_up=scipy.array(xchunk_corr_up)
221 ychunk_corr_up=scipy.array(ychunk_corr_up)
224 #STEP 2: actually do the fit
226 #Find furthest point of chunk and add it a bit; the fit must converge
228 xchunk_high=max(xchunk_corr_up)
229 xchunk_high+=(xchunk_high/10)
231 #Here are the linearized start parameters for the WLC.
232 #[lambd=1/Lo , pii=1/P]
234 p0=[(1/xchunk_high),(1/(3.5e-10))]
235 p0_plfix=[(1/xchunk_high)]
238 fixme: remove these comments after testing
240 def dist(px,py,linex,liney):
241 distancesx=scipy.array([(px-x)**2 for x in linex])
242 minindex=np.argmin(distancesx)
243 print minindex, px, linex[0], linex[-1]
244 return (py-liney[minindex])**2
250 return (np.exp(2*z)+1)/(np.exp(2*z)-1)
252 def x_fjc(params,f,T=T):
254 fjc function for ODR fitting
260 #x=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
261 x=(1/lambd)*(coth(f*(1/pii)/therm) - (therm*pii)/f)
264 def x_fjc_plfix(params,f,pl_value=pl_value,T=T):
266 fjc function for ODR fitting
272 #y=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
273 x=(1/lambd)*(coth(f*(1/pii)/therm) - (therm*pii)/f)
277 realdata=scipy.odr.RealData(ychunk_corr_up,xchunk_corr_up)
279 model=scipy.odr.Model(x_fjc_plfix)
280 o = scipy.odr.ODR(realdata, model, p0_plfix)
282 model=scipy.odr.Model(x_fjc)
283 o = scipy.odr.ODR(realdata, model, p0)
285 o.set_job(fit_type=2)
287 fit_out=[(1/i) for i in out.beta]
289 #Calculate fit errors from output standard deviations.
290 #We must propagate the error because we fit the *inverse* parameters!
291 #The error = (error of the inverse)*(value**2)
293 for sd,value in zip(out.sd_beta, fit_out):
294 err_real=sd*(value**2)
295 fit_errors.append(err_real)
297 def fjc_eval(y,params,pl_value,T):
299 Evaluates the WLC function
309 Kb=(1.38065e-23) #boltzmann constant
310 therm=Kb*T #so we have thermal energy
311 #return ( (therm*pii/4.0) * (((1-(x*lambd))**-2.0) - 1 + (4.0*x*lambd)) )
312 return (1/lambd)*(coth(y*(1/pii)/therm) - (therm*pii)/y)
316 #STEP 3: plotting the fit
317 #obtain domain to plot the fit - from contact point to last_index plus 20 points
318 thule_index=last_index+10
319 if thule_index > len(xvector): #for rare cases in which we fit something at the END of whole curve.
320 thule_index = len(xvector)
321 #reverse etc. the domain
322 ychunk=yvector[clicked_points[0].index:thule_index]
325 y_evalchunk=np.linspace(min(ychunk),max(ychunk),100)
327 #Empty y-chunk. It happens whenever we set the contact point after a recognized peak,
328 #or other buggy situations. Kludge to live with it now...
329 ychunk=yvector[:thule_index]
330 y_evalchunk=np.linspace(min(ychunk),max(ychunk),100)
332 yfit_down=[-y for y in y_evalchunk]
333 yfit_corr_down=[y+clicked_points[0].graph_coords[1] for y in yfit_down]
334 yfit_corr_down=scipy.array(yfit_corr_down)
336 #the fitted curve: reflip, re-uncorrect
337 xfit=fjc_eval(yfit_corr_down, out.beta, pl_value,T)
340 xfit_chunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xfit]
342 #xfit_chunk_corr_up=scipy.array(xfit_chunk_corr_up)
343 #deltay=yfit_down[0]-yvector[clicked_points[0].index]
345 #This is a terrible, terrible kludge to find the point where it should normalize (and from where it should plot)
347 for index in scipy.arange(1,len(xfit_chunk_corr_up),1):
348 xxxdists.append((clicked_points[0].graph_coords[0]-xfit_chunk_corr_up[index])**2)
349 normalize_index=xxxdists.index(min(xxxdists))
352 deltay=yfit_down[normalize_index]-clicked_points[0].graph_coords[1]
353 yfit_corr_down=[y-deltay for y in yfit_down]
356 #calculate fit quality
357 #creates dense y vector
358 yqeval=np.linspace(np.min(ychunk_corr_up)/2,np.max(ychunk_corr_up)*2,10*len(ychunk_corr_up))
359 #corresponding fitted x
360 xqeval=fjc_eval(yqeval,out.beta,pl_value,T)
363 for qindex in np.arange(0,len(ychunk_corr_up)):
364 qsum+=dist(xchunk_corr_up[qindex],ychunk_corr_up[qindex],xqeval,yqeval)
365 qstd=np.sqrt(qsum/len(ychunk_corr_up))
369 return fit_out, yfit_corr_down[normalize_index+1:], xfit_chunk_corr_up[normalize_index+1:], fit_errors, qstd
371 return fit_out, yfit_corr_down[normalize_index+1:], xfit_chunk_corr_up[normalize_index+1:], None, qstd
373 def efjc_fit(self,clicked_points,xvector,yvector, pl_value, T=293.0, return_errors=False):
375 Extended Freely-jointed chain function
376 ref: F Oesterhelt, M Rief and H E Gaub, New Journal of Physics 1 (1999) 6.1–6.11
377 Please note that 2-parameter fitting (length and kl) usually does not converge, use fixed kl
380 clicked_points[0] is the contact point (calculated or hand-clicked)
381 clicked_points[1] and [2] are edges of chunk
384 #Fixed parameters from reference
385 Kb=(1.38065e-2) #in pN.nm
387 Lh=0.280 #helical, nm
391 #STEP 1: Prepare the vectors to apply the fit.
393 #indexes of the selected chunk
394 first_index=min(clicked_points[1].index, clicked_points[2].index)
395 last_index=max(clicked_points[1].index, clicked_points[2].index)
397 #getting the chunk and reverting it
398 xchunk,ychunk=xvector[first_index:last_index],yvector[first_index:last_index]
401 #put contact point at zero and flip around the contact point (the fit wants a positive growth for extension and force)
402 xchunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xchunk]
403 ychunk_corr_up=[-(y-clicked_points[0].graph_coords[1]) for y in ychunk]
407 xchunk_corr_up=scipy.array(xchunk_corr_up)
408 ychunk_corr_up=scipy.array(ychunk_corr_up)
410 xchunk_corr_up_nm=xchunk_corr_up*1e9
411 ychunk_corr_up_pn=ychunk_corr_up*1e12
414 #STEP 2: actually do the fit
416 #Find furthest point of chunk and add it a bit; the fit must converge
418 xchunk_high=max(xchunk_corr_up_nm)
419 xchunk_high+=(xchunk_high/10.0)
421 #Here are the linearized start parameters for the WLC.
423 #max number of monomers (all helical)for a contour length xchunk_high
424 excessNs=xchunk_high/(Lp)
425 p0=[excessNs,(1.0/(0.7))]
426 p0_plfix=[(excessNs)]
428 def dist(px,py,linex,liney):
429 distancesx=scipy.array([(px-x)**2 for x in linex])
430 minindex=np.argmin(distancesx)
431 return (py-liney[minindex])**2
434 dG0=12.4242 #3kt in pN.nm
435 dL=0.078 #planar-helical
440 Lh=0.280 #helical, nm
445 return Lp/(np.exp(dG/therm)+1)+Lh/(np.exp(-dG/therm)+1)
451 return 1.0/np.tanh(z)
453 def x_efjc(params,f,T=T,Ks=Ks):
455 efjc function for ODR fitting
464 x=Ns*Lfactor(f)*(coth(beta)-1.0/beta)+Ns*f/Ks
467 def x_efjc_plfix(params,f,kl_value=pl_value,T=T,Ks=Ks):
469 efjc function for ODR fitting
478 x=Ns*Lfactor(f)*(coth(beta)-1.0/beta)+Ns*f/Ks
482 realdata=scipy.odr.RealData(ychunk_corr_up_pn,xchunk_corr_up_nm)
484 model=scipy.odr.Model(x_efjc_plfix)
485 o = scipy.odr.ODR(realdata, model, p0_plfix)
487 print 'WARNING eFJC fit with free pl sometimes does not converge'
488 model=scipy.odr.Model(x_efjc)
489 o = scipy.odr.ODR(realdata, model, p0)
491 o.set_job(fit_type=2)
498 kfit=1e-9/out.beta[1]
506 #Calculate fit errors from output standard deviations.
508 fit_errors.append(out.sd_beta[0]*Lp*1e-9)
510 fit_errors.append(1e9*out.sd_beta[1]*kfit**2)
514 def efjc_eval(y,params,pl_value,T=T,Lfactor=Lfactor,Ks=Ks):
516 Evaluates the eFJC function
526 Kb=(1.38065e-2) #boltzmann constant
527 therm=Kb*T #so we have thermal energy
530 x= Ns*Lfactor(y)*(coth(beta)-1.0/beta)+Ns*y/Ks
534 #STEP 3: plotting the fit
535 #obtain domain to plot the fit - from contact point to last_index plus 20 points
536 thule_index=last_index+10
537 if thule_index > len(xvector): #for rare cases in which we fit something at the END of whole curve.
538 thule_index = len(xvector)
539 #reverse etc. the domain
540 ychunk=yvector[clicked_points[0].index:thule_index]
543 y_evalchunk=np.linspace(min(ychunk),max(ychunk),100)
545 #Empty y-chunk. It happens whenever we set the contact point after a recognized peak,
546 #or other buggy situations. Kludge to live with it now...
547 ychunk=yvector[:thule_index]
548 y_evalchunk=np.linspace(min(ychunk),max(ychunk),100)
550 yfit_down=[-y for y in y_evalchunk]
551 yfit_corr_down=[y+clicked_points[0].graph_coords[1] for y in yfit_down]
552 yfit_corr_down=scipy.array(yfit_corr_down)
554 #the fitted curve: reflip, re-uncorrect
555 xfit=efjc_eval(1e12*yfit_corr_down, out.beta, pl_value,T)*1e-9
558 xfit_chunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xfit]
560 #xfit_chunk_corr_up=scipy.array(xfit_chunk_corr_up)
561 #deltay=yfit_down[0]-yvector[clicked_points[0].index]
563 #This is a terrible, terrible kludge to find the point where it should normalize (and from where it should plot)
565 for index in scipy.arange(1,len(xfit_chunk_corr_up),1):
566 xxxdists.append((clicked_points[0].graph_coords[0]-xfit_chunk_corr_up[index])**2)
567 normalize_index=xxxdists.index(min(xxxdists))
570 deltay=yfit_down[normalize_index]-clicked_points[0].graph_coords[1]
571 yfit_corr_down=[y-deltay for y in yfit_down]
573 #calculate fit quality
574 #creates dense y vector
575 yqeval=np.linspace(np.min(ychunk_corr_up_pn)/2,np.max(ychunk_corr_up_pn)*2,10*len(ychunk_corr_up_pn))
576 #corresponding fitted x
577 xqeval=efjc_eval(yqeval,out.beta,pl_value)
580 for qindex in np.arange(0,len(ychunk_corr_up_pn)):
581 qsum+=dist(xchunk_corr_up_nm[qindex],ychunk_corr_up_pn[qindex],xqeval,yqeval)
582 qstd=1e-12*np.sqrt(qsum/len(ychunk_corr_up_pn))
585 return fit_out, yfit_corr_down[normalize_index+1:], xfit_chunk_corr_up[normalize_index+1:], fit_errors, qstd
587 return fit_out, yfit_corr_down[normalize_index+1:], xfit_chunk_corr_up[normalize_index+1:], None, qstd
592 def do_wlc(self,args):
601 def do_fjc(self,args):
610 def do_fit(self,args):
614 Fits an entropic elasticity function to a given chunk of the curve.
616 First you have to click a contact point.
617 Then you have to click the two edges of the data you want to fit.
619 Fit quality compares the distance to the fit with the thermal noise (a good fit should be close to 1)
621 The fit function depends on the fit_function variable. You can set it with the command
622 "set fit_function wlc" or "set fit_function fjc" depending on the function you prefer.
624 For WLC, the function is the simple polynomial worm-like chain as proposed by
625 C.Bustamante, J.F.Marko, E.D.Siggia and S.Smith (Science. 1994
626 Sep 9;265(5178):1599-600.)
629 C.Ray and B.B. Akhremitchev; http://www.chem.duke.edu/~boris/research/force_spectroscopy/fit_efjc.pdf
632 F Oesterhelt, M Rief and H E Gaub, New Journal of Physics 1 (1999) 6.1–6.11 (section 4.2)
633 NOTE: use fixed pl for better results.
636 pl=[value] : Use a fixed persistent length (WLC) or Kuhn length (FJC) for the fit. If pl is not given,
637 the fit will be a 2-variable
638 fit. DO NOT put spaces between 'pl', '=' and the value.
639 The value must be in nanometers.
641 t=[value] : Use a user-defined temperature. The value must be in
642 kelvins; by default it is 293 K.
643 DO NOT put spaces between 't', '=' and the value.
645 noauto : allows for clicking the contact point by
646 hand (otherwise it is automatically estimated) the first time.
647 If subsequent measurements are made, the same contact point
650 reclick : redefines by hand the contact point, if noauto has been used before
651 but the user is unsatisfied of the previously choosen contact point.
653 Syntax: fit [pl=(value)] [t=value] [noauto]
656 T=self.config['temperature']
657 for arg in args.split():
658 #look for a persistent length argument.
660 pl_expression=arg.split('=')
661 pl_value=float(pl_expression[1]) #actual value
662 #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
663 if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
664 t_expression=arg.split('=')
665 T=float(t_expression[1])
667 #use the currently displayed plot for the fit
668 displayed_plot=self._get_displayed_plot()
670 #handle contact point arguments correctly
671 if 'reclick' in args.split():
672 print 'Click contact point'
673 contact_point=self._measure_N_points(N=1, whatset=1)[0]
674 contact_point_index=contact_point.index
675 self.wlccontact_point=contact_point
676 self.wlccontact_index=contact_point.index
677 self.wlccurrent=self.current.path
678 elif 'noauto' in args.split():
679 if self.wlccontact_index==None or self.wlccurrent != self.current.path:
680 print 'Click contact point'
681 contact_point=self._measure_N_points(N=1, whatset=1)[0]
682 contact_point_index=contact_point.index
683 self.wlccontact_point=contact_point
684 self.wlccontact_index=contact_point.index
685 self.wlccurrent=self.current.path
687 contact_point=self.wlccontact_point
688 contact_point_index=self.wlccontact_index
690 cindex=self.find_contact_point()
691 contact_point=ClickedPoint()
692 contact_point.absolute_coords=displayed_plot.vectors[1][0][cindex], displayed_plot.vectors[1][1][cindex]
693 contact_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
694 contact_point.is_marker=True
696 print 'Click edges of chunk'
697 points=self._measure_N_points(N=2, whatset=1)
698 points=[contact_point]+points
701 if self.config['fit_function']=='wlc':
702 params, yfit, xfit, fit_errors,qstd = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T, return_errors=True )
703 name_of_charlength='Persistent length'
704 elif self.config['fit_function']=='fjc':
705 params, yfit, xfit, fit_errors,qstd = self.fjc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T, return_errors=True )
706 name_of_charlength='Kuhn length'
707 elif self.config['fit_function']=='efjc':
708 params, yfit, xfit, fit_errors,qstd = self.efjc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T, return_errors=True )
709 name_of_charlength='Kuhn length (e)'
711 print 'No recognized fit function defined!'
712 print 'Set your fit function to wlc, fjc or efjc.'
716 print 'Fit not possible. Probably wrong interval -did you click two *different* points?'
719 #FIXME: print "Kuhn length" for FJC
720 print 'Fit function:',self.config['fit_function']
721 print 'Contour length: ',params[0]*(1.0e+9),' nm'
722 to_dump='contour '+self.current.path+' '+str(params[0]*(1.0e+9))+' nm'
723 self.outlet.push(to_dump)
724 if len(params)==2: #if we did choose 2-value fit
725 print name_of_charlength+': ',params[1]*(1.0e+9),' nm'
726 to_dump='persistent '+self.current.path+' '+str(params[1]*(1.0e+9))+' nm'
727 self.outlet.push(to_dump)
730 fit_nm=[i*(10**9) for i in fit_errors]
731 print 'Standard deviation (contour length)', fit_nm[0]
733 print 'Standard deviation ('+name_of_charlength+')', fit_nm[1]
735 print 'Fit quality: '+str(qstd/np.std(displayed_plot.vectors[1][1][-20:-1]))
738 #add the clicked points in the final PlotObject
739 clickvector_x, clickvector_y=[], []
741 clickvector_x.append(item.graph_coords[0])
742 clickvector_y.append(item.graph_coords[1])
744 #create a custom PlotObject to gracefully plot the fit along the curves
746 fitplot=copy.deepcopy(displayed_plot)
747 fitplot.add_set(xfit,yfit)
748 fitplot.add_set(clickvector_x,clickvector_y)
750 #FIXME: this colour/styles stuff must be solved at the root!
751 if fitplot.styles==[]:
752 fitplot.styles=[None,None,None,'scatter']
754 fitplot.styles+=[None,'scatter']
756 if fitplot.colors==[]:
757 fitplot.colors=[None,None,None,None]
759 fitplot.colors+=[None,None]
761 self._send_plot([fitplot])
768 def find_contact_point(self,plot=False):
770 Finds the contact point on the curve.
772 The current algorithm (thanks to Francesco Musiani, francesco.musiani@unibo.it and Massimo Sandal) is:
773 - take care of the PicoForce trigger bug - exclude retraction portions with too high standard deviation
774 - fit the second half of the retraction curve to a line
775 - if the fit is not almost horizontal, take a smaller chunk and repeat
776 - otherwise, we have something horizontal
777 - so take the average of horizontal points and use it as a baseline
779 Then, start from the rise of the retraction curve and look at the first point below the
782 FIXME: should be moved, probably to generalvclamp.py
788 outplot=self.subtract_curves(1)
789 xret=outplot.vectors[1][0]
790 ydiff=outplot.vectors[1][1]
792 xext=plot.vectors[0][0]
793 yext=plot.vectors[0][1]
794 xret2=plot.vectors[1][0]
795 yret=plot.vectors[1][1]
797 #taking care of the picoforce trigger bug: we exclude portions of the curve that have too much
798 #standard deviation. yes, a lot of magic is here.
800 monlength=len(xret)-int(len(xret)/20)
803 monchunk=scipy.array(ydiff[monlength:finalength])
804 if abs(np.std(monchunk)) < 2e-10:
806 else: #move away from the monster
807 monlength-=int(len(xret)/50)
808 finalength-=int(len(xret)/50)
811 #take half of the thing
812 endlength=int(len(xret)/2)
817 xchunk=yext[endlength:monlength]
818 ychunk=yext[endlength:monlength]
819 regr=scipy.stats.linregress(xchunk,ychunk)[0:2]
820 #we stop if we found an almost-horizontal fit or if we're going too short...
821 #FIXME: 0.1 and 6 here are "magic numbers" (although reasonable)
822 if (abs(regr[1]) > 0.1) and ( endlength < len(xret)-int(len(xret)/6) ) :
828 ymean=np.mean(ychunk) #baseline
833 #find the first point below the calculated baseline
839 #The algorithm didn't find anything below the baseline! It should NEVER happen
847 def find_contact_point2(self, debug=False):
849 TO BE DEVELOPED IN THE FUTURE
850 Finds the contact point on the curve.
852 FIXME: should be moved, probably to generalvclamp.py
855 #raw_plot=self.current.curve.default_plots()[0]
856 raw_plot=self.plots[0]
857 '''xext=self.plots[0].vectors[0][0]
858 yext=self.plots[0].vectors[0][1]
859 xret2=self.plots[0].vectors[1][0]
860 yret=self.plots[0].vectors[1][1]
862 xext=raw_plot.vectors[0][0]
863 yext=raw_plot.vectors[0][1]
864 xret2=raw_plot.vectors[1][0]
865 yret=raw_plot.vectors[1][1]
867 first_point=[xext[0], yext[0]]
868 last_point=[xext[-1], yext[-1]]
870 #regr=scipy.polyfit(first_point, last_point,1)[0:2]
871 diffx=abs(first_point[0]-last_point[0])
872 diffy=abs(first_point[1]-last_point[1])
874 #using polyfit results in numerical errors. good old algebra.
876 b=first_point[1]-(a*first_point[0])
877 baseline=scipy.polyval((a,b), xext)
879 ysub=[item-basitem for item,basitem in zip(yext,baseline)]
881 contact=ysub.index(min(ysub))
883 return xext,ysub,contact
885 #now, exploit a ClickedPoint instance to calculate index...
887 dummy.absolute_coords=(x_intercept,y_intercept)
888 dummy.find_graph_coords(xret2,yret)
891 return dummy.index, regr, regr_contact
897 def x_do_contact(self,args):
899 DEBUG COMMAND to be activated in the future
901 xext,ysub,contact=self.find_contact_point2(debug=True)
903 contact_plot=self.plots[0]
904 contact_plot.add_set(xext,ysub)
905 contact_plot.add_set([xext[contact]],[self.plots[0].vectors[0][1][contact]])
906 #contact_plot.add_set([first_point[0]],[first_point[1]])
907 #contact_plot.add_set([last_point[0]],[last_point[1]])
908 contact_plot.styles=[None,None,None,'scatter']
909 self._send_plot([contact_plot])
913 index,regr,regr_contact=self.find_contact_point2(debug=True)
916 raw_plot=self.current.curve.default_plots()[0]
917 xret=raw_plot.vectors[0][0]
918 #nc_line=[(item*regr[0])+regr[1] for item in x_nc]
919 nc_line=scipy.polyval(regr,xret)
920 c_line=scipy.polyval(regr_contact,xret)
923 contact_plot=self.current.curve.default_plots()[0]
924 contact_plot.add_set(xret, nc_line)
925 contact_plot.add_set(xret, c_line)
926 contact_plot.styles=[None,None,None,None]
927 #contact_plot.styles.append(None)
928 contact_plot.destination=1
929 self._send_plot([contact_plot])