X-Git-Url: http://git.tremily.us/?p=calibcant.git;a=blobdiff_plain;f=calibcant%2Fanalyze.py;h=1033dcc3b692865f9b87f56d7018fad882dc7ad5;hp=7ed75645f0b3d0b1028a9b4dfbf24f574ba1938d;hb=904f98c200a40abb6f86304c3ae912a01b9bac11;hpb=d102f5faf7064c208cea650249f4e28322d7f0b0 diff --git a/calibcant/analyze.py b/calibcant/analyze.py old mode 100755 new mode 100644 index 7ed7564..1033dcc --- a/calibcant/analyze.py +++ b/calibcant/analyze.py @@ -1,293 +1,524 @@ -#!/usr/bin/python -# # calibcant - tools for thermally calibrating AFM cantilevers # -# Copyright (C) 2007,2008, William Trevor King +# Copyright (C) 2008-2012 W. Trevor King # -# This program 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. +# This file is part of calibcant. # -# This program 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. +# 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. # -# You should have received a copy of the GNU General Public License -# along with this program; if not, write to the Free Software -# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA -# 02111-1307, USA. +# 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. # -# The author may be contacted at on the Internet, or -# write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St., -# Philadelphia PA 19104, USA. +# You should have received a copy of the GNU General Public License along with +# calibcant. If not, see . -""" -Separate the more general calib_analyze() from the other calib_*() -functions in calibcant. Also provide a command line interface -for analyzing data acquired through other workflows. - -The relevent physical quantities are : - Vzp_out Output z-piezo voltage (what we generate) - Vzp Applied z-piezo voltage (after external ZPGAIN) - Zp The z-piezo position - Zcant The cantilever vertical deflection - Vphoto The photodiode vertical deflection voltage (what we measure) - Fcant The force on the cantilever - T The temperature of the cantilever and surrounding solution - (another thing we measure) - k_b Boltzmann's constant - -Which are related by the parameters : - zpGain Vzp_out / Vzp - zpSensitivity Zp / Vzp - photoSensitivity Vphoto / Zcant - k_cant Fcant / Zcant +"""Calculate `k` from arrays of bumps, temperatures, and vibrations. + +Separate the more general `calib_analyze()` from the other `calib_*()` +functions in calibcant. + +The relevent physical quantities are : + Vzp_out Output z-piezo voltage (what we generate) + Vzp Applied z-piezo voltage (after external ZPGAIN) + Zp The z-piezo position + Zcant The cantilever vertical deflection + Vphoto The photodiode vertical deflection voltage (what we measure) + Fcant The force on the cantilever + T The temperature of the cantilever and surrounding solution + (another thing we measure) + k_b Boltzmann's constant + +Which are related by the parameters: + zp_gain Vzp_out / Vzp + zp_sensitivity Zp / Vzp + photo_sensitivity Vphoto / Zcant + k_cant Fcant / Zcant + + +>>> import os +>>> import tempfile +>>> import numpy +>>> 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 = 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)) + +>>> k,k_s = calib_analyze(bumps=bumps, temperatures=temperatures, +... vibrations=vibrations) +>>> (k, k_s) # doctest: +ELLIPSIS +(0.0493..., 0.00248...) + +Most of the error in this example comes from uncertainty in the +photodiode sensitivity (bumps). + +>>> k_s/k # doctest: +ELLIPSIS +0.0503... +>>> bumps.std()/bumps.mean() # doctest: +ELLIPSIS +0.0251... +>>> temperatures.std()/temperatures.mean() # doctest: +ELLIPSIS +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 numpy -from splittable_kwargs import splittableKwargsFunction, \ - make_splittable_kwargs_function -import data_logger +import h5py as _h5py +import numpy as _numpy +try: + from scipy.constants import Boltzmann as _kB # in J/K +except ImportError: + from scipy.constants import Bolzmann as _kB # in J/K +# Bolzmann -> Boltzmann patch submitted: +# http://projects.scipy.org/scipy/ticket/1417 +# Fixed in scipy commit 4716d91, Apr 2, 2011, during work after v0.9.0rc5. -import common # common module for the calibcant package -import config # config module for the calibcant package -import T_analyze # T_analyze module for the calibcant package +try: + import matplotlib as _matplotlib + import matplotlib.pyplot as _matplotlib_pyplot + import time as _time # for timestamping lines on plots +except (ImportError, RuntimeError), e: + _matplotlib = None + _matplotlib_import_error = e -kb = 1.3806504e-23 # Boltzmann's constant in J/K +from h5config.storage.hdf5 import HDF5_Storage as _HDF5_Storage +from h5config.storage.hdf5 import h5_create_group as _h5_create_group -#@splittableKwargsFunction((calib_plot, 'bumps', 'Ts', 'vibs')) -# Some of the child functions aren't yet defined, so postpone -# make-splittable until later in the module. -def calib_analyze(bumps, Ts, vibs, **kwargs) : - """ - Analyze data from get_calibration_data() - return (k, k_s, !!!a2_r, T_r, one_o_Vp2_r) - Inputs (all are arrays of recorded data) : - bumps measured (V_photodiode / nm_tip) proportionality constant - Ts measured temperature (K) - vibs measured V_photodiode variance in free solution (V**2) - Outputs : - k cantilever spring constant (in N/m, or equivalently nN/nm) - k_s standard deviation in our estimate of k - !!!a2_r relative error in a**2 - !!!T_r relative error in T - !!!one_o_Vp2_r relative error in 1/Vphotodiode_variance - Notes : - We're assuming vib is mostly from thermal cantilever vibrations - (and then only from vibrations in the single vertical degree of freedom), - and not from other noise sources. - The various relative errors are returned to help you gauge the - largest source of random error in your measurement of k. - If one of them is small, don't bother repeating that measurment too often. - If one is large, try repeating that measurement more. - Remember that you need enough samples to have a valid error estimate in - the first place, and that none of this addresses any systematic errors. +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): + """Analyze data from `get_calibration_data()` + + Inputs (all are arrays of recorded data): + bumps measured (V_photodiode / nm_tip) proportionality constant + temperatures measured temperature (K) + vibrations measured V_photodiode variance in free solution (V**2) + Outputs: + k cantilever spring constant (in N/m, or equivalently nN/nm) + k_s standard deviation in our estimate of k + + Notes: + + We're assuming vib is mostly from thermal cantilever vibrations + (and then only from vibrations in the single vertical degree of + freedom), and not from other noise sources. + + If the error is large, check the relative errors + (`x.std()/x.mean()`)of your input arrays. If one of them is + small, don't bother repeating that measurment too often. If one + is large, try repeating that measurement more. Remember that you + need enough samples to have a valid error estimate in the first + place, and that none of this addresses any systematic errors. """ - calib_plot_kwargs, = calib_analyze._splitargs(calib_analyze, kwargs) - photoSensitivity2 = bumps**2 - one_o_Vphoto2 = 1/vibs - - photoSensitivity2_m = photoSensitivity2.mean() - T_m = Ts.mean() - one_o_Vphoto2_m = one_o_Vphoto2.mean() - # Vphoto / photoSensitivity = x - # k = kb T / = kb T photoSensitiviy**2 * (1e9nm/m)**2 / - # + ps_m = bumps.mean() # ps for photo-sensitivity + ps_s = bumps.std() + T_m = temperatures.mean() + T_s = temperatures.std() + v2_m = vibrations.mean() # average voltage variance + v2_s = vibrations.std() + + # Vphoto / photo_sensitivity = x + # k = kB T / = kB T photo_sensitivity**2 / Vphoto_var # - # units, photoSensitivity = Vphoto/(Zcant in nm), - # so Vphoto/photoSensitivity = Zcant in nm - # so k = J/K * K / nm^2 * (1e9nm/m)**2 = N/m - k = kb * T_m * photoSensitivity2_m * one_o_Vphoto2_m * 1e18 - - # propogation of errors !!! - # first, get standard deviations - photoSensitivity2_s = photoSensitivity2.std() - T_s = Ts.std() - one_o_Vphoto2_s = one_o_Vphoto2.std() - # !!!!now, get relative errors - photoSensitivity2_r = photoSensitivity2_s / photoSensitivity2_m - T_r = T_s / T_m - one_o_Vphoto2_r = one_o_Vphoto2_s / one_o_Vphoto2_m - - k_s = k*(photoSensitivity2_r**2 + T_r**2 + one_o_Vphoto2_r**2)**0.5 - - calib_plot(bumps, Ts, vibs, **calib_plot_kwargs) - - return (k, k_s, - photoSensitivity2_m, photoSensitivity2_s, - T_m, T_s, one_o_Vphoto2_m, one_o_Vphoto2_s) - -@splittableKwargsFunction() -def string_errors(k_m, k_s, - photoSensitivity2_m, photoSensitivity2_s, - T_m, T_s, - one_o_Vphoto2_m, one_o_Vphoto2_s) : - k_r = k_s / k_m - photoSensitivity2_r = photoSensitivity2_s / photoSensitivity2_m - T_r = T_s / T_m - one_o_Vphoto2_r = one_o_Vphoto2_s / one_o_Vphoto2_m - string = "Variable (units) : mean +/- std. dev. (relative error)\n" - string += "Cantilever k (N/m) : %g +/- %g (%g)\n" \ - % (k_m, k_s, k_r) - string += "photoSensitivity**2 (V/nm)**2 : %g +/- %g (%g)\n" \ - % (photoSensitivity2_m, photoSensitivity2_s, photoSensitivity2_r) - string += "T (K) : %g +/- %g (%g)\n" \ - % (T_m, T_s, T_r) - string += "1/Vphoto**2 (1/V)**2 : %g +/- %g (%g)" \ - % (one_o_Vphoto2_m, one_o_Vphoto2_s, one_o_Vphoto2_r) - return string - -@splittableKwargsFunction() -def calib_save(bumps, Ts, vibs, log_dir=None) : - """ - Save a dictonary with the bump, T, and vib data. - """ - if log_dir != None : - data = {'bump':bumps, 'T':Ts, 'vib':vibs} - log = data_logger.data_log(log_dir, noclobber_logsubdir=False, - log_name="calib") - log.write_dict_of_arrays(data) - -def calib_load(datafile) : - "Load the dictionary data, using data_logger.date_load()" - dl = data_logger.data_load() - data = dl.read_dict_of_arrays(path) - return (data['bump'], data['T'], data['vib']) - -def calib_save_analysis(k, k_s, - photoSensitivity2_m, photoSensitivity2_s, - T_m, T_s, one_o_Vphoto2_m, one_o_Vphoto2_s, - log_dir=None) : - if log_dir != None : - log = data_logger.data_log(log_dir, noclobber_logsubdir=False, - log_name="calib_analysis_text") - log.write_binary(string_errors(k, k_s, - photoSensitivity2_m, photoSensitivity2_s, - T_m, T_s, - one_o_Vphoto2_m, one_o_Vphoto2_s) - +'\n') - -@splittableKwargsFunction() -def calib_plot(bumps, Ts, vibs, plotVerbose=False) : - if plotVerbose or config.PYLAB_VERBOSE : - common._import_pylab() - common._pylab.figure(config.BASE_FIGNUM+4) - common._pylab.subplot(311) - common._pylab.plot(bumps, 'g.-') - common._pylab.title('Photodiode sensitivity (V/nm)') - common._pylab.subplot(312) - common._pylab.plot(Ts, 'r.-') - common._pylab.title('Temperature (K)') - common._pylab.subplot(313) - common._pylab.plot(vibs, 'b.-') - common._pylab.title('Thermal deflection variance (Volts**2)') - common._flush_plot() - -make_splittable_kwargs_function(calib_analyze, - (calib_plot, 'bumps', 'Ts', 'vibs')) - -@splittableKwargsFunction((calib_analyze, 'bumps', 'Ts', 'vibs')) -def calib_load_analyze_tweaks(bump_tweaks, vib_tweaks, T_tweaks=None) : - raise NotImplementedError - a = read_tweaked_bumps(bump_tweaks) - vib = V_photo_variance_from_file(vib_tweaks) - if T_tweaks == None: - pass - return analyze_calibration_data(a, T, vib, log_dir=log_dir) - -# commandline interface functions -import scipy.io, sys - -def array_from_string(string): - ret = [] - for num in string.split(',') : - ret.append(float(num)) - assert len(ret) > 0 - return numpy.array(ret) - -def read_data(ifile): - "ifile can be a filename string or open (seekable) file object" - unlabeled_data=scipy.io.read_array(file) - return unlabeled_data - -def get_array(string, filename, name) : - "get an array from supplied command line options" - if string != None : - array = array_from_string(string) - elif filename != None : - array = read_data(filename) - else : - raise Exception, "no %s information given" % (name) - return array - -if __name__ == '__main__' : - # command line interface - from optparse import OptionParser - - usage_string = ('%prog \n' - '2008, W. Trevor King.\n' - '\n' - 'Takes arrays of Vphotodiode sensitivity (V/nm), Temperature (K), \n' - 'and Vibration variance (V**2) as comma seperated lists.\n' - 'Returns the cantilever spring constant (pN/nm).\n' - 'for example:\n' - ' $ %prog -b 0.02,0.03,0.025 -t 298.2,300.1 -v 6e-9,5.5e-9\n' - ) - parser = OptionParser(usage=usage_string, version='%prog '+common.VERSION) - parser.add_option('-b', '--bump-string', dest='bump_string', - help='comma seperated photodiode sensitivities (V/nm)', - type='string', metavar='BUMPS') - parser.add_option('-t', '--temp-string', dest='temp_string', - help='comma seperated temperatures (K)', - type='string', metavar='TEMPS') - parser.add_option('-v', '--vib-string', dest='vib_string', - help='comma seperated vibration variances (V**2)', - type='string', metavar='VIBS') - parser.add_option('-B', '--bump-file', dest='bump_file', - help='comma seperated photodiode sensitivities (V/nm)', - type='string', metavar='BUMPFILE') - parser.add_option('-T', '--temp-file', dest='temp_file', - help='comma seperated temperatures (K)', - type='string', metavar='TEMPFILE') - parser.add_option('-V', '--vib-file', dest='vib_file', - help='comma seperated vibration variances (V**2)', - type='string', metavar='VIBFILE') - parser.add_option('-C', '--celsius', dest='celsius', - help='Use Celsius input temperatures instead of Kelvin (defaul %default)\n', - action='store_true', default=False) - parser.add_option('-o', '--output-file', dest='ofilename', - help='write output to FILE (default stdout)', - type='string', metavar='FILE') - parser.add_option('-p', '--plot-inputs', dest='plot', - help='plot the input arrays to check their distribution', - action='store_true', default=False) - parser.add_option('', '--verbose', dest='verbose', action='store_true', - help='print lots of debugging information', - default=False) - - options,args = parser.parse_args() - parser.destroy() - - bumps = get_array(options.bump_string, options.bump_file, "bump") - temps = get_array(options.temp_string, options.temp_file, "temp") - vibs = get_array(options.vib_string, options.vib_file, "vib") - - if options.plot : - calib_plot(bumps, temps, vibs) - - if options.celsius : - for i in range(len(temps)) : - temps[i] = T_analyze.C_to_K(temps[i]) - - km,ks,ps2m,ps2s,Tm,Ts,ooVp2m,ooVp2s = calib_analyze(bumps, temps, vibs) - - if options.verbose : - print >> sys.stderr, \ - string_errors(km,ks,ps2m,ps2s,Tm,Ts,ooVp2m,ooVp2s) - - if options.ofilename != None : - print >> file(options.ofilename, 'w'), km - else : - print km + # units, photo_sensitivity = Vphoto/(Zcant in m), + # so Vphoto/photo_sensitivity = Zcant in m + # so k = J/K * K / m^2 = J / m^2 = N/m + k = _kB * T_m * ps_m**2 / v2_m + + # propogation of errors + # dk/dT = k/T + dk_T = k/T_m * T_s + # dk/dps = 2k/ps + dk_ps = 2*k/ps_m * ps_s + # dk/dv2 = -k/v2 + dk_v = -k/v2_m * v2_s + + k_s = _numpy.sqrt(dk_T**2 + dk_ps**2 + dk_v**2) + + _LOG.info('variable (units) : ' + 'mean +/- std. dev. (relative error)') + _LOG.info('cantilever k (N/m) : %g +/- %g (%g)' % (k, k_s, k_s/k)) + _LOG.info('photo sensitivity (V/m) : %g +/- %g (%g)' + % (ps_m, ps_s, ps_s/ps_m)) + _LOG.info('T (K) : %g +/- %g (%g)' + % (T_m, T_s, T_s/T_m)) + _LOG.info('vibration variance (V^2) : %g +/- %g (%g)' + % (v2_m, v2_s, v2_s/v2_m)) + + if _package_config['matplotlib']: + calib_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): + if not _matplotlib: + raise _matplotlib_import_error + figure = _matplotlib_pyplot.figure() + + bump_axes = figure.add_subplot(3, 1, 1) + T_axes = figure.add_subplot(3, 1, 2) + vib_axes = figure.add_subplot(3, 1, 3) + + timestamp = _time.strftime('%H%M%S') + bump_axes.set_title('cantilever calibration %s' % timestamp) + + bump_axes.plot(bumps, 'g.-') + bump_axes.set_ylabel('photodiode sensitivity (V/m)') + T_axes.plot(temperatures, 'r.-') + T_axes.set_ylabel('temperature (K)') + vib_axes.plot(vibrations, 'b.-') + vib_axes.set_ylabel('thermal deflection variance (V^2)') + + if hasattr(figure, 'show'): + 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/') + 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']): + bump_group = '%sbump/%d/' % (group, i) + (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_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)) + 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 + 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 + 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_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'] + 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'] + 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))