hooke/plugin/fit.py

   1 '''
   2 FIT
   3
   4 Force spectroscopy curves basic fitting plugin.
   5 Licensed under the GNU GPL version 2
   6
   7 Non-standard Dependencies:
   8 procplots.py (plot processing plugin)
   9 '''
  10 from libhooke import WX_GOOD, ClickedPoint
  11 import wxversion
  12 wxversion.select(WX_GOOD)
  13 #from wx import PostEvent
  14 #from wx.lib.newevent import NewEvent
  15 import scipy
  16 import scipy.odr
  17 import numpy as np
  18 import copy
  19 import Queue
  20
  21 global measure_wlc
  22 global EVT_MEASURE_WLC
  23
  24 #measure_wlc, EVT_MEASURE_WLC = NewEvent()
  25
  26 global events_from_fit
  27 events_from_fit=Queue.Queue() #GUI ---> CLI COMMUNICATION
  28
  29
  30 class fitCommands(object):
  31
  32     def _plug_init(self):
  33         self.wlccurrent=None
  34         self.wlccontact_point=None
  35         self.wlccontact_index=None
  36
  37     def wlc_fit(self,clicked_points,xvector,yvector, pl_value, T=293, return_errors=False):
  38         '''
  39         Worm-like chain model fitting.
  40         The function is the simple polynomial worm-like chain as proposed by C.Bustamante, J.F.Marko, E.D.Siggia
  41         and S.Smith (Science. 1994 Sep 9;265(5178):1599-600.)
  42         '''
  43
  44         '''
  45         clicked_points[0] = contact point (calculated or hand-clicked)
  46         clicked_points[1] and [2] are edges of chunk
  47         '''
  48
  49         #STEP 1: Prepare the vectors to apply the fit.
  50
  51
  52         if pl_value is not None:
  53             pl_value=pl_value/(10**9)
  54
  55         #indexes of the selected chunk
  56         first_index=min(clicked_points[1].index, clicked_points[2].index)
  57         last_index=max(clicked_points[1].index, clicked_points[2].index)
  58
  59         #getting the chunk and reverting it
  60         xchunk,ychunk=xvector[first_index:last_index],yvector[first_index:last_index]
  61         xchunk.reverse()
  62         ychunk.reverse()
  63         #put contact point at zero and flip around the contact point (the fit wants a positive growth for extension and force)
  64         xchunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xchunk]
  65         ychunk_corr_up=[-(y-clicked_points[0].graph_coords[1]) for y in ychunk]
  66
  67         #make them arrays
  68         xchunk_corr_up=scipy.array(xchunk_corr_up)
  69         ychunk_corr_up=scipy.array(ychunk_corr_up)
  70
  71
  72         #STEP 2: actually do the fit
  73
  74         #Find furthest point of chunk and add it a bit; the fit must converge
  75         #from an excess!
  76         xchunk_high=max(xchunk_corr_up)
  77         xchunk_high+=(xchunk_high/10)
  78
  79         #Here are the linearized start parameters for the WLC.
  80         #[lambd=1/Lo , pii=1/P]
  81
  82         p0=[(1/xchunk_high),(1/(3.5e-10))]
  83         p0_plfix=[(1/xchunk_high)]
  84         '''
  85         ODR STUFF
  86         fixme: remove these comments after testing
  87         '''
  88
  89
  90         def f_wlc(params,x,T=T):
  91             '''
  92             wlc function for ODR fitting
  93             '''
  94             lambd,pii=params
  95             Kb=(1.38065e-23)
  96             therm=Kb*T
  97             y=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
  98             return y
  99
 100         def f_wlc_plfix(params,x,pl_value=pl_value,T=T):
 101             '''
 102             wlc function for ODR fitting
 103             '''
 104             lambd=params
 105             pii=1/pl_value
 106             Kb=(1.38065e-23)
 107             therm=Kb*T
 108             y=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
 109             return y
 110
 111         #make the ODR fit
 112         realdata=scipy.odr.RealData(xchunk_corr_up,ychunk_corr_up)
 113         if pl_value:
 114             model=scipy.odr.Model(f_wlc_plfix)
 115             o = scipy.odr.ODR(realdata, model, p0_plfix)
 116         else:
 117             model=scipy.odr.Model(f_wlc)
 118             o = scipy.odr.ODR(realdata, model, p0)
 119
 120         o.set_job(fit_type=2)
 121         out=o.run()
 122         fit_out=[(1/i) for i in out.beta]
 123
 124         #Calculate fit errors from output standard deviations.
 125         #We must propagate the error because we fit the *inverse* parameters!
 126         #The error = (error of the inverse)*(value**2)
 127         fit_errors=[]
 128         for sd,value in zip(out.sd_beta, fit_out):
 129             err_real=sd*(value**2)
 130             fit_errors.append(err_real)
 131
 132         def wlc_eval(x,params,pl_value,T):
 133             '''
 134             Evaluates the WLC function
 135             '''
 136             if not pl_value:
 137                 lambd, pii = params
 138             else:
 139                 lambd = params
 140
 141             if pl_value:
 142                 pii=1/pl_value
 143
 144             Kb=(1.38065e-23) #boltzmann constant
 145             therm=Kb*T #so we have thermal energy
 146
 147             return ( (therm*pii/4.0) * (((1-(x*lambd))**-2.0) - 1 + (4.0*x*lambd)) )
 148
 149         #STEP 3: plotting the fit
 150
 151         #obtain domain to plot the fit - from contact point to last_index plus 20 points
 152         thule_index=last_index+10
 153         if thule_index > len(xvector): #for rare cases in which we fit something at the END of whole curve.
 154             thule_index = len(xvector)
 155         #reverse etc. the domain
 156         xfit_chunk=xvector[clicked_points[0].index:thule_index]
 157         xfit_chunk.reverse()
 158         xfit_chunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xfit_chunk]
 159         xfit_chunk_corr_up=scipy.array(xfit_chunk_corr_up)
 160
 161         #the fitted curve: reflip, re-uncorrect
 162         yfit=wlc_eval(xfit_chunk_corr_up, out.beta, pl_value,T)
 163         yfit_down=[-y for y in yfit]
 164         yfit_corr_down=[y+clicked_points[0].graph_coords[1] for y in yfit_down]
 165
 166         if return_errors:
 167             return fit_out, yfit_corr_down, xfit_chunk, fit_errors
 168         else:
 169             return fit_out, yfit_corr_down, xfit_chunk, None
 170
 171     def fjc_fit(self,clicked_points,xvector,yvector, pl_value, T=293, return_errors=False):
 172         '''
 173         Freely-jointed chain function
 174         ref: C.Ray and B.B. Akhremitchev; http://www.chem.duke.edu/~boris/research/force_spectroscopy/fit_efjc.pdf
 175         '''
 176         '''
 177         clicked_points[0] is the contact point (calculated or hand-clicked)
 178         clicked_points[1] and [2] are edges of chunk
 179         '''
 180
 181         #STEP 1: Prepare the vectors to apply the fit.
 182         if pl_value is not None:
 183             pl_value=pl_value/(10**9)
 184
 185         #indexes of the selected chunk
 186         first_index=min(clicked_points[1].index, clicked_points[2].index)
 187         last_index=max(clicked_points[1].index, clicked_points[2].index)
 188
 189         #getting the chunk and reverting it
 190         xchunk,ychunk=xvector[first_index:last_index],yvector[first_index:last_index]
 191         xchunk.reverse()
 192         ychunk.reverse()
 193         #put contact point at zero and flip around the contact point (the fit wants a positive growth for extension and force)
 194         xchunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xchunk]
 195         ychunk_corr_up=[-(y-clicked_points[0].graph_coords[1]) for y in ychunk]
 196
 197
 198         #make them arrays
 199         xchunk_corr_up=scipy.array(xchunk_corr_up)
 200         ychunk_corr_up=scipy.array(ychunk_corr_up)
 201
 202
 203         #STEP 2: actually do the fit
 204
 205         #Find furthest point of chunk and add it a bit; the fit must converge
 206         #from an excess!
 207         xchunk_high=max(xchunk_corr_up)
 208         xchunk_high+=(xchunk_high/10)
 209
 210         #Here are the linearized start parameters for the WLC.
 211         #[lambd=1/Lo , pii=1/P]
 212
 213         p0=[(1/xchunk_high),(1/(3.5e-10))]
 214         p0_plfix=[(1/xchunk_high)]
 215         '''
 216         ODR STUFF
 217         fixme: remove these comments after testing
 218         '''
 219         def coth(z):
 220             '''
 221             hyperbolic cotangent
 222             '''
 223             return (np.exp(2*z)+1)/(np.exp(2*z)-1)
 224
 225         def x_fjc(params,f,T=T):
 226             '''
 227             fjc function for ODR fitting
 228             '''
 229             lambd,pii=params
 230             Kb=(1.38065e-23)
 231             therm=Kb*T
 232
 233             #x=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
 234             x=(1/lambd)*(coth(f*(1/pii)/therm) - (therm*pii)/f)
 235             return x
 236
 237         def x_fjc_plfix(params,f,pl_value=pl_value,T=T):
 238             '''
 239             fjc function for ODR fitting
 240             '''
 241             lambd=params
 242             pii=1/pl_value
 243             Kb=(1.38065e-23)
 244             therm=Kb*T
 245             #y=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
 246             x=(1/lambd)*(coth(f*(1/pii)/therm) - (therm*pii)/f)
 247             return x
 248
 249         #make the ODR fit
 250         realdata=scipy.odr.RealData(ychunk_corr_up,xchunk_corr_up)
 251         if pl_value:
 252             model=scipy.odr.Model(x_fjc_plfix)
 253             o = scipy.odr.ODR(realdata, model, p0_plfix)
 254         else:
 255             model=scipy.odr.Model(x_fjc)
 256             o = scipy.odr.ODR(realdata, model, p0)
 257
 258         o.set_job(fit_type=2)
 259         out=o.run()
 260         fit_out=[(1/i) for i in out.beta]
 261
 262         #Calculate fit errors from output standard deviations.
 263         #We must propagate the error because we fit the *inverse* parameters!
 264         #The error = (error of the inverse)*(value**2)
 265         fit_errors=[]
 266         for sd,value in zip(out.sd_beta, fit_out):
 267             err_real=sd*(value**2)
 268             fit_errors.append(err_real)
 269
 270         def fjc_eval(y,params,pl_value,T):
 271             '''
 272             Evaluates the WLC function
 273             '''
 274             if not pl_value:
 275                 lambd, pii = params
 276             else:
 277                 lambd = params
 278
 279             if pl_value:
 280                 pii=1/pl_value
 281
 282             Kb=(1.38065e-23) #boltzmann constant
 283             therm=Kb*T #so we have thermal energy
 284             #return ( (therm*pii/4.0) * (((1-(x*lambd))**-2.0) - 1 + (4.0*x*lambd)) )
 285             return (1/lambd)*(coth(y*(1/pii)/therm) - (therm*pii)/y)
 286
 287
 288         #STEP 3: plotting the fit
 289         #obtain domain to plot the fit - from contact point to last_index plus 20 points
 290         thule_index=last_index+10
 291         if thule_index > len(xvector): #for rare cases in which we fit something at the END of whole curve.
 292             thule_index = len(xvector)
 293         #reverse etc. the domain
 294         ychunk=yvector[clicked_points[0].index:thule_index]
 295
 296         if len(ychunk)>0:
 297             y_evalchunk=np.linspace(min(ychunk),max(ychunk),100)
 298         else:
 299             #Empty y-chunk. It happens whenever we set the contact point after a recognized peak,
 300             #or other buggy situations. Kludge to live with it now...
 301             ychunk=yvector[:thule_index]
 302             y_evalchunk=np.linspace(min(ychunk),max(ychunk),100)
 303
 304         yfit_down=[-y for y in y_evalchunk]
 305         yfit_corr_down=[y+clicked_points[0].graph_coords[1] for y in yfit_down]
 306         yfit_corr_down=scipy.array(yfit_corr_down)
 307
 308         #the fitted curve: reflip, re-uncorrect
 309         xfit=fjc_eval(yfit_corr_down, out.beta, pl_value,T)
 310         xfit=list(xfit)
 311         xfit.reverse()
 312         xfit_chunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xfit]
 313
 314         #xfit_chunk_corr_up=scipy.array(xfit_chunk_corr_up)
 315         #deltay=yfit_down[0]-yvector[clicked_points[0].index]
 316
 317         #This is a terrible, terrible kludge to find the point where it should normalize (and from where it should plot)
 318         xxxdists=[]
 319         for index in scipy.arange(1,len(xfit_chunk_corr_up),1):
 320             xxxdists.append((clicked_points[0].graph_coords[0]-xfit_chunk_corr_up[index])**2)
 321         normalize_index=xxxdists.index(min(xxxdists))
 322         #End of kludge
 323
 324         deltay=yfit_down[normalize_index]-clicked_points[0].graph_coords[1]
 325         yfit_corr_down=[y-deltay for y in yfit_down]
 326
 327         if return_errors:
 328             return fit_out, yfit_corr_down[normalize_index+1:], xfit_chunk_corr_up[normalize_index+1:], fit_errors
 329         else:
 330             return fit_out, yfit_corr_down[normalize_index+1:], xfit_chunk_corr_up[normalize_index+1:], None
 331
 332
 333     def do_wlc(self,args):
 334         '''
 335         WLC
 336         (fit.py plugin)
 337
 338         See the fit command
 339         '''
 340         self.do_fit(args)
 341
 342     def do_fjc(self,args):
 343         '''
 344         FJC
 345         (fit.py plugin)
 346
 347         See the fit command
 348         '''
 349         self.do_fit(args)
 350
 351     def do_fit(self,args):
 352         '''
 353         FIT
 354         (fit.py plugin)
 355         Fits an entropic elasticity function to a given chunk of the curve.
 356
 357         First you have to click a contact point.
 358         Then you have to click the two edges of the data you want to fit.
 359
 360         The fit function depends on the fit_function variable. You can set it with the command
 361         "set fit_function wlc" or  "set fit_function fjc" depending on the function you prefer.
 362
 363         For WLC, the function is the simple polynomial worm-like chain as proposed by
 364         C.Bustamante, J.F.Marko, E.D.Siggia and S.Smith (Science. 1994
 365         Sep 9;265(5178):1599-600.)
 366
 367         For FJC, ref:
 368         C.Ray and B.B. Akhremitchev; http://www.chem.duke.edu/~boris/research/force_spectroscopy/fit_efjc.pdf
 369
 370         Arguments:
 371         pl=[value] : Use a fixed persistent length (WLC) or Kuhn length (FJC) for the fit. If pl is not given,
 372                      the fit will be a 2-variable
 373                      fit. DO NOT put spaces between 'pl', '=' and the value.
 374                      The value must be in nanometers.
 375
 376         t=[value] : Use a user-defined temperature. The value must be in
 377                     kelvins; by default it is 293 K.
 378                     DO NOT put spaces between 't', '=' and the value.
 379
 380         noauto : allows for clicking the contact point by
 381                  hand (otherwise it is automatically estimated) the first time.
 382                  If subsequent measurements are made, the same contact point
 383                  clicked is used
 384
 385         reclick : redefines by hand the contact point, if noauto has been used before
 386                   but the user is unsatisfied of the previously choosen contact point.
 387         ---------
 388         Syntax: fit [pl=(value)] [t=value] [noauto]
 389         '''
 390         pl_value=None
 391         T=self.config['temperature']
 392         for arg in args.split():
 393             #look for a persistent length argument.
 394             if 'pl=' in arg:
 395                 pl_expression=arg.split('=')
 396                 pl_value=float(pl_expression[1]) #actual value
 397             #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
 398             if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
 399                 t_expression=arg.split('=')
 400                 T=float(t_expression[1])
 401
 402         #use the currently displayed plot for the fit
 403         displayed_plot=self._get_displayed_plot()
 404
 405         #handle contact point arguments correctly
 406         if 'reclick' in args.split():
 407             print 'Click contact point'
 408             contact_point=self._measure_N_points(N=1, whatset=1)[0]
 409             contact_point_index=contact_point.index
 410             self.wlccontact_point=contact_point
 411             self.wlccontact_index=contact_point.index
 412             self.wlccurrent=self.current.path
 413         elif 'noauto' in args.split():
 414             if self.wlccontact_index==None or self.wlccurrent != self.current.path:
 415                 print 'Click contact point'
 416                 contact_point=self._measure_N_points(N=1, whatset=1)[0]
 417                 contact_point_index=contact_point.index
 418                 self.wlccontact_point=contact_point
 419                 self.wlccontact_index=contact_point.index
 420                 self.wlccurrent=self.current.path
 421             else:
 422                 contact_point=self.wlccontact_point
 423                 contact_point_index=self.wlccontact_index
 424         else:
 425             cindex=self.find_contact_point()
 426             contact_point=ClickedPoint()
 427             contact_point.absolute_coords=displayed_plot.vectors[1][0][cindex], displayed_plot.vectors[1][1][cindex]
 428             contact_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
 429             contact_point.is_marker=True
 430
 431         print 'Click edges of chunk'
 432         points=self._measure_N_points(N=2, whatset=1)
 433         points=[contact_point]+points
 434         try:
 435             if self.config['fit_function']=='wlc':
 436                 params, yfit, xfit, fit_errors = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T, return_errors=True )
 437                 name_of_charlength='Persistent length'
 438             elif self.config['fit_function']=='fjc':
 439                 params, yfit, xfit, fit_errors = self.fjc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T, return_errors=True )
 440                 name_of_charlength='Kuhn length'
 441             else:
 442                 print 'No recognized fit function defined!'
 443                 print 'Set your fit function to wlc or fjc.'
 444                 return
 445
 446         except:
 447             print 'Fit not possible. Probably wrong interval -did you click two *different* points?'
 448             return
 449
 450         #FIXME: print "Kuhn length" for FJC
 451         print 'Fit function:',self.config['fit_function']
 452         print 'Contour length: ',params[0]*(1.0e+9),' nm'
 453         to_dump='contour '+self.current.path+' '+str(params[0]*(1.0e+9))+' nm'
 454         self.outlet.push(to_dump)
 455         if len(params)==2: #if we did choose 2-value fit
 456             print name_of_charlength+': ',params[1]*(1.0e+9),' nm'
 457             to_dump='persistent '+self.current.path+' '+str(params[1]*(1.0e+9))+' nm'
 458             self.outlet.push(to_dump)
 459
 460         if fit_errors:
 461             fit_nm=[i*(10**9) for i in fit_errors]
 462             print 'Standard deviation (contour length)', fit_nm[0]
 463             if len(fit_nm)>1:
 464                 print 'Standard deviation ('+name_of_charlength+')', fit_nm[1]
 465
 466
 467         #add the clicked points in the final PlotObject
 468         clickvector_x, clickvector_y=[], []
 469         for item in points:
 470             clickvector_x.append(item.graph_coords[0])
 471             clickvector_y.append(item.graph_coords[1])
 472
 473         #create a custom PlotObject to gracefully plot the fit along the curves
 474
 475         fitplot=copy.deepcopy(displayed_plot)
 476         fitplot.add_set(xfit,yfit)
 477         fitplot.add_set(clickvector_x,clickvector_y)
 478
 479         #FIXME: this colour/styles stuff must be solved at the root!
 480         if fitplot.styles==[]:
 481             fitplot.styles=[None,None,None,'scatter']
 482         else:
 483             fitplot.styles+=[None,'scatter']
 484
 485         if fitplot.colors==[]:
 486             fitplot.colors=[None,None,None,None]
 487         else:
 488             fitplot.colors+=[None,None]
 489
 490         self._send_plot([fitplot])
 491
 492
 493
 494     #----------
 495
 496
 497     def find_contact_point(self,plot=False):
 498         '''
 499         Finds the contact point on the curve.
 500
 501         The current algorithm (thanks to Francesco Musiani, francesco.musiani@unibo.it and Massimo Sandal) is:
 502         - take care of the PicoForce trigger bug - exclude retraction portions with too high standard deviation
 503         - fit the second half of the retraction curve to a line
 504         - if the fit is not almost horizontal, take a smaller chunk and repeat
 505         - otherwise, we have something horizontal
 506         - so take the average of horizontal points and use it as a baseline
 507
 508         Then, start from the rise of the retraction curve and look at the first point below the
 509         baseline.
 510
 511         FIXME: should be moved, probably to generalvclamp.py
 512         '''
 513
 514         if not plot:
 515             plot=self.plots[0]
 516
 517         outplot=self.subtract_curves(1)
 518         xret=outplot.vectors[1][0]
 519         ydiff=outplot.vectors[1][1]
 520
 521         xext=plot.vectors[0][0]
 522         yext=plot.vectors[0][1]
 523         xret2=plot.vectors[1][0]
 524         yret=plot.vectors[1][1]
 525
 526         #taking care of the picoforce trigger bug: we exclude portions of the curve that have too much
 527         #standard deviation. yes, a lot of magic is here.
 528         monster=True
 529         monlength=len(xret)-int(len(xret)/20)
 530         finalength=len(xret)
 531         while monster:
 532             monchunk=scipy.array(ydiff[monlength:finalength])
 533             if abs(np.std(monchunk)) < 2e-10:
 534                 monster=False
 535             else: #move away from the monster
 536                 monlength-=int(len(xret)/50)
 537                 finalength-=int(len(xret)/50)
 538
 539
 540         #take half of the thing
 541         endlength=int(len(xret)/2)
 542
 543         ok=False
 544
 545         while not ok:
 546             xchunk=yext[endlength:monlength]
 547             ychunk=yext[endlength:monlength]
 548             regr=scipy.stats.linregress(xchunk,ychunk)[0:2]
 549             #we stop if we found an almost-horizontal fit or if we're going too short...
 550             #FIXME: 0.1 and 6 here are "magic numbers" (although reasonable)
 551             if (abs(regr[1]) > 0.1) and ( endlength < len(xret)-int(len(xret)/6) ) :
 552                 endlength+=10
 553             else:
 554                 ok=True
 555
 556
 557         ymean=np.mean(ychunk) #baseline
 558
 559         index=0
 560         point = ymean+1
 561
 562         #find the first point below the calculated baseline
 563         while point > ymean:
 564             try:
 565                 point=yret[index]
 566                 index+=1
 567             except IndexError:
 568                 #The algorithm didn't find anything below the baseline! It should NEVER happen
 569                 index=0
 570                 return index
 571
 572         return index
 573
 574
 575
 576     def find_contact_point2(self, debug=False):
 577         '''
 578         TO BE DEVELOPED IN THE FUTURE
 579         Finds the contact point on the curve.
 580
 581         FIXME: should be moved, probably to generalvclamp.py
 582         '''
 583
 584         #raw_plot=self.current.curve.default_plots()[0]
 585         raw_plot=self.plots[0]
 586         '''xext=self.plots[0].vectors[0][0]
 587         yext=self.plots[0].vectors[0][1]
 588         xret2=self.plots[0].vectors[1][0]
 589         yret=self.plots[0].vectors[1][1]
 590         '''
 591         xext=raw_plot.vectors[0][0]
 592         yext=raw_plot.vectors[0][1]
 593         xret2=raw_plot.vectors[1][0]
 594         yret=raw_plot.vectors[1][1]
 595
 596         first_point=[xext[0], yext[0]]
 597         last_point=[xext[-1], yext[-1]]
 598
 599         #regr=scipy.polyfit(first_point, last_point,1)[0:2]
 600         diffx=abs(first_point[0]-last_point[0])
 601         diffy=abs(first_point[1]-last_point[1])
 602
 603         #using polyfit results in numerical errors. good old algebra.
 604         a=diffy/diffx
 605         b=first_point[1]-(a*first_point[0])
 606         baseline=scipy.polyval((a,b), xext)
 607
 608         ysub=[item-basitem for item,basitem in zip(yext,baseline)]
 609
 610         contact=ysub.index(min(ysub))
 611
 612         return xext,ysub,contact
 613
 614         #now, exploit a ClickedPoint instance to calculate index...
 615         dummy=ClickedPoint()
 616         dummy.absolute_coords=(x_intercept,y_intercept)
 617         dummy.find_graph_coords(xret2,yret)
 618
 619         if debug:
 620             return dummy.index, regr, regr_contact
 621         else:
 622             return dummy.index
 623
 624
 625
 626     def x_do_contact(self,args):
 627         '''
 628         DEBUG COMMAND to be activated in the future
 629         '''
 630         xext,ysub,contact=self.find_contact_point2(debug=True)
 631
 632         contact_plot=self.plots[0]
 633         contact_plot.add_set(xext,ysub)
 634         contact_plot.add_set([xext[contact]],[self.plots[0].vectors[0][1][contact]])
 635         #contact_plot.add_set([first_point[0]],[first_point[1]])
 636         #contact_plot.add_set([last_point[0]],[last_point[1]])
 637         contact_plot.styles=[None,None,None,'scatter']
 638         self._send_plot([contact_plot])
 639         return
 640
 641
 642         index,regr,regr_contact=self.find_contact_point2(debug=True)
 643         print regr
 644         print regr_contact
 645         raw_plot=self.current.curve.default_plots()[0]
 646         xret=raw_plot.vectors[0][0]
 647         #nc_line=[(item*regr[0])+regr[1] for item in x_nc]
 648         nc_line=scipy.polyval(regr,xret)
 649         c_line=scipy.polyval(regr_contact,xret)
 650
 651
 652         contact_plot=self.current.curve.default_plots()[0]
 653         contact_plot.add_set(xret, nc_line)
 654         contact_plot.add_set(xret, c_line)
 655         contact_plot.styles=[None,None,None,None]
 656         #contact_plot.styles.append(None)
 657         contact_plot.destination=1
 658         self._send_plot([contact_plot])