hooke/plugin/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
  15 import wxversion
  16 wxversion.select(WX_GOOD)
  17 #from wx import PostEvent
  18 #from wx.lib.newevent import NewEvent
  19 import scipy
  20 import scipy.odr
  21 import scipy.stats
  22 import numpy as np
  23 import copy
  24 import Queue
  25
  26 global measure_wlc
  27 global EVT_MEASURE_WLC
  28
  29 #measure_wlc, EVT_MEASURE_WLC = NewEvent()
  30
  31 global events_from_fit
  32 events_from_fit=Queue.Queue() #GUI ---> CLI COMMUNICATION
  33
  34
  35 class fitCommands(object):
  36
  37     def _plug_init(self):
  38         self.wlccurrent=None
  39         self.wlccontact_point=None
  40         self.wlccontact_index=None
  41
  42     def dist2fit(self):
  43         '''Calculates the average distance from data to fit, scaled by the standard deviation
  44         of the free cantilever area (thermal noise)
  45         '''
  46
  47
  48
  49     def wlc_fit(self,clicked_points,xvector,yvector, pl_value, T=293, return_errors=False):
  50         '''
  51         Worm-like chain model fitting.
  52         The function is the simple polynomial worm-like chain as proposed by C.Bustamante, J.F.Marko, E.D.Siggia
  53         and S.Smith (Science. 1994 Sep 9;265(5178):1599-600.)
  54         '''
  55
  56         '''
  57         clicked_points[0] = contact point (calculated or hand-clicked)
  58         clicked_points[1] and [2] are edges of chunk
  59         '''
  60
  61         #STEP 1: Prepare the vectors to apply the fit.
  62         if pl_value is not None:
  63             pl_value=pl_value/(10**9)
  64
  65         #indexes of the selected chunk
  66         first_index=min(clicked_points[1].index, clicked_points[2].index)
  67         last_index=max(clicked_points[1].index, clicked_points[2].index)
  68
  69         #getting the chunk and reverting it
  70         xchunk,ychunk=xvector[first_index:last_index],yvector[first_index:last_index]
  71         xchunk.reverse()
  72         ychunk.reverse()
  73         #put contact point at zero and flip around the contact point (the fit wants a positive growth for extension and force)
  74         xchunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xchunk]
  75         ychunk_corr_up=[-(y-clicked_points[0].graph_coords[1]) for y in ychunk]
  76
  77         #make them arrays
  78         xchunk_corr_up=scipy.array(xchunk_corr_up)
  79         ychunk_corr_up=scipy.array(ychunk_corr_up)
  80
  81
  82         #STEP 2: actually do the fit
  83
  84         #Find furthest point of chunk and add it a bit; the fit must converge
  85         #from an excess!
  86         xchunk_high=max(xchunk_corr_up)
  87         xchunk_high+=(xchunk_high/10)
  88
  89         #Here are the linearized start parameters for the WLC.
  90         #[lambd=1/Lo , pii=1/P]
  91
  92         p0=[(1/xchunk_high),(1/(3.5e-10))]
  93         p0_plfix=[(1/xchunk_high)]
  94         '''
  95         ODR STUFF
  96         fixme: remove these comments after testing
  97         '''
  98         def dist(px,py,linex,liney):
  99             distancesx=scipy.array([(px-x)**2 for x in linex])
 100             minindex=np.argmin(distancesx)
 101             print px, linex[0], linex[-1]
 102             return (py-liney[minindex])**2
 103
 104         def f_wlc(params,x,T=T):
 105             '''
 106             wlc function for ODR fitting
 107             '''
 108             lambd,pii=params
 109             Kb=(1.38065e-23)
 110             therm=Kb*T
 111             y=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
 112             return y
 113
 114         def f_wlc_plfix(params,x,pl_value=pl_value,T=T):
 115             '''
 116             wlc function for ODR fitting
 117             '''
 118             lambd=params
 119             pii=1/pl_value
 120             Kb=(1.38065e-23)
 121             therm=Kb*T
 122             y=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
 123             return y
 124
 125         #make the ODR fit
 126         realdata=scipy.odr.RealData(xchunk_corr_up,ychunk_corr_up)
 127         if pl_value:
 128             model=scipy.odr.Model(f_wlc_plfix)
 129             o = scipy.odr.ODR(realdata, model, p0_plfix)
 130         else:
 131             model=scipy.odr.Model(f_wlc)
 132             o = scipy.odr.ODR(realdata, model, p0)
 133
 134         o.set_job(fit_type=2)
 135         out=o.run()
 136         fit_out=[(1/i) for i in out.beta]
 137
 138         #Calculate fit errors from output standard deviations.
 139         #We must propagate the error because we fit the *inverse* parameters!
 140         #The error = (error of the inverse)*(value**2)
 141         fit_errors=[]
 142         for sd,value in zip(out.sd_beta, fit_out):
 143             err_real=sd*(value**2)
 144             fit_errors.append(err_real)
 145
 146         def wlc_eval(x,params,pl_value,T):
 147             '''
 148             Evaluates the WLC function
 149             '''
 150             if not pl_value:
 151                 lambd, pii = params
 152             else:
 153                 lambd = params
 154
 155             if pl_value:
 156                 pii=1/pl_value
 157
 158             Kb=(1.38065e-23) #boltzmann constant
 159             therm=Kb*T #so we have thermal energy
 160
 161             return ( (therm*pii/4.0) * (((1-(x*lambd))**-2.0) - 1 + (4.0*x*lambd)) )
 162
 163
 164         #STEP 3: plotting the fit
 165
 166         #obtain domain to plot the fit - from contact point to last_index plus 20 points
 167         thule_index=last_index+10
 168         if thule_index > len(xvector): #for rare cases in which we fit something at the END of whole curve.
 169             thule_index = len(xvector)
 170         #reverse etc. the domain
 171         xfit_chunk=xvector[clicked_points[0].index:thule_index]
 172         xfit_chunk.reverse()
 173         xfit_chunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xfit_chunk]
 174         xfit_chunk_corr_up=scipy.array(xfit_chunk_corr_up)
 175
 176         #the fitted curve: reflip, re-uncorrect
 177         yfit=wlc_eval(xfit_chunk_corr_up, out.beta, pl_value,T)
 178         yfit_down=[-y for y in yfit]
 179         yfit_corr_down=[y+clicked_points[0].graph_coords[1] for y in yfit_down]
 180
 181
 182         #calculate fit quality
 183         qsum=0
 184         yqeval=wlc_eval(xchunk_corr_up,out.beta,pl_value,T)
 185         #we need to cut the extra from thuleindex
 186         for qindex in np.arange(0,len(yqeval)):
 187             qsum+=(yqeval[qindex]-ychunk_corr_up[qindex])**2
 188         qstd=np.sqrt(qsum/len(ychunk_corr_up))
 189
 190
 191         if return_errors:
 192             return fit_out, yfit_corr_down, xfit_chunk, fit_errors, qstd
 193         else:
 194             return fit_out, yfit_corr_down, xfit_chunk, None, qstd
 195
 196     def fjc_fit(self,clicked_points,xvector,yvector, pl_value, T=293, return_errors=False):
 197         '''
 198         Freely-jointed chain function
 199         ref: C.Ray and B.B. Akhremitchev; http://www.chem.duke.edu/~boris/research/force_spectroscopy/fit_efjc.pdf
 200         '''
 201         '''
 202         clicked_points[0] is the contact point (calculated or hand-clicked)
 203         clicked_points[1] and [2] are edges of chunk
 204         '''
 205
 206         #STEP 1: Prepare the vectors to apply the fit.
 207         if pl_value is not None:
 208             pl_value=pl_value/(10**9)
 209
 210         #indexes of the selected chunk
 211         first_index=min(clicked_points[1].index, clicked_points[2].index)
 212         last_index=max(clicked_points[1].index, clicked_points[2].index)
 213
 214         #getting the chunk and reverting it
 215         xchunk,ychunk=xvector[first_index:last_index],yvector[first_index:last_index]
 216         xchunk.reverse()
 217         ychunk.reverse()
 218         #put contact point at zero and flip around the contact point (the fit wants a positive growth for extension and force)
 219         xchunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xchunk]
 220         ychunk_corr_up=[-(y-clicked_points[0].graph_coords[1]) for y in ychunk]
 221
 222
 223         #make them arrays
 224         xchunk_corr_up=scipy.array(xchunk_corr_up)
 225         ychunk_corr_up=scipy.array(ychunk_corr_up)
 226
 227
 228         #STEP 2: actually do the fit
 229
 230         #Find furthest point of chunk and add it a bit; the fit must converge
 231         #from an excess!
 232         xchunk_high=max(xchunk_corr_up)
 233         xchunk_high+=(xchunk_high/10)
 234
 235         #Here are the linearized start parameters for the WLC.
 236         #[lambd=1/Lo , pii=1/P]
 237
 238         p0=[(1/xchunk_high),(1/(3.5e-10))]
 239         p0_plfix=[(1/xchunk_high)]
 240         '''
 241         ODR STUFF
 242         fixme: remove these comments after testing
 243         '''
 244         def dist(px,py,linex,liney):
 245             distancesx=scipy.array([(px-x)**2 for x in linex])
 246             minindex=np.argmin(distancesx)
 247             print minindex, px, linex[0], linex[-1]
 248             return (py-liney[minindex])**2
 249
 250         def coth(z):
 251             '''
 252             hyperbolic cotangent
 253             '''
 254             return (np.exp(2*z)+1)/(np.exp(2*z)-1)
 255
 256         def x_fjc(params,f,T=T):
 257             '''
 258             fjc function for ODR fitting
 259             '''
 260             lambd,pii=params
 261             Kb=(1.38065e-23)
 262             therm=Kb*T
 263
 264             #x=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
 265             x=(1/lambd)*(coth(f*(1/pii)/therm) - (therm*pii)/f)
 266             return x
 267
 268         def x_fjc_plfix(params,f,pl_value=pl_value,T=T):
 269             '''
 270             fjc function for ODR fitting
 271             '''
 272             lambd=params
 273             pii=1/pl_value
 274             Kb=(1.38065e-23)
 275             therm=Kb*T
 276             #y=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
 277             x=(1/lambd)*(coth(f*(1/pii)/therm) - (therm*pii)/f)
 278             return x
 279
 280         #make the ODR fit
 281         realdata=scipy.odr.RealData(ychunk_corr_up,xchunk_corr_up)
 282         if pl_value:
 283             model=scipy.odr.Model(x_fjc_plfix)
 284             o = scipy.odr.ODR(realdata, model, p0_plfix)
 285         else:
 286             model=scipy.odr.Model(x_fjc)
 287             o = scipy.odr.ODR(realdata, model, p0)
 288
 289         o.set_job(fit_type=2)
 290         out=o.run()
 291         fit_out=[(1/i) for i in out.beta]
 292
 293         #Calculate fit errors from output standard deviations.
 294         #We must propagate the error because we fit the *inverse* parameters!
 295         #The error = (error of the inverse)*(value**2)
 296         fit_errors=[]
 297         for sd,value in zip(out.sd_beta, fit_out):
 298             err_real=sd*(value**2)
 299             fit_errors.append(err_real)
 300
 301         def fjc_eval(y,params,pl_value,T):
 302             '''
 303             Evaluates the WLC function
 304             '''
 305             if not pl_value:
 306                 lambd, pii = params
 307             else:
 308                 lambd = params
 309
 310             if pl_value:
 311                 pii=1/pl_value
 312
 313             Kb=(1.38065e-23) #boltzmann constant
 314             therm=Kb*T #so we have thermal energy
 315             #return ( (therm*pii/4.0) * (((1-(x*lambd))**-2.0) - 1 + (4.0*x*lambd)) )
 316             return (1/lambd)*(coth(y*(1/pii)/therm) - (therm*pii)/y)
 317
 318
 319
 320         #STEP 3: plotting the fit
 321         #obtain domain to plot the fit - from contact point to last_index plus 20 points
 322         thule_index=last_index+10
 323         if thule_index > len(xvector): #for rare cases in which we fit something at the END of whole curve.
 324             thule_index = len(xvector)
 325         #reverse etc. the domain
 326         ychunk=yvector[clicked_points[0].index:thule_index]
 327
 328         if len(ychunk)>0:
 329             y_evalchunk=np.linspace(min(ychunk),max(ychunk),100)
 330         else:
 331             #Empty y-chunk. It happens whenever we set the contact point after a recognized peak,
 332             #or other buggy situations. Kludge to live with it now...
 333             ychunk=yvector[:thule_index]
 334             y_evalchunk=np.linspace(min(ychunk),max(ychunk),100)
 335
 336         yfit_down=[-y for y in y_evalchunk]
 337         yfit_corr_down=[y+clicked_points[0].graph_coords[1] for y in yfit_down]
 338         yfit_corr_down=scipy.array(yfit_corr_down)
 339
 340         #the fitted curve: reflip, re-uncorrect
 341         xfit=fjc_eval(yfit_corr_down, out.beta, pl_value,T)
 342         xfit=list(xfit)
 343         xfit.reverse()
 344         xfit_chunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xfit]
 345
 346         #xfit_chunk_corr_up=scipy.array(xfit_chunk_corr_up)
 347         #deltay=yfit_down[0]-yvector[clicked_points[0].index]
 348
 349         #This is a terrible, terrible kludge to find the point where it should normalize (and from where it should plot)
 350         xxxdists=[]
 351         for index in scipy.arange(1,len(xfit_chunk_corr_up),1):
 352             xxxdists.append((clicked_points[0].graph_coords[0]-xfit_chunk_corr_up[index])**2)
 353         normalize_index=xxxdists.index(min(xxxdists))
 354         #End of kludge
 355
 356         deltay=yfit_down[normalize_index]-clicked_points[0].graph_coords[1]
 357         yfit_corr_down=[y-deltay for y in yfit_down]
 358
 359
 360         #calculate fit quality
 361         #creates dense y vector
 362         yqeval=np.linspace(np.min(ychunk_corr_up)/2,np.max(ychunk_corr_up)*2,10*len(ychunk_corr_up))
 363         #corresponding fitted x
 364         xqeval=fjc_eval(yqeval,out.beta,pl_value,T)
 365
 366         qsum=0
 367         for qindex in np.arange(0,len(ychunk_corr_up)):
 368             qsum+=dist(xchunk_corr_up[qindex],ychunk_corr_up[qindex],xqeval,yqeval)
 369         qstd=np.sqrt(qsum/len(ychunk_corr_up))
 370
 371
 372         if return_errors:
 373             return fit_out, yfit_corr_down[normalize_index+1:], xfit_chunk_corr_up[normalize_index+1:], fit_errors, qstd
 374         else:
 375             return fit_out, yfit_corr_down[normalize_index+1:], xfit_chunk_corr_up[normalize_index+1:], None, qstd
 376
 377     def efjc_fit(self,clicked_points,xvector,yvector, pl_value, T=293.0, return_errors=False):
 378         '''
 379         Extended Freely-jointed chain function
 380         ref: F Oesterhelt, M Rief and H E Gaub, New Journal of Physics 1 (1999) 6.1–6.11
 381         Please note that 2-parameter fitting (length and kl) usually does not converge, use fixed kl
 382         '''
 383         '''
 384         clicked_points[0] is the contact point (calculated or hand-clicked)
 385         clicked_points[1] and [2] are edges of chunk
 386
 387         '''
 388         #Fixed parameters from reference
 389         Kb=(1.38065e-2) #in pN.nm
 390         Lp=0.358 #planar, nm
 391         Lh=0.280 #helical, nm
 392         Ks=150e3  #pN/nm
 393
 394
 395         #STEP 1: Prepare the vectors to apply the fit.
 396
 397         #indexes of the selected chunk
 398         first_index=min(clicked_points[1].index, clicked_points[2].index)
 399         last_index=max(clicked_points[1].index, clicked_points[2].index)
 400
 401         #getting the chunk and reverting it
 402         xchunk,ychunk=xvector[first_index:last_index],yvector[first_index:last_index]
 403         xchunk.reverse()
 404         ychunk.reverse()
 405         #put contact point at zero and flip around the contact point (the fit wants a positive growth for extension and force)
 406         xchunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xchunk]
 407         ychunk_corr_up=[-(y-clicked_points[0].graph_coords[1]) for y in ychunk]
 408
 409
 410         #make them arrays
 411         xchunk_corr_up=scipy.array(xchunk_corr_up)
 412         ychunk_corr_up=scipy.array(ychunk_corr_up)
 413
 414         xchunk_corr_up_nm=xchunk_corr_up*1e9
 415         ychunk_corr_up_pn=ychunk_corr_up*1e12
 416
 417
 418         #STEP 2: actually do the fit
 419
 420         #Find furthest point of chunk and add it a bit; the fit must converge
 421         #from an excess!
 422         xchunk_high=max(xchunk_corr_up_nm)
 423         xchunk_high+=(xchunk_high/10.0)
 424
 425         #Here are the linearized start parameters for the WLC.
 426         #[Ns , pii=1/P]
 427         #max number of monomers (all helical)for a contour length xchunk_high
 428         excessNs=xchunk_high/(Lp)
 429         p0=[excessNs,(1.0/(0.7))]
 430         p0_plfix=[(excessNs)]
 431
 432         def dist(px,py,linex,liney):
 433             distancesx=scipy.array([(px-x)**2 for x in linex])
 434             minindex=np.argmin(distancesx)
 435             return (py-liney[minindex])**2
 436
 437         def deltaG(f):
 438             dG0=12.4242 #3kt in pN.nm
 439             dL=0.078 #planar-helical
 440             return dG0-f*dL
 441
 442         def Lfactor(f,T=T):
 443             Lp=0.358 #planar, nm
 444             Lh=0.280 #helical, nm
 445             Kb=(1.38065e-2)
 446             therm=Kb*T
 447             dG=deltaG(f)
 448
 449             return Lp/(np.exp(dG/therm)+1)+Lh/(np.exp(-dG/therm)+1)
 450
 451         def coth(z):
 452             '''
 453             hyperbolic cotangent
 454             '''
 455             return 1.0/np.tanh(z)
 456
 457         def x_efjc(params,f,T=T,Ks=Ks):
 458             '''
 459             efjc function for ODR fitting
 460             '''
 461
 462             Ns=params[0]
 463             invkl=params[1]
 464             Kb=(1.38065e-2)
 465             therm=Kb*T
 466             beta=(f/therm)/invkl
 467
 468             x=Ns*Lfactor(f)*(coth(beta)-1.0/beta)+Ns*f/Ks
 469             return x
 470
 471         def x_efjc_plfix(params,f,kl_value=pl_value,T=T,Ks=Ks):
 472             '''
 473             efjc function for ODR fitting
 474             '''
 475
 476             Ns=params
 477             invkl=1.0/kl_value
 478             Kb=(1.38065e-2)
 479             therm=Kb*T
 480             beta=(f/therm)/invkl
 481
 482             x=Ns*Lfactor(f)*(coth(beta)-1.0/beta)+Ns*f/Ks
 483             return x
 484
 485         #make the ODR fit
 486         realdata=scipy.odr.RealData(ychunk_corr_up_pn,xchunk_corr_up_nm)
 487         if pl_value:
 488             model=scipy.odr.Model(x_efjc_plfix)
 489             o = scipy.odr.ODR(realdata, model, p0_plfix)
 490         else:
 491             print 'WARNING eFJC fit with free pl sometimes does not converge'
 492             model=scipy.odr.Model(x_efjc)
 493             o = scipy.odr.ODR(realdata, model, p0)
 494
 495         o.set_job(fit_type=2)
 496         out=o.run()
 497
 498
 499         Ns=out.beta[0]
 500         Lc=Ns*Lp*1e-9
 501         if len(out.beta)>1:
 502             kfit=1e-9/out.beta[1]
 503             kfitnm=1/out.beta[1]
 504         else:
 505             kfit=1e-9*pl_value
 506             kfitnm=pl_value
 507
 508         fit_out=[Lc, kfit]
 509
 510         #Calculate fit errors from output standard deviations.
 511         fit_errors=[]
 512         fit_errors.append(out.sd_beta[0]*Lp*1e-9)
 513         if len(out.beta)>1:
 514             fit_errors.append(1e9*out.sd_beta[1]*kfit**2)
 515
 516
 517
 518         def efjc_eval(y,params,pl_value,T=T,Lfactor=Lfactor,Ks=Ks):
 519             '''
 520             Evaluates the eFJC function
 521             '''
 522             if not pl_value:
 523                 Ns, invkl = params
 524             else:
 525                 Ns = params
 526
 527             if pl_value:
 528                 invkl=1.0/pl_value
 529
 530             Kb=(1.38065e-2) #boltzmann constant
 531             therm=Kb*T #so we have thermal energy
 532             beta=(y/therm)/invkl
 533
 534             x= Ns*Lfactor(y)*(coth(beta)-1.0/beta)+Ns*y/Ks
 535
 536             return x
 537
 538         #STEP 3: plotting the fit
 539         #obtain domain to plot the fit - from contact point to last_index plus 20 points
 540         thule_index=last_index+10
 541         if thule_index > len(xvector): #for rare cases in which we fit something at the END of whole curve.
 542             thule_index = len(xvector)
 543         #reverse etc. the domain
 544         ychunk=yvector[clicked_points[0].index:thule_index]
 545
 546         if len(ychunk)>0:
 547             y_evalchunk=np.linspace(min(ychunk),max(ychunk),100)
 548         else:
 549             #Empty y-chunk. It happens whenever we set the contact point after a recognized peak,
 550             #or other buggy situations. Kludge to live with it now...
 551             ychunk=yvector[:thule_index]
 552             y_evalchunk=np.linspace(min(ychunk),max(ychunk),100)
 553
 554         yfit_down=[-y for y in y_evalchunk]
 555         yfit_corr_down=[y+clicked_points[0].graph_coords[1] for y in yfit_down]
 556         yfit_corr_down=scipy.array(yfit_corr_down)
 557
 558         #the fitted curve: reflip, re-uncorrect
 559         xfit=efjc_eval(1e12*yfit_corr_down, out.beta, pl_value,T)*1e-9
 560         xfit=list(xfit)
 561         xfit.reverse()
 562         xfit_chunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xfit]
 563
 564         #xfit_chunk_corr_up=scipy.array(xfit_chunk_corr_up)
 565         #deltay=yfit_down[0]-yvector[clicked_points[0].index]
 566
 567         #This is a terrible, terrible kludge to find the point where it should normalize (and from where it should plot)
 568         xxxdists=[]
 569         for index in scipy.arange(1,len(xfit_chunk_corr_up),1):
 570             xxxdists.append((clicked_points[0].graph_coords[0]-xfit_chunk_corr_up[index])**2)
 571         normalize_index=xxxdists.index(min(xxxdists))
 572         #End of kludge
 573
 574         deltay=yfit_down[normalize_index]-clicked_points[0].graph_coords[1]
 575         yfit_corr_down=[y-deltay for y in yfit_down]
 576
 577         #calculate fit quality
 578         #creates dense y vector
 579         yqeval=np.linspace(np.min(ychunk_corr_up_pn)/2,np.max(ychunk_corr_up_pn)*2,10*len(ychunk_corr_up_pn))
 580         #corresponding fitted x
 581         xqeval=efjc_eval(yqeval,out.beta,pl_value)
 582
 583         qsum=0
 584         for qindex in np.arange(0,len(ychunk_corr_up_pn)):
 585             qsum+=dist(xchunk_corr_up_nm[qindex],ychunk_corr_up_pn[qindex],xqeval,yqeval)
 586         qstd=1e-12*np.sqrt(qsum/len(ychunk_corr_up_pn))
 587
 588         if return_errors:
 589             return fit_out, yfit_corr_down[normalize_index+1:], xfit_chunk_corr_up[normalize_index+1:], fit_errors, qstd
 590         else:
 591             return fit_out, yfit_corr_down[normalize_index+1:], xfit_chunk_corr_up[normalize_index+1:], None, qstd
 592
 593
 594
 595
 596     def do_wlc(self,args):
 597         '''
 598         WLC
 599         (fit.py plugin)
 600
 601         See the fit command
 602         '''
 603         self.do_fit(args)
 604
 605     def do_fjc(self,args):
 606         '''
 607         FJC
 608         (fit.py plugin)
 609
 610         See the fit command
 611         '''
 612         self.do_fit(args)
 613
 614     def do_fit(self,args):
 615         '''
 616         FIT
 617         (fit.py plugin)
 618         Fits an entropic elasticity function to a given chunk of the curve.
 619
 620         First you have to click a contact point.
 621         Then you have to click the two edges of the data you want to fit.
 622
 623         Fit quality compares the distance to the fit with the thermal noise (a good fit should be close to 1)
 624
 625         The fit function depends on the fit_function variable. You can set it with the command
 626         "set fit_function wlc" or  "set fit_function fjc" depending on the function you prefer.
 627
 628         For WLC, the function is the simple polynomial worm-like chain as proposed by
 629         C.Bustamante, J.F.Marko, E.D.Siggia and S.Smith (Science. 1994
 630         Sep 9;265(5178):1599-600.)
 631
 632         For FJC, ref:
 633         C.Ray and B.B. Akhremitchev; http://www.chem.duke.edu/~boris/research/force_spectroscopy/fit_efjc.pdf
 634
 635         For eFJC, ref:
 636         F Oesterhelt, M Rief and H E Gaub, New Journal of Physics 1 (1999) 6.1–6.11 (section 4.2)
 637         NOTE: use fixed pl for better results.
 638
 639         Arguments:
 640         pl=[value] : Use a fixed persistent length (WLC) or Kuhn length (FJC) for the fit. If pl is not given,
 641                      the fit will be a 2-variable
 642                      fit. DO NOT put spaces between 'pl', '=' and the value.
 643                      The value must be in nanometers.
 644
 645         t=[value] : Use a user-defined temperature. The value must be in
 646                     kelvins; by default it is 293 K.
 647                     DO NOT put spaces between 't', '=' and the value.
 648
 649         noauto : allows for clicking the contact point by
 650                  hand (otherwise it is automatically estimated) the first time.
 651                  If subsequent measurements are made, the same contact point
 652                  clicked is used
 653
 654         reclick : redefines by hand the contact point, if noauto has been used before
 655                   but the user is unsatisfied of the previously choosen contact point.
 656         ---------
 657         Syntax: fit [pl=(value)] [t=value] [noauto]
 658         '''
 659         pl_value=None
 660         T=self.config['temperature']
 661         for arg in args.split():
 662             #look for a persistent length argument.
 663             if 'pl=' in arg:
 664                 pl_expression=arg.split('=')
 665                 pl_value=float(pl_expression[1]) #actual value
 666             #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
 667             if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
 668                 t_expression=arg.split('=')
 669                 T=float(t_expression[1])
 670
 671         #use the currently displayed plot for the fit
 672         displayed_plot=self._get_displayed_plot()
 673
 674         #handle contact point arguments correctly
 675         if 'reclick' in args.split():
 676             print 'Click contact point'
 677             contact_point=self._measure_N_points(N=1, whatset=1)[0]
 678             contact_point_index=contact_point.index
 679             self.wlccontact_point=contact_point
 680             self.wlccontact_index=contact_point.index
 681             self.wlccurrent=self.current.path
 682         elif 'noauto' in args.split():
 683             if self.wlccontact_index==None or self.wlccurrent != self.current.path:
 684                 print 'Click contact point'
 685                 contact_point=self._measure_N_points(N=1, whatset=1)[0]
 686                 contact_point_index=contact_point.index
 687                 self.wlccontact_point=contact_point
 688                 self.wlccontact_index=contact_point.index
 689                 self.wlccurrent=self.current.path
 690             else:
 691                 contact_point=self.wlccontact_point
 692                 contact_point_index=self.wlccontact_index
 693         else:
 694             cindex=self.find_contact_point()
 695             contact_point=ClickedPoint()
 696             contact_point.absolute_coords=displayed_plot.vectors[1][0][cindex], displayed_plot.vectors[1][1][cindex]
 697             contact_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
 698             contact_point.is_marker=True
 699
 700         print 'Click edges of chunk'
 701         points=self._measure_N_points(N=2, whatset=1)
 702         points=[contact_point]+points
 703
 704         try:
 705             if self.config['fit_function']=='wlc':
 706                 params, yfit, xfit, fit_errors,qstd = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T, return_errors=True )
 707                 name_of_charlength='Persistent length'
 708             elif self.config['fit_function']=='fjc':
 709                 params, yfit, xfit, fit_errors,qstd = self.fjc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T, return_errors=True )
 710                 name_of_charlength='Kuhn length'
 711             elif self.config['fit_function']=='efjc':
 712                 params, yfit, xfit, fit_errors,qstd = self.efjc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T, return_errors=True )
 713                 name_of_charlength='Kuhn length (e)'
 714             else:
 715                 print 'No recognized fit function defined!'
 716                 print 'Set your fit function to wlc, fjc or efjc.'
 717                 return
 718
 719         except:
 720             print 'Fit not possible. Probably wrong interval -did you click two *different* points?'
 721             return
 722
 723         #FIXME: print "Kuhn length" for FJC
 724         print 'Fit function:',self.config['fit_function']
 725         print 'Contour length: ',params[0]*(1.0e+9),' nm'
 726         to_dump='contour '+self.current.path+' '+str(params[0]*(1.0e+9))+' nm'
 727         self.outlet.push(to_dump)
 728         if len(params)==2: #if we did choose 2-value fit
 729             print name_of_charlength+': ',params[1]*(1.0e+9),' nm'
 730             to_dump='persistent '+self.current.path+' '+str(params[1]*(1.0e+9))+' nm'
 731             self.outlet.push(to_dump)
 732
 733         if fit_errors:
 734             fit_nm=[i*(10**9) for i in fit_errors]
 735             print 'Standard deviation (contour length)', fit_nm[0]
 736             if len(fit_nm)>1:
 737                 print 'Standard deviation ('+name_of_charlength+')', fit_nm[1]
 738
 739         print 'Fit quality: '+str(qstd/np.std(displayed_plot.vectors[1][1][-20:-1]))
 740
 741
 742         #add the clicked points in the final PlotObject
 743         clickvector_x, clickvector_y=[], []
 744         for item in points:
 745             clickvector_x.append(item.graph_coords[0])
 746             clickvector_y.append(item.graph_coords[1])
 747
 748         #create a custom PlotObject to gracefully plot the fit along the curves
 749
 750         fitplot=copy.deepcopy(displayed_plot)
 751         fitplot.add_set(xfit,yfit)
 752         fitplot.add_set(clickvector_x,clickvector_y)
 753
 754         #FIXME: this colour/styles stuff must be solved at the root!
 755         if fitplot.styles==[]:
 756             fitplot.styles=[None,None,None,'scatter']
 757         else:
 758             fitplot.styles+=[None,'scatter']
 759
 760         if fitplot.colors==[]:
 761             fitplot.colors=[None,None,None,None]
 762         else:
 763             fitplot.colors+=[None,None]
 764
 765         self._send_plot([fitplot])
 766
 767
 768
 769     #----------
 770
 771
 772     def find_contact_point(self,plot=False):
 773         '''
 774         Finds the contact point on the curve.
 775
 776         The current algorithm (thanks to Francesco Musiani, francesco.musiani@unibo.it and Massimo Sandal) is:
 777         - take care of the PicoForce trigger bug - exclude retraction portions with too high standard deviation
 778         - fit the second half of the retraction curve to a line
 779         - if the fit is not almost horizontal, take a smaller chunk and repeat
 780         - otherwise, we have something horizontal
 781         - so take the average of horizontal points and use it as a baseline
 782
 783         Then, start from the rise of the retraction curve and look at the first point below the
 784         baseline.
 785
 786         FIXME: should be moved, probably to generalvclamp.py
 787         '''
 788
 789         if not plot:
 790             plot=self.plots[0]
 791
 792         outplot=self.subtract_curves(1)
 793         xret=outplot.vectors[1][0]
 794         ydiff=outplot.vectors[1][1]
 795
 796         xext=plot.vectors[0][0]
 797         yext=plot.vectors[0][1]
 798         xret2=plot.vectors[1][0]
 799         yret=plot.vectors[1][1]
 800
 801         #taking care of the picoforce trigger bug: we exclude portions of the curve that have too much
 802         #standard deviation. yes, a lot of magic is here.
 803         monster=True
 804         monlength=len(xret)-int(len(xret)/20)
 805         finalength=len(xret)
 806         while monster:
 807             monchunk=scipy.array(ydiff[monlength:finalength])
 808             if abs(np.std(monchunk)) < 2e-10:
 809                 monster=False
 810             else: #move away from the monster
 811                 monlength-=int(len(xret)/50)
 812                 finalength-=int(len(xret)/50)
 813
 814
 815         #take half of the thing
 816         endlength=int(len(xret)/2)
 817
 818         ok=False
 819
 820         while not ok:
 821             xchunk=yext[endlength:monlength]
 822             ychunk=yext[endlength:monlength]
 823             regr=scipy.stats.linregress(xchunk,ychunk)[0:2]
 824             #we stop if we found an almost-horizontal fit or if we're going too short...
 825             #FIXME: 0.1 and 6 here are "magic numbers" (although reasonable)
 826             if (abs(regr[1]) > 0.1) and ( endlength < len(xret)-int(len(xret)/6) ) :
 827                 endlength+=10
 828             else:
 829                 ok=True
 830
 831
 832         ymean=np.mean(ychunk) #baseline
 833
 834         index=0
 835         point = ymean+1
 836
 837         #find the first point below the calculated baseline
 838         while point > ymean:
 839             try:
 840                 point=yret[index]
 841                 index+=1
 842             except IndexError:
 843                 #The algorithm didn't find anything below the baseline! It should NEVER happen
 844                 index=0
 845                 return index
 846
 847         return index
 848
 849
 850
 851     def find_contact_point2(self, debug=False):
 852         '''
 853         TO BE DEVELOPED IN THE FUTURE
 854         Finds the contact point on the curve.
 855
 856         FIXME: should be moved, probably to generalvclamp.py
 857         '''
 858
 859         #raw_plot=self.current.curve.default_plots()[0]
 860         raw_plot=self.plots[0]
 861         '''xext=self.plots[0].vectors[0][0]
 862         yext=self.plots[0].vectors[0][1]
 863         xret2=self.plots[0].vectors[1][0]
 864         yret=self.plots[0].vectors[1][1]
 865         '''
 866         xext=raw_plot.vectors[0][0]
 867         yext=raw_plot.vectors[0][1]
 868         xret2=raw_plot.vectors[1][0]
 869         yret=raw_plot.vectors[1][1]
 870
 871         first_point=[xext[0], yext[0]]
 872         last_point=[xext[-1], yext[-1]]
 873
 874         #regr=scipy.polyfit(first_point, last_point,1)[0:2]
 875         diffx=abs(first_point[0]-last_point[0])
 876         diffy=abs(first_point[1]-last_point[1])
 877
 878         #using polyfit results in numerical errors. good old algebra.
 879         a=diffy/diffx
 880         b=first_point[1]-(a*first_point[0])
 881         baseline=scipy.polyval((a,b), xext)
 882
 883         ysub=[item-basitem for item,basitem in zip(yext,baseline)]
 884
 885         contact=ysub.index(min(ysub))
 886
 887         return xext,ysub,contact
 888
 889         #now, exploit a ClickedPoint instance to calculate index...
 890         dummy=ClickedPoint()
 891         dummy.absolute_coords=(x_intercept,y_intercept)
 892         dummy.find_graph_coords(xret2,yret)
 893
 894         if debug:
 895             return dummy.index, regr, regr_contact
 896         else:
 897             return dummy.index
 898
 899
 900
 901     def x_do_contact(self,args):
 902         '''
 903         DEBUG COMMAND to be activated in the future
 904         '''
 905         xext,ysub,contact=self.find_contact_point2(debug=True)
 906
 907         contact_plot=self.plots[0]
 908         contact_plot.add_set(xext,ysub)
 909         contact_plot.add_set([xext[contact]],[self.plots[0].vectors[0][1][contact]])
 910         #contact_plot.add_set([first_point[0]],[first_point[1]])
 911         #contact_plot.add_set([last_point[0]],[last_point[1]])
 912         contact_plot.styles=[None,None,None,'scatter']
 913         self._send_plot([contact_plot])
 914         return
 915
 916
 917         index,regr,regr_contact=self.find_contact_point2(debug=True)
 918         print regr
 919         print regr_contact
 920         raw_plot=self.current.curve.default_plots()[0]
 921         xret=raw_plot.vectors[0][0]
 922         #nc_line=[(item*regr[0])+regr[1] for item in x_nc]
 923         nc_line=scipy.polyval(regr,xret)
 924         c_line=scipy.polyval(regr_contact,xret)
 925
 926
 927         contact_plot=self.current.curve.default_plots()[0]
 928         contact_plot.add_set(xret, nc_line)
 929         contact_plot.add_set(xret, c_line)
 930         contact_plot.styles=[None,None,None,None]
 931         #contact_plot.styles.append(None)
 932         contact_plot.destination=1
 933         self._send_plot([contact_plot])