fit.py

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