hooke/plugin/fit.py

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