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 from igor.binarywave import load as _loadibw
  38
  39 from .. import curve as curve
  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 = _loadibw(path)
  72         blocks,info = self._translate_ibw(data)
  73         return (blocks, info)
  74
  75     def _translate_ibw(self, data):
  76         if data['version'] != 5:
  77             raise NotImplementedError(
  78                 'IBW version {} (< 5) not supported'.format(data['version']))
  79             # We need version 5 for multidimensional arrays.
  80
  81         # Parse the note into a dictionary
  82         note = {}
  83         for line in data['wave']['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         data['wave']['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':data,
 103             'time':note['Seconds'],
 104             'spring constant (N/m)':float(note['SpringConstant']),
 105             'temperature (K)':self._temperature(note),
 106             }
 107
 108         # Extract data blocks
 109         blocks = []
 110         indexes = [int(i) for i in note['Indexes'].split(',')]
 111         assert indexes[0] == 0, indexes
 112         for i,start in enumerate(indexes[:-1]):
 113             stop = indexes[i+1]
 114             blocks.append(
 115                 self._scale_block(
 116                     data['wave']['wData'][start:stop+1,:], info, i))
 117
 118         return (blocks, info)
 119
 120     def _scale_block(self, data, info, index):
 121         """Convert the block from its native format to a `numpy.float`
 122         array in SI units.
 123         """
 124         # MFP3D's native data dimensions match Hooke's (<point>, <column>) layout.
 125         shape = 3
 126         # raw column indices
 127         columns = info['raw info']['wave']['labels'][1][1:]
 128         # Depending on your MFP3D version:
 129         #   VerDate 80501.0207: ['Raw', 'Defl', 'LVDT', 'Time']
 130         #   VerDate 80501.041:  ['Raw', 'Defl', 'LVDT']
 131         if 'Time' in columns:
 132             n_col = 3
 133         else:
 134             n_col = 2
 135         ret = curve.Data(
 136             shape=(data.shape[0], n_col),
 137             dtype=numpy.float,
 138             info=copy.deepcopy(info)
 139             )
 140
 141         version = info['raw info']['wave']['note']['VerDate']
 142         if version == '80501.041':
 143             name = ['approach', 'retract', 'pause'][index]
 144         elif version == '80501.0207':
 145             name = ['approach', 'pause', 'retract'][index]
 146         else:
 147             raise NotImplementedError()
 148         ret.info['name'] = name
 149         ret.info['raw data'] = data # store the raw data
 150
 151         z_rcol = columns.index('LVDT')
 152         d_rcol = columns.index('Defl')
 153
 154         # scaled column indices
 155         ret.info['columns'] = ['z piezo (m)', 'deflection (m)']
 156         z_scol = ret.info['columns'].index('z piezo (m)')
 157         d_scol = ret.info['columns'].index('deflection (m)')
 158
 159         # Leading '-' because increasing voltage extends the piezo,
 160         # moving the tip towards the surface (positive indentation),
 161         # but it makes more sense to me to have it increase away from
 162         # the surface (positive separation).
 163         ret[:,z_scol] = -data[:,z_rcol].astype(ret.dtype)
 164
 165         # Leading '-' because deflection voltage increases as the tip
 166         # moves away from the surface, but it makes more sense to me
 167         # to have it increase as it moves toward the surface (positive
 168         # tension on the protein chain).
 169         ret[:,d_scol] = -data[:,d_rcol]
 170
 171         if 'Time' in columns:
 172             ret.info['columns'].append('time (s)')
 173             t_rcol = columns.index('Time')
 174             t_scol = ret.info['columns'].index('time (s)')
 175             ret[:,t_scol] = data[:,t_rcol]
 176
 177         return ret
 178
 179     def _temperature(self, note):
 180         # I'm not sure which field we should be using here.  Options are:
 181         #   StartHeadTemp
 182         #   StartScannerTemp
 183         #   StartBioHeaterTemp
 184         #   EndScannerTemp
 185         #   EndHeadTemp
 186         # I imagine the 'Start*Temp' fields were measured at
 187         # 'StartTempSeconds' at the beginning of a series of curves,
 188         # while our particular curve was initiated at 'Seconds'.
 189         #   python -c "from hooke.hooke import Hooke;
 190         #              h=Hooke();
 191         #              h.run_command('load playlist',
 192         #                  {'input':'test/data/vclamp_mfp3d/playlist'});
 193         #              x = [(int(c.info['raw info']['wave']['note']['Seconds'])
 194         #                    - int(c.info['raw info']['wave']['note']['StartTempSeconds']))
 195         #                   for c in h.playlists.current().items()];
 196         #              print 'average', float(sum(x))/len(x);
 197         #              print 'range', min(x), max(x);
 198         #              print x"
 199         # For the Line*Point*.ibw series, the difference increases slowly
 200         #   46, 46, 47, 47, 48, 49, 49, 50, 50, 51, 51, 52, 52, 53, 53, 54,...
 201         # However, for the Image*.ibw series, the difference increase
 202         # is much faster:
 203         #   21, 38, 145, 150, 171, 181
 204         # This makes the 'Start*Temp' fields less and less relevant as
 205         # the experiment continues.  Still, I suppose it's better than
 206         # nothing.
 207         #
 208         # The 'Thermal' fields seem to be related to cantilever calibration.
 209         celsius = unicode(note['StartHeadTemp'], 'latin-1')
 210         if celsius.endswith(u' \u00b0C'):
 211             number = celsius.split(None, 1)[0]
 212             return float(number) + 273.15  # Convert to Kelvin.
 213         else:
 214             raise NotImplementedError(
 215                 'unkown temperature format: %s' % repr(celsius))