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 experiment as experiment
  41 from . import Driver as Driver
  42 from .igorbinarywave import loadibw
  43
  44
  45 __version__='0.0.0.20100604'
  46
  47
  48 class MFP3DDriver (Driver):
  49     """Handle Asylum Research's MFP3D data format.
  50     """
  51     def __init__(self):
  52         super(MFP3DDriver, self).__init__(name='mfp3d')
  53
  54     def is_me(self, path):
  55         """Look for identifying fields in the IBW note.
  56         """
  57         if os.path.isdir(path):
  58             return False
  59         if not path.endswith('.ibw'):
  60             return False
  61         targets = ['Version:', 'XOPVersion:', 'ForceNote:']
  62         found = [False]*len(targets)
  63         for line in open(path, 'rU'):
  64             for i,ft in enumerate(zip(found, targets)):
  65                 f,t = ft
  66                 if f == False and line.startswith(t):
  67                     found[i] = True
  68         if min(found) == True:
  69             return True
  70         return False
  71
  72     def read(self, path, info=None):
  73         data,bin_info,wave_info = loadibw(path)
  74         blocks,info = self._translate_ibw(data, bin_info, wave_info)
  75         info['filetype'] = self.name
  76         info['experiment'] = experiment.VelocityClamp()
  77         return (blocks, info)
  78
  79     def _translate_ibw(self, data, bin_info, wave_info):
  80         if bin_info['version'] != 5:
  81             raise NotImplementedError('IBW version %d (< 5) not supported'
  82                                       % bin_info['version'])
  83             # We need version 5 for multidimensional arrays.
  84
  85         # Parse the note into a dictionary
  86         note = {}
  87         for line in bin_info['note'].split('\r'):
  88             fields = [x.strip() for x in line.split(':', 1)]
  89             key = fields[0]
  90             if len(fields) == 2:
  91                 value = fields[1]
  92             else:
  93                 value = None
  94             note[key] = value
  95         bin_info['note'] = note
  96
  97         # Ensure a valid MFP3D file version.
  98         if note['VerDate'] not in ['80501.041', '80501.0207']:
  99             raise Exception(note['VerDate'])
 100             raise NotImplementedError(
 101                 '%s file version %s not supported (yet!)\n%s'
 102                 % (self.name, note['VerDate'], pprint.pformat(note)))
 103
 104         # Parse known parameters into standard Hooke format.
 105         info = {
 106             'raw info':{'bin':bin_info,
 107                         'wave':wave_info},
 108             'time':note['Seconds'],
 109             'spring constant (N/m)':float(note['SpringConstant']),
 110             'temperature (K)':self._temperature(note),
 111             }
 112
 113         # Extract data blocks
 114         blocks = []
 115         indexes = [int(i) for i in note['Indexes'].split(',')]
 116         assert indexes[0] == 0, indexes
 117         for i,start in enumerate(indexes[:-1]):
 118             stop = indexes[i+1]
 119             blocks.append(self._scale_block(data[start:stop+1,:], info, i))
 120
 121         return (blocks, info)
 122
 123     def _scale_block(self, data, info, index):
 124         """Convert the block from its native format to a `numpy.float`
 125         array in SI units.
 126         """
 127         # MFP3D's native data dimensions match Hooke's (<point>, <column>) layout.
 128         shape = 3
 129         # raw column indices
 130         columns = info['raw info']['bin']['dimLabels'][1]
 131         # Depending on your MFP3D version:
 132         #   VerDate 80501.0207: ['Raw', 'Defl', 'LVDT', 'Time']
 133         #   VerDate 80501.041:  ['Raw', 'Defl', 'LVDT']
 134         if 'Time' in columns:
 135             n_col = 3
 136         else:
 137             n_col = 2
 138         ret = curve.Data(
 139             shape=(data.shape[0], n_col),
 140             dtype=numpy.float,
 141             info=copy.deepcopy(info)
 142             )
 143
 144         version = info['raw info']['bin']['note']['VerDate']
 145         if version == '80501.041':
 146             name = ['approach', 'retract', 'pause'][index]
 147         elif version == '80501.0207':
 148             name = ['approach', 'pause', 'retract'][index]
 149         else:
 150             raise NotImplementedError()
 151         ret.info['name'] = name
 152         ret.info['raw data'] = data # store the raw data
 153
 154         z_rcol = columns.index('LVDT')
 155         d_rcol = columns.index('Defl')
 156
 157         # scaled column indices
 158         ret.info['columns'] = ['z piezo (m)', 'deflection (m)']
 159         z_scol = ret.info['columns'].index('z piezo (m)')
 160         d_scol = ret.info['columns'].index('deflection (m)')
 161
 162         # Leading '-' because increasing voltage extends the piezo,
 163         # moving the tip towards the surface (positive indentation),
 164         # but it makes more sense to me to have it increase away from
 165         # the surface (positive separation).
 166         ret[:,z_scol] = -data[:,z_rcol].astype(ret.dtype)
 167
 168         # Leading '-' because deflection voltage increases as the tip
 169         # moves away from the surface, but it makes more sense to me
 170         # to have it increase as it moves toward the surface (positive
 171         # tension on the protein chain).
 172         ret[:,d_scol] = -data[:,d_rcol]
 173
 174         if 'Time' in columns:
 175             ret.info['columns'].append('time (s)')
 176             t_rcol = columns.index('Time')
 177             t_scol = ret.info['columns'].index('time (s)')
 178             ret[:,t_scol] = data[:,t_rcol]
 179
 180         return ret
 181
 182     def _temperature(self, note):
 183         # I'm not sure which field we should be using here.  Options are:
 184         #   StartHeadTemp
 185         #   StartScannerTemp
 186         #   StartBioHeaterTemp
 187         #   EndScannerTemp
 188         #   EndHeadTemp
 189         # I imagine the 'Start*Temp' fields were measured at
 190         # 'StartTempSeconds' at the beginning of a series of curves,
 191         # while our particular curve was initiated at 'Seconds'.
 192         #   python -c "from hooke.hooke import Hooke;
 193         #              h=Hooke();
 194         #              h.run_command('load playlist',
 195         #                  {'input':'test/data/vclamp_mfp3d/playlist'});
 196         #              x = [(int(c.info['raw info']['bin']['note']['Seconds'])
 197         #                    - int(c.info['raw info']['bin']['note']['StartTempSeconds']))
 198         #                   for c in h.playlists.current().items()];
 199         #              print 'average', float(sum(x))/len(x);
 200         #              print 'range', min(x), max(x);
 201         #              print x"
 202         # For the Line*Point*.ibw series, the difference increases slowly
 203         #   46, 46, 47, 47, 48, 49, 49, 50, 50, 51, 51, 52, 52, 53, 53, 54,...
 204         # However, for the Image*.ibw series, the difference increase
 205         # is much faster:
 206         #   21, 38, 145, 150, 171, 181
 207         # This makes the 'Start*Temp' fields less and less relevant as
 208         # the experiment continues.  Still, I suppose it's better than
 209         # nothing.
 210         #
 211         # The 'Thermal' fields seem to be related to cantilever calibration.
 212         celsius = unicode(note['StartHeadTemp'], 'latin-1')
 213         if celsius.endswith(u' \u00b0C'):
 214             number = celsius.split(None, 1)[0]
 215             return float(number) + 273.15  # Convert to Kelvin.
 216         else:
 217             raise NotImplementedError(
 218                 'unkown temperature format: %s' % repr(celsius))