hooke/plugin/fit.py

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