autopeak.py

   1 #!/usr/bin/env python
   2
   3 from libhooke import WX_GOOD, ClickedPoint
   4 import wxversion
   5 wxversion.select(WX_GOOD)
   6 from wx import PostEvent
   7 import numpy as np
   8 import scipy as sp
   9 import copy
  10 import os.path
  11 import time
  12
  13 import warnings
  14 warnings.simplefilter('ignore',np.RankWarning)
  15
  16
  17 class autopeakCommands:
  18
  19     def do_autopeak(self,args):
  20         '''
  21         AUTOPEAK
  22         (autopeak.py)
  23         Automatically performs a number of analyses on the peaks of the given curve.
  24         Currently it automatically:
  25         - fits peaks with WLC function
  26         - measures peak maximum forces with a baseline
  27         - measures slope in proximity of peak maximum
  28         Requires flatten plotmanipulator , fit.py plugin , flatfilts.py plugin with convfilt
  29
  30         Syntax:
  31         autopeak [rebase] [pl=value] [t=value] [noauto] [reclick]
  32
  33         rebase : Re-asks baseline interval
  34
  35         pl=[value] : Use a fixed persistent length for the fit. If pl is not given,
  36                      the fit will be a 2-variable
  37                      fit. DO NOT put spaces between 'pl', '=' and the value.
  38                      The value must be in meters.
  39                      Scientific notation like 0.35e-9 is fine.
  40
  41         t=[value] : Use a user-defined temperature. The value must be in
  42                     kelvins; by default it is 293 K.
  43                     DO NOT put spaces between 't', '=' and the value.
  44
  45         noauto : allows for clicking the contact point by
  46                  hand (otherwise it is automatically estimated) the first time.
  47                  If subsequent measurements are made, the same contact point
  48                  clicked the first time is used
  49
  50         reclick : redefines by hand the contact point, if noauto has been used before
  51                   but the user is unsatisfied of the previously choosen contact point.
  52
  53         usepoints : fit interval by number of points instead than by nanometers
  54
  55         noflatten : does not use the "flatten" plot manipulator
  56
  57         When you first issue the command, it will ask for the filename. If you are giving the filename
  58         of an existing file, autopeak will resume it and append measurements to it. If you are giving
  59         a new filename, it will create the file and append to it until you close Hooke.
  60
  61
  62         Useful variables (to set with SET command):
  63         ---
  64         temperature= temperature of the system for wlc fit (in K)
  65
  66         auto_slope_span = number of points on which measure the slope, for slope
  67
  68         auto_fit_nm = number of nm to fit before the peak maximum, for WLC (if usepoints false)
  69         auto_fit_points = number of points to fit before the peak maximum, for WLC (if usepoints true)
  70
  71         baseline_clicks = -1: no baseline, f=0 at the contact point (whether hand-picked or automatically found)
  72                            0: automatic baseline
  73                            1: decide baseline with a single click and length defined in auto_left_baseline
  74                            2: let user click points of baseline
  75         auto_left_baseline = length in nm to use as baseline from the right point (if baseline_clicks=0 , 1)
  76         auto_right_baseline = distance in nm of peak-most baseline point from last peak (if baseline_clicks = 0)
  77         '''
  78
  79         #MACROS.
  80         #FIXME: to move outside function
  81         def fit_interval_nm(start_index,plot,nm,backwards):
  82             '''
  83             Calculates the number of points to fit, given a fit interval in nm
  84             start_index: index of point
  85             plot: plot to use
  86             backwards: if true, finds a point backwards.
  87             '''
  88             whatset=1 #FIXME: should be decidable
  89             x_vect=plot.vectors[1][0]
  90
  91             c=0
  92             i=start_index
  93             start=x_vect[start_index]
  94             maxlen=len(x_vect)
  95             while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
  96                 if i==0 or i==maxlen-1: #we reached boundaries of vector!
  97                     return c
  98
  99                 if backwards:
 100                     i-=1
 101                 else:
 102                     i+=1
 103                 c+=1
 104             return c
 105
 106         def pickup_contact_point():
 107             '''macro to pick up the contact point by clicking'''
 108             contact_point=self._measure_N_points(N=1, whatset=1)[0]
 109             contact_point_index=contact_point.index
 110             self.wlccontact_point=contact_point
 111             self.wlccontact_index=contact_point.index
 112             self.wlccurrent=self.current.path
 113             return contact_point, contact_point_index
 114
 115         def find_current_peaks(noflatten):
 116             #Find peaks.
 117             defplot=self.current.curve.default_plots()[0]
 118             if not noflatten:
 119                 flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
 120                 defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
 121             peak_location,peak_size=self.has_peaks(defplot, self.convfilt_config['mindeviation'])
 122             return peak_location, peak_size
 123
 124         #default fit etc. variables
 125         pl_value=None
 126         T=self.config['temperature']
 127
 128         slope_span=int(self.config['auto_slope_span'])
 129         delta_force=10
 130         rebase=False #if true=we select rebase
 131         noflatten=False #if true=we avoid flattening
 132
 133         #initialize output data vectors
 134         c_lengths=[]
 135         p_lengths=[]
 136         sigma_c_lengths=[]
 137         sigma_p_lengths=[]
 138         forces=[]
 139         slopes=[]
 140
 141         #pick up plot
 142         displayed_plot=self._get_displayed_plot(0)
 143
 144         #COMMAND LINE PARSING
 145         #--Using points instead of nm interval
 146         if 'usepoints' in args.split():
 147             fit_points=int(self.config['auto_fit_points'])
 148             usepoints=True
 149         else:
 150             fit_points=None
 151             usepoints=False
 152         #--Recalculate baseline
 153         if 'rebase' in args or (self.basecurrent != self.current.path):
 154             rebase=True
 155
 156         if 'noflatten' in args:
 157             noflatten=True
 158
 159         #--Custom persistent length / custom temperature
 160         for arg in args.split():
 161             #look for a persistent length argument.
 162             if 'pl=' in arg:
 163                 pl_expression=arg.split('=')
 164                 pl_value=float(pl_expression[1]) #actual value
 165             #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
 166             if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
 167                 t_expression=arg.split('=')
 168                 T=float(t_expression[1])
 169         #--Contact point arguments
 170         if 'reclick' in args.split():
 171             print 'Click contact point'
 172             contact_point, contact_point_index = pickup_contact_point()
 173         elif 'noauto' in args.split():
 174             if self.wlccontact_index==None or self.wlccurrent != self.current.path:
 175                 print 'Click contact point'
 176                 contact_point , contact_point_index = pickup_contact_point()
 177             else:
 178                 contact_point=self.wlccontact_point
 179                 contact_point_index=self.wlccontact_index
 180         else:
 181             #Automatically find contact point
 182             cindex=self.find_contact_point()
 183             contact_point=self._clickize(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], cindex)
 184         #--END COMMAND LINE PARSING--
 185
 186
 187         peak_location, peak_size = find_current_peaks(noflatten)
 188
 189         if len(peak_location) == 0:
 190             print 'No peaks to fit.'
 191             return
 192
 193         fitplot=copy.deepcopy(displayed_plot)
 194
 195         #Pick up force baseline
 196         if rebase:
 197             clicks=self.config['baseline_clicks']
 198             if clicks==0:
 199                 self.basepoints=[]
 200                 base_index_0=peak_location[-1]+fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
 201                 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_0))
 202                 base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
 203                 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 204             elif clicks>0:
 205                 print 'Select baseline'
 206                 if clicks==1:
 207                     self.basepoints=self._measure_N_points(N=1, whatset=whatset)
 208                     base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
 209                     self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 210                 else:
 211                     self.basepoints=self._measure_N_points(N=2, whatset=whatset)
 212
 213             self.basecurrent=self.current.path
 214
 215         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 216         boundaries.sort()
 217         to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
 218         avg=np.mean(to_average)
 219
 220         clicks=self.config['baseline_clicks']
 221         if clicks==-1:
 222             try:
 223                 avg=displayed_plot.vectors[1][1][contact_point_index]
 224             except:
 225                 avg=displayed_plot.vectors[1][1][cindex]
 226
 227         for peak in peak_location:
 228             #WLC FITTING
 229             #define fit interval
 230             if not usepoints:
 231                 fit_points=fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
 232             peak_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak)
 233             other_fit_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak-fit_points)
 234
 235             #points for the fit
 236             points=[contact_point, peak_point, other_fit_point]
 237
 238             if abs(peak_point.index-other_fit_point.index) < 2:
 239                 continue
 240
 241             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)
 242
 243
 244             #Measure forces
 245             delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
 246             y=min(delta_to_measure)
 247             #save force values (pN)
 248             #Measure slopes
 249             slope=self.linefit_between(peak-slope_span,peak)[0]
 250
 251
 252             #check fitted data and, if right, add peak to the measurement
 253             #FIXME: code duplication
 254             if len(params)==1: #if we did choose 1-value fit
 255                 p_lengths.append(pl_value)
 256                 c_lengths.append(params[0]*(1.0e+9))
 257                 sigma_p_lengths.append(0)
 258                 sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
 259                 forces.append(abs(y-avg)*(1.0e+12))
 260                 slopes.append(slope)
 261                 #Add WLC fit lines to plot
 262                 fitplot.add_set(xfit,yfit)
 263                 if len(fitplot.styles)==0:
 264                     fitplot.styles=[]
 265                     fitplot.colors=[]
 266                 else:
 267                     fitplot.styles.append(None)
 268                     fitplot.colors.append(None)
 269             else: #2-value fit
 270                 p_leng=params[1]*(1.0e+9)
 271                 #check if persistent length makes sense. otherwise, discard peak.
 272                 if p_leng>self.config['auto_min_p'] and p_leng<self.config['auto_max_p']:
 273                     p_lengths.append(p_leng)
 274                     c_lengths.append(params[0]*(1.0e+9))
 275                     sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
 276                     sigma_p_lengths.append(fit_errors[1]*(1.0e+9))
 277                     forces.append(abs(y-avg)*(1.0e+12))
 278                     slopes.append(slope)
 279
 280                     #Add WLC fit lines to plot
 281                     fitplot.add_set(xfit,yfit)
 282                     if len(fitplot.styles)==0:
 283                         fitplot.styles=[]
 284                         fitplot.colors=[]
 285                     else:
 286                         fitplot.styles.append(None)
 287                         fitplot.colors.append(None)
 288                 else:
 289                     pass
 290
 291
 292         #add basepoints to fitplot
 293         fitplot.add_set([self.basepoints[0].graph_coords[0],self.basepoints[1].graph_coords[0]],[self.basepoints[0].graph_coords[1],self.basepoints[1].graph_coords[1]])
 294         fitplot.styles.append('scatter')
 295         fitplot.colors.append(None)
 296
 297         #Show wlc fits and peak locations
 298         self._send_plot([fitplot])
 299         #self.do_peaks('')
 300
 301         print 'Measurements for all peaks detected:'
 302         print 'contour (nm)', c_lengths
 303         print 'sigma contour (nm)',sigma_c_lengths
 304         print 'p (nm)',p_lengths
 305         print 'sigma p (nm)',sigma_p_lengths
 306         print 'forces (pN)',forces
 307         print 'slopes (N/m)',slopes
 308
 309         #Ask the user what peaks to ignore from analysis.
 310         print 'Peaks to ignore (0,1...n from contact point,return to take all)'
 311         print 'N to discard measurement'
 312         exclude_raw=raw_input('Input:')
 313         if exclude_raw=='N':
 314             print 'Discarded.'
 315             return
 316         if not exclude_raw=='':
 317             exclude=exclude_raw.split(',')
 318             try:
 319                 exclude=[int(item) for item in exclude]
 320                 for i in exclude:
 321                     c_lengths[i]=None
 322                     p_lengths[i]=None
 323                     forces[i]=None
 324                     slopes[i]=None
 325                     sigma_c_lengths[i]=None
 326                     sigma_p_lengths[i]=None
 327             except:
 328                  print 'Bad input, taking all...'
 329         #Clean data vectors from ignored peaks
 330         #FIXME:code duplication
 331         c_lengths=[item for item in c_lengths if item != None]
 332         p_lengths=[item for item in p_lengths if item != None]
 333         forces=[item for item in forces if item != None]
 334         slopes=[item for item in slopes if item != None]
 335         sigma_c_lengths=[item for item in sigma_c_lengths if item != None]
 336         sigma_p_lengths=[item for item in sigma_p_lengths if item != None]
 337
 338         print 'Measurements for chosen peaks:'
 339         print 'contour (nm)',c_lengths
 340         print 'sigma contour (nm)',sigma_c_lengths
 341         print 'p (nm)',p_lengths
 342         print 'sigma p (nm)',sigma_p_lengths
 343         print 'forces (pN)',forces
 344         print 'slopes (N/m)',slopes
 345
 346         #Save file info
 347         if self.autofile=='':
 348             self.autofile=raw_input('Autopeak filename? (return to ignore) ')
 349             if self.autofile=='':
 350                 print 'Not saved.'
 351                 return
 352
 353         if not os.path.exists(self.autofile):
 354             f=open(self.autofile,'w+')
 355             f.write('Analysis started '+time.asctime()+'\n')
 356             f.write('----------------------------------------\n')
 357             f.write('; Contour length (nm)  ;  Persistence length (nm) ;  Max.Force (pN)  ;  Slope (N/m) ;  Sigma contour (nm) ; Sigma persistence (nm)\n')
 358             f.close()
 359
 360         print 'Saving...'
 361         f=open(self.autofile,'a+')
 362
 363         f.write(self.current.path+'\n')
 364         for i in range(len(c_lengths)):
 365             f.write(' ; '+str(c_lengths[i])+' ; '+str(p_lengths[i])+' ; '+str(forces[i])+' ; '+str(slopes[i])+' ; '+str(sigma_c_lengths[i])+' ; '+str(sigma_p_lengths[i])+'\n')
 366
 367         f.close()
 368         self.do_note('autopeak')
 369