fit.py

   1 #!/usr/bin/env python
   2
   3 '''
   4 FIT
   5
   6 Force spectroscopy curves basic fitting plugin.
   7 Licensed under the GNU GPL version 2
   8
   9 Non-standard Dependencies:
  10 procplots.py (plot processing plugin)
  11 '''
  12 from libhooke import WX_GOOD, ClickedPoint
  13 import wxversion
  14 wxversion.select(WX_GOOD)
  15 #from wx import PostEvent
  16 #from wx.lib.newevent import NewEvent
  17 import scipy
  18 import scipy.odr
  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:
  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         #print clicked_points[0].graph_coords[1]
 298         y_evalchunk=np.linspace(min(ychunk),max(ychunk),100)
 299
 300         yfit_down=[-y for y in y_evalchunk]
 301         yfit_corr_down=[y+clicked_points[0].graph_coords[1] for y in yfit_down]
 302         yfit_corr_down=scipy.array(yfit_corr_down)
 303
 304         #the fitted curve: reflip, re-uncorrect
 305         xfit=fjc_eval(yfit_corr_down, out.beta, pl_value,T)
 306         xfit=list(xfit)
 307         xfit.reverse()
 308         xfit_chunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xfit]
 309
 310         #xfit_chunk_corr_up=scipy.array(xfit_chunk_corr_up)
 311         #deltay=yfit_down[0]-yvector[clicked_points[0].index]
 312
 313         #This is a terrible, terrible kludge to find the point where it should normalize (and from where it should plot)
 314         xxxdists=[]
 315         for index in scipy.arange(1,len(xfit_chunk_corr_up),1):
 316             xxxdists.append((clicked_points[0].graph_coords[0]-xfit_chunk_corr_up[index])**2)
 317         normalize_index=xxxdists.index(min(xxxdists))
 318         #End of kludge
 319
 320         deltay=yfit_down[normalize_index]-clicked_points[0].graph_coords[1]
 321         yfit_corr_down=[y-deltay for y in yfit_down]
 322
 323         if return_errors:
 324             return fit_out, yfit_corr_down[normalize_index+1:], xfit_chunk_corr_up[normalize_index+1:], fit_errors
 325         else:
 326             return fit_out, yfit_corr_down[normalize_index+1:], xfit_chunk_corr_up[normalize_index+1:], None
 327
 328
 329     def do_wlc(self,args):
 330         '''
 331         WLC
 332         (fit.py plugin)
 333
 334         See the fit command
 335         '''
 336         self.do_fit(args)
 337
 338     def do_fjc(self,args):
 339         '''
 340         FJC
 341         (fit.py plugin)
 342
 343         See the fit command
 344         '''
 345         self.do_fit(args)
 346
 347     def do_fit(self,args):
 348         '''
 349         FIT
 350         (fit.py plugin)
 351         Fits an entropic elasticity function to a given chunk of the curve.
 352
 353         First you have to click a contact point.
 354         Then you have to click the two edges of the data you want to fit.
 355
 356         The fit function depends on the fit_function variable. You can set it with the command
 357         "set fit_function wlc" or  "set fit_function fjc" depending on the function you prefer.
 358
 359         For WLC, the function is the simple polynomial worm-like chain as proposed by
 360         C.Bustamante, J.F.Marko, E.D.Siggia and S.Smith (Science. 1994
 361         Sep 9;265(5178):1599-600.)
 362
 363         For FJC, ref:
 364         C.Ray and B.B. Akhremitchev; http://www.chem.duke.edu/~boris/research/force_spectroscopy/fit_efjc.pdf
 365
 366         Arguments:
 367         pl=[value] : Use a fixed persistent length (WLC) or Kuhn length (FJC) for the fit. If pl is not given,
 368                      the fit will be a 2-variable
 369                      fit. DO NOT put spaces between 'pl', '=' and the value.
 370                      The value must be in nanometers.
 371
 372         t=[value] : Use a user-defined temperature. The value must be in
 373                     kelvins; by default it is 293 K.
 374                     DO NOT put spaces between 't', '=' and the value.
 375
 376         noauto : allows for clicking the contact point by
 377                  hand (otherwise it is automatically estimated) the first time.
 378                  If subsequent measurements are made, the same contact point
 379                  clicked is used
 380
 381         reclick : redefines by hand the contact point, if noauto has been used before
 382                   but the user is unsatisfied of the previously choosen contact point.
 383         ---------
 384         Syntax: fit [pl=(value)] [t=value] [noauto]
 385         '''
 386         pl_value=None
 387         T=self.config['temperature']
 388         for arg in args.split():
 389             #look for a persistent length argument.
 390             if 'pl=' in arg:
 391                 pl_expression=arg.split('=')
 392                 pl_value=float(pl_expression[1]) #actual value
 393             #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
 394             if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
 395                 t_expression=arg.split('=')
 396                 T=float(t_expression[1])
 397
 398         #use the currently displayed plot for the fit
 399         displayed_plot=self._get_displayed_plot()
 400
 401         #handle contact point arguments correctly
 402         if 'reclick' in args.split():
 403             print 'Click contact point'
 404             contact_point=self._measure_N_points(N=1, whatset=1)[0]
 405             contact_point_index=contact_point.index
 406             self.wlccontact_point=contact_point
 407             self.wlccontact_index=contact_point.index
 408             self.wlccurrent=self.current.path
 409         elif 'noauto' in args.split():
 410             if self.wlccontact_index==None or self.wlccurrent != self.current.path:
 411                 print 'Click contact point'
 412                 contact_point=self._measure_N_points(N=1, whatset=1)[0]
 413                 contact_point_index=contact_point.index
 414                 self.wlccontact_point=contact_point
 415                 self.wlccontact_index=contact_point.index
 416                 self.wlccurrent=self.current.path
 417             else:
 418                 contact_point=self.wlccontact_point
 419                 contact_point_index=self.wlccontact_index
 420         else:
 421             cindex=self.find_contact_point()
 422             contact_point=ClickedPoint()
 423             contact_point.absolute_coords=displayed_plot.vectors[1][0][cindex], displayed_plot.vectors[1][1][cindex]
 424             contact_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
 425             contact_point.is_marker=True
 426
 427         print 'Click edges of chunk'
 428         points=self._measure_N_points(N=2, whatset=1)
 429         points=[contact_point]+points
 430         try:
 431             if self.config['fit_function']=='wlc':
 432                 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 )
 433                 name_of_charlength='Persistent length'
 434             elif self.config['fit_function']=='fjc':
 435                 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 )
 436                 name_of_charlength='Kuhn length'
 437             else:
 438                 print 'No recognized fit function defined!'
 439                 print 'Set your fit function to wlc or fjc.'
 440                 return
 441
 442         except:
 443             print 'Fit not possible. Probably wrong interval -did you click two *different* points?'
 444             return
 445
 446         #FIXME: print "Kuhn length" for FJC
 447         print 'Fit function:',self.config['fit_function']
 448         print 'Contour length: ',params[0]*(1.0e+9),' nm'
 449         to_dump='contour '+self.current.path+' '+str(params[0]*(1.0e+9))+' nm'
 450         self.outlet.push(to_dump)
 451         if len(params)==2: #if we did choose 2-value fit
 452             print name_of_charlength+': ',params[1]*(1.0e+9),' nm'
 453             to_dump='persistent '+self.current.path+' '+str(params[1]*(1.0e+9))+' nm'
 454             self.outlet.push(to_dump)
 455
 456         if fit_errors:
 457             fit_nm=[i*(10**9) for i in fit_errors]
 458             print 'Standard deviation (contour length)', fit_nm[0]
 459             if len(fit_nm)>1:
 460                 print 'Standard deviation ('+name_of_charlength+')', fit_nm[1]
 461
 462
 463         #add the clicked points in the final PlotObject
 464         clickvector_x, clickvector_y=[], []
 465         for item in points:
 466             clickvector_x.append(item.graph_coords[0])
 467             clickvector_y.append(item.graph_coords[1])
 468
 469         #create a custom PlotObject to gracefully plot the fit along the curves
 470
 471         fitplot=copy.deepcopy(displayed_plot)
 472         fitplot.add_set(xfit,yfit)
 473         fitplot.add_set(clickvector_x,clickvector_y)
 474
 475         #FIXME: this colour/styles stuff must be solved at the root!
 476         if fitplot.styles==[]:
 477             fitplot.styles=[None,None,None,'scatter']
 478         else:
 479             fitplot.styles+=[None,'scatter']
 480
 481         if fitplot.colors==[]:
 482             fitplot.colors=[None,None,None,None]
 483         else:
 484             fitplot.colors+=[None,None]
 485
 486         self._send_plot([fitplot])
 487
 488
 489
 490     #----------
 491
 492
 493     def find_contact_point(self,plot=False):
 494         '''
 495         Finds the contact point on the curve.
 496
 497         The current algorithm (thanks to Francesco Musiani, francesco.musiani@unibo.it and Massimo Sandal) is:
 498         - take care of the PicoForce trigger bug - exclude retraction portions with too high standard deviation
 499         - fit the second half of the retraction curve to a line
 500         - if the fit is not almost horizontal, take a smaller chunk and repeat
 501         - otherwise, we have something horizontal
 502         - so take the average of horizontal points and use it as a baseline
 503
 504         Then, start from the rise of the retraction curve and look at the first point below the
 505         baseline.
 506
 507         FIXME: should be moved, probably to generalvclamp.py
 508         '''
 509
 510         if not plot:
 511             plot=self.plots[0]
 512
 513         outplot=self.subtract_curves(1)
 514         xret=outplot.vectors[1][0]
 515         ydiff=outplot.vectors[1][1]
 516
 517         xext=plot.vectors[0][0]
 518         yext=plot.vectors[0][1]
 519         xret2=plot.vectors[1][0]
 520         yret=plot.vectors[1][1]
 521
 522         #taking care of the picoforce trigger bug: we exclude portions of the curve that have too much
 523         #standard deviation. yes, a lot of magic is here.
 524         monster=True
 525         monlength=len(xret)-int(len(xret)/20)
 526         finalength=len(xret)
 527         while monster:
 528             monchunk=scipy.array(ydiff[monlength:finalength])
 529             if abs(np.std(monchunk)) < 2e-10:
 530                 monster=False
 531             else: #move away from the monster
 532                 monlength-=int(len(xret)/50)
 533                 finalength-=int(len(xret)/50)
 534
 535
 536         #take half of the thing
 537         endlength=int(len(xret)/2)
 538
 539         ok=False
 540
 541         while not ok:
 542             xchunk=yext[endlength:monlength]
 543             ychunk=yext[endlength:monlength]
 544             regr=scipy.stats.linregress(xchunk,ychunk)[0:2]
 545             #we stop if we found an almost-horizontal fit or if we're going too short...
 546             #FIXME: 0.1 and 6 here are "magic numbers" (although reasonable)
 547             if (abs(regr[1]) > 0.1) and ( endlength < len(xret)-int(len(xret)/6) ) :
 548                 endlength+=10
 549             else:
 550                 ok=True
 551
 552
 553         ymean=np.mean(ychunk) #baseline
 554
 555         index=0
 556         point = ymean+1
 557
 558         #find the first point below the calculated baseline
 559         while point > ymean:
 560             try:
 561                 point=yret[index]
 562                 index+=1
 563             except IndexError:
 564                 #The algorithm didn't find anything below the baseline! It should NEVER happen
 565                 index=0
 566                 return index
 567
 568         return index
 569
 570
 571
 572     def find_contact_point2(self, debug=False):
 573         '''
 574         TO BE DEVELOPED IN THE FUTURE
 575         Finds the contact point on the curve.
 576
 577         FIXME: should be moved, probably to generalvclamp.py
 578         '''
 579
 580         #raw_plot=self.current.curve.default_plots()[0]
 581         raw_plot=self.plots[0]
 582         '''xext=self.plots[0].vectors[0][0]
 583         yext=self.plots[0].vectors[0][1]
 584         xret2=self.plots[0].vectors[1][0]
 585         yret=self.plots[0].vectors[1][1]
 586         '''
 587         xext=raw_plot.vectors[0][0]
 588         yext=raw_plot.vectors[0][1]
 589         xret2=raw_plot.vectors[1][0]
 590         yret=raw_plot.vectors[1][1]
 591
 592         first_point=[xext[0], yext[0]]
 593         last_point=[xext[-1], yext[-1]]
 594
 595         #regr=scipy.polyfit(first_point, last_point,1)[0:2]
 596         diffx=abs(first_point[0]-last_point[0])
 597         diffy=abs(first_point[1]-last_point[1])
 598
 599         #using polyfit results in numerical errors. good old algebra.
 600         a=diffy/diffx
 601         b=first_point[1]-(a*first_point[0])
 602         baseline=scipy.polyval((a,b), xext)
 603
 604         ysub=[item-basitem for item,basitem in zip(yext,baseline)]
 605
 606         contact=ysub.index(min(ysub))
 607
 608         return xext,ysub,contact
 609
 610         #now, exploit a ClickedPoint instance to calculate index...
 611         dummy=ClickedPoint()
 612         dummy.absolute_coords=(x_intercept,y_intercept)
 613         dummy.find_graph_coords(xret2,yret)
 614
 615         if debug:
 616             return dummy.index, regr, regr_contact
 617         else:
 618             return dummy.index
 619
 620
 621
 622     def x_do_contact(self,args):
 623         '''
 624         DEBUG COMMAND to be activated in the future
 625         '''
 626         xext,ysub,contact=self.find_contact_point2(debug=True)
 627
 628         contact_plot=self.plots[0]
 629         contact_plot.add_set(xext,ysub)
 630         contact_plot.add_set([xext[contact]],[self.plots[0].vectors[0][1][contact]])
 631         #contact_plot.add_set([first_point[0]],[first_point[1]])
 632         #contact_plot.add_set([last_point[0]],[last_point[1]])
 633         contact_plot.styles=[None,None,None,'scatter']
 634         self._send_plot([contact_plot])
 635         return
 636
 637
 638         index,regr,regr_contact=self.find_contact_point2(debug=True)
 639         print regr
 640         print regr_contact
 641         raw_plot=self.current.curve.default_plots()[0]
 642         xret=raw_plot.vectors[0][0]
 643         #nc_line=[(item*regr[0])+regr[1] for item in x_nc]
 644         nc_line=scipy.polyval(regr,xret)
 645         c_line=scipy.polyval(regr_contact,xret)
 646
 647
 648         contact_plot=self.current.curve.default_plots()[0]
 649         contact_plot.add_set(xret, nc_line)
 650         contact_plot.add_set(xret, c_line)
 651         contact_plot.styles=[None,None,None,None]
 652         #contact_plot.styles.append(None)
 653         contact_plot.destination=1
 654         self._send_plot([contact_plot])
 655