hooke/plugin/jumpstat.py

   1 # Copyright (C) 2010 Fabrizio Benedetti
   2 #                    W. Trevor King <wking@drexel.edu>
   3 #
   4 # This file is part of Hooke.
   5 #
   6 # Hooke is free software: you can redistribute it and/or
   7 # modify it under the terms of the GNU Lesser General Public
   8 # License as published by the Free Software Foundation, either
   9 # version 3 of the License, or (at your option) any later version.
  10 #
  11 # Hooke is distributed in the hope that it will be useful,
  12 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 # GNU Lesser General Public License for more details.
  15 #
  16 # You should have received a copy of the GNU Lesser General Public
  17 # License along with Hooke.  If not, see
  18 # <http://www.gnu.org/licenses/>.
  19
  20 from ..libhooke import WX_GOOD, ClickedPoint
  21 import wxversion
  22 wxversion.select(WX_GOOD)
  23 from wx import PostEvent
  24 import numpy as np
  25 import scipy as sp
  26 import copy
  27 import os.path
  28 import time
  29 import sys
  30 import warnings
  31 warnings.simplefilter('ignore',np.RankWarning)
  32
  33
  34 class jumpstatCommands():
  35
  36     def do_jumpstat(self,args):
  37      '''
  38      JUMPSTAT
  39      jumpstat.py
  40      Based on the convolution recognition automatically give:
  41      - the delta distance between the peaks,
  42      - the delta-force from the top of the peaks and subsequent relaxation,
  43      - the delta-force from the top of the peaks and the baseline
  44      The command allow also to remove the unwanted peaks that can be due to interference.
  45      When you first issue the command, it will ask for the filename. If you are giving the filename
  46      of an existing file, autopeak will resume it and append measurements to it. If you are giving
  47      a new filename, it will create the file and append to it until you close Hooke.
  48      You can also define a minimun deviation of the peaks.
  49
  50      Syntax:
  51      jumpstat [deviation]
  52      deviation = number of times the convolution signal is above the noise absolute deviation.
  53      '''
  54
  55
  56      #finding the max and the minimum positions for all the peaks
  57      noflatten=False
  58      #we use if else only to avoid a "bad input" message from find_current_peaks
  59      if (len(args)==0):
  60              max_peaks_location, peak_size=self.find_current_peaks(noflatten)
  61              min_peaks_location, pks2=self.find_current_peaks(noflatten, True, False)
  62      else:
  63              max_peaks_location, peak_size=self.find_current_peaks(noflatten, args)
  64              min_peaks_location, pks2=self.find_current_peaks(noflatten, args, False)
  65
  66
  67      #print "max_peaks_location:  "+str(len(max_peaks_location))
  68      #print "min_peaks_location:  "+str(len(min_peaks_location))
  69
  70      #if no peaks, we have nothing to plot. exit.
  71      if len(max_peaks_location)==0:
  72             print "No peaks on this curve."
  73             return
  74
  75      if len(max_peaks_location)!=len(min_peaks_location):
  76             print "Something went wrong in peaks recognition, number of minima is different from number of maxima. Exiting."
  77             return
  78
  79      #otherwise, we plot the peak locations.
  80      xplotted_ret=self.plots[0].vectors[1][0]
  81      yplotted_ret=self.plots[0].vectors[1][1]
  82      xgood=[xplotted_ret[index] for index in max_peaks_location]
  83      ygood=[yplotted_ret[index] for index in max_peaks_location]
  84
  85      xafter=[xplotted_ret[index] for index in min_peaks_location]
  86      yafter=[yplotted_ret[index] for index in min_peaks_location]
  87
  88      recplot=self._get_displayed_plot()
  89      recplot2=self._get_displayed_plot()
  90      recplot.vectors.append([xgood,ygood])
  91      recplot2.vectors.append([xafter,yafter])
  92
  93      if recplot.styles==[]:
  94          recplot.styles=[None,None,'scatter']
  95          recplot.colors=[None,None,None]
  96      else:
  97          recplot.styles+=['scatter']
  98          recplot.colors+=[None]
  99
 100      if recplot2.styles==[]:
 101          recplot2.styles=[None,None,None]
 102          recplot2.colors=[None,'1.0',None]
 103      else:
 104          recplot2.styles+=['scatter']
 105          recplot2.colors+=['0.5']
 106
 107      self._send_plot([recplot])
 108      self._send_plot([recplot2])
 109
 110
 111      #finding the baseline
 112      self.basepoints=self.baseline_points(max_peaks_location, recplot)
 113      boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 114      boundaries.sort()
 115      to_average=recplot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
 116      avg=np.mean(to_average)
 117
 118
 119      dist=[]
 120      jumpforce=[]
 121      force=[]
 122
 123      #we calculate the distance vector
 124      for g in range(len(max_peaks_location)-1):
 125        dist.append((10**9)*(xplotted_ret[max_peaks_location[g]]-xplotted_ret[max_peaks_location[g+1]]))
 126      print "Distance values for the peaks in nm:"
 127      print dist
 128
 129      #the jump-force vector
 130      for g in range(len(max_peaks_location)):
 131       jumpforce.append((10**12)     *(yplotted_ret[min_peaks_location[g]] -yplotted_ret[max_peaks_location[g]])   )
 132      print "Force values for the jumps of the peaks in pN:"
 133      print jumpforce
 134
 135      #the force from baseline vector
 136      for g in range(len(max_peaks_location)):
 137       force.append((10**12)*(avg-yplotted_ret[max_peaks_location[g]]))
 138      print "Force values for the peaks in pN:"
 139      print force
 140
 141
 142
 143      #Now ask for the peaks that we don't want
 144      print 'Peaks to ignore (0,1...n from contact point,return to take all)'
 145      print 'N to discard measurement'
 146      exclude_raw=raw_input('Input:')
 147      if exclude_raw=='N':
 148         print 'Discarded.'
 149         return
 150
 151      if not exclude_raw=='':
 152         exclude=exclude_raw.split(',')
 153         #we convert in numbers the input
 154         try:
 155             exclude=[int(item) for item in exclude]
 156         except:
 157             print 'Bad input, taking nothing.'
 158             return
 159
 160 #    we remove the peaks that we don't want from the list, we need a counter beacause if we remove
 161 #    a peaks the other peaks in the list are shifted by one at each step
 162         count=0
 163         for a in exclude:
 164           if (a==0):
 165              max_peaks_location=max_peaks_location[1:]
 166              min_peaks_location=min_peaks_location[1:]
 167           else:
 168              new_a=a-count
 169              max_peaks_location=  max_peaks_location[0:new_a]+max_peaks_location[new_a+1:]
 170              min_peaks_location=  min_peaks_location[0:new_a]+min_peaks_location[new_a+1:]
 171              peak_size=            peak_size[0:new_a]+peak_size[new_a+1:]
 172           count+=1
 173
 174
 175      #print "max_peaks_location:  "+str(len(max_peaks_location))
 176      #print "min_peaks_location:  "+str(len(min_peaks_location))
 177
 178
 179      dist=[]
 180      jumpforce=[]
 181      force=[]
 182      #we recalculate the distances and the forces after the removing of the unwanted peaks
 183      for g in range(len(max_peaks_location)-1):
 184          dist.append(xplotted_ret[max_peaks_location[g]]-xplotted_ret[max_peaks_location[g+1]])
 185      for g in range(len(max_peaks_location)):
 186          jumpforce.append( yplotted_ret[min_peaks_location[g]] - yplotted_ret[max_peaks_location[g]]   )
 187      for g in range(len(max_peaks_location)):
 188          force.append(avg  -  yplotted_ret[max_peaks_location[g]])
 189
 190
 191
 192
 193
 194         #Save file info
 195      if self.autofile=='':
 196             self.autofile=raw_input('Jumpstat filename? (return to ignore) ')
 197             if self.autofile=='':
 198                 print 'Not saved.'
 199                 return
 200
 201      if not os.path.exists(self.autofile):
 202             f=open(self.autofile,'w+')
 203             f.write('Analysis started '+time.asctime()+'\n')
 204             f.write('----------------------------------------\n')
 205             f.write('; Delta Distance length (m);  Jump Force pN;  Standard Force pN\n')
 206             f.write(self.current.path+'\n')
 207             for k in range(len(dist)):
 208                f.write(";")
 209                f.write(str(dist[k])+";"+str(jumpforce[k])+";"+str(force[k])+"\n"   )
 210             f.write("\n")
 211             f.close()
 212
 213      else:
 214             f=open(self.autofile,'a+')
 215             f.write(self.current.path+'\n')
 216             for k in range(len(dist)):
 217               f.write(";")
 218               f.write(str(dist[k])+";"+str(jumpforce[k])+";"+str(force[k])+"\n"   )
 219             f.write("\n")
 220             f.close()
 221
 222      print 'Saving...'