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