plugins/fit.py

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