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         When you first issue the command, it will ask for the filename. If you are giving the filename
  56         of an existing file, autopeak will resume it and append measurements to it. If you are giving
  57         a new filename, it will create the file and append to it until you close Hooke.
  58
  59
  60         Useful variables (to set with SET command):
  61         ---
  62         temperature= temperature of the system for wlc fit (in K)
  63
  64         auto_slope_span = number of points on which measure the slope, for slope
  65
  66         auto_fit_nm = number of nm to fit before the peak maximum, for WLC (if usepoints false)
  67         auto_fit_points = number of points to fit before the peak maximum, for WLC (if usepoints true)
  68
  69         baseline_clicks = 0: automatic baseline
  70                           1: decide baseline with a single click and length defined in auto_left_baseline
  71                           2: let user click points of baseline
  72         auto_left_baseline = length in nm to use as baseline from the right point (if baseline_clicks=0 , 1)
  73         auto_right_baseline = distance in nm of peak-most baseline point from last peak (if baseline_clicks = 0)
  74         '''
  75
  76         #MACROS.
  77         #FIXME: to move outside function
  78         def fit_interval_nm(start_index,plot,nm,backwards):
  79             '''
  80             Calculates the number of points to fit, given a fit interval in nm
  81             start_index: index of point
  82             plot: plot to use
  83             backwards: if true, finds a point backwards.
  84             '''
  85             x_vect=plot.vectors[1][0]
  86
  87             c=0
  88             i=start_index
  89             start=x_vect[start_index]
  90             maxlen=len(x_vect)
  91             while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
  92                 if i==0 or i==maxlen-1: #we reached boundaries of vector!
  93                     return c
  94
  95                 if backwards:
  96                     i-=1
  97                 else:
  98                     i+=1
  99                 c+=1
 100             return c
 101
 102         def pickup_contact_point():
 103             '''macro to pick up the contact point by clicking'''
 104             contact_point=self._measure_N_points(N=1, whatset=1)[0]
 105             contact_point_index=contact_point.index
 106             self.wlccontact_point=contact_point
 107             self.wlccontact_index=contact_point.index
 108             self.wlccurrent=self.current.path
 109             return contact_point, contact_point_index
 110
 111         def find_current_peaks():
 112             #Find peaks.
 113             defplot=self.current.curve.default_plots()[0]
 114             flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
 115             defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
 116             peak_location,peak_size=self.has_peaks(defplot, self.convfilt_config['mindeviation'])
 117             return peak_location, peak_size
 118
 119
 120
 121
 122         #default fit etc. variables
 123         pl_value=None
 124         T=self.config['temperature']
 125
 126         slope_span=int(self.config['auto_slope_span'])
 127         delta_force=10
 128         rebase=False #if true=we select rebase
 129
 130         #initialize output data vectors
 131         c_lengths=[]
 132         p_lengths=[]
 133         forces=[]
 134         slopes=[]
 135
 136         #pick up plot
 137         displayed_plot=self._get_displayed_plot(0)
 138
 139         #COMMAND LINE PARSING
 140         #--Using points instead of nm interval
 141         if 'usepoints' in args.split():
 142             fit_points=int(self.config['auto_fit_points'])
 143             usepoints=True
 144         else:
 145             fit_points=None
 146             usepoints=False
 147         #--Recalculate baseline
 148         if 'rebase' in args or (self.basecurrent != self.current.path):
 149             rebase=True
 150
 151         #--Custom persistent length / custom temperature
 152         for arg in args.split():
 153             #look for a persistent length argument.
 154             if 'pl=' in arg:
 155                 pl_expression=arg.split('=')
 156                 pl_value=float(pl_expression[1]) #actual value
 157             #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
 158             if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
 159                 t_expression=arg.split('=')
 160                 T=float(t_expression[1])
 161         #--Contact point arguments
 162         if 'reclick' in args.split():
 163             print 'Click contact point'
 164             contact_point, contact_point_index = pickup_contact_point()
 165         elif 'noauto' in args.split():
 166             if self.wlccontact_index==None or self.wlccurrent != self.current.path:
 167                 print 'Click contact point'
 168                 contact_point , contact_point_index = pickup_contact_point()
 169             else:
 170                 contact_point=self.wlccontact_point
 171                 contact_point_index=self.wlccontact_index
 172         else:
 173             #Automatically find contact point
 174             cindex=self.find_contact_point()
 175             contact_point=self._clickize(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], cindex)
 176             '''
 177             contact_point=ClickedPoint()
 178             contact_point.absolute_coords=displayed_plot.vectors[1][0][cindex], displayed_plot.vectors[1][1][cindex]
 179             contact_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
 180             contact_point.is_marker=True
 181             '''
 182         #--END COMMAND LINE PARSING--
 183
 184
 185         peak_location, peak_size = find_current_peaks()
 186
 187         fitplot=copy.deepcopy(displayed_plot)
 188
 189         #Pick up force baseline
 190         if rebase:
 191             clicks=self.config['baseline_clicks']
 192             if clicks==0:
 193                 self.basepoints=[]
 194                 base_index_0=peak_location[-1]+fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
 195                 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_0))
 196                 base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
 197                 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 198             elif clicks>0:
 199                 print 'Select baseline'
 200                 if clicks==1:
 201                     self.basepoints=self._measure_N_points(N=1, whatset=whatset)
 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                 else:
 205                     self.basepoints=self._measure_N_points(N=2, whatset=whatset)
 206
 207             self.basecurrent=self.current.path
 208
 209         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 210         boundaries.sort()
 211         to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
 212         avg=np.mean(to_average)
 213
 214
 215         for peak in peak_location:
 216             #WLC FITTING
 217             #define fit interval
 218             if not usepoints:
 219                 fit_points=fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
 220             peak_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak)
 221             other_fit_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak-fit_points)
 222
 223             #points for the fit
 224             points=[contact_point, peak_point, other_fit_point]
 225
 226             if abs(peak_point.index-other_fit_point.index) < 2:
 227                 continue
 228
 229             params, yfit, xfit = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], pl_value, T)
 230
 231             #save wlc values (nm)
 232             c_lengths.append(params[0]*(1.0e+9))
 233             if len(params)==2: #if we did choose 2-value fit
 234                 p_lengths.append(params[1]*(1.0e+9))
 235             else:
 236                 p_lengths.append(pl_value)
 237             #Add WLC fit lines to plot
 238             fitplot.add_set(xfit,yfit)
 239
 240             if len(fitplot.styles)==0:
 241                 fitplot.styles=[]
 242             else:
 243                 fitplot.styles.append(None)
 244
 245             #Measure forces
 246             delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
 247             y=min(delta_to_measure)
 248             #save force values (pN)
 249             forces.append(abs(y-avg)*(1.0e+12))
 250
 251             #Measure slopes
 252             slope=self.linefit_between(peak-slope_span,peak)[0]
 253             slopes.append(slope)
 254
 255         #add basepoints to fitplot
 256         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]])
 257         fitplot.styles.append('scatter')
 258
 259         #Show wlc fits and peak locations
 260         self._send_plot([fitplot])
 261         self.do_peaks('')
 262
 263
 264         #Ask the user what peaks to ignore from analysis.
 265         print 'Peaks to ignore (0,1...n from contact point,return to take all)'
 266         print 'N to discard measurement'
 267         exclude_raw=raw_input('Input:')
 268         if exclude_raw=='N':
 269             print 'Discarded.'
 270             return
 271         if not exclude_raw=='':
 272             exclude=exclude_raw.split(',')
 273             try:
 274                 exclude=[int(item) for item in exclude]
 275                 for i in exclude:
 276                     c_lengths[i]=None
 277                     p_lengths[i]=None
 278                     forces[i]=None
 279                     slopes[i]=None
 280             except:
 281                  print 'Bad input, taking all...'
 282         #Clean data vectors from ignored peaks
 283         c_lengths=[item for item in c_lengths if item != None]
 284         p_lengths=[item for item in p_lengths if item != None]
 285         forces=[item for item in forces if item != None]
 286         slopes=[item for item in slopes if item != None]
 287         print 'contour (nm)',c_lengths
 288         print 'p (nm)',p_lengths
 289         print 'forces (pN)',forces
 290         print 'slopes (N/m)',slopes
 291
 292         #Save file info
 293         if self.autofile=='':
 294             self.autofile=raw_input('Autopeak filename? (return to ignore) ')
 295             if self.autofile=='':
 296                 print 'Not saved.'
 297                 return
 298
 299         if not os.path.exists(self.autofile):
 300             f=open(self.autofile,'w+')
 301             f.write('Analysis started '+time.asctime()+'\n')
 302             f.write('----------------------------------------\n')
 303             f.write('; Contour length (nm)  ;  Persistence length (nm) ;  Max.Force (pN)  ;  Slope (N/m) \n')
 304             f.close()
 305
 306         print 'Saving...'
 307         f=open(self.autofile,'a+')
 308
 309         f.write(self.current.path+'\n')
 310         for i in range(len(c_lengths)):
 311             f.write(' ; '+str(c_lengths[i])+' ; '+str(p_lengths[i])+' ; '+str(forces[i])+' ; '+str(slopes[i])+'\n')
 312
 313         f.close()
 314         self.do_note('autopeak')
 315