2 # -*- coding: utf-8 -*-
7 Force spectroscopy curves basic fitting plugin.
8 Licensed under the GNU GPL version 2
10 Non-standard Dependencies:
11 procplots.py (plot processing plugin)
13 from ..libhooke import WX_GOOD, ClickedPoint
16 wxversion.select(WX_GOOD)
17 #from wx import PostEvent
18 #from wx.lib.newevent import NewEvent
27 global EVT_MEASURE_WLC
29 #measure_wlc, EVT_MEASURE_WLC = NewEvent()
31 global events_from_fit
32 events_from_fit=Queue.Queue() #GUI ---> CLI COMMUNICATION
35 class fitCommands(object):
39 self.wlccontact_point=None
40 self.wlccontact_index=None
43 '''Calculates the average distance from data to fit, scaled by the standard deviation
44 of the free cantilever area (thermal noise)
49 def wlc_fit(self,clicked_points,xvector,yvector, pl_value, T=293, return_errors=False):
51 Worm-like chain model fitting.
52 The function is the simple polynomial worm-like chain as proposed by C.Bustamante, J.F.Marko, E.D.Siggia
53 and S.Smith (Science. 1994 Sep 9;265(5178):1599-600.)
57 clicked_points[0] = contact point (calculated or hand-clicked)
58 clicked_points[1] and [2] are edges of chunk
61 #STEP 1: Prepare the vectors to apply the fit.
62 if pl_value is not None:
63 pl_value=pl_value/(10**9)
65 #indexes of the selected chunk
66 first_index=min(clicked_points[1].index, clicked_points[2].index)
67 last_index=max(clicked_points[1].index, clicked_points[2].index)
69 #getting the chunk and reverting it
70 xchunk,ychunk=xvector[first_index:last_index],yvector[first_index:last_index]
73 #put contact point at zero and flip around the contact point (the fit wants a positive growth for extension and force)
74 xchunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xchunk]
75 ychunk_corr_up=[-(y-clicked_points[0].graph_coords[1]) for y in ychunk]
78 xchunk_corr_up=scipy.array(xchunk_corr_up)
79 ychunk_corr_up=scipy.array(ychunk_corr_up)
82 #STEP 2: actually do the fit
84 #Find furthest point of chunk and add it a bit; the fit must converge
86 xchunk_high=max(xchunk_corr_up)
87 xchunk_high+=(xchunk_high/10)
89 #Here are the linearized start parameters for the WLC.
90 #[lambd=1/Lo , pii=1/P]
92 p0=[(1/xchunk_high),(1/(3.5e-10))]
93 p0_plfix=[(1/xchunk_high)]
96 fixme: remove these comments after testing
98 def dist(px,py,linex,liney):
99 distancesx=scipy.array([(px-x)**2 for x in linex])
100 minindex=np.argmin(distancesx)
101 print px, linex[0], linex[-1]
102 return (py-liney[minindex])**2
104 def f_wlc(params,x,T=T):
106 wlc function for ODR fitting
111 y=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
114 def f_wlc_plfix(params,x,pl_value=pl_value,T=T):
116 wlc function for ODR fitting
122 y=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
126 realdata=scipy.odr.RealData(xchunk_corr_up,ychunk_corr_up)
128 model=scipy.odr.Model(f_wlc_plfix)
129 o = scipy.odr.ODR(realdata, model, p0_plfix)
131 model=scipy.odr.Model(f_wlc)
132 o = scipy.odr.ODR(realdata, model, p0)
134 o.set_job(fit_type=2)
136 fit_out=[(1/i) for i in out.beta]
138 #Calculate fit errors from output standard deviations.
139 #We must propagate the error because we fit the *inverse* parameters!
140 #The error = (error of the inverse)*(value**2)
142 for sd,value in zip(out.sd_beta, fit_out):
143 err_real=sd*(value**2)
144 fit_errors.append(err_real)
146 def wlc_eval(x,params,pl_value,T):
148 Evaluates the WLC function
158 Kb=(1.38065e-23) #boltzmann constant
159 therm=Kb*T #so we have thermal energy
161 return ( (therm*pii/4.0) * (((1-(x*lambd))**-2.0) - 1 + (4.0*x*lambd)) )
164 #STEP 3: plotting the fit
166 #obtain domain to plot the fit - from contact point to last_index plus 20 points
167 thule_index=last_index+10
168 if thule_index > len(xvector): #for rare cases in which we fit something at the END of whole curve.
169 thule_index = len(xvector)
170 #reverse etc. the domain
171 xfit_chunk=xvector[clicked_points[0].index:thule_index]
173 xfit_chunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xfit_chunk]
174 xfit_chunk_corr_up=scipy.array(xfit_chunk_corr_up)
176 #the fitted curve: reflip, re-uncorrect
177 yfit=wlc_eval(xfit_chunk_corr_up, out.beta, pl_value,T)
178 yfit_down=[-y for y in yfit]
179 yfit_corr_down=[y+clicked_points[0].graph_coords[1] for y in yfit_down]
182 #calculate fit quality
184 yqeval=wlc_eval(xchunk_corr_up,out.beta,pl_value,T)
185 #we need to cut the extra from thuleindex
186 for qindex in np.arange(0,len(yqeval)):
187 qsum+=(yqeval[qindex]-ychunk_corr_up[qindex])**2
188 qstd=np.sqrt(qsum/len(ychunk_corr_up))
192 return fit_out, yfit_corr_down, xfit_chunk, fit_errors, qstd
194 return fit_out, yfit_corr_down, xfit_chunk, None, qstd
196 def fjc_fit(self,clicked_points,xvector,yvector, pl_value, T=293, return_errors=False):
198 Freely-jointed chain function
199 ref: C.Ray and B.B. Akhremitchev; http://www.chem.duke.edu/~boris/research/force_spectroscopy/fit_efjc.pdf
202 clicked_points[0] is the contact point (calculated or hand-clicked)
203 clicked_points[1] and [2] are edges of chunk
206 #STEP 1: Prepare the vectors to apply the fit.
207 if pl_value is not None:
208 pl_value=pl_value/(10**9)
210 #indexes of the selected chunk
211 first_index=min(clicked_points[1].index, clicked_points[2].index)
212 last_index=max(clicked_points[1].index, clicked_points[2].index)
214 #getting the chunk and reverting it
215 xchunk,ychunk=xvector[first_index:last_index],yvector[first_index:last_index]
218 #put contact point at zero and flip around the contact point (the fit wants a positive growth for extension and force)
219 xchunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xchunk]
220 ychunk_corr_up=[-(y-clicked_points[0].graph_coords[1]) for y in ychunk]
224 xchunk_corr_up=scipy.array(xchunk_corr_up)
225 ychunk_corr_up=scipy.array(ychunk_corr_up)
228 #STEP 2: actually do the fit
230 #Find furthest point of chunk and add it a bit; the fit must converge
232 xchunk_high=max(xchunk_corr_up)
233 xchunk_high+=(xchunk_high/10)
235 #Here are the linearized start parameters for the WLC.
236 #[lambd=1/Lo , pii=1/P]
238 p0=[(1/xchunk_high),(1/(3.5e-10))]
239 p0_plfix=[(1/xchunk_high)]
242 fixme: remove these comments after testing
244 def dist(px,py,linex,liney):
245 distancesx=scipy.array([(px-x)**2 for x in linex])
246 minindex=np.argmin(distancesx)
247 print minindex, px, linex[0], linex[-1]
248 return (py-liney[minindex])**2
254 return (np.exp(2*z)+1)/(np.exp(2*z)-1)
256 def x_fjc(params,f,T=T):
258 fjc function for ODR fitting
264 #x=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
265 x=(1/lambd)*(coth(f*(1/pii)/therm) - (therm*pii)/f)
268 def x_fjc_plfix(params,f,pl_value=pl_value,T=T):
270 fjc function for ODR fitting
276 #y=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
277 x=(1/lambd)*(coth(f*(1/pii)/therm) - (therm*pii)/f)
281 realdata=scipy.odr.RealData(ychunk_corr_up,xchunk_corr_up)
283 model=scipy.odr.Model(x_fjc_plfix)
284 o = scipy.odr.ODR(realdata, model, p0_plfix)
286 model=scipy.odr.Model(x_fjc)
287 o = scipy.odr.ODR(realdata, model, p0)
289 o.set_job(fit_type=2)
291 fit_out=[(1/i) for i in out.beta]
293 #Calculate fit errors from output standard deviations.
294 #We must propagate the error because we fit the *inverse* parameters!
295 #The error = (error of the inverse)*(value**2)
297 for sd,value in zip(out.sd_beta, fit_out):
298 err_real=sd*(value**2)
299 fit_errors.append(err_real)
301 def fjc_eval(y,params,pl_value,T):
303 Evaluates the WLC function
313 Kb=(1.38065e-23) #boltzmann constant
314 therm=Kb*T #so we have thermal energy
315 #return ( (therm*pii/4.0) * (((1-(x*lambd))**-2.0) - 1 + (4.0*x*lambd)) )
316 return (1/lambd)*(coth(y*(1/pii)/therm) - (therm*pii)/y)
320 #STEP 3: plotting the fit
321 #obtain domain to plot the fit - from contact point to last_index plus 20 points
322 thule_index=last_index+10
323 if thule_index > len(xvector): #for rare cases in which we fit something at the END of whole curve.
324 thule_index = len(xvector)
325 #reverse etc. the domain
326 ychunk=yvector[clicked_points[0].index:thule_index]
329 y_evalchunk=np.linspace(min(ychunk),max(ychunk),100)
331 #Empty y-chunk. It happens whenever we set the contact point after a recognized peak,
332 #or other buggy situations. Kludge to live with it now...
333 ychunk=yvector[:thule_index]
334 y_evalchunk=np.linspace(min(ychunk),max(ychunk),100)
336 yfit_down=[-y for y in y_evalchunk]
337 yfit_corr_down=[y+clicked_points[0].graph_coords[1] for y in yfit_down]
338 yfit_corr_down=scipy.array(yfit_corr_down)
340 #the fitted curve: reflip, re-uncorrect
341 xfit=fjc_eval(yfit_corr_down, out.beta, pl_value,T)
344 xfit_chunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xfit]
346 #xfit_chunk_corr_up=scipy.array(xfit_chunk_corr_up)
347 #deltay=yfit_down[0]-yvector[clicked_points[0].index]
349 #This is a terrible, terrible kludge to find the point where it should normalize (and from where it should plot)
351 for index in scipy.arange(1,len(xfit_chunk_corr_up),1):
352 xxxdists.append((clicked_points[0].graph_coords[0]-xfit_chunk_corr_up[index])**2)
353 normalize_index=xxxdists.index(min(xxxdists))
356 deltay=yfit_down[normalize_index]-clicked_points[0].graph_coords[1]
357 yfit_corr_down=[y-deltay for y in yfit_down]
360 #calculate fit quality
361 #creates dense y vector
362 yqeval=np.linspace(np.min(ychunk_corr_up)/2,np.max(ychunk_corr_up)*2,10*len(ychunk_corr_up))
363 #corresponding fitted x
364 xqeval=fjc_eval(yqeval,out.beta,pl_value,T)
367 for qindex in np.arange(0,len(ychunk_corr_up)):
368 qsum+=dist(xchunk_corr_up[qindex],ychunk_corr_up[qindex],xqeval,yqeval)
369 qstd=np.sqrt(qsum/len(ychunk_corr_up))
373 return fit_out, yfit_corr_down[normalize_index+1:], xfit_chunk_corr_up[normalize_index+1:], fit_errors, qstd
375 return fit_out, yfit_corr_down[normalize_index+1:], xfit_chunk_corr_up[normalize_index+1:], None, qstd
377 def efjc_fit(self,clicked_points,xvector,yvector, pl_value, T=293.0, return_errors=False):
379 Extended Freely-jointed chain function
380 ref: F Oesterhelt, M Rief and H E Gaub, New Journal of Physics 1 (1999) 6.1–6.11
381 Please note that 2-parameter fitting (length and kl) usually does not converge, use fixed kl
384 clicked_points[0] is the contact point (calculated or hand-clicked)
385 clicked_points[1] and [2] are edges of chunk
388 #Fixed parameters from reference
389 Kb=(1.38065e-2) #in pN.nm
391 Lh=0.280 #helical, nm
395 #STEP 1: Prepare the vectors to apply the fit.
397 #indexes of the selected chunk
398 first_index=min(clicked_points[1].index, clicked_points[2].index)
399 last_index=max(clicked_points[1].index, clicked_points[2].index)
401 #getting the chunk and reverting it
402 xchunk,ychunk=xvector[first_index:last_index],yvector[first_index:last_index]
405 #put contact point at zero and flip around the contact point (the fit wants a positive growth for extension and force)
406 xchunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xchunk]
407 ychunk_corr_up=[-(y-clicked_points[0].graph_coords[1]) for y in ychunk]
411 xchunk_corr_up=scipy.array(xchunk_corr_up)
412 ychunk_corr_up=scipy.array(ychunk_corr_up)
414 xchunk_corr_up_nm=xchunk_corr_up*1e9
415 ychunk_corr_up_pn=ychunk_corr_up*1e12
418 #STEP 2: actually do the fit
420 #Find furthest point of chunk and add it a bit; the fit must converge
422 xchunk_high=max(xchunk_corr_up_nm)
423 xchunk_high+=(xchunk_high/10.0)
425 #Here are the linearized start parameters for the WLC.
427 #max number of monomers (all helical)for a contour length xchunk_high
428 excessNs=xchunk_high/(Lp)
429 p0=[excessNs,(1.0/(0.7))]
430 p0_plfix=[(excessNs)]
432 def dist(px,py,linex,liney):
433 distancesx=scipy.array([(px-x)**2 for x in linex])
434 minindex=np.argmin(distancesx)
435 return (py-liney[minindex])**2
438 dG0=12.4242 #3kt in pN.nm
439 dL=0.078 #planar-helical
444 Lh=0.280 #helical, nm
449 return Lp/(np.exp(dG/therm)+1)+Lh/(np.exp(-dG/therm)+1)
455 return 1.0/np.tanh(z)
457 def x_efjc(params,f,T=T,Ks=Ks):
459 efjc function for ODR fitting
468 x=Ns*Lfactor(f)*(coth(beta)-1.0/beta)+Ns*f/Ks
471 def x_efjc_plfix(params,f,kl_value=pl_value,T=T,Ks=Ks):
473 efjc function for ODR fitting
482 x=Ns*Lfactor(f)*(coth(beta)-1.0/beta)+Ns*f/Ks
486 realdata=scipy.odr.RealData(ychunk_corr_up_pn,xchunk_corr_up_nm)
488 model=scipy.odr.Model(x_efjc_plfix)
489 o = scipy.odr.ODR(realdata, model, p0_plfix)
491 print 'WARNING eFJC fit with free pl sometimes does not converge'
492 model=scipy.odr.Model(x_efjc)
493 o = scipy.odr.ODR(realdata, model, p0)
495 o.set_job(fit_type=2)
502 kfit=1e-9/out.beta[1]
510 #Calculate fit errors from output standard deviations.
512 fit_errors.append(out.sd_beta[0]*Lp*1e-9)
514 fit_errors.append(1e9*out.sd_beta[1]*kfit**2)
518 def efjc_eval(y,params,pl_value,T=T,Lfactor=Lfactor,Ks=Ks):
520 Evaluates the eFJC function
530 Kb=(1.38065e-2) #boltzmann constant
531 therm=Kb*T #so we have thermal energy
534 x= Ns*Lfactor(y)*(coth(beta)-1.0/beta)+Ns*y/Ks
538 #STEP 3: plotting the fit
539 #obtain domain to plot the fit - from contact point to last_index plus 20 points
540 thule_index=last_index+10
541 if thule_index > len(xvector): #for rare cases in which we fit something at the END of whole curve.
542 thule_index = len(xvector)
543 #reverse etc. the domain
544 ychunk=yvector[clicked_points[0].index:thule_index]
547 y_evalchunk=np.linspace(min(ychunk),max(ychunk),100)
549 #Empty y-chunk. It happens whenever we set the contact point after a recognized peak,
550 #or other buggy situations. Kludge to live with it now...
551 ychunk=yvector[:thule_index]
552 y_evalchunk=np.linspace(min(ychunk),max(ychunk),100)
554 yfit_down=[-y for y in y_evalchunk]
555 yfit_corr_down=[y+clicked_points[0].graph_coords[1] for y in yfit_down]
556 yfit_corr_down=scipy.array(yfit_corr_down)
558 #the fitted curve: reflip, re-uncorrect
559 xfit=efjc_eval(1e12*yfit_corr_down, out.beta, pl_value,T)*1e-9
562 xfit_chunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xfit]
564 #xfit_chunk_corr_up=scipy.array(xfit_chunk_corr_up)
565 #deltay=yfit_down[0]-yvector[clicked_points[0].index]
567 #This is a terrible, terrible kludge to find the point where it should normalize (and from where it should plot)
569 for index in scipy.arange(1,len(xfit_chunk_corr_up),1):
570 xxxdists.append((clicked_points[0].graph_coords[0]-xfit_chunk_corr_up[index])**2)
571 normalize_index=xxxdists.index(min(xxxdists))
574 deltay=yfit_down[normalize_index]-clicked_points[0].graph_coords[1]
575 yfit_corr_down=[y-deltay for y in yfit_down]
577 #calculate fit quality
578 #creates dense y vector
579 yqeval=np.linspace(np.min(ychunk_corr_up_pn)/2,np.max(ychunk_corr_up_pn)*2,10*len(ychunk_corr_up_pn))
580 #corresponding fitted x
581 xqeval=efjc_eval(yqeval,out.beta,pl_value)
584 for qindex in np.arange(0,len(ychunk_corr_up_pn)):
585 qsum+=dist(xchunk_corr_up_nm[qindex],ychunk_corr_up_pn[qindex],xqeval,yqeval)
586 qstd=1e-12*np.sqrt(qsum/len(ychunk_corr_up_pn))
589 return fit_out, yfit_corr_down[normalize_index+1:], xfit_chunk_corr_up[normalize_index+1:], fit_errors, qstd
591 return fit_out, yfit_corr_down[normalize_index+1:], xfit_chunk_corr_up[normalize_index+1:], None, qstd
596 def do_wlc(self,args):
605 def do_fjc(self,args):
614 def do_fit(self,args):
618 Fits an entropic elasticity function to a given chunk of the curve.
620 First you have to click a contact point.
621 Then you have to click the two edges of the data you want to fit.
623 Fit quality compares the distance to the fit with the thermal noise (a good fit should be close to 1)
625 The fit function depends on the fit_function variable. You can set it with the command
626 "set fit_function wlc" or "set fit_function fjc" depending on the function you prefer.
628 For WLC, the function is the simple polynomial worm-like chain as proposed by
629 C.Bustamante, J.F.Marko, E.D.Siggia and S.Smith (Science. 1994
630 Sep 9;265(5178):1599-600.)
633 C.Ray and B.B. Akhremitchev; http://www.chem.duke.edu/~boris/research/force_spectroscopy/fit_efjc.pdf
636 F Oesterhelt, M Rief and H E Gaub, New Journal of Physics 1 (1999) 6.1–6.11 (section 4.2)
637 NOTE: use fixed pl for better results.
640 pl=[value] : Use a fixed persistent length (WLC) or Kuhn length (FJC) for the fit. If pl is not given,
641 the fit will be a 2-variable
642 fit. DO NOT put spaces between 'pl', '=' and the value.
643 The value must be in nanometers.
645 t=[value] : Use a user-defined temperature. The value must be in
646 kelvins; by default it is 293 K.
647 DO NOT put spaces between 't', '=' and the value.
649 noauto : allows for clicking the contact point by
650 hand (otherwise it is automatically estimated) the first time.
651 If subsequent measurements are made, the same contact point
654 reclick : redefines by hand the contact point, if noauto has been used before
655 but the user is unsatisfied of the previously choosen contact point.
657 Syntax: fit [pl=(value)] [t=value] [noauto]
660 T=self.config['temperature']
661 for arg in args.split():
662 #look for a persistent length argument.
664 pl_expression=arg.split('=')
665 pl_value=float(pl_expression[1]) #actual value
666 #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
667 if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
668 t_expression=arg.split('=')
669 T=float(t_expression[1])
671 #use the currently displayed plot for the fit
672 displayed_plot=self._get_displayed_plot()
674 #handle contact point arguments correctly
675 if 'reclick' in args.split():
676 print 'Click contact point'
677 contact_point=self._measure_N_points(N=1, whatset=1)[0]
678 contact_point_index=contact_point.index
679 self.wlccontact_point=contact_point
680 self.wlccontact_index=contact_point.index
681 self.wlccurrent=self.current.path
682 elif 'noauto' in args.split():
683 if self.wlccontact_index==None or self.wlccurrent != self.current.path:
684 print 'Click contact point'
685 contact_point=self._measure_N_points(N=1, whatset=1)[0]
686 contact_point_index=contact_point.index
687 self.wlccontact_point=contact_point
688 self.wlccontact_index=contact_point.index
689 self.wlccurrent=self.current.path
691 contact_point=self.wlccontact_point
692 contact_point_index=self.wlccontact_index
694 cindex=self.find_contact_point()
695 contact_point=ClickedPoint()
696 contact_point.absolute_coords=displayed_plot.vectors[1][0][cindex], displayed_plot.vectors[1][1][cindex]
697 contact_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
698 contact_point.is_marker=True
700 print 'Click edges of chunk'
701 points=self._measure_N_points(N=2, whatset=1)
702 points=[contact_point]+points
705 if self.config['fit_function']=='wlc':
706 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 )
707 name_of_charlength='Persistent length'
708 elif self.config['fit_function']=='fjc':
709 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 )
710 name_of_charlength='Kuhn length'
711 elif self.config['fit_function']=='efjc':
712 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 )
713 name_of_charlength='Kuhn length (e)'
715 print 'No recognized fit function defined!'
716 print 'Set your fit function to wlc, fjc or efjc.'
720 print 'Fit not possible. Probably wrong interval -did you click two *different* points?'
723 #FIXME: print "Kuhn length" for FJC
724 print 'Fit function:',self.config['fit_function']
725 print 'Contour length: ',params[0]*(1.0e+9),' nm'
726 to_dump='contour '+self.current.path+' '+str(params[0]*(1.0e+9))+' nm'
727 self.outlet.push(to_dump)
728 if len(params)==2: #if we did choose 2-value fit
729 print name_of_charlength+': ',params[1]*(1.0e+9),' nm'
730 to_dump='persistent '+self.current.path+' '+str(params[1]*(1.0e+9))+' nm'
731 self.outlet.push(to_dump)
734 fit_nm=[i*(10**9) for i in fit_errors]
735 print 'Standard deviation (contour length)', fit_nm[0]
737 print 'Standard deviation ('+name_of_charlength+')', fit_nm[1]
739 print 'Fit quality: '+str(qstd/np.std(displayed_plot.vectors[1][1][-20:-1]))
742 #add the clicked points in the final PlotObject
743 clickvector_x, clickvector_y=[], []
745 clickvector_x.append(item.graph_coords[0])
746 clickvector_y.append(item.graph_coords[1])
748 #create a custom PlotObject to gracefully plot the fit along the curves
750 fitplot=copy.deepcopy(displayed_plot)
751 fitplot.add_set(xfit,yfit)
752 fitplot.add_set(clickvector_x,clickvector_y)
754 #FIXME: this colour/styles stuff must be solved at the root!
755 if fitplot.styles==[]:
756 fitplot.styles=[None,None,None,'scatter']
758 fitplot.styles+=[None,'scatter']
760 if fitplot.colors==[]:
761 fitplot.colors=[None,None,None,None]
763 fitplot.colors+=[None,None]
765 self._send_plot([fitplot])
772 def find_contact_point(self,plot=False):
774 Finds the contact point on the curve.
776 The current algorithm (thanks to Francesco Musiani, francesco.musiani@unibo.it and Massimo Sandal) is:
777 - take care of the PicoForce trigger bug - exclude retraction portions with too high standard deviation
778 - fit the second half of the retraction curve to a line
779 - if the fit is not almost horizontal, take a smaller chunk and repeat
780 - otherwise, we have something horizontal
781 - so take the average of horizontal points and use it as a baseline
783 Then, start from the rise of the retraction curve and look at the first point below the
786 FIXME: should be moved, probably to generalvclamp.py
792 outplot=self.subtract_curves(1)
793 xret=outplot.vectors[1][0]
794 ydiff=outplot.vectors[1][1]
796 xext=plot.vectors[0][0]
797 yext=plot.vectors[0][1]
798 xret2=plot.vectors[1][0]
799 yret=plot.vectors[1][1]
801 #taking care of the picoforce trigger bug: we exclude portions of the curve that have too much
802 #standard deviation. yes, a lot of magic is here.
804 monlength=len(xret)-int(len(xret)/20)
807 monchunk=scipy.array(ydiff[monlength:finalength])
808 if abs(np.std(monchunk)) < 2e-10:
810 else: #move away from the monster
811 monlength-=int(len(xret)/50)
812 finalength-=int(len(xret)/50)
815 #take half of the thing
816 endlength=int(len(xret)/2)
821 xchunk=yext[endlength:monlength]
822 ychunk=yext[endlength:monlength]
823 regr=scipy.stats.linregress(xchunk,ychunk)[0:2]
824 #we stop if we found an almost-horizontal fit or if we're going too short...
825 #FIXME: 0.1 and 6 here are "magic numbers" (although reasonable)
826 if (abs(regr[1]) > 0.1) and ( endlength < len(xret)-int(len(xret)/6) ) :
832 ymean=np.mean(ychunk) #baseline
837 #find the first point below the calculated baseline
843 #The algorithm didn't find anything below the baseline! It should NEVER happen
851 def find_contact_point2(self, debug=False):
853 TO BE DEVELOPED IN THE FUTURE
854 Finds the contact point on the curve.
856 FIXME: should be moved, probably to generalvclamp.py
859 #raw_plot=self.current.curve.default_plots()[0]
860 raw_plot=self.plots[0]
861 '''xext=self.plots[0].vectors[0][0]
862 yext=self.plots[0].vectors[0][1]
863 xret2=self.plots[0].vectors[1][0]
864 yret=self.plots[0].vectors[1][1]
866 xext=raw_plot.vectors[0][0]
867 yext=raw_plot.vectors[0][1]
868 xret2=raw_plot.vectors[1][0]
869 yret=raw_plot.vectors[1][1]
871 first_point=[xext[0], yext[0]]
872 last_point=[xext[-1], yext[-1]]
874 #regr=scipy.polyfit(first_point, last_point,1)[0:2]
875 diffx=abs(first_point[0]-last_point[0])
876 diffy=abs(first_point[1]-last_point[1])
878 #using polyfit results in numerical errors. good old algebra.
880 b=first_point[1]-(a*first_point[0])
881 baseline=scipy.polyval((a,b), xext)
883 ysub=[item-basitem for item,basitem in zip(yext,baseline)]
885 contact=ysub.index(min(ysub))
887 return xext,ysub,contact
889 #now, exploit a ClickedPoint instance to calculate index...
891 dummy.absolute_coords=(x_intercept,y_intercept)
892 dummy.find_graph_coords(xret2,yret)
895 return dummy.index, regr, regr_contact
901 def x_do_contact(self,args):
903 DEBUG COMMAND to be activated in the future
905 xext,ysub,contact=self.find_contact_point2(debug=True)
907 contact_plot=self.plots[0]
908 contact_plot.add_set(xext,ysub)
909 contact_plot.add_set([xext[contact]],[self.plots[0].vectors[0][1][contact]])
910 #contact_plot.add_set([first_point[0]],[first_point[1]])
911 #contact_plot.add_set([last_point[0]],[last_point[1]])
912 contact_plot.styles=[None,None,None,'scatter']
913 self._send_plot([contact_plot])
917 index,regr,regr_contact=self.find_contact_point2(debug=True)
920 raw_plot=self.current.curve.default_plots()[0]
921 xret=raw_plot.vectors[0][0]
922 #nc_line=[(item*regr[0])+regr[1] for item in x_nc]
923 nc_line=scipy.polyval(regr,xret)
924 c_line=scipy.polyval(regr_contact,xret)
927 contact_plot=self.current.curve.default_plots()[0]
928 contact_plot.add_set(xret, nc_line)
929 contact_plot.add_set(xret, c_line)
930 contact_plot.styles=[None,None,None,None]
931 #contact_plot.styles.append(None)
932 contact_plot.destination=1
933 self._send_plot([contact_plot])