Only calculate relative errors if a previous value exists.
[calibcant.git] / calibcant / analyze.py
index 430ffe9c4e1165d867cebf2193e773bf9ac93002..6a840004a6aee3426c9f5a9af0e8adb1cbe4403f 100644 (file)
@@ -1,22 +1,20 @@
 # calibcant - tools for thermally calibrating AFM cantilevers
 #
-# Copyright (C) 2008-2011 W. Trevor King <wking@drexel.edu>
+# Copyright (C) 2008-2012 W. Trevor King <wking@drexel.edu>
 #
 # This file is part of calibcant.
 #
-# calibcant is free software: you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation, either
-# version 3 of the License, or (at your option) any later version.
+# calibcant is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
 #
-# calibcant is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU Lesser General Public License for more details.
+# calibcant is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+# A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
 #
-# You should have received a copy of the GNU Lesser General Public
-# License along with calibcant.  If not, see
-# <http://www.gnu.org/licenses/>.
+# You should have received a copy of the GNU General Public License along with
+# calibcant.  If not, see <http://www.gnu.org/licenses/>.
 
 """Calculate `k` from arrays of bumps, temperatures, and vibrations.
 
@@ -44,13 +42,14 @@ Which are related by the parameters:
 >>> import os
 >>> import tempfile
 >>> import numpy
->>> from pypiezo.config import pprint_HDF5
->>> from .config import HDF5_CalibrationConfig
+>>> from h5config.storage.hdf5 import HDF5_Storage, pprint_HDF5
+>>> from .config import CalibrationConfig
 
 >>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
 >>> os.close(fd)
 
->>> calibration_config = HDF5_CalibrationConfig(filename, '/calib/config/')
+>>> calibration_config = CalibrationConfig(storage=HDF5_Storage(
+...         filename=filename, group='/calib/config/'))
 >>> bumps = numpy.array((15.9e6, 16.9e6, 16.3e6))
 >>> temperatures = numpy.array((295, 295.2, 294.8))
 >>> vibrations = numpy.array((2.20e-5, 2.22e-5, 2.21e-5))
@@ -79,15 +78,21 @@ photodiode sensitivity (bumps).
 /
   /calib
     /calib/config
-      <HDF5 dataset "num-bumps": shape (), type "|S2">
+      <HDF5 dataset "bump": shape (), type "|S1">
+<BLANKLINE>
+      <HDF5 dataset "num-bumps": shape (), type "<i4">
         10
-      <HDF5 dataset "num-temperatures": shape (), type "|S2">
+      <HDF5 dataset "num-temperatures": shape (), type "<i4">
         10
-      <HDF5 dataset "num-vibrations": shape (), type "|S2">
+      <HDF5 dataset "num-vibrations": shape (), type "<i4">
         20
-      <HDF5 dataset "temperature-sleep": shape (), type "|S1">
+      <HDF5 dataset "temperature": shape (), type "|S1">
+<BLANKLINE>
+      <HDF5 dataset "temperature-sleep": shape (), type "<i4">
         1
-      <HDF5 dataset "vibration-spacing": shape (), type "|S5">
+      <HDF5 dataset "vibration": shape (), type "|S1">
+<BLANKLINE>
+      <HDF5 dataset "vibration-spacing": shape (), type "<f8">
         5e-05
     /calib/processed
       /calib/processed/spring-constant
@@ -140,14 +145,15 @@ except (ImportError, RuntimeError), e:
     _matplotlib = None
     _matplotlib_import_error = e
 
-from pypiezo.config import h5_create_group as _h5_create_group
+from h5config.storage.hdf5 import HDF5_Storage as _HDF5_Storage
+from h5config.storage.hdf5 import h5_create_group as _h5_create_group
 
 from . import LOG as _LOG
-from . import base_config as _base_config
+from . import package_config as _package_config
 from .bump_analyze import bump_analyze as _bump_analyze
 from .bump_analyze import bump_load as _bump_load
 from .bump_analyze import bump_save as _bump_save
-from .config import HDF5_CalibrationConfig as _HDF5_CalibrationConfig
+from .config import CalibrationConfig as _CalibrationConfig
 from .T_analyze import T_analyze as _temperature_analyze
 from .T_analyze import T_load as _temperature_load
 from .T_analyze import T_save as _temperature_save
@@ -215,7 +221,7 @@ def calib_analyze(bumps, temperatures, vibrations):
     _LOG.info('vibration variance (V^2) : %g +/- %g (%g)'
               % (v2_m, v2_s, v2_s/v2_m))
 
-    if _base_config['matplotlib']:
+    if _package_config['matplotlib']:
         calib_plot(bumps, temperatures, vibrations)
 
     return (k, k_s)
@@ -226,7 +232,8 @@ def calib_save(filename, group='/', bumps=None, temperatures=None,
         cwg = _h5_create_group(f, group)
         if calibration_config is not None:
             config_cwg = _h5_create_group(cwg, 'config')
-            calibration_config.save(group=config_cwg)
+            storage = _HDF5_Storage()
+            storage.save(config=calibration_config, group=config_cwg)
         if bumps is not None:
             try:
                 del cwg['raw/photodiode-sensitivity/data']
@@ -296,15 +303,16 @@ def calib_load(filename, group='/'):
         except KeyError:
             pass
         try:
-            k = f[group+'processed/spring-constant/data'][...]
+            k = float(f[group+'processed/spring-constant/data'][...])
         except KeyError:
             pass
         try:
-            k_s = f[group+'processed/spring-constant/standard-deviation'][...]
+            k_s = float(
+                f[group+'processed/spring-constant/standard-deviation'][...])
         except KeyError:
             pass
-    calibration_config = _HDF5_CalibrationConfig(
-        filename=filename, group=group+'config/')
+    calibration_config = _CalibrationConfig(storage=_HDF5_Storage(
+            filename=filename, group=group+'config/'))
     calibration_config.load()
     return (bumps, temperatures, vibrations, calibration_config, k, k_s)
 
@@ -327,7 +335,8 @@ def calib_plot(bumps, temperatures, vibrations):
     vib_axes.plot(vibrations, 'b.-')
     vib_axes.set_ylabel('thermal deflection variance (V^2)')
 
-    figure.show()
+    if hasattr(figure, 'show'):
+        figure.show()
 
 
 def calib_load_all(filename, group='/'):
@@ -337,13 +346,12 @@ def calib_load_all(filename, group='/'):
         filename, group+'calibration/')
     bump_details = []
     for i in range(calibration_config['num-bumps']):
-        (raw_bump,bump_config,z_channel_config,z_axis_config,
-         deflection_channel_config,processed_bump) = _bump_load(
+        (raw_bump,bump_config,z_axis_config,deflection_channel_config,
+         processed_bump) = _bump_load(
             filename=filename, group='%sbump/%d/' % (group, i))
         bump_details.append({
                 'raw_bump': raw_bump,
                 'bump_config': bump_config,
-                'z_channel_config': z_channel_config,
                 'z_axis_config': z_axis_config,
                 'deflection_channel_config': deflection_channel_config,
                 'processed_bump': processed_bump,
@@ -387,28 +395,43 @@ def calib_analyze_all(filename, group='/', maximum_relative_error=1e-5,
     assert group.endswith('/'), group
     bumps,temperatures,vibrations,calibration_config,k,k_s = calib_load(
         filename, group+'calibration/')
+    if bumps is None:
+        bumps = _numpy.zeros(
+            (calibration_config['num-bumps'],), dtype=float)
+    if temperatures is None:
+        temperatures = _numpy.zeros(
+            (calibration_config['num-temperatures'],), dtype=float)
+    if vibrations is None:
+        vibrations = _numpy.zeros(
+            (calibration_config['num-vibrations'],), dtype=float)
     changed_bump = changed_temperature = changed_vibration = False
     for i in range(calibration_config['num-bumps']):
+        _changed_bump = False
         bump_group = '%sbump/%d/' % (group, i)
-        (raw_bump,bump_config,z_channel_config,z_axis_config,
+        (raw_bump,bump_config,z_axis_config,
          deflection_channel_config,processed_bump) = _bump_load(
             filename=filename, group=bump_group)
         sensitivity = _bump_analyze(
             data=raw_bump, bump_config=bump_config,
-            z_channel_config=z_channel_config,
             z_axis_config=z_axis_config,
             deflection_channel_config=deflection_channel_config)
         bumps[i] = sensitivity
-        rel_error = abs(sensitivity - processed_bump)/processed_bump
-        if rel_error > maximum_relative_error:
-            _LOG.warn(("new analysis doesn't match for bump %d: %g -> %g "
-                       "(difference: %g, relative error: %g)")
-                      % (i, processed_bump, sensitivity,
-                         sensitivity-processed_bump, rel_error))
+        if processed_bump is None:
+            _changed_bump = True            
+            _LOG.warn('new analysis for bump %d: %g' % (i, sensitivity))
+        else:
+            rel_error = abs(sensitivity - processed_bump)/processed_bump
+            if rel_error > maximum_relative_error:
+                _changed_bump = True
+                _LOG.warn(("new analysis doesn't match for bump %d: %g -> %g "
+                           "(difference: %g, relative error: %g)")
+                          % (i, processed_bump, sensitivity,
+                             sensitivity-processed_bump, rel_error))
+        if _changed_bump and not dry_run:
             changed_bump = True
-            if not dry_run:
-                _bump_save(filename, bump_group, processed_bump=sensitivity)
+            _bump_save(filename, bump_group, processed_bump=sensitivity)
     for i in range(calibration_config['num-temperatures']):
+        _changed_temperature = False
         temperature_group = '%stemperature/%d/' % (group, i)
         (raw_temperature,temperature_config,processed_temperature
          ) = _temperature_load(
@@ -416,19 +439,25 @@ def calib_analyze_all(filename, group='/', maximum_relative_error=1e-5,
         temperature = _temperature_analyze(
             raw_temperature, temperature_config)
         temperatures[i] = temperature
-        rel_error = abs(temperature - processed_temperature
-                        )/processed_temperature
-        if rel_error > maximum_relative_error:
-            _LOG.warn(("new analysis doesn't match for temperature %d: "
-                       "%g -> %g (difference: %g, relative error: %g)")
-                      % (i, processed_temperature, temperature,
-                         temperature-processed_temperature, rel_error))
+        if processed_temperature is None:
+            _changed_temperature = True            
+            _LOG.warn('new analysis for temperature %d: %g' % (i, temperature))
+        else:
+            rel_error = abs(temperature - processed_temperature
+                            )/processed_temperature
+            if rel_error > maximum_relative_error:
+                _changed_temperature = True
+                _LOG.warn(("new analysis doesn't match for temperature %d: "
+                           "%g -> %g (difference: %g, relative error: %g)")
+                          % (i, processed_temperature, temperature,
+                             temperature-processed_temperature, rel_error))
+        if _changed_temperature and not dry_run:
             changed_temperature = True
-            if not dry_run:
-                _temperature_save(
-                    filename, temperature_group,
-                    processed_T=temperature)
+            _temperature_save(
+                filename, temperature_group,
+                processed_T=temperature)
     for i in range(calibration_config['num-vibrations']):
+        _changed_vibration = False
         vibration_group = '%svibration/%d/' % (group, i)
         (raw_vibration,vibration_config,deflection_channel_config,
          processed_vibration) = _vibration_load(
@@ -437,16 +466,21 @@ def calib_analyze_all(filename, group='/', maximum_relative_error=1e-5,
             deflection=raw_vibration, vibration_config=vibration_config,
             deflection_channel_config=deflection_channel_config)
         vibrations[i] = variance
-        rel_error = abs(variance - processed_vibration)/processed_vibration
-        if rel_error > maximum_relative_error:
-            _LOG.warn(("new analysis doesn't match for vibration %d: %g -> %g "
-                       "(difference: %g, relative error: %g)")
-                      % (i, processed_vibration, variance,
-                         variance-processed_vibration, rel_error))
+        if processed_vibration is None:
+            _changed_vibration = True
+            _LOG.warn('new analysis for vibration %d: %g' % (i, variance))
+        else:
+            rel_error = abs(variance - processed_vibration)/processed_vibration
+            if rel_error > maximum_relative_error:
+                _changed_vibration = True
+                _LOG.warn(("new analysis doesn't match for vibration %d: "
+                           "%g -> %g (difference: %g, relative error: %g)")
+                          % (i, processed_vibration, variance,
+                             variance-processed_vibration, rel_error))
+        if _changed_vibration and not dry_run:
             changed_vibration = True
-            if not dry_run:
-                _vibration_save(
-                    filename, vibration_group, processed_vibration=variance)
+            _vibration_save(
+                filename, vibration_group, processed_vibration=variance)
 
     calib_group = '%scalibration/' % group
 
@@ -459,20 +493,32 @@ def calib_analyze_all(filename, group='/', maximum_relative_error=1e-5,
 
     new_k,new_k_s = calib_analyze(
         bumps=bumps, temperatures=temperatures, vibrations=vibrations)
-    rel_error = abs(new_k-k)/k
-    if rel_error > maximum_relative_error:
-        _LOG.warn(("new analysis doesn't match for k: %g -> %g "
-                   "(difference: %g, relative error: %g)")
-                  % (k, new_k, new_k-k, rel_error))
-        if not dry_run:
-            calib_save(filename, calib_group, k=new_k)
-    rel_error = abs(new_k_s-k_s)/k_s
-    if rel_error > maximum_relative_error:
-        _LOG.warn(("new analysis doesn't match for k_s: %g -> %g "
-                   "(difference: %g, relative error: %g)")
-                  % (k_s, new_k_s, new_k_s-k_s, rel_error))
-        if not dry_run:
-            calib_save(filename, calib_group, k_s=new_k_s)
+    new_calib_k = False
+    if k is None:
+        new_calib_k = True
+        _LOG.warn('new analysis for k: %g' % new_k)
+    else:
+        rel_error = abs(new_k-k)/k
+        if rel_error > maximum_relative_error:
+            new_calib_k = True
+            _LOG.warn(("new analysis doesn't match for k: %g -> %g "
+                       "(difference: %g, relative error: %g)")
+                      % (k, new_k, new_k-k, rel_error)) 
+    if new_calib_k and not dry_run:
+        calib_save(filename, calib_group, k=new_k)
+    new_calib_k_s = False
+    if k_s is None:
+        new_calib_k_s = True
+        _LOG.warn('new analysis for k_s: %g' % new_k_s)
+    else:
+        rel_error = abs(new_k_s-k_s)/k_s
+        if rel_error > maximum_relative_error:
+            new_calib_k_s = True
+            _LOG.warn(("new analysis doesn't match for k_s: %g -> %g "
+                       "(difference: %g, relative error: %g)")
+                      % (k_s, new_k_s, new_k_s-k_s, rel_error))
+    if new_calib_k_s and not dry_run:
+        calib_save(filename, calib_group, k_s=new_k_s)
     return (new_k, new_k_s)
 
 def calib_plot_all(bumps, bump_details, temperatures, temperature_details,
@@ -482,7 +528,6 @@ def calib_plot_all(bumps, bump_details, temperatures, temperature_details,
     for i,bump in enumerate(bump_details):
         sensitivity = _bump_analyze(
             data=bump['raw_bump'], bump_config=bump['bump_config'],
-            z_channel_config=bump['z_channel_config'],
             z_axis_config=bump['z_axis_config'],
             deflection_channel_config=bump['deflection_channel_config'],
             plot=True)