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