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