From c48aef65ea19fd0fb8fb4afaa94f8a4e401c945c Mon Sep 17 00:00:00 2001
From: "fabrizio.benedetti.82" <devnull@localhost>
Date: Wed, 27 Jan 2010 08:02:51 +0000
Subject: [PATCH] Added two plugin: jumpstat.py and multidistance.py. Change to
 the function has_peaks. Moved some basics routines in the new file
 curvetools.py.

---
 autopeak.py      |  68 +---------------
 curvetools.py    |  92 +++++++++++++++++++++
 flatfilts.py     |  28 +++++--
 hooke.conf       |   2 +
 jumpstat.py      | 204 +++++++++++++++++++++++++++++++++++++++++++++++
 multidistance.py |  19 +----
 6 files changed, 324 insertions(+), 89 deletions(-)
 create mode 100755 curvetools.py
 create mode 100755 jumpstat.py

diff --git a/autopeak.py b/autopeak.py
index 8a7baaf..49f2646 100644
--- a/autopeak.py
+++ b/autopeak.py
@@ -82,52 +82,7 @@ class autopeakCommands:
         auto_min_p ; auto_max_p = Minimum and maximum persistence length (if using WLC) or Kuhn length (if using FJC)
                                 outside of which the peak is automatically discarded (in nm)
         '''
-        
-        #MACROS.
-        #FIXME: to move outside function
-        def fit_interval_nm(start_index,plot,nm,backwards):
-            '''
-            Calculates the number of points to fit, given a fit interval in nm
-            start_index: index of point
-            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
-            i=start_index
-            start=x_vect[start_index]
-            maxlen=len(x_vect)
-            while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
-                if i==0 or i==maxlen-1: #we reached boundaries of vector!
-                    return c
-                
-                if backwards:
-                    i-=1
-                else:
-                    i+=1
-                c+=1
-            return c
-                
-        def pickup_contact_point():
-            '''macro to pick up the contact point by clicking'''
-            contact_point=self._measure_N_points(N=1, whatset=1)[0]
-            contact_point_index=contact_point.index
-            self.wlccontact_point=contact_point
-            self.wlccontact_index=contact_point.index
-            self.wlccurrent=self.current.path
-            return contact_point, contact_point_index
-        
-        def find_current_peaks(noflatten):
-            #Find peaks.
-            defplot=self.current.curve.default_plots()[0]
-            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']
@@ -191,7 +146,7 @@ class autopeakCommands:
         #--END COMMAND LINE PARSING--
         
         
-        peak_location, peak_size = find_current_peaks(noflatten)
+        peak_location, peak_size = self.find_current_peaks(noflatten)
         
         if len(peak_location) == 0:
             print 'No peaks to fit.'
@@ -201,23 +156,7 @@ class autopeakCommands:
         
         #Pick up force baseline
         if rebase:
-            clicks=self.config['baseline_clicks']
-            if clicks==0:
-                self.basepoints=[]
-                base_index_0=peak_location[-1]+fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
-                self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_0))
-                base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
-                self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
-            elif clicks>0:
-                print 'Select baseline'
-                if clicks==1:
-                    self.basepoints=self._measure_N_points(N=1, whatset=whatset)
-                    base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
-                    self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
-                else:
-                    self.basepoints=self._measure_N_points(N=2, whatset=whatset)
-            
-            self.basecurrent=self.current.path
+            self.basepoints=self.baseline_points(peak_location, displayed_plot)
         
         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
         boundaries.sort()
@@ -235,7 +174,7 @@ class autopeakCommands:
             #WLC FITTING
             #define fit interval
             if not usepoints:
-                fit_points=fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
+                fit_points=self.fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
             peak_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak)
             other_fit_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak-fit_points)
             
@@ -311,7 +250,6 @@ class autopeakCommands:
         
         #Show wlc fits and peak locations
         self._send_plot([fitplot])
-        #self.do_peaks('')
         
         print 'Using fit function: ',self.config['fit_function']
         print 'Measurements for all peaks detected:'
diff --git a/curvetools.py b/curvetools.py
new file mode 100755
index 0000000..51a8fd4
--- /dev/null
+++ b/curvetools.py
@@ -0,0 +1,92 @@
+# -*- coding: utf-8 -*-
+from libhooke import WX_GOOD, ClickedPoint
+import wxversion
+wxversion.select(WX_GOOD)
+from wx import PostEvent
+import numpy as np
+import scipy as sp
+import copy
+import os.path
+import time
+
+
+class curvetoolsCommands:
+
+      def fit_interval_nm(self,start_index,plot,nm,backwards):
+	  '''
+	  Calculates the number of points to fit, given a fit interval in nm
+	  start_index: index of point
+	  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
+	  i=start_index
+	  start=x_vect[start_index]
+	  maxlen=len(x_vect)
+	  while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
+	      if i==0 or i==maxlen-1: #we reached boundaries of vector!
+		  return c
+	      
+	      if backwards:
+		  i-=1
+	      else:
+		  i+=1
+	      c+=1
+	  return c
+
+
+
+      def find_current_peaks(self,noflatten, a=True, maxpeak=True):
+	    #Find peaks.
+	    if a==True:
+		  a=self.convfilt_config['mindeviation']
+	    try:
+		  abs_devs=float(a)
+	    except:
+		  print "Bad input, using default."
+		  abs_devs=self.convfilt_config['mindeviation']
+
+	    defplot=self.current.curve.default_plots()[0]
+	    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
+	    pk_location,peak_size=self.has_peaks(defplot, abs_devs, maxpeak)
+	    return pk_location, peak_size
+
+
+      def pickup_contact_point(self,N=1,whatset=1):
+	    '''macro to pick up the contact point by clicking'''
+	    contact_point=self._measure_N_points(N=1, whatset=1)[0]
+	    contact_point_index=contact_point.index
+	    self.wlccontact_point=contact_point
+	    self.wlccontact_index=contact_point.index
+	    self.wlccurrent=self.current.path
+	    return contact_point, contact_point_index
+
+
+
+      def baseline_points(self,peak_location, displayed_plot):
+            clicks=self.config['baseline_clicks']
+            if clicks==0:
+                self.basepoints=[]
+                base_index_0=peak_location[-1]+self.fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
+                self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_0))
+                base_index_1=self.basepoints[0].index+self.fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
+                self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
+            elif clicks>0:
+                print 'Select baseline'
+                if clicks==1:
+                    self.basepoints=self._measure_N_points(N=1, whatset=1)
+                    base_index_1=self.basepoints[0].index+self.fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
+                    self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
+                else:
+                    self.basepoints=self._measure_N_points(N=2, whatset=1)
+            
+            self.basecurrent=self.current.path
+            return self.basepoints
+
+
+
diff --git a/flatfilts.py b/flatfilts.py
index 3e9f9c8..0843481 100755
--- a/flatfilts.py
+++ b/flatfilts.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python
+# -*- coding: utf-8 -*-
 
 '''
 FLATFILTS
@@ -147,7 +148,7 @@ class flatfiltsCommands:
     #-----Convolution-based peak recognition and filtering.
     #Requires the libpeakspot.py library
     
-    def has_peaks(self, plot, abs_devs=None):
+    def has_peaks(self, plot, abs_devs=None, maxpeak=True, window=10):
         '''
         Finds peak position in a force curve.
         FIXME: should be moved in libpeakspot.py
@@ -176,12 +177,24 @@ class flatfiltsCommands:
         above=lps.abovenoise(convoluted,noise_level,cut_index,abs_devs)     
         peak_location,peak_size=lps.find_peaks(above,seedouble=self.convfilt_config['seedouble'])
                 
-        #take the maximum
+        #take the minimum or the maximum of a peak
         for i in range(len(peak_location)):
             peak=peak_location[i]
-            maxpk=min(yret[peak-10:peak+10])
-            index_maxpk=yret[peak-10:peak+10].index(maxpk)+(peak-10)
-            peak_location[i]=index_maxpk
+            valpk=min(yret[peak-window:peak+window])  #maximum in force (near the unfolding point)
+            index_pk=yret[peak-window:peak+window].index(valpk)+(peak-window)            
+
+            if maxpeak==False:
+               valpk=max(yret[peak:peak+window]) #minimum in force, near the baseline
+               index_pk=yret[peak:peak+window].index(valpk)+(peak)
+
+#  Let's explain that for the minimum.  Immaging that we know that there is a peak at position/region 100 and you have found its y-value,
+#  Now you look in the array, from 100-10 to 100+10  (if the window is 10).
+#  This "100-10 to 100+10" is substancially a new array with its index. In this array you have 20
+#  elements, so the index of your y-value will be 10.
+#  Now to find the index in the TOTAL array you have to add the "position" of the "region" (that in this case
+#  correspond to 100) and also substract the window size ---> (+100-10)
+
+            peak_location[i]=index_pk
             
         return peak_location,peak_size
     
@@ -319,7 +332,10 @@ class flatfiltsCommands:
                 item.peak_size=peak_size
                 item.curve=None #empty the item object, to further avoid memory leak
                 notflat_list.append(item)
-        
+
+            for i in range(1000):
+		 k=0
+
         #Warn that no flattening had been done.
         if not ('flatten' in self.config['plotmanips']):
             print 'Flatten manipulator was not found. Processing was done without flattening.'
diff --git a/hooke.conf b/hooke.conf
index bda5613..73e47bc 100755
--- a/hooke.conf
+++ b/hooke.conf
@@ -23,6 +23,7 @@ This section defines which plugins have to be loaded by Hooke.
     -->
 <plugins>
     <fit/>
+    <curvetools/>
     <procplots/>
     <flatfilts/>
     <generalclamp/>
@@ -37,6 +38,7 @@ This section defines which plugins have to be loaded by Hooke.
     <pcluster/>
     <generaltccd/>
     <multidistance/>
+    <jumpstat/>
 </plugins>
 
 <!--
diff --git a/jumpstat.py b/jumpstat.py
new file mode 100755
index 0000000..07f6e05
--- /dev/null
+++ b/jumpstat.py
@@ -0,0 +1,204 @@
+# -*- coding: utf-8 -*-
+from libhooke import WX_GOOD, ClickedPoint
+import wxversion
+wxversion.select(WX_GOOD)
+from wx import PostEvent
+import numpy as np
+import scipy as sp
+import copy
+import os.path
+import time
+import sys
+import warnings
+warnings.simplefilter('ignore',np.RankWarning)
+
+
+class jumpstatCommands():
+    
+    def do_jumpstat(self,args):
+     '''
+     JUMPSTAT
+     jumpstat.py
+     Based on the convolution recognition automatically give:
+     - the delta distance between the peaks,
+     - the delta-force from the top of the peaks and subsequent relaxation,
+     - the delta-force from the top of the peaks and the baseline
+     The command allow also to remove the unwanted peaks that can be due to interference.
+     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.
+     You can also define a minimun deviation of the peaks.
+     
+     Syntax:
+     jumpstat [deviation]
+     deviation = number of times the convolution signal is above the noise absolute deviation.
+     '''
+
+
+     #finding the max and the minimum positions for all the peaks
+     noflatten=False
+     #we use if else only to avoid a "bad input" message from find_current_peaks
+     if (len(args)==0):
+             max_peaks_location, peak_size=self.find_current_peaks(noflatten)
+	     min_peaks_location, pks2=self.find_current_peaks(noflatten, True, False)
+     else:
+	     max_peaks_location, peak_size=self.find_current_peaks(noflatten, args)
+	     min_peaks_location, pks2=self.find_current_peaks(noflatten, args, False)
+
+
+     #print "max_peaks_location:  "+str(len(max_peaks_location))
+     #print "min_peaks_location:  "+str(len(min_peaks_location))
+
+     #if no peaks, we have nothing to plot. exit.
+     if len(max_peaks_location)==0:
+            print "No peaks on this curve."
+            return
+
+     if len(max_peaks_location)!=len(min_peaks_location):
+            print "Something went wrong in peaks recognition, number of minima is different from number of maxima. Exiting."
+            return
+
+     #otherwise, we plot the peak locations.
+     xplotted_ret=self.plots[0].vectors[1][0]
+     yplotted_ret=self.plots[0].vectors[1][1]
+     xgood=[xplotted_ret[index] for index in max_peaks_location]
+     ygood=[yplotted_ret[index] for index in max_peaks_location]
+
+     xafter=[xplotted_ret[index] for index in min_peaks_location]
+     yafter=[yplotted_ret[index] for index in min_peaks_location]
+
+     recplot=self._get_displayed_plot()
+     recplot2=self._get_displayed_plot()
+     recplot.vectors.append([xgood,ygood])
+     recplot2.vectors.append([xafter,yafter])
+
+     if recplot.styles==[]:
+         recplot.styles=[None,None,'scatter']
+         recplot.colors=[None,None,None]
+     else:
+         recplot.styles+=['scatter']
+         recplot.colors+=[None]
+
+     if recplot2.styles==[]:
+         recplot2.styles=[None,None,None]
+         recplot2.colors=[None,'1.0',None]
+     else:
+         recplot2.styles+=['scatter']
+         recplot2.colors+=['0.5']
+
+     self._send_plot([recplot])
+     self._send_plot([recplot2])
+
+
+     #finding the baseline
+     self.basepoints=self.baseline_points(max_peaks_location, recplot)    
+     boundaries=[self.basepoints[0].index, self.basepoints[1].index]
+     boundaries.sort()
+     to_average=recplot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
+     avg=np.mean(to_average)
+
+
+     dist=[]
+     jumpforce=[]
+     force=[]
+
+     #we calculate the distance vector
+     for g in range(len(max_peaks_location)-1):
+       dist.append((10**9)*(xplotted_ret[max_peaks_location[g]]-xplotted_ret[max_peaks_location[g+1]]))
+     print "Distance values for the peaks in nm:"
+     print dist
+
+     #the jump-force vector
+     for g in range(len(max_peaks_location)):
+      jumpforce.append((10**12)     *(yplotted_ret[min_peaks_location[g]] -yplotted_ret[max_peaks_location[g]])   )
+     print "Force values for the jumps of the peaks in pN:"
+     print jumpforce
+
+     #the force from baseline vector
+     for g in range(len(max_peaks_location)):
+      force.append((10**12)*(avg-yplotted_ret[max_peaks_location[g]]))
+     print "Force values for the peaks in pN:"
+     print force
+     
+
+
+     #Now ask for the peaks that we don't want
+     print 'Peaks to ignore (0,1...n from contact point,return to take all)'
+     print 'N to discard measurement'
+     exclude_raw=raw_input('Input:')
+     if exclude_raw=='N':
+        print 'Discarded.'
+        return
+     
+     if not exclude_raw=='':
+        exclude=exclude_raw.split(',')
+	#we convert in numbers the input
+        try:
+            exclude=[int(item) for item in exclude]
+        except:
+            print 'Bad input, taking nothing.'
+	    return
+
+#    we remove the peaks that we don't want from the list, we need a counter beacause if we remove
+#    a peaks the other peaks in the list are shifted by one at each step
+        count=0
+        for a in exclude:
+          if (a==0):
+             max_peaks_location=max_peaks_location[1:]
+	     min_peaks_location=min_peaks_location[1:]
+          else:
+             new_a=a-count
+             max_peaks_location=  max_peaks_location[0:new_a]+max_peaks_location[new_a+1:]
+             min_peaks_location=  min_peaks_location[0:new_a]+min_peaks_location[new_a+1:]
+             peak_size=            peak_size[0:new_a]+peak_size[new_a+1:]
+          count+=1
+
+
+     #print "max_peaks_location:  "+str(len(max_peaks_location))
+     #print "min_peaks_location:  "+str(len(min_peaks_location))
+
+
+     dist=[]
+     jumpforce=[]
+     force=[]
+     #we recalculate the distances and the forces after the removing of the unwanted peaks
+     for g in range(len(max_peaks_location)-1):
+         dist.append(xplotted_ret[max_peaks_location[g]]-xplotted_ret[max_peaks_location[g+1]])
+     for g in range(len(max_peaks_location)):
+         jumpforce.append( yplotted_ret[min_peaks_location[g]] - yplotted_ret[max_peaks_location[g]]   )
+     for g in range(len(max_peaks_location)):
+         force.append(avg  -  yplotted_ret[max_peaks_location[g]])
+
+
+
+
+
+        #Save file info
+     if self.autofile=='':
+            self.autofile=raw_input('Jumpstat filename? (return to ignore) ')
+            if self.autofile=='':
+                print 'Not saved.'
+                return
+
+     if not os.path.exists(self.autofile):
+            f=open(self.autofile,'w+')
+            f.write('Analysis started '+time.asctime()+'\n')
+            f.write('----------------------------------------\n')
+            f.write('; Delta Distance length (m);  Jump Force pN;  Standard Force pN\n')
+	    f.write(self.current.path+'\n')
+	    for k in range(len(dist)):
+	       f.write(";")
+	       f.write(str(dist[k])+";"+str(jumpforce[k])+";"+str(force[k])+"\n"   )
+            f.write("\n")
+            f.close()
+            
+     else:
+	    f=open(self.autofile,'a+')
+	    f.write(self.current.path+'\n')
+	    for k in range(len(dist)):
+	      f.write(";")
+	      f.write(str(dist[k])+";"+str(jumpforce[k])+";"+str(force[k])+"\n"   )
+	    f.write("\n")
+	    f.close()
+ 
+     print 'Saving...'
\ No newline at end of file
diff --git a/multidistance.py b/multidistance.py
index c93900f..903ba4f 100644
--- a/multidistance.py
+++ b/multidistance.py
@@ -31,26 +31,9 @@ class multidistanceCommands:
      deviation = number of times the convolution signal is above the noise absolute deviation.
      '''
 
-     def find_current_peaks(noflatten, a):
-            #Find peaks.
-	    if len(a)==0:
-                 a=self.convfilt_config['mindeviation']
-            try:
-                 abs_devs=float(a)
-            except:
-		 print "Bad input, using default."
-                 abs_devs=self.convfilt_config['mindeviation']
-
-            defplot=self.current.curve.default_plots()[0]
-            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
-            pk_loc,peak_size=self.has_peaks(defplot, abs_devs)
-            return pk_loc, peak_size
-
       
      noflatten=False
-     peaks_location, peak_size=find_current_peaks(noflatten, args)
+     peaks_location, peak_size=self.find_current_peaks(noflatten)
      
      #if no peaks, we have nothing to plot. exit.
      if len(peaks_location)==0:
-- 
2.26.2