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