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