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         params, yfit, xfit = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T)
 236         try:
 237             params, yfit, xfit = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T)
 238         except:
 239             print 'Fit not possible. Probably wrong interval -did you click two *different* points?'
 240             return
 241
 242         print 'Contour length: ',params[0]*(1.0e+9),' nm'
 243         if len(params)==2: #if we did choose 2-value fit
 244             print 'Persistent length: ',params[1]*(1.0e+9),' nm'
 245
 246         #add the clicked points in the final PlotObject
 247         clickvector_x, clickvector_y=[], []
 248         for item in points:
 249             clickvector_x.append(item.graph_coords[0])
 250             clickvector_y.append(item.graph_coords[1])
 251
 252         #create a custom PlotObject to gracefully plot the fit along the curves
 253
 254         fitplot=copy.deepcopy(displayed_plot)
 255         fitplot.add_set(xfit,yfit)
 256         fitplot.add_set(clickvector_x,clickvector_y)
 257
 258         if fitplot.styles==[]:
 259             fitplot.styles=[None,None,None,'scatter']
 260         else:
 261             fitplot.styles+=[None,'scatter']
 262
 263         self._send_plot([fitplot])
 264
 265     def find_contact_point(self):
 266         '''
 267         Finds the contact point on the curve.
 268
 269         The current algorithm (thanks to Francesco Musiani, francesco.musiani@unibo.it and Massimo Sandal) is:
 270         - take care of the PicoForce trigger bug - exclude retraction portions with too high standard deviation
 271         - fit the second half of the retraction curve to a line
 272         - if the fit is not almost horizontal, take a smaller chunk and repeat
 273         - otherwise, we have something horizontal
 274         - so take the average of horizontal points and use it as a baseline
 275
 276         Then, start from the rise of the retraction curve and look at the first point below the
 277         baseline.
 278
 279         FIXME: should be moved, probably to generalvclamp.py
 280         '''
 281         outplot=self.subtract_curves(1)
 282         xret=outplot.vectors[1][0]
 283         ydiff=outplot.vectors[1][1]
 284
 285         xext=self.plots[0].vectors[0][0]
 286         yext=self.plots[0].vectors[0][1]
 287         xret2=self.plots[0].vectors[1][0]
 288         yret=self.plots[0].vectors[1][1]
 289
 290         #taking care of the picoforce trigger bug: we exclude portions of the curve that have too much
 291         #standard deviation. yes, a lot of magic is here.
 292         monster=True
 293         monlength=len(xret)-int(len(xret)/20)
 294         finalength=len(xret)
 295         while monster:
 296             monchunk=scipy.array(ydiff[monlength:finalength])
 297             if abs(scipy.stats.std(monchunk)) < 2e-10:
 298                 monster=False
 299             else: #move away from the monster
 300                 monlength-=int(len(xret)/50)
 301                 finalength-=int(len(xret)/50)
 302
 303
 304         #take half of the thing
 305         endlength=int(len(xret)/2)
 306
 307         ok=False
 308
 309         while not ok:
 310             xchunk=yext[endlength:monlength]
 311             ychunk=yext[endlength:monlength]
 312             regr=scipy.stats.linregress(xchunk,ychunk)[0:2]
 313             #we stop if we found an almost-horizontal fit or if we're going too short...
 314             #FIXME: 0.1 and 6 here are "magic numbers" (although reasonable)
 315             if (abs(regr[1]) > 0.1) and ( endlength < len(xret)-int(len(xret)/6) ) :
 316                 endlength+=10
 317             else:
 318                 ok=True
 319
 320
 321         ymean=scipy.mean(ychunk) #baseline
 322
 323         index=0
 324         point = ymean+1
 325
 326         #find the first point below the calculated baseline
 327         while point > ymean:
 328             try:
 329                 point=yret[index]
 330                 index+=1
 331             except IndexError:
 332                 #The algorithm didn't find anything below the baseline! It should NEVER happen
 333                 index=0
 334                 return index
 335
 336         return index
 337
 338
 339
 340     def find_contact_point2(self, debug=False):
 341         '''
 342         TO BE DEVELOPED IN THE FUTURE
 343         Finds the contact point on the curve.
 344
 345         FIXME: should be moved, probably to generalvclamp.py
 346         '''
 347         outplot=self.subtract_curves(1)
 348         xret=outplot.vectors[1][0]
 349         ydiff=outplot.vectors[1][1]
 350
 351         raw_plot=self.current.curve.default_plots()[0]
 352
 353         '''xext=self.plots[0].vectors[0][0]
 354         yext=self.plots[0].vectors[0][1]
 355         xret2=self.plots[0].vectors[1][0]
 356         yret=self.plots[0].vectors[1][1]
 357         '''
 358         xext=raw_plot.vectors[0][0]
 359         yext=raw_plot.vectors[0][1]
 360         xret2=raw_plot.vectors[1][0]
 361         yret=raw_plot.vectors[1][1]
 362         #taking care of the picoforce trigger bug: we exclude portions of the curve that have too much
 363         #standard deviation. yes, a lot of magic is here.
 364
 365         monlength=len(xext)
 366         #take half of the thing
 367         endlength=int(len(xext)/2)
 368
 369         ok=False
 370         xchunk=xext[endlength:monlength]
 371         ychunk=yext[endlength:monlength]
 372         regr=scipy.polyfit(xchunk,ychunk,1)[0:2]
 373
 374         '''
 375         while not ok:
 376             xchunk=yext[endlength:monlength]
 377             ychunk=yext[endlength:monlength]
 378             print len(xchunk)
 379             #regr=scipy.stats.linregress(xchunk,ychunk)[0:2]
 380             regr=scipy.polyfit(xchunk,ychunk,1)[0:2]
 381             #we stop if we found an almost-horizontal fit or if we're going too short...
 382             #FIXME: 0.1 and 6 here are "magic numbers" (although reasonable)
 383             if (abs(regr[1]) > 0.1) and ( endlength < len(xret)-int(len(xret)/6) ) :
 384                 endlength+=10
 385             else:
 386                 ok=True
 387         '''
 388
 389         '''
 390         ymean=scipy.mean(ychunk) #baseline
 391
 392         index=0
 393         point = ymean+1
 394
 395         #find the first point below the calculated baseline
 396         while point > ymean:
 397             try:
 398                 point=yret[index]
 399                 index+=1
 400             except IndexError:
 401                 #The algorithm didn't find anything below the baseline! It should NEVER happen
 402                 index=0
 403         '''
 404         #fit the contact line
 405         n_contact=100
 406         x_contact_chunk=xret2[0:n_contact]
 407         y_contact_chunk=yret[0:n_contact]
 408
 409         regr_contact=scipy.polyfit(x_contact_chunk,y_contact_chunk,1)[0:2]
 410         x_intercept=(regr_contact[1]-regr[1])/(regr[0]-regr_contact[0])
 411         y_intercept=(regr_contact[0]*x_intercept + regr_contact[1])
 412
 413         #now, exploit a ClickedPoint instance to calculate index...
 414         dummy=ClickedPoint()
 415         dummy.absolute_coords=(x_intercept,y_intercept)
 416         dummy.find_graph_coords(xret2,yret)
 417
 418         if debug:
 419             return dummy.index, regr, regr_contact
 420         else:
 421             return dummy.index
 422
 423
 424     def x_do_contact(self,args):
 425         '''
 426         DEBUG COMMAND to be activated in the future
 427         '''
 428         index,regr,regr_contact=self.find_contact_point2(debug=True)
 429         print regr
 430         print regr_contact
 431         raw_plot=self.current.curve.default_plots()[0]
 432         xret=raw_plot.vectors[0][0]
 433         #nc_line=[(item*regr[0])+regr[1] for item in x_nc]
 434         nc_line=scipy.polyval(regr,xret)
 435         c_line=scipy.polyval(regr_contact,xret)
 436
 437
 438         contact_plot=self.current.curve.default_plots()[0]
 439         contact_plot.add_set(xret, nc_line)
 440         contact_plot.add_set(xret, c_line)
 441         contact_plot.styles=[None,None,None,None]
 442         #contact_plot.styles.append(None)
 443         contact_plot.destination=1
 444         self._send_plot([contact_plot])
 445