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 . import Driver as Driver
  41 from .igorbinarywave import loadibw
  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))