hooke/driver/mfp3d.py

   1 # Copyright (C) 2008-2010 A. Seeholzer
   2 #                         Alberto Gomez-Casado
   3 #                         Richard Naud <richard.naud@epfl.ch>
   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
  10 # under the terms of the GNU Lesser General Public License as
  11 # published by the Free Software Foundation, either version 3 of the
  12 # License, or (at your option) any later version.
  13 #
  14 # Hooke is distributed in the hope that it will be useful, but WITHOUT
  15 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  16 # or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
  17 # Public License for more details.
  18 #
  19 # You should have received a copy of the GNU Lesser General Public
  20 # License along with Hooke.  If not, see
  21 # <http://www.gnu.org/licenses/>.
  22
  23 """Driver for MFP-3D files.
  24
  25 This driver reads IGOR binary waves.
  26
  27 AUTHORS:
  28 Matlab version: Richard Naud August 2008 (http://lcn.epfl.ch/~naud/)
  29 Python port: A. Seeholzer October 2008
  30 Hooke submission: Rolf Schmidt, Alberto Gomez-Casado 2009
  31 """
  32
  33 import copy
  34 import os.path
  35 import pprint
  36
  37 import numpy
  38
  39 from .. import curve as curve
  40 from ..util.igorbinarywave import loadibw
  41 from . import Driver as Driver
  42
  43
  44 __version__='0.0.0.20100604'
  45
  46
  47 class MFP3DDriver (Driver):
  48     """Handle Asylum Research's MFP3D data format.
  49     """
  50     def __init__(self):
  51         super(MFP3DDriver, self).__init__(name='mfp3d')
  52
  53     def is_me(self, path):
  54         """Look for identifying fields in the IBW note.
  55         """
  56         if os.path.isdir(path):
  57             return False
  58         if not path.endswith('.ibw'):
  59             return False
  60         targets = ['Version:', 'XOPVersion:', 'ForceNote:']
  61         found = [False]*len(targets)
  62         for line in open(path, 'rU'):
  63             for i,ft in enumerate(zip(found, targets)):
  64                 f,t = ft
  65                 if f == False and line.startswith(t):
  66                     found[i] = True
  67         if min(found) == True:
  68             return True
  69         return False
  70
  71     def read(self, path, info=None):
  72         data,bin_info,wave_info = loadibw(path)
  73         blocks,info = self._translate_ibw(data, bin_info, wave_info)
  74         return (blocks, info)
  75
  76     def _translate_ibw(self, data, bin_info, wave_info):
  77         if bin_info['version'] != 5:
  78             raise NotImplementedError('IBW version %d (< 5) not supported'
  79                                       % bin_info['version'])
  80             # We need version 5 for multidimensional arrays.
  81
  82         # Parse the note into a dictionary
  83         note = {}
  84         for line in bin_info['note'].split('\r'):
  85             fields = [x.strip() for x in line.split(':', 1)]
  86             key = fields[0]
  87             if len(fields) == 2:
  88                 value = fields[1]
  89             else:
  90                 value = None
  91             note[key] = value
  92         bin_info['note'] = note
  93
  94         # Ensure a valid MFP3D file version.
  95         if note['VerDate'] not in ['80501.041', '80501.0207']:
  96             raise Exception(note['VerDate'])
  97             raise NotImplementedError(
  98                 '%s file version %s not supported (yet!)\n%s'
  99                 % (self.name, note['VerDate'], pprint.pformat(note)))
 100
 101         # Parse known parameters into standard Hooke format.
 102         info = {
 103             'raw info':{'bin':bin_info,
 104                         'wave':wave_info},
 105             'time':note['Seconds'],
 106             'spring constant (N/m)':float(note['SpringConstant']),
 107             'temperature (K)':self._temperature(note),
 108             }
 109
 110         # Extract data blocks
 111         blocks = []
 112         indexes = [int(i) for i in note['Indexes'].split(',')]
 113         assert indexes[0] == 0, indexes
 114         for i,start in enumerate(indexes[:-1]):
 115             stop = indexes[i+1]
 116             blocks.append(self._scale_block(data[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']['bin']['dimLabels'][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']['bin']['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']['bin']['note']['Seconds'])
 194         #                    - int(c.info['raw info']['bin']['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))