fit.py

   1 #!/usr/bin/env python
   2
   3 '''
   4 FIT
   5
   6 Force spectroscopy curves basic fitting plugin.
   7 Licensed under the GNU GPL version 2
   8
   9 Non-standard Dependencies:
  10 procplots.py (plot processing plugin)
  11 '''
  12 from libhooke import WX_GOOD, ClickedPoint
  13 import wxversion
  14 wxversion.select(WX_GOOD)
  15 #from wx import PostEvent
  16 #from wx.lib.newevent import NewEvent
  17 import scipy
  18 import scipy.odr
  19 import numpy as np
  20 import copy
  21 import Queue
  22
  23 global measure_wlc
  24 global EVT_MEASURE_WLC
  25
  26 #measure_wlc, EVT_MEASURE_WLC = NewEvent()
  27
  28 global events_from_fit
  29 events_from_fit=Queue.Queue() #GUI ---> CLI COMMUNICATION
  30
  31
  32 class fitCommands:
  33
  34     def _plug_init(self):
  35         self.wlccurrent=None
  36         self.wlccontact_point=None
  37         self.wlccontact_index=None
  38
  39     def wlc_fit(self,clicked_points,xvector,yvector, pl_value, T=293, return_errors=False):
  40         '''
  41         Worm-like chain model fitting.
  42         The function is the simple polynomial worm-like chain as proposed by C.Bustamante, J.F.Marko, E.D.Siggia
  43         and S.Smith (Science. 1994 Sep 9;265(5178):1599-600.)
  44         '''
  45
  46         '''clicked_points[0] = contact point (calculated or hand-clicked)
  47         clicked_points[1] and [2] are edges of chunk'''
  48
  49         #STEP 1: Prepare the vectors to apply the fit.
  50
  51         if pl_value is not None:
  52             pl_value=pl_value/(10**9)
  53
  54         #indexes of the selected chunk
  55         first_index=min(clicked_points[1].index, clicked_points[2].index)
  56         last_index=max(clicked_points[1].index, clicked_points[2].index)
  57
  58         #getting the chunk and reverting it
  59         xchunk,ychunk=xvector[first_index:last_index],yvector[first_index:last_index]
  60         xchunk.reverse()
  61         ychunk.reverse()
  62         #put contact point at zero and flip around the contact point (the fit wants a positive growth for extension and force)
  63         xchunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xchunk]
  64         ychunk_corr_up=[-(y-clicked_points[0].graph_coords[1]) for y in ychunk]
  65
  66         #make them arrays
  67         xchunk_corr_up=scipy.array(xchunk_corr_up)
  68         ychunk_corr_up=scipy.array(ychunk_corr_up)
  69
  70
  71         #STEP 2: actually do the fit
  72
  73         #Find furthest point of chunk and add it a bit; the fit must converge
  74         #from an excess!
  75         xchunk_high=max(xchunk_corr_up)
  76         xchunk_high+=(xchunk_high/10)
  77
  78         #Here are the linearized start parameters for the WLC.
  79         #[lambd=1/Lo , pii=1/P]
  80
  81         p0=[(1/xchunk_high),(1/(3.5e-10))]
  82         p0_plfix=[(1/xchunk_high)]
  83         '''
  84         ODR STUFF
  85         fixme: remove these comments after testing
  86         '''
  87
  88
  89         def f_wlc(params,x,T=T):
  90             '''
  91             wlc function for ODR fitting
  92             '''
  93             lambd,pii=params
  94             Kb=(1.38065e-23)
  95             therm=Kb*T
  96             y=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
  97             return y
  98
  99         def f_wlc_plfix(params,x,pl_value=pl_value,T=T):
 100             '''
 101             wlc function for ODR fitting
 102             '''
 103             lambd=params
 104             pii=1/pl_value
 105             Kb=(1.38065e-23)
 106             therm=Kb*T
 107             y=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
 108             return y
 109
 110         #make the ODR fit
 111         realdata=scipy.odr.RealData(xchunk_corr_up,ychunk_corr_up)
 112         if pl_value:
 113             model=scipy.odr.Model(f_wlc_plfix)
 114             o = scipy.odr.ODR(realdata, model, p0_plfix)
 115         else:
 116             model=scipy.odr.Model(f_wlc)
 117             o = scipy.odr.ODR(realdata, model, p0)
 118
 119         o.set_job(fit_type=2)
 120         out=o.run()
 121         fit_out=[(1/i) for i in out.beta]
 122
 123         #Calculate fit errors from output standard deviations.
 124         #We must propagate the error because we fit the *inverse* parameters!
 125         #The error = (error of the inverse)*(value**2)
 126         fit_errors=[]
 127         for sd,value in zip(out.sd_beta, fit_out):
 128             err_real=sd*(value**2)
 129             fit_errors.append(err_real)
 130
 131         def wlc_eval(x,params,pl_value,T):
 132             '''
 133             Evaluates the WLC function
 134             '''
 135             if not pl_value:
 136                 lambd, pii = params
 137             else:
 138                 lambd = params
 139
 140             if pl_value:
 141                 pii=1/pl_value
 142
 143             Kb=(1.38065e-23) #boltzmann constant
 144             therm=Kb*T #so we have thermal energy
 145
 146             return ( (therm*pii/4.0) * (((1-(x*lambd))**-2.0) - 1 + (4.0*x*lambd)) )
 147
 148         #STEP 3: plotting the fit
 149
 150         #obtain domain to plot the fit - from contact point to last_index plus 20 points
 151         thule_index=last_index+10
 152         if thule_index > len(xvector): #for rare cases in which we fit something at the END of whole curve.
 153             thule_index = len(xvector)
 154         #reverse etc. the domain
 155         xfit_chunk=xvector[clicked_points[0].index:thule_index]
 156         xfit_chunk.reverse()
 157         xfit_chunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xfit_chunk]
 158         xfit_chunk_corr_up=scipy.array(xfit_chunk_corr_up)
 159
 160         #the fitted curve: reflip, re-uncorrect
 161         yfit=wlc_eval(xfit_chunk_corr_up, out.beta, pl_value,T)
 162         yfit_down=[-y for y in yfit]
 163         yfit_corr_down=[y+clicked_points[0].graph_coords[1] for y in yfit_down]
 164
 165         if return_errors:
 166             return fit_out, yfit_corr_down, xfit_chunk, fit_errors
 167         else:
 168             return fit_out, yfit_corr_down, xfit_chunk, None
 169
 170
 171     def do_wlc(self,args):
 172         '''
 173         WLC
 174         (fit plugin)
 175         Fits a worm-like chain entropic rise to a given chunk of the curve.
 176
 177         First you have to click a contact point.
 178         Then you have to click the two edges of the data you want to fit.
 179         The function is the simple polynomial worm-like chain as proposed by
 180         C.Bustamante, J.F.Marko, E.D.Siggia and S.Smith (Science. 1994
 181         Sep 9;265(5178):1599-600.)
 182
 183         Arguments:
 184         pl=[value] : Use a fixed persistent length for the fit. If pl is not given,
 185                      the fit will be a 2-variable
 186                      fit. DO NOT put spaces between 'pl', '=' and the value.
 187                      The value must be in nanometers.
 188
 189         t=[value] : Use a user-defined temperature. The value must be in
 190                     kelvins; by default it is 293 K.
 191                     DO NOT put spaces between 't', '=' and the value.
 192
 193         noauto : allows for clicking the contact point by
 194                  hand (otherwise it is automatically estimated) the first time.
 195                  If subsequent measurements are made, the same contact point
 196                  clicked is used
 197
 198         reclick : redefines by hand the contact point, if noauto has been used before
 199                   but the user is unsatisfied of the previously choosen contact point.
 200         ---------
 201         Syntax: wlc [pl=(value)] [t=value] [noauto]
 202         '''
 203         pl_value=None
 204         T=self.config['temperature']
 205         for arg in args.split():
 206             #look for a persistent length argument.
 207             if 'pl=' in arg:
 208                 pl_expression=arg.split('=')
 209                 pl_value=float(pl_expression[1]) #actual value
 210             #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
 211             if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
 212                 t_expression=arg.split('=')
 213                 T=float(t_expression[1])
 214
 215         #use the currently displayed plot for the fit
 216         displayed_plot=self._get_displayed_plot()
 217
 218         #handle contact point arguments correctly
 219         if 'reclick' in args.split():
 220             print 'Click contact point'
 221             contact_point=self._measure_N_points(N=1, whatset=1)[0]
 222             contact_point_index=contact_point.index
 223             self.wlccontact_point=contact_point
 224             self.wlccontact_index=contact_point.index
 225             self.wlccurrent=self.current.path
 226         elif 'noauto' in args.split():
 227             if self.wlccontact_index==None or self.wlccurrent != self.current.path:
 228                 print 'Click contact point'
 229                 contact_point=self._measure_N_points(N=1, whatset=1)[0]
 230                 contact_point_index=contact_point.index
 231                 self.wlccontact_point=contact_point
 232                 self.wlccontact_index=contact_point.index
 233                 self.wlccurrent=self.current.path
 234             else:
 235                 contact_point=self.wlccontact_point
 236                 contact_point_index=self.wlccontact_index
 237         else:
 238             cindex=self.find_contact_point()
 239             contact_point=ClickedPoint()
 240             contact_point.absolute_coords=displayed_plot.vectors[1][0][cindex], displayed_plot.vectors[1][1][cindex]
 241             contact_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
 242             contact_point.is_marker=True
 243
 244         print 'Click edges of chunk'
 245         points=self._measure_N_points(N=2, whatset=1)
 246         points=[contact_point]+points
 247         try:
 248             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 )
 249         except:
 250             print 'Fit not possible. Probably wrong interval -did you click two *different* points?'
 251             return
 252
 253         print 'Contour length: ',params[0]*(1.0e+9),' nm'
 254         to_dump='contour '+self.current.path+' '+str(params[0]*(1.0e+9))+' nm'
 255         self.outlet.push(to_dump)
 256         if len(params)==2: #if we did choose 2-value fit
 257             print 'Persistent length: ',params[1]*(1.0e+9),' nm'
 258             to_dump='persistent '+self.current.path+' '+str(params[1]*(1.0e+9))+' nm'
 259             self.outlet.push(to_dump)
 260
 261         if fit_errors:
 262             fit_nm=[i*(10**9) for i in fit_errors]
 263             print 'Standard deviation (contour length)', fit_nm[0]
 264             if len(fit_nm)>1:
 265                 print 'Standard deviation (persistent length)', fit_nm[1]
 266
 267
 268         #add the clicked points in the final PlotObject
 269         clickvector_x, clickvector_y=[], []
 270         for item in points:
 271             clickvector_x.append(item.graph_coords[0])
 272             clickvector_y.append(item.graph_coords[1])
 273
 274         #create a custom PlotObject to gracefully plot the fit along the curves
 275
 276         fitplot=copy.deepcopy(displayed_plot)
 277         fitplot.add_set(xfit,yfit)
 278         fitplot.add_set(clickvector_x,clickvector_y)
 279
 280         if fitplot.styles==[]:
 281             fitplot.styles=[None,None,None,'scatter']
 282         else:
 283             fitplot.styles+=[None,'scatter']
 284
 285         self._send_plot([fitplot])
 286
 287     def find_contact_point(self,plot=False):
 288         '''
 289         Finds the contact point on the curve.
 290
 291         The current algorithm (thanks to Francesco Musiani, francesco.musiani@unibo.it and Massimo Sandal) is:
 292         - take care of the PicoForce trigger bug - exclude retraction portions with too high standard deviation
 293         - fit the second half of the retraction curve to a line
 294         - if the fit is not almost horizontal, take a smaller chunk and repeat
 295         - otherwise, we have something horizontal
 296         - so take the average of horizontal points and use it as a baseline
 297
 298         Then, start from the rise of the retraction curve and look at the first point below the
 299         baseline.
 300
 301         FIXME: should be moved, probably to generalvclamp.py
 302         '''
 303
 304         if not plot:
 305             plot=self.plots[0]
 306
 307         outplot=self.subtract_curves(1)
 308         xret=outplot.vectors[1][0]
 309         ydiff=outplot.vectors[1][1]
 310
 311         xext=plot.vectors[0][0]
 312         yext=plot.vectors[0][1]
 313         xret2=plot.vectors[1][0]
 314         yret=plot.vectors[1][1]
 315
 316         #taking care of the picoforce trigger bug: we exclude portions of the curve that have too much
 317         #standard deviation. yes, a lot of magic is here.
 318         monster=True
 319         monlength=len(xret)-int(len(xret)/20)
 320         finalength=len(xret)
 321         while monster:
 322             monchunk=scipy.array(ydiff[monlength:finalength])
 323             if abs(scipy.stats.std(monchunk)) < 2e-10:
 324                 monster=False
 325             else: #move away from the monster
 326                 monlength-=int(len(xret)/50)
 327                 finalength-=int(len(xret)/50)
 328
 329
 330         #take half of the thing
 331         endlength=int(len(xret)/2)
 332
 333         ok=False
 334
 335         while not ok:
 336             xchunk=yext[endlength:monlength]
 337             ychunk=yext[endlength:monlength]
 338             regr=scipy.stats.linregress(xchunk,ychunk)[0:2]
 339             #we stop if we found an almost-horizontal fit or if we're going too short...
 340             #FIXME: 0.1 and 6 here are "magic numbers" (although reasonable)
 341             if (abs(regr[1]) > 0.1) and ( endlength < len(xret)-int(len(xret)/6) ) :
 342                 endlength+=10
 343             else:
 344                 ok=True
 345
 346
 347         ymean=scipy.mean(ychunk) #baseline
 348
 349         index=0
 350         point = ymean+1
 351
 352         #find the first point below the calculated baseline
 353         while point > ymean:
 354             try:
 355                 point=yret[index]
 356                 index+=1
 357             except IndexError:
 358                 #The algorithm didn't find anything below the baseline! It should NEVER happen
 359                 index=0
 360                 return index
 361
 362         return index
 363
 364
 365
 366     def find_contact_point2(self, debug=False):
 367         '''
 368         TO BE DEVELOPED IN THE FUTURE
 369         Finds the contact point on the curve.
 370
 371         FIXME: should be moved, probably to generalvclamp.py
 372         '''
 373
 374         #raw_plot=self.current.curve.default_plots()[0]
 375         raw_plot=self.plots[0]
 376         '''xext=self.plots[0].vectors[0][0]
 377         yext=self.plots[0].vectors[0][1]
 378         xret2=self.plots[0].vectors[1][0]
 379         yret=self.plots[0].vectors[1][1]
 380         '''
 381         xext=raw_plot.vectors[0][0]
 382         yext=raw_plot.vectors[0][1]
 383         xret2=raw_plot.vectors[1][0]
 384         yret=raw_plot.vectors[1][1]
 385
 386         first_point=[xext[0], yext[0]]
 387         last_point=[xext[-1], yext[-1]]
 388
 389         #regr=scipy.polyfit(first_point, last_point,1)[0:2]
 390         diffx=abs(first_point[0]-last_point[0])
 391         diffy=abs(first_point[1]-last_point[1])
 392
 393         #using polyfit results in numerical errors. good old algebra.
 394         a=diffy/diffx
 395         b=first_point[1]-(a*first_point[0])
 396         baseline=scipy.polyval((a,b), xext)
 397
 398         ysub=[item-basitem for item,basitem in zip(yext,baseline)]
 399
 400         contact=ysub.index(min(ysub))
 401
 402         return xext,ysub,contact
 403
 404         #now, exploit a ClickedPoint instance to calculate index...
 405         dummy=ClickedPoint()
 406         dummy.absolute_coords=(x_intercept,y_intercept)
 407         dummy.find_graph_coords(xret2,yret)
 408
 409         if debug:
 410             return dummy.index, regr, regr_contact
 411         else:
 412             return dummy.index
 413
 414
 415
 416     def x_do_contact(self,args):
 417         '''
 418         DEBUG COMMAND to be activated in the future
 419         '''
 420         xext,ysub,contact=self.find_contact_point2(debug=True)
 421
 422         contact_plot=self.plots[0]
 423         contact_plot.add_set(xext,ysub)
 424         contact_plot.add_set([xext[contact]],[self.plots[0].vectors[0][1][contact]])
 425         #contact_plot.add_set([first_point[0]],[first_point[1]])
 426         #contact_plot.add_set([last_point[0]],[last_point[1]])
 427         contact_plot.styles=[None,None,None,'scatter']
 428         self._send_plot([contact_plot])
 429         return
 430
 431
 432         index,regr,regr_contact=self.find_contact_point2(debug=True)
 433         print regr
 434         print regr_contact
 435         raw_plot=self.current.curve.default_plots()[0]
 436         xret=raw_plot.vectors[0][0]
 437         #nc_line=[(item*regr[0])+regr[1] for item in x_nc]
 438         nc_line=scipy.polyval(regr,xret)
 439         c_line=scipy.polyval(regr_contact,xret)
 440
 441
 442         contact_plot=self.current.curve.default_plots()[0]
 443         contact_plot.add_set(xret, nc_line)
 444         contact_plot.add_set(xret, c_line)
 445         contact_plot.styles=[None,None,None,None]
 446         #contact_plot.styles.append(None)
 447         contact_plot.destination=1
 448         self._send_plot([contact_plot])
 449