hooke/driver/mfp3d.py

   1 # Copyright (C) 2008-2012 A. Seeholzer
   2 #                         Alberto Gomez-Casado <a.gomezcasado@tnw.utwente.nl>
   3 #                         Richard Naud
   4 #                         Rolf Schmidt <rschmidt@alcor.concordia.ca>
   5 #                         W. Trevor King <wking@drexel.edu>
   6 #
   7 # This file is part of Hooke.
   8 #
   9 # Hooke is free software: you can redistribute it and/or modify it under the
  10 # terms of the GNU Lesser General Public License as published by the Free
  11 # Software Foundation, either version 3 of the License, or (at your option) any
  12 # later version.
  13 #
  14 # Hooke is distributed in the hope that it will be useful, but WITHOUT ANY
  15 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  16 # A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
  17 # details.
  18 #
  19 # You should have received a copy of the GNU Lesser General Public License
  20 # along with Hooke.  If not, see <http://www.gnu.org/licenses/>.
  21
  22 """Driver for MFP-3D files.
  23
  24 This driver reads IGOR binary waves.
  25
  26 AUTHORS:
  27 Matlab version: Richard Naud August 2008 (http://lcn.epfl.ch/~naud/)
  28 Python port: A. Seeholzer October 2008
  29 Hooke submission: Rolf Schmidt, Alberto Gomez-Casado 2009
  30 """
  31
  32 import copy
  33 import os.path
  34 import pprint
  35
  36 import numpy
  37
  38 from .. import curve as curve
  39 from ..util.igorbinarywave import loadibw
  40 from . import Driver as Driver
  41
  42
  43 __version__='0.0.0.20100604'
  44
  45
  46 class MFP3DDriver (Driver):
  47     """Handle Asylum Research's MFP3D data format.
  48     """
  49     def __init__(self):
  50         super(MFP3DDriver, self).__init__(name='mfp3d')
  51
  52     def is_me(self, path):
  53         """Look for identifying fields in the IBW note.
  54         """
  55         if os.path.isdir(path):
  56             return False
  57         if not path.endswith('.ibw'):
  58             return False
  59         targets = ['Version:', 'XOPVersion:', 'ForceNote:']
  60         found = [False]*len(targets)
  61         for line in open(path, 'rU'):
  62             for i,ft in enumerate(zip(found, targets)):
  63                 f,t = ft
  64                 if f == False and line.startswith(t):
  65                     found[i] = True
  66         if min(found) == True:
  67             return True
  68         return False
  69
  70     def read(self, path, info=None):
  71         data,bin_info,wave_info = loadibw(path)
  72         blocks,info = self._translate_ibw(data, bin_info, wave_info)
  73         return (blocks, info)
  74
  75     def _translate_ibw(self, data, bin_info, wave_info):
  76         if bin_info['version'] != 5:
  77             raise NotImplementedError('IBW version %d (< 5) not supported'
  78                                       % bin_info['version'])
  79             # We need version 5 for multidimensional arrays.
  80
  81         # Parse the note into a dictionary
  82         note = {}
  83         for line in bin_info['note'].split('\r'):
  84             fields = [x.strip() for x in line.split(':', 1)]
  85             key = fields[0]
  86             if len(fields) == 2:
  87                 value = fields[1]
  88             else:
  89                 value = None
  90             note[key] = value
  91         bin_info['note'] = note
  92
  93         # Ensure a valid MFP3D file version.
  94         if note['VerDate'] not in ['80501.041', '80501.0207']:
  95             raise Exception(note['VerDate'])
  96             raise NotImplementedError(
  97                 '%s file version %s not supported (yet!)\n%s'
  98                 % (self.name, note['VerDate'], pprint.pformat(note)))
  99
 100         # Parse known parameters into standard Hooke format.
 101         info = {
 102             'raw info':{'bin':bin_info,
 103                         'wave':wave_info},
 104             'time':note['Seconds'],
 105             'spring constant (N/m)':float(note['SpringConstant']),
 106             'temperature (K)':self._temperature(note),
 107             }
 108
 109         # Extract data blocks
 110         blocks = []
 111         indexes = [int(i) for i in note['Indexes'].split(',')]
 112         assert indexes[0] == 0, indexes
 113         for i,start in enumerate(indexes[:-1]):
 114             stop = indexes[i+1]
 115             blocks.append(self._scale_block(data[start:stop+1,:], info, i))
 116
 117         return (blocks, info)
 118
 119     def _scale_block(self, data, info, index):
 120         """Convert the block from its native format to a `numpy.float`
 121         array in SI units.
 122         """
 123         # MFP3D's native data dimensions match Hooke's (<point>, <column>) layout.
 124         shape = 3
 125         # raw column indices
 126         columns = info['raw info']['bin']['dimLabels'][1]
 127         # Depending on your MFP3D version:
 128         #   VerDate 80501.0207: ['Raw', 'Defl', 'LVDT', 'Time']
 129         #   VerDate 80501.041:  ['Raw', 'Defl', 'LVDT']
 130         if 'Time' in columns:
 131             n_col = 3
 132         else:
 133             n_col = 2
 134         ret = curve.Data(
 135             shape=(data.shape[0], n_col),
 136             dtype=numpy.float,
 137             info=copy.deepcopy(info)
 138             )
 139
 140         version = info['raw info']['bin']['note']['VerDate']
 141         if version == '80501.041':
 142             name = ['approach', 'retract', 'pause'][index]
 143         elif version == '80501.0207':
 144             name = ['approach', 'pause', 'retract'][index]
 145         else:
 146             raise NotImplementedError()
 147         ret.info['name'] = name
 148         ret.info['raw data'] = data # store the raw data
 149
 150         z_rcol = columns.index('LVDT')
 151         d_rcol = columns.index('Defl')
 152
 153         # scaled column indices
 154         ret.info['columns'] = ['z piezo (m)', 'deflection (m)']
 155         z_scol = ret.info['columns'].index('z piezo (m)')
 156         d_scol = ret.info['columns'].index('deflection (m)')
 157
 158         # Leading '-' because increasing voltage extends the piezo,
 159         # moving the tip towards the surface (positive indentation),
 160         # but it makes more sense to me to have it increase away from
 161         # the surface (positive separation).
 162         ret[:,z_scol] = -data[:,z_rcol].astype(ret.dtype)
 163
 164         # Leading '-' because deflection voltage increases as the tip
 165         # moves away from the surface, but it makes more sense to me
 166         # to have it increase as it moves toward the surface (positive
 167         # tension on the protein chain).
 168         ret[:,d_scol] = -data[:,d_rcol]
 169
 170         if 'Time' in columns:
 171             ret.info['columns'].append('time (s)')
 172             t_rcol = columns.index('Time')
 173             t_scol = ret.info['columns'].index('time (s)')
 174             ret[:,t_scol] = data[:,t_rcol]
 175
 176         return ret
 177
 178     def _temperature(self, note):
 179         # I'm not sure which field we should be using here.  Options are:
 180         #   StartHeadTemp
 181         #   StartScannerTemp
 182         #   StartBioHeaterTemp
 183         #   EndScannerTemp
 184         #   EndHeadTemp
 185         # I imagine the 'Start*Temp' fields were measured at
 186         # 'StartTempSeconds' at the beginning of a series of curves,
 187         # while our particular curve was initiated at 'Seconds'.
 188         #   python -c "from hooke.hooke import Hooke;
 189         #              h=Hooke();
 190         #              h.run_command('load playlist',
 191         #                  {'input':'test/data/vclamp_mfp3d/playlist'});
 192         #              x = [(int(c.info['raw info']['bin']['note']['Seconds'])
 193         #                    - int(c.info['raw info']['bin']['note']['StartTempSeconds']))
 194         #                   for c in h.playlists.current().items()];
 195         #              print 'average', float(sum(x))/len(x);
 196         #              print 'range', min(x), max(x);
 197         #              print x"
 198         # For the Line*Point*.ibw series, the difference increases slowly
 199         #   46, 46, 47, 47, 48, 49, 49, 50, 50, 51, 51, 52, 52, 53, 53, 54,...
 200         # However, for the Image*.ibw series, the difference increase
 201         # is much faster:
 202         #   21, 38, 145, 150, 171, 181
 203         # This makes the 'Start*Temp' fields less and less relevant as
 204         # the experiment continues.  Still, I suppose it's better than
 205         # nothing.
 206         #
 207         # The 'Thermal' fields seem to be related to cantilever calibration.
 208         celsius = unicode(note['StartHeadTemp'], 'latin-1')
 209         if celsius.endswith(u' \u00b0C'):
 210             number = celsius.split(None, 1)[0]
 211             return float(number) + 273.15  # Convert to Kelvin.
 212         else:
 213             raise NotImplementedError(
 214                 'unkown temperature format: %s' % repr(celsius))