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