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