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