updated with density limit
[hooke.git] / autopeak.py
index 795189d3aa9d59f8854b2dba9b0b04fbacae4064..ffacd826bc9756c899f47db79a46c7b92bdc0bad 100644 (file)
@@ -52,6 +52,8 @@ class autopeakCommands:
         
         usepoints : fit interval by number of points instead than by nanometers
         
+        noflatten : does not use the "flatten" plot manipulator
+        
         When you first issue the command, it will ask for the filename. If you are giving the filename
         of an existing file, autopeak will resume it and append measurements to it. If you are giving
         a new filename, it will create the file and append to it until you close Hooke.
@@ -66,9 +68,10 @@ class autopeakCommands:
         auto_fit_nm = number of nm to fit before the peak maximum, for WLC (if usepoints false)
         auto_fit_points = number of points to fit before the peak maximum, for WLC (if usepoints true)
         
-        baseline_clicks = 0: automatic baseline
-                          1: decide baseline with a single click and length defined in auto_left_baseline
-                          2: let user click points of baseline
+        baseline_clicks = -1: no baseline, f=0 at the contact point (whether hand-picked or automatically found)
+                           0: automatic baseline
+                           1: decide baseline with a single click and length defined in auto_left_baseline
+                           2: let user click points of baseline
         auto_left_baseline = length in nm to use as baseline from the right point (if baseline_clicks=0 , 1)
         auto_right_baseline = distance in nm of peak-most baseline point from last peak (if baseline_clicks = 0)
         '''
@@ -82,6 +85,7 @@ class autopeakCommands:
             plot: plot to use
             backwards: if true, finds a point backwards.
             '''
+            whatset=1 #FIXME: should be decidable
             x_vect=plot.vectors[1][0]
             
             c=0
@@ -108,17 +112,15 @@ class autopeakCommands:
             self.wlccurrent=self.current.path
             return contact_point, contact_point_index
         
-        def find_current_peaks():
+        def find_current_peaks(noflatten):
             #Find peaks.
             defplot=self.current.curve.default_plots()[0]
-            flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
-            defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
+            if not noflatten:
+                flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
+                defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
             peak_location,peak_size=self.has_peaks(defplot, self.convfilt_config['mindeviation'])
             return peak_location, peak_size
-        
-        
-                
-                
+    
         #default fit etc. variables
         pl_value=None
         T=self.config['temperature']
@@ -126,10 +128,13 @@ class autopeakCommands:
         slope_span=int(self.config['auto_slope_span'])
         delta_force=10
         rebase=False #if true=we select rebase
+        noflatten=False #if true=we avoid flattening
         
         #initialize output data vectors
         c_lengths=[]
         p_lengths=[]
+        sigma_c_lengths=[]
+        sigma_p_lengths=[]
         forces=[]
         slopes=[]
         
@@ -148,6 +153,9 @@ class autopeakCommands:
         if 'rebase' in args or (self.basecurrent != self.current.path):
             rebase=True 
         
+        if 'noflatten' in args:
+            noflatten=True
+        
         #--Custom persistent length / custom temperature
         for arg in args.split():
             #look for a persistent length argument.
@@ -176,7 +184,11 @@ class autopeakCommands:
         #--END COMMAND LINE PARSING--
         
         
-        peak_location, peak_size = find_current_peaks()
+        peak_location, peak_size = find_current_peaks(noflatten)
+        
+        if len(peak_location) == 0:
+            print 'No peaks to fit.'
+            return
         
         fitplot=copy.deepcopy(displayed_plot)
         
@@ -205,6 +217,12 @@ class autopeakCommands:
         to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
         avg=np.mean(to_average)
         
+        clicks=self.config['baseline_clicks']
+        if clicks==-1:
+            try:
+                avg=displayed_plot.vectors[1][1][contact_point_index]
+            except:
+                avg=displayed_plot.vectors[1][1][cindex]
         
         for peak in peak_location:
             #WLC FITTING
@@ -220,7 +238,7 @@ class autopeakCommands:
             if abs(peak_point.index-other_fit_point.index) < 2:
                 continue
             
-            params, yfit, xfit = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], pl_value, T)
+            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)
             
                 
             #Measure forces
@@ -234,22 +252,28 @@ class autopeakCommands:
             #check fitted data and, if right, add peak to the measurement
             #FIXME: code duplication
             if len(params)==1: #if we did choose 1-value fit
-                p_lengths.append(params[1]*(1.0e+9))
+                p_lengths.append(pl_value)
                 c_lengths.append(params[0]*(1.0e+9))
+                sigma_p_lengths.append(0)
+                sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
                 forces.append(abs(y-avg)*(1.0e+12))
                 slopes.append(slope)     
                 #Add WLC fit lines to plot
                 fitplot.add_set(xfit,yfit)
                 if len(fitplot.styles)==0:
                     fitplot.styles=[]
+                    fitplot.colors=[]
                 else:
                     fitplot.styles.append(None)
-            else:
+                    fitplot.colors.append(None)
+            else: #2-value fit
                 p_leng=params[1]*(1.0e+9)
                 #check if persistent length makes sense. otherwise, discard peak.
                 if p_leng>self.config['auto_min_p'] and p_leng<self.config['auto_max_p']:
                     p_lengths.append(p_leng)       
                     c_lengths.append(params[0]*(1.0e+9))
+                    sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
+                    sigma_p_lengths.append(fit_errors[1]*(1.0e+9))
                     forces.append(abs(y-avg)*(1.0e+12))
                     slopes.append(slope)     
                     
@@ -257,8 +281,10 @@ class autopeakCommands:
                     fitplot.add_set(xfit,yfit)
                     if len(fitplot.styles)==0:
                         fitplot.styles=[]
+                        fitplot.colors=[]
                     else:
                         fitplot.styles.append(None)
+                        fitplot.colors.append(None)
                 else:
                     pass
  
@@ -266,12 +292,19 @@ class autopeakCommands:
         #add basepoints to fitplot
         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]]) 
         fitplot.styles.append('scatter')
-        
+        fitplot.colors.append(None)
         
         #Show wlc fits and peak locations
         self._send_plot([fitplot])
         #self.do_peaks('')
         
+        print 'Measurements for all peaks detected:'
+        print 'contour (nm)', c_lengths
+        print 'sigma contour (nm)',sigma_c_lengths
+        print 'p (nm)',p_lengths
+        print 'sigma p (nm)',sigma_p_lengths
+        print 'forces (pN)',forces
+        print 'slopes (N/m)',slopes
         
         #Ask the user what peaks to ignore from analysis.
         print 'Peaks to ignore (0,1...n from contact point,return to take all)'
@@ -289,15 +322,24 @@ class autopeakCommands:
                     p_lengths[i]=None
                     forces[i]=None
                     slopes[i]=None
+                    sigma_c_lengths[i]=None
+                    sigma_p_lengths[i]=None
             except:
                  print 'Bad input, taking all...'
         #Clean data vectors from ignored peaks        
+        #FIXME:code duplication
         c_lengths=[item for item in c_lengths if item != None]
         p_lengths=[item for item in p_lengths if item != None]
         forces=[item for item in forces if item != None]
         slopes=[item for item in slopes if item != None]    
+        sigma_c_lengths=[item for item in sigma_c_lengths if item != None]    
+        sigma_p_lengths=[item for item in sigma_p_lengths if item != None]    
+        
+        print 'Measurements for chosen peaks:'
         print 'contour (nm)',c_lengths
+        print 'sigma contour (nm)',sigma_c_lengths
         print 'p (nm)',p_lengths
+        print 'sigma p (nm)',sigma_p_lengths
         print 'forces (pN)',forces
         print 'slopes (N/m)',slopes
         
@@ -312,7 +354,7 @@ class autopeakCommands:
             f=open(self.autofile,'w+')
             f.write('Analysis started '+time.asctime()+'\n')
             f.write('----------------------------------------\n')
-            f.write('; Contour length (nm)  ;  Persistence length (nm) ;  Max.Force (pN)  ;  Slope (N/m) \n')
+            f.write('; Contour length (nm)  ;  Persistence length (nm) ;  Max.Force (pN)  ;  Slope (N/m) ;  Sigma contour (nm) ; Sigma persistence (nm)\n')
             f.close()
             
         print 'Saving...'
@@ -320,7 +362,7 @@ class autopeakCommands:
         
         f.write(self.current.path+'\n')
         for i in range(len(c_lengths)):
-            f.write(' ; '+str(c_lengths[i])+' ; '+str(p_lengths[i])+' ; '+str(forces[i])+' ; '+str(slopes[i])+'\n')
+            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')
             
         f.close()
         self.do_note('autopeak')