X-Git-Url: http://git.tremily.us/?p=calibcant.git;a=blobdiff_plain;f=calibcant%2Fanalyze.py;h=67d6908418321369350094e39bf225af04d10b29;hp=dcfc7471304426f2b58df0b23e182f9860859e9f;hb=e78ae8cd9fec9fe163b35b371807dc481b84abf2;hpb=277131b4942c353fbecaf30d944e1bd4360511ef diff --git a/calibcant/analyze.py b/calibcant/analyze.py index dcfc747..67d6908 100644 --- a/calibcant/analyze.py +++ b/calibcant/analyze.py @@ -1,26 +1,24 @@ # calibcant - tools for thermally calibrating AFM cantilevers # -# Copyright (C) 2008-2011 W. Trevor King +# Copyright (C) 2008-2012 W. Trevor King # # 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 -# . +# You should have received a copy of the GNU General Public License along with +# calibcant. If not, see . """Calculate `k` from arrays of bumps, temperatures, and vibrations. -Separate the more general `calib_analyze()` from the other `calib_*()` +Separate the more general `analyze()` from the other calibration functions in calibcant. The relevent physical quantities are : @@ -41,22 +39,15 @@ Which are related by the parameters: k_cant Fcant / Zcant ->>> import os ->>> import tempfile >>> import numpy ->>> from h5config.storage.hdf5 import HDF5_Storage, pprint_HDF5 ->>> from .config import CalibrationConfig +>>> from .config import CalibrateConfig ->>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-') ->>> os.close(fd) - ->>> calibration_config = CalibrationConfig(storage=HDF5_Storage( -... filename=filename, group='/calib/config/')) +>>> config = CalibrateConfig() >>> 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)) ->>> k,k_s = calib_analyze(bumps=bumps, temperatures=temperatures, +>>> k,k_s = analyze(bumps=bumps, temperatures=temperatures, ... vibrations=vibrations) >>> (k, k_s) # doctest: +ELLIPSIS (0.0493..., 0.00248...) @@ -72,61 +63,6 @@ photodiode sensitivity (bumps). 0.000553... >>> vibrations.std()/vibrations.mean() # doctest: +ELLIPSIS 0.00369... - ->>> calib_save(filename=filename, group='/calib/', -... bumps=bumps, temperatures=temperatures, vibrations=vibrations, -... calibration_config=calibration_config, k=k, k_s=k_s) ->>> pprint_HDF5(filename) # doctest: +ELLIPSIS, +REPORT_UDIFF -/ - /calib - /calib/config - - - - 10 - - 10 - - 20 - - - - 1 - - - - 5e-05 - /calib/processed - /calib/processed/spring-constant - - 0.0493... - - 0.00248... - - N/m - /calib/raw - /calib/raw/photodiode-sensitivity - - [ 15900000. 16900000. 16300000.] - - V/m - /calib/raw/temperature - - [ 295. 295.2 294.8] - - K - /calib/raw/thermal-vibration-variance - - [ 2.20...e-05 2.220...e-05 2.210...e-05] - - V^2 - ->>> bumps,temperatures,vibrations,calibration_config,k,k_s = calib_load( -... filename=filename, group='/calib/') ->>> (k, k_s) # doctest: +ELLIPSIS -(0.0493..., 0.00248...) - ->>> os.remove(filename) """ import h5py as _h5py @@ -147,24 +83,21 @@ except (ImportError, RuntimeError), e: _matplotlib = None _matplotlib_import_error = e -from h5config.storage.hdf5 import HDF5_Storage as _HDF5_Storage from h5config.storage.hdf5 import h5_create_group as _h5_create_group +from pypiezo.base import get_axis_name as _get_axis_name from . import LOG as _LOG 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 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 -from .vib_analyze import vib_analyze as _vibration_analyze -from .vib_analyze import vib_load as _vibration_load -from .vib_analyze import vib_save as _vibration_save - - -def calib_analyze(bumps, temperatures, vibrations): + +from .bump_analyze import analyze as _bump_analyze +from .bump_analyze import save as _bump_save +from .temperature_analyze import analyze as _temperature_analyze +from .temperature_analyze import save as _temperature_save +from .vibration_analyze import analyze as _vibration_analyze +from .vibration_analyze import save as _vibration_save + + +def analyze(bumps, temperatures, vibrations): """Analyze data from `get_calibration_data()` Inputs (all are arrays of recorded data): @@ -224,101 +157,12 @@ def calib_analyze(bumps, temperatures, vibrations): % (v2_m, v2_s, v2_s/v2_m)) if _package_config['matplotlib']: - calib_plot(bumps, temperatures, vibrations) + plot(bumps, temperatures, vibrations) return (k, k_s) -def calib_save(filename, group='/', bumps=None, temperatures=None, - vibrations=None, calibration_config=None, k=None, k_s=None): - with _h5py.File(filename, 'a') as f: - cwg = _h5_create_group(f, group) - if calibration_config is not None: - config_cwg = _h5_create_group(cwg, 'config') - storage = _HDF5_Storage() - storage.save(config=calibration_config, group=config_cwg) - if bumps is not None: - try: - del cwg['raw/photodiode-sensitivity/data'] - except KeyError: - pass - try: - del cwg['raw/photodiode-sensitivity/units'] - except KeyError: - pass - cwg['raw/photodiode-sensitivity/data'] = bumps - cwg['raw/photodiode-sensitivity/units'] = 'V/m' - if temperatures is not None: - try: - del cwg['raw/temperature/data'] - except KeyError: - pass - try: - del cwg['raw/temperature/units'] - except KeyError: - pass - cwg['raw/temperature/data'] = temperatures - cwg['raw/temperature/units'] = 'K' - if vibrations is not None: - try: - del cwg['raw/thermal-vibration-variance/data'] - except KeyError: - pass - try: - del cwg['raw/thermal-vibration-variance/units'] - except KeyError: - pass - cwg['raw/thermal-vibration-variance/data'] = vibrations - cwg['raw/thermal-vibration-variance/units'] = 'V^2' - if k is not None: - try: - del cwg['processed/spring-constant/data'] - except KeyError: - pass - try: - del cwg['processed/spring-constant/units'] - except KeyError: - pass - cwg['processed/spring-constant/data'] = k - cwg['processed/spring-constant/units'] = 'N/m' - if k_s is not None: - try: - del cwg['processed/spring-constant/standard-deviation'] - except KeyError: - pass - cwg['processed/spring-constant/standard-deviation'] = k_s - -def calib_load(filename, group='/'): - assert group.endswith('/') - bumps = temperatures = vibrations = k = k_s = None - configs = [] - with _h5py.File(filename, 'a') as f: - try: - bumps = f[group+'raw/photodiode-sensitivity/data'][...] - except KeyError: - pass - try: - temperatures = f[group+'raw/temperature/data'][...] - except KeyError: - pass - try: - vibrations = f[group+'raw/thermal-vibration-variance/data'][...] - except KeyError: - pass - try: - k = float(f[group+'processed/spring-constant/data'][...]) - except KeyError: - pass - try: - k_s = float( - f[group+'processed/spring-constant/standard-deviation'][...]) - except KeyError: - pass - calibration_config = _CalibrationConfig(storage=_HDF5_Storage( - filename=filename, group=group+'config/')) - calibration_config.load() - return (bumps, temperatures, vibrations, calibration_config, k, k_s) - -def calib_plot(bumps, temperatures, vibrations): + +def plot(bumps, temperatures, vibrations): if not _matplotlib: raise _matplotlib_import_error figure = _matplotlib_pyplot.figure() @@ -337,182 +181,170 @@ def calib_plot(bumps, temperatures, vibrations): vib_axes.plot(vibrations, 'b.-') vib_axes.set_ylabel('thermal deflection variance (V^2)') - figure.show() - - -def calib_load_all(filename, group='/'): - "Load all data from a `calib()` run." - assert group.endswith('/'), group - bumps,temperatures,vibrations,calibration_config,k,k_s = calib_load( - filename, group+'calibration/') - bump_details = [] - for i in range(calibration_config['num-bumps']): - (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_axis_config': z_axis_config, - 'deflection_channel_config': deflection_channel_config, - 'processed_bump': processed_bump, - }) - temperature_details = [] - for i in range(calibration_config['num-temperatures']): - (raw_temperature,temperature_config,processed_temperature - ) = _temperature_load( - filename=filename, group='%stemperature/%d/' % (group, i)) - temperature_details.append({ - 'raw_temperature': raw_temperature, - 'temperature_config': temperature_config, - 'processed_temperature': processed_temperature, - }) - vibration_details = [] - for i in range(calibration_config['num-vibrations']): - (raw_vibration,vibration_config,deflection_channel_config, - processed_vibration) = _vibration_load( - filename=filename, group='%svibration/%d/' % (group, i)) - vibration_details.append({ - 'raw_vibration': raw_vibration, - 'vibration_config': vibration_config, - 'deflection_channel_config': deflection_channel_config, - 'processed_vibration': processed_vibration, - }) - return { - 'bumps': bumps, - 'bump_details': bump_details, - 'temperatures': temperatures, - 'temperature_details': temperature_details, - 'vibrations': vibrations, - 'vibration_details': vibration_details, - 'calibration_config': calibration_config, - 'k': k, - 'k_s': k_s, - } - -def calib_analyze_all(filename, group='/', maximum_relative_error=1e-5, - dry_run=False): - "(Re)analyze all data from a `calib()` run." - assert group.endswith('/'), group - bumps,temperatures,vibrations,calibration_config,k,k_s = calib_load( - filename, group+'calibration/') - changed_bump = changed_temperature = changed_vibration = False - for i in range(calibration_config['num-bumps']): - bump_group = '%sbump/%d/' % (group, i) - (raw_bump,bump_config,z_channel_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 hasattr(figure, 'show'): + figure.show() + return figure +_plot = plot # alternative name for use inside analyze_all() + + +def analyze_all(config, data, raw_data, maximum_relative_error=1e-5, + filename=None, group=None, plot=False, dry_run=False): + "(Re)analyze (and possibly plot) all data from a `calib()` run." + if not data.get('bump', None): + data['bump'] = _numpy.zeros((config['num-bumps'],), dtype=float) + if not data.get('temperature', None): + data['temperature'] = _numpy.zeros( + (config['num-temperatures'],), dtype=float) + if not data.get('vibrations', None): + data['vibration'] = _numpy.zeros( + (config['num-vibrations'],), dtype=float) + axis_config = config['afm']['piezo'].select_config( + setting_name='axes', + attribute_value=config['afm']['main-axis'], + get_attribute=_get_axis_name) + input_config = config['afm']['piezo'].select_config( + setting_name='inputs', attribute_value='deflection') + bumps_changed = temperatures_changed = vibrations_changed = False + if not isinstance(group, _h5py.Group) and not dry_run: + f = _h5py.File(filename, mode) + group = _h5_create_group(f, group) + else: + f = None + try: + for i,bump in enumerate(raw_data['bump']): + data['bump'][i],changed = check_bump( + index=i, bump=bump, z_axis_config=axis_config, + deflection_channel_config=input_config, plot=plot, + maximum_relative_error=maximum_relative_error) + if changed and not dry_run: + bumps_changed = True + bump_group = _h5_create_group(group, 'bump/{}'.format(i)) + _bump_save(group=bump_group, processed=data['bump'][i]) + for i,temperature in enumerate(raw_data['temperature']): + data['temperature'][i],changed = check_temperature( + index=i, temperature=temperature, + maximum_relative_error=maximum_relative_error) + if changed and not dry_run: + temperatures_changed = True + temperature_group = _h5_create_group( + group, 'temperature/{}'.format(i)) + _temperature_save( + group=temerature_group, processed=data['temperature'][i]) + for i,vibration in enumerate(raw_data['vibration']): + data['vibration'][i],changed = check_vibration( + index=i, vibration=vibration, + deflection_channel_config=input_config, plot=plot, + maximum_relative_error=maximum_relative_error) + if changed and not dry_run: + vibrations_changed = True + vibration_group = _h5_create_group( + group, 'vibration/{}'.format(i)) + _vibration_save( + group=vibration_group, processed=data['vibration']) + k,k_s,changed = check_calibration( + k=data['processed']['spring_constant'], + k_s=data['processed']['spring_constant_deviation'], + bumps=data['bump'], + temperatures=data['temperature'], vibrations=data['vibration'], + maximum_relative_error=maximum_relative_error) + if (changed or bumps_changed or temperatures_changed or + vibrations_changed) and not dry_run: + calibration_group = _h5_create_group(group, 'calibration') + if bumps_changed: + calib_save(group=calibration_group, bump=data['bump']) + if temperatures_changed: + calib_save( + group=calibration_group, temperature=data['temperature']) + if vibrations_changed: + calib_save( + group=calibration_group, vibration=data['vibration']) + if changed: + calib_save(group=calibration_group, k=k, k_s=k_s) + finally: + if f: + f.close() + if plot: + _plot(bumps=data['raw']['bump'], + temperatures=data['raw']['temperature'], + vibrations=data['raw']['vibration']) + return (k, k_s) + +def check_bump(index, bump, maximum_relative_error, **kwargs): + changed = False + sensitivity = _bump_analyze( + config=bump['config']['bump'], data=bump['raw'], **kwargs) + if bump.get('processed', None) is None: + changed = True + _LOG.warn('new analysis for bump {}: {}'.format(index, sensitivity)) + else: + rel_error = abs(sensitivity - bump['processed'])/bump['processed'] 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)) - changed_bump = True - if not dry_run: - _bump_save(filename, bump_group, processed_bump=sensitivity) - for i in range(calibration_config['num-temperatures']): - temperature_group = '%stemperature/%d/' % (group, i) - (raw_temperature,temperature_config,processed_temperature - ) = _temperature_load( - filename=filename, group=temperature_group) - temperature = _temperature_analyze( - raw_temperature, temperature_config) - temperatures[i] = temperature - rel_error = abs(temperature - processed_temperature - )/processed_temperature + changed = True + _LOG.warn(("new analysis doesn't match for bump {}: {} -> {} " + "(difference: {}, relative error: {})").format( + index, bump['processed'], sensitivity, + sensitivity-bump['processed'], rel_error)) + return (sensitivity, changed) + +def check_temperature(index, temperature, maximum_relative_error, **kwargs): + changed = False + temp = _temperature_analyze( + config=temperature['config']['temperature'], + temperature=temperature['raw'], **kwargs) + if temperature.get('processed', None) is None: + changed = True + _LOG.warn('new analysis for temperature {}: {}'.format(index, temp)) + else: + rel_error = abs(temp - temperature['processed'] + )/temperature['processed'] 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)) - changed_temperature = True - if not dry_run: - _temperature_save( - filename, temperature_group, - processed_T=temperature) - for i in range(calibration_config['num-vibrations']): - vibration_group = '%svibration/%d/' % (group, i) - (raw_vibration,vibration_config,deflection_channel_config, - processed_vibration) = _vibration_load( - filename=filename, group=vibration_group) - variance = _vibration_analyze( - deflection=raw_vibration, vibration_config=vibration_config, - deflection_channel_config=deflection_channel_config) - vibrations[i] = variance - rel_error = abs(variance - processed_vibration)/processed_vibration + changed = True + _LOG.warn(("new analysis doesn't match for temperature " + "{} -> {} (difference: {}, relative error: {})" + ).format( + index, temperature['processed'], temp, + temp-temperature['processed'], rel_error)) + return (temp, changed) + +def check_vibration(index, vibration, maximum_relative_error, **kwargs): + changed = False + variance = _vibration_analyze( + config=vibration['config']['vibration'], + deflection=vibration['raw'], **kwargs) + if vibration.get('processed', None) is None: + changed = True + _LOG.warn('new analysis for temperature {}: {}'.format( + index, variance)) + else: + rel_error = abs(variance-vibration['processed'])/vibration['processed'] 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)) - changed_vibration = True - if not dry_run: - _vibration_save( - filename, vibration_group, processed_vibration=variance) - - calib_group = '%scalibration/' % group - - if changed_bump and not dry_run: - calib_save(filename, calib_group, bumps=bumps) - if changed_temperature and not dry_run: - calib_save(filename, calib_group, temperatures=temperatures) - if changed_vibration and not dry_run: - calib_save(filename, calib_group, vibrations=vibrations) - - 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) - return (new_k, new_k_s) - -def calib_plot_all(bumps, bump_details, temperatures, temperature_details, - vibrations, vibration_details, calibration_config, k, k_s, - maximum_relative_error=1e-5): - calib_plot(bumps, temperatures, vibrations) - 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) - rel_error = abs(sensitivity - bump['processed_bump'] - )/bump['processed_bump'] + _LOG.warn(("new analysis doesn't match for vibration {}: {} != {} " + "(difference: {}, relative error: {})").format( + index, variance, vibration['processed'], + variance-vibration['processed'], rel_error)) + return (variance, changed) + +def check_calibration(k, k_s, maximum_relative_error, **kwargs): + changed = False + new_k,new_k_s = analyze(**kwargs) + if k is None: + changed = True + _LOG.warn('new analysis for the spring constant: {}'.format(new_k)) + else: + rel_error = abs(new_k-k)/k if rel_error > maximum_relative_error: - _LOG.warn(("new analysis doesn't match for bump %d: %g != %g " - "(difference: %g, relative error: %g)") - % (i, sensitivity, bump['processed_bump'], - sensitivity-bump['processed_bump'], rel_error)) - # nothing interesting to plot for temperatures... - for i,vibration in enumerate(vibration_details): - variance = _vibration_analyze( - deflection=vibration['raw_vibration'], - vibration_config=vibration['vibration_config'], - deflection_channel_config=vibration['deflection_channel_config'], - plot=True) - rel_error = abs(variance - vibration['processed_vibration'] - )/vibration['processed_vibration'] + _LOG.warn(("new analysis doesn't match for the spring constant: " + "{} != {} (difference: {}, relative error: {})").format( + new_k, k, new_k-k, rel_error)) + if k_s is None: + changed = True + _LOG.warn('new analysis for the spring constant deviation: {}'.format( + new_k_s)) + else: + rel_error = abs(new_k-k)/k if rel_error > maximum_relative_error: - _LOG.warn(("new analysis doesn't match for vibration %d: %g != %g " - "(difference: %g, relative error: %g)") - % (i, variance, vibration['processed_vibration'], - variance-vibration['processed_vibration'], rel_error)) + _LOG.warn( + ("new analysis doesn't match for the spring constant deviation" + ": {} != {} (difference: {}, relative error: {})").format( + new_k_s, k_s, new_k_s-k_s, rel_error)) + return (new_k, new_k_s, changed)