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,plot=False):
 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
 305         if not plot:
 306             plot=self.plots[0]
 307
 308         outplot=self.subtract_curves(1)
 309         xret=outplot.vectors[1][0]
 310         ydiff=outplot.vectors[1][1]
 311
 312         xext=plot.vectors[0][0]
 313         yext=plot.vectors[0][1]
 314         xret2=plot.vectors[1][0]
 315         yret=plot.vectors[1][1]
 316
 317         #taking care of the picoforce trigger bug: we exclude portions of the curve that have too much
 318         #standard deviation. yes, a lot of magic is here.
 319         monster=True
 320         monlength=len(xret)-int(len(xret)/20)
 321         finalength=len(xret)
 322         while monster:
 323             monchunk=scipy.array(ydiff[monlength:finalength])
 324             if abs(scipy.stats.std(monchunk)) < 2e-10:
 325                 monster=False
 326             else: #move away from the monster
 327                 monlength-=int(len(xret)/50)
 328                 finalength-=int(len(xret)/50)
 329
 330
 331         #take half of the thing
 332         endlength=int(len(xret)/2)
 333
 334         ok=False
 335
 336         while not ok:
 337             xchunk=yext[endlength:monlength]
 338             ychunk=yext[endlength:monlength]
 339             regr=scipy.stats.linregress(xchunk,ychunk)[0:2]
 340             #we stop if we found an almost-horizontal fit or if we're going too short...
 341             #FIXME: 0.1 and 6 here are "magic numbers" (although reasonable)
 342             if (abs(regr[1]) > 0.1) and ( endlength < len(xret)-int(len(xret)/6) ) :
 343                 endlength+=10
 344             else:
 345                 ok=True
 346
 347
 348         ymean=scipy.mean(ychunk) #baseline
 349
 350         index=0
 351         point = ymean+1
 352
 353         #find the first point below the calculated baseline
 354         while point > ymean:
 355             try:
 356                 point=yret[index]
 357                 index+=1
 358             except IndexError:
 359                 #The algorithm didn't find anything below the baseline! It should NEVER happen
 360                 index=0
 361                 return index
 362
 363         return index
 364
 365
 366
 367     def find_contact_point2(self, debug=False):
 368         '''
 369         TO BE DEVELOPED IN THE FUTURE
 370         Finds the contact point on the curve.
 371
 372         FIXME: should be moved, probably to generalvclamp.py
 373         '''
 374
 375         #raw_plot=self.current.curve.default_plots()[0]
 376         raw_plot=self.plots[0]
 377         '''xext=self.plots[0].vectors[0][0]
 378         yext=self.plots[0].vectors[0][1]
 379         xret2=self.plots[0].vectors[1][0]
 380         yret=self.plots[0].vectors[1][1]
 381         '''
 382         xext=raw_plot.vectors[0][0]
 383         yext=raw_plot.vectors[0][1]
 384         xret2=raw_plot.vectors[1][0]
 385         yret=raw_plot.vectors[1][1]
 386
 387         first_point=[xext[0], yext[0]]
 388         last_point=[xext[-1], yext[-1]]
 389
 390         #regr=scipy.polyfit(first_point, last_point,1)[0:2]
 391         diffx=abs(first_point[0]-last_point[0])
 392         diffy=abs(first_point[1]-last_point[1])
 393
 394         #using polyfit results in numerical errors. good old algebra.
 395         a=diffy/diffx
 396         b=first_point[1]-(a*first_point[0])
 397         baseline=scipy.polyval((a,b), xext)
 398
 399         ysub=[item-basitem for item,basitem in zip(yext,baseline)]
 400
 401         contact=ysub.index(min(ysub))
 402
 403         return xext,ysub,contact
 404
 405         #now, exploit a ClickedPoint instance to calculate index...
 406         dummy=ClickedPoint()
 407         dummy.absolute_coords=(x_intercept,y_intercept)
 408         dummy.find_graph_coords(xret2,yret)
 409
 410         if debug:
 411             return dummy.index, regr, regr_contact
 412         else:
 413             return dummy.index
 414
 415
 416
 417     def x_do_contact(self,args):
 418         '''
 419         DEBUG COMMAND to be activated in the future
 420         '''
 421         xext,ysub,contact=self.find_contact_point2(debug=True)
 422
 423         contact_plot=self.plots[0]
 424         contact_plot.add_set(xext,ysub)
 425         contact_plot.add_set([xext[contact]],[self.plots[0].vectors[0][1][contact]])
 426         #contact_plot.add_set([first_point[0]],[first_point[1]])
 427         #contact_plot.add_set([last_point[0]],[last_point[1]])
 428         contact_plot.styles=[None,None,None,'scatter']
 429         self._send_plot([contact_plot])
 430         return
 431
 432
 433         index,regr,regr_contact=self.find_contact_point2(debug=True)
 434         print regr
 435         print regr_contact
 436         raw_plot=self.current.curve.default_plots()[0]
 437         xret=raw_plot.vectors[0][0]
 438         #nc_line=[(item*regr[0])+regr[1] for item in x_nc]
 439         nc_line=scipy.polyval(regr,xret)
 440         c_line=scipy.polyval(regr_contact,xret)
 441
 442
 443         contact_plot=self.current.curve.default_plots()[0]
 444         contact_plot.add_set(xret, nc_line)
 445         contact_plot.add_set(xret, c_line)
 446         contact_plot.styles=[None,None,None,None]
 447         #contact_plot.styles.append(None)
 448         contact_plot.destination=1
 449         self._send_plot([contact_plot])
 450