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             #T=293 #temperature FIXME:should be user-modifiable!
 145             therm=Kb*T #so we have thermal energy
 146
 147             return ( (therm*pii/4.0) * (((1-(x*lambd))**-2.0) - 1 + (4.0*x*lambd)) )
 148
 149         #STEP 3: plotting the fit
 150
 151         #obtain domain to plot the fit - from contact point to last_index plus 20 points
 152         thule_index=last_index+10
 153         if thule_index > len(xvector): #for rare cases in which we fit something at the END of whole curve.
 154             thule_index = len(xvector)
 155         #reverse etc. the domain
 156         xfit_chunk=xvector[clicked_points[0].index:thule_index]
 157         xfit_chunk.reverse()
 158         xfit_chunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xfit_chunk]
 159         xfit_chunk_corr_up=scipy.array(xfit_chunk_corr_up)
 160
 161         #the fitted curve: reflip, re-uncorrect
 162         yfit=wlc_eval(xfit_chunk_corr_up, out.beta, pl_value,T)
 163         yfit_down=[-y for y in yfit]
 164         yfit_corr_down=[y+clicked_points[0].graph_coords[1] for y in yfit_down]
 165
 166         if return_errors:
 167             return fit_out, yfit_corr_down, xfit_chunk, fit_errors
 168         else:
 169             return fit_out, yfit_corr_down, xfit_chunk, None
 170
 171
 172     def do_wlc(self,args):
 173         '''
 174         WLC
 175         (fit plugin)
 176         Fits a worm-like chain entropic rise to a given chunk of the curve.
 177
 178         First you have to click a contact point.
 179         Then you have to click the two edges of the data you want to fit.
 180         The function is the simple polynomial worm-like chain as proposed by
 181         C.Bustamante, J.F.Marko, E.D.Siggia and S.Smith (Science. 1994
 182         Sep 9;265(5178):1599-600.)
 183
 184         Arguments:
 185         pl=[value] : Use a fixed persistent length for the fit. If pl is not given,
 186                      the fit will be a 2-variable
 187                      fit. DO NOT put spaces between 'pl', '=' and the value.
 188                      The value must be in nanometers.
 189
 190         t=[value] : Use a user-defined temperature. The value must be in
 191                     kelvins; by default it is 293 K.
 192                     DO NOT put spaces between 't', '=' and the value.
 193
 194         noauto : allows for clicking the contact point by
 195                  hand (otherwise it is automatically estimated) the first time.
 196                  If subsequent measurements are made, the same contact point
 197                  clicked is used
 198
 199         reclick : redefines by hand the contact point, if noauto has been used before
 200                   but the user is unsatisfied of the previously choosen contact point.
 201         ---------
 202         Syntax: wlc [pl=(value)] [t=value] [noauto]
 203         '''
 204         pl_value=None
 205         T=self.config['temperature']
 206         for arg in args.split():
 207             #look for a persistent length argument.
 208             if 'pl=' in arg:
 209                 pl_expression=arg.split('=')
 210                 pl_value=float(pl_expression[1]) #actual value
 211             #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
 212             if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
 213                 t_expression=arg.split('=')
 214                 T=float(t_expression[1])
 215
 216         #use the currently displayed plot for the fit
 217         displayed_plot=self._get_displayed_plot()
 218
 219         #handle contact point arguments correctly
 220         if 'reclick' in args.split():
 221             print 'Click contact point'
 222             contact_point=self._measure_N_points(N=1, whatset=1)[0]
 223             contact_point_index=contact_point.index
 224             self.wlccontact_point=contact_point
 225             self.wlccontact_index=contact_point.index
 226             self.wlccurrent=self.current.path
 227         elif 'noauto' in args.split():
 228             if self.wlccontact_index==None or self.wlccurrent != self.current.path:
 229                 print 'Click contact point'
 230                 contact_point=self._measure_N_points(N=1, whatset=1)[0]
 231                 contact_point_index=contact_point.index
 232                 self.wlccontact_point=contact_point
 233                 self.wlccontact_index=contact_point.index
 234                 self.wlccurrent=self.current.path
 235             else:
 236                 contact_point=self.wlccontact_point
 237                 contact_point_index=self.wlccontact_index
 238         else:
 239             cindex=self.find_contact_point()
 240             contact_point=ClickedPoint()
 241             contact_point.absolute_coords=displayed_plot.vectors[1][0][cindex], displayed_plot.vectors[1][1][cindex]
 242             contact_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
 243             contact_point.is_marker=True
 244
 245         print 'Click edges of chunk'
 246         points=self._measure_N_points(N=2, whatset=1)
 247         points=[contact_point]+points
 248         try:
 249             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 )
 250         except:
 251             print 'Fit not possible. Probably wrong interval -did you click two *different* points?'
 252             return
 253
 254         print 'Contour length: ',params[0]*(1.0e+9),' nm'
 255         to_dump='contour '+self.current.path+' '+str(params[0]*(1.0e+9))+' nm'
 256         self.outlet.push(to_dump)
 257         if len(params)==2: #if we did choose 2-value fit
 258             print 'Persistent length: ',params[1]*(1.0e+9),' nm'
 259             to_dump='persistent '+self.current.path+' '+str(params[1]*(1.0e+9))+' nm'
 260             self.outlet.push(to_dump)
 261
 262         if fit_errors:
 263             fit_nm=[i*(10**9) for i in fit_errors]
 264             print 'Standard deviation (contour length)', fit_nm[0]
 265             if len(fit_nm)>1:
 266                 print 'Standard deviation (persistent length)', fit_nm[1]
 267
 268
 269         #add the clicked points in the final PlotObject
 270         clickvector_x, clickvector_y=[], []
 271         for item in points:
 272             clickvector_x.append(item.graph_coords[0])
 273             clickvector_y.append(item.graph_coords[1])
 274
 275         #create a custom PlotObject to gracefully plot the fit along the curves
 276
 277         fitplot=copy.deepcopy(displayed_plot)
 278         fitplot.add_set(xfit,yfit)
 279         fitplot.add_set(clickvector_x,clickvector_y)
 280
 281         if fitplot.styles==[]:
 282             fitplot.styles=[None,None,None,'scatter']
 283         else:
 284             fitplot.styles+=[None,'scatter']
 285
 286         self._send_plot([fitplot])
 287
 288     def find_contact_point(self):
 289         '''
 290         Finds the contact point on the curve.
 291
 292         The current algorithm (thanks to Francesco Musiani, francesco.musiani@unibo.it and Massimo Sandal) is:
 293         - take care of the PicoForce trigger bug - exclude retraction portions with too high standard deviation
 294         - fit the second half of the retraction curve to a line
 295         - if the fit is not almost horizontal, take a smaller chunk and repeat
 296         - otherwise, we have something horizontal
 297         - so take the average of horizontal points and use it as a baseline
 298
 299         Then, start from the rise of the retraction curve and look at the first point below the
 300         baseline.
 301
 302         FIXME: should be moved, probably to generalvclamp.py
 303         '''
 304         outplot=self.subtract_curves(1)
 305         xret=outplot.vectors[1][0]
 306         ydiff=outplot.vectors[1][1]
 307
 308         xext=self.plots[0].vectors[0][0]
 309         yext=self.plots[0].vectors[0][1]
 310         xret2=self.plots[0].vectors[1][0]
 311         yret=self.plots[0].vectors[1][1]
 312
 313         #taking care of the picoforce trigger bug: we exclude portions of the curve that have too much
 314         #standard deviation. yes, a lot of magic is here.
 315         monster=True
 316         monlength=len(xret)-int(len(xret)/20)
 317         finalength=len(xret)
 318         while monster:
 319             monchunk=scipy.array(ydiff[monlength:finalength])
 320             if abs(scipy.stats.std(monchunk)) < 2e-10:
 321                 monster=False
 322             else: #move away from the monster
 323                 monlength-=int(len(xret)/50)
 324                 finalength-=int(len(xret)/50)
 325
 326
 327         #take half of the thing
 328         endlength=int(len(xret)/2)
 329
 330         ok=False
 331
 332         while not ok:
 333             xchunk=yext[endlength:monlength]
 334             ychunk=yext[endlength:monlength]
 335             regr=scipy.stats.linregress(xchunk,ychunk)[0:2]
 336             #we stop if we found an almost-horizontal fit or if we're going too short...
 337             #FIXME: 0.1 and 6 here are "magic numbers" (although reasonable)
 338             if (abs(regr[1]) > 0.1) and ( endlength < len(xret)-int(len(xret)/6) ) :
 339                 endlength+=10
 340             else:
 341                 ok=True
 342
 343
 344         ymean=scipy.mean(ychunk) #baseline
 345
 346         index=0
 347         point = ymean+1
 348
 349         #find the first point below the calculated baseline
 350         while point > ymean:
 351             try:
 352                 point=yret[index]
 353                 index+=1
 354             except IndexError:
 355                 #The algorithm didn't find anything below the baseline! It should NEVER happen
 356                 index=0
 357                 return index
 358
 359         return index
 360
 361
 362
 363     def find_contact_point2(self, debug=False):
 364         '''
 365         TO BE DEVELOPED IN THE FUTURE
 366         Finds the contact point on the curve.
 367
 368         FIXME: should be moved, probably to generalvclamp.py
 369         '''
 370
 371         #raw_plot=self.current.curve.default_plots()[0]
 372         raw_plot=self.plots[0]
 373         '''xext=self.plots[0].vectors[0][0]
 374         yext=self.plots[0].vectors[0][1]
 375         xret2=self.plots[0].vectors[1][0]
 376         yret=self.plots[0].vectors[1][1]
 377         '''
 378         xext=raw_plot.vectors[0][0]
 379         yext=raw_plot.vectors[0][1]
 380         xret2=raw_plot.vectors[1][0]
 381         yret=raw_plot.vectors[1][1]
 382
 383         first_point=[xext[0], yext[0]]
 384         last_point=[xext[-1], yext[-1]]
 385
 386         #regr=scipy.polyfit(first_point, last_point,1)[0:2]
 387         diffx=abs(first_point[0]-last_point[0])
 388         diffy=abs(first_point[1]-last_point[1])
 389
 390         #using polyfit results in numerical errors. good old algebra.
 391         a=diffy/diffx
 392         b=first_point[1]-(a*first_point[0])
 393         baseline=scipy.polyval((a,b), xext)
 394
 395         ysub=[item-basitem for item,basitem in zip(yext,baseline)]
 396
 397         contact=ysub.index(min(ysub))
 398
 399         return xext,ysub,contact
 400
 401         #now, exploit a ClickedPoint instance to calculate index...
 402         dummy=ClickedPoint()
 403         dummy.absolute_coords=(x_intercept,y_intercept)
 404         dummy.find_graph_coords(xret2,yret)
 405
 406         if debug:
 407             return dummy.index, regr, regr_contact
 408         else:
 409             return dummy.index
 410
 411
 412
 413     def x_do_contact(self,args):
 414         '''
 415         DEBUG COMMAND to be activated in the future
 416         '''
 417         xext,ysub,contact=self.find_contact_point2(debug=True)
 418
 419         contact_plot=self.plots[0]
 420         contact_plot.add_set(xext,ysub)
 421         contact_plot.add_set([xext[contact]],[self.plots[0].vectors[0][1][contact]])
 422         #contact_plot.add_set([first_point[0]],[first_point[1]])
 423         #contact_plot.add_set([last_point[0]],[last_point[1]])
 424         contact_plot.styles=[None,None,None,'scatter']
 425         self._send_plot([contact_plot])
 426         return
 427
 428
 429         index,regr,regr_contact=self.find_contact_point2(debug=True)
 430         print regr
 431         print regr_contact
 432         raw_plot=self.current.curve.default_plots()[0]
 433         xret=raw_plot.vectors[0][0]
 434         #nc_line=[(item*regr[0])+regr[1] for item in x_nc]
 435         nc_line=scipy.polyval(regr,xret)
 436         c_line=scipy.polyval(regr_contact,xret)
 437
 438
 439         contact_plot=self.current.curve.default_plots()[0]
 440         contact_plot.add_set(xret, nc_line)
 441         contact_plot.add_set(xret, c_line)
 442         contact_plot.styles=[None,None,None,None]
 443         #contact_plot.styles.append(None)
 444         contact_plot.destination=1
 445         self._send_plot([contact_plot])
 446