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