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