calibcant/vibration_analyze.py

   1 # calibcant - tools for thermally calibrating AFM cantilevers
   2 #
   3 # Copyright (C) 2008-2012 W. Trevor King <wking@drexel.edu>
   4 #
   5 # This file is part of calibcant.
   6 #
   7 # calibcant is free software: you can redistribute it and/or modify it under
   8 # the terms of the GNU General Public License as published by the Free Software
   9 # Foundation, either version 3 of the License, or (at your option) any later
  10 # version.
  11 #
  12 # calibcant is distributed in the hope that it will be useful, but WITHOUT ANY
  13 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  14 # A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
  15 #
  16 # You should have received a copy of the GNU General Public License along with
  17 # calibcant.  If not, see <http://www.gnu.org/licenses/>.
  18
  19 """Thermal vibration analysis.
  20
  21 Separate the more general `analyze()` from the other `vibration`
  22 functions in calibcant.
  23
  24 The relevent physical quantities are :
  25   Vphoto   The photodiode vertical deflection voltage (what we measure)
  26
  27 >>> import os
  28 >>> from pprint import pprint
  29 >>> import random
  30 >>> import tempfile
  31 >>> import numpy
  32 >>> from .config import VibrationConfig
  33 >>> from h5config.storage.hdf5 import pprint_HDF5
  34 >>> from pypiezo.test import get_piezo_config
  35 >>> from pypiezo.base import convert_volts_to_bits
  36
  37 >>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
  38 >>> os.close(fd)
  39
  40 >>> piezo_config = get_piezo_config()
  41 >>> config = VibrationConfig()
  42 >>> config['frequency'] = 50e3
  43
  44 We'll be generating a test vibration time series with the following
  45 parameters.  Make sure these are all floats to avoid accidental
  46 integer division in later steps.
  47
  48 >>> m = 5e-11       # kg
  49 >>> gamma = 1.6e-6  # N*s/m
  50 >>> k = 0.05        # N/m
  51 >>> T = 1/config['frequency']
  52 >>> T  # doctest: +ELLIPSIS
  53 2...e-05
  54 >>> N = int(2**15)  # count
  55 >>> F_sigma = 1e3   # N
  56
  57 where `T` is the sampling period, `N` is the number of samples, and
  58 `F_sigma` is the standard deviation of the white-noise external force.
  59 Note that the resonant frequency is less than the Nyquist frequency so
  60 we don't have to worry too much about aliasing.
  61
  62 >>> w0 = numpy.sqrt(k/m)
  63 >>> f0 = w0/(2*numpy.pi)
  64 >>> f0  # doctest: +ELLIPSIS
  65 5032.9...
  66 >>> f_nyquist = config['frequency']/2
  67 >>> f_nyquist  # doctest: +ELLIPSIS
  68 25000.0...
  69
  70 The damping ratio is
  71
  72 >>> damping = gamma / (2*m*w0)
  73 >>> damping  # doctest: +ELLIPSIS
  74 0.505...
  75
  76 The quality factor is
  77
  78 >>> Q = m*w0 / gamma
  79 >>> Q  # doctest: +ELLIPSIS
  80 0.988...
  81 >>> (1 / (2*damping)) / Q  # doctest: +ELLIPSIS
  82 1.000...
  83
  84 We expect the white-noise power spectral density (PSD) to be a flat
  85 line at
  86
  87 >>> F0 = F_sigma**2 * 2 * T
  88
  89 because the integral from `0` `1/2T` should be `F_sigma**2`.
  90
  91 The expected time series PSD parameters are
  92
  93 >>> A = f0
  94 >>> B = gamma/(m*2*numpy.pi)
  95 >>> C = F0/(m**2*(2*numpy.pi)**4)
  96
  97 Simulate a time series with the proper PSD using center-differencing.
  98
  99   m\ddot{x} + \gamma \dot{x} + kx = F
 100
 101   m \frac{x_{i+1} - 2x_i + x_{i-1}}{T**2}
 102     + \gamma \frac{x_{i+1}-x_{i-1}}{T}
 103     + kx_i = F_i
 104
 105   a x_{i+1} + b x_{i} + c x_{i-1} = F_i
 106
 107 where `T` is the sampling period, `i=t/T` is the measurement index,
 108 `a=m/T**2+gamma/2T`, `b=-2m/T**2+k`, and `c=m/T**2-gamma/2T`.
 109 Rearranging and shifting to `j=i-1`
 110
 111   x_j = \frac{F_{i-1} - bx_{i-1} - cx_{i-2}}{a}
 112
 113 >>> a = m/T**2 + gamma/(2*T)
 114 >>> b = -2*m/T**2 + k
 115 >>> c = m/T**2 - gamma/(2*T)
 116 >>> x = numpy.zeros((N+2,), dtype=numpy.float)  # two extra initial points
 117 >>> F = numpy.zeros((N,), dtype=numpy.float)
 118 >>> for i in range(2, x.size):
 119 ...     Fp = random.gauss(mu=0, sigma=F_sigma)
 120 ...     F[i-2] = Fp
 121 ...     xp = x[i-1]
 122 ...     xpp = x[i-2]
 123 ...     x[i] = (Fp - b*xp - c*xpp)/a
 124 >>> x = x[2:]  # drop the initial points
 125
 126 Convert the simulated data to bits.
 127
 128 >>> deflection = x
 129 >>> deflection_bits = convert_volts_to_bits(
 130 ...     piezo_config.select_config('inputs', 'deflection'), x)
 131
 132 Analyze the simulated data.
 133
 134 >>> naive = analyze_naive(deflection)
 135 >>> naive  # doctest: +SKIP
 136 136871517.43486854
 137 >>> abs(naive / 136.9e6 - 1) < 0.1
 138 True
 139
 140 >>> processed = analyze(
 141 ...     deflection_bits, config,
 142 ...     piezo_config.select_config('inputs', 'deflection'))
 143 >>> processed  # doctest: +SKIP
 144 136457906.25574699
 145
 146 >>> plot(deflection=deflection_bits, config=config)
 147 >>> save(filename=filename, group='/vibration/',
 148 ...     raw=deflection_bits, config=config,
 149 ...     deflection_channel_config=piezo_config.select_config(
 150 ...         'inputs', 'deflection'),
 151 ...     processed=processed)
 152
 153 >>> pprint_HDF5(filename)  # doctest: +ELLIPSIS, +REPORT_UDIFF
 154 /
 155   /vibration
 156     /vibration/config
 157       /vibration/config/deflection
 158         <HDF5 dataset "analog-reference": shape (), type "|S6">
 159           ground
 160         <HDF5 dataset "channel": shape (), type "<i4">
 161           0
 162         <HDF5 dataset "conversion-coefficients": shape (2,), type "<i4">
 163           [0 1]
 164         <HDF5 dataset "conversion-origin": shape (), type "<i4">
 165           0
 166         <HDF5 dataset "device": shape (), type "|S12">
 167           /dev/comedi0
 168         <HDF5 dataset "inverse-conversion-coefficients": shape (2,), type "<i4">
 169           [0 1]
 170         <HDF5 dataset "inverse-conversion-origin": shape (), type "<i4">
 171           0
 172         <HDF5 dataset "maxdata": shape (), type "<i4">
 173           100
 174         <HDF5 dataset "name": shape (), type "|S10">
 175           deflection
 176         <HDF5 dataset "range": shape (), type "<i4">
 177           0
 178         <HDF5 dataset "subdevice": shape (), type "<i4">
 179           -1
 180       /vibration/config/vibration
 181         <HDF5 dataset "chunk-size": shape (), type "<i4">
 182           2048
 183         <HDF5 dataset "frequency": shape (), type "<f8">
 184           50000.0
 185         <HDF5 dataset "maximum-fit-frequency": shape (), type "<f8">
 186           25000.0
 187         <HDF5 dataset "minimum-fit-frequency": shape (), type "<f8">
 188           500.0
 189         <HDF5 dataset "model": shape (), type "|S12">
 190           Breit-Wigner
 191         <HDF5 dataset "overlap": shape (), type "|b1">
 192           False
 193         <HDF5 dataset "sample-time": shape (), type "<i4">
 194           1
 195         <HDF5 dataset "window": shape (), type "|S4">
 196           Hann
 197     /vibration/processed
 198       <HDF5 dataset "data": shape (), type "<f8">
 199         ...
 200       <HDF5 dataset "units": shape (), type "|S6">
 201         V^2/Hz
 202     /vibration/raw
 203       <HDF5 dataset "data": shape (32768,), type "<f8">
 204         [...]
 205       <HDF5 dataset "units": shape (), type "|S4">
 206         bits
 207
 208 >>> data = load(filename=filename, group='/vibration/')
 209
 210 >>> pprint(data)  # doctest: +ELLIPSIS
 211 {'config': {'vibration': <InputChannelConfig ...>},
 212  'processed': ...,
 213  'raw': array([...])}
 214 >>> data['processed']  # doctest: +SKIP
 215 136457906.25574699
 216 >>> abs(data['processed'] / 136.5e6 - 1) < 0.1
 217 True
 218
 219 >>> os.remove(filename)
 220 """
 221
 222 import os as _os
 223 import time as _time
 224
 225 import h5py as _h5py
 226 import numpy as _numpy
 227 from scipy.optimize import leastsq as _leastsq
 228
 229 try:
 230     import matplotlib as _matplotlib
 231     import matplotlib.pyplot as _matplotlib_pyplot
 232     import time as _time  # for timestamping lines on plots
 233 except (ImportError, RuntimeError), e:
 234     _matplotlib = None
 235     _matplotlib_import_error = e
 236
 237 from h5config.storage.hdf5 import HDF5_Storage as _HDF5_Storage
 238 from h5config.storage.hdf5 import h5_create_group as _h5_create_group
 239 import FFT_tools as _FFT_tools
 240 from pypiezo.base import convert_bits_to_volts as _convert_bits_to_volts
 241 from pypiezo.config import InputChannelConfig as _InputChannelConfig
 242
 243 from . import LOG as _LOG
 244 from . import package_config as _package_config
 245 from .config import Variance as _Variance
 246 from .config import BreitWigner as _BreitWigner
 247 from .config import OffsetBreitWigner as _OffsetBreitWigner
 248 from .config import VibrationConfig as _VibrationConfig
 249 from .util import SaveSpec as _SaveSpec
 250 from .util import save as _save
 251 from .util import load as _load
 252
 253
 254 def analyze_naive(deflection):
 255     """Calculate the deflection variance in Volts**2.
 256
 257     This method is simple and easy to understand, but it highly
 258     succeptible to noise, drift, etc.
 259
 260     Inputs
 261       deflection : numpy array with deflection timeseries in Volts.
 262     """
 263     std = deflection.std()
 264     var = std**2
 265     _LOG.debug('naive deflection variance: %g V**2' % var)
 266     return var
 267
 268 def analyze(deflection, config, deflection_channel_config,
 269             plot=False):
 270     """Calculate the deflection variance in Volts**2.
 271
 272     Improves on `analyze_naive()` by first converting `Vphoto(t)`
 273     to frequency space, and fitting a Breit-Wigner in the relevant
 274     frequency range (see cantilever_calib.pdf for derivation).
 275     However, there may be cases where the fit is thrown off by noise
 276     spikes in frequency space.  To protect from errors, the fitted
 277     variance is compared to the naive variance (where all noise is
 278     included), and the minimum variance is returned.
 279
 280     Inputs:
 281       deflection        Vphoto deflection input in bits.
 282       config            `.config.VibrationConfig` instance
 283       deflection_channel_config
 284                         deflection `pypiezo.config.ChannelConfig` instance
 285       plot              boolean overriding matplotlib config setting.
 286
 287     The conversion to frequency space generates an average power
 288     spectrum by breaking the data into windowed chunks and averaging
 289     the power spectrums for the chunks together. See
 290     `FFT_tools.unitary_avg_power_spectrum()` for details.
 291     """
 292     # convert the data from bits to volts
 293     deflection_v = _convert_bits_to_volts(
 294         deflection_channel_config, deflection)
 295     mean = deflection_v.mean()
 296     _LOG.debug('average thermal deflection (Volts): %g' % mean)
 297
 298     naive_variance = analyze_naive(deflection_v)
 299     if config['model'] == _Variance:
 300         return naive_variance
 301
 302     # Compute the average power spectral density per unit time (in uV**2/Hz)
 303     _LOG.debug('compute the averaged power spectral density in uV**2/Hz')
 304     freq_axis,power = _FFT_tools.unitary_avg_power_spectrum(
 305         (deflection_v - mean)*1e6, config['frequency'],
 306         config['chunk-size'], config['overlap'],
 307         config['window'])
 308
 309     A,B,C,D = fit_psd(
 310         freq_axis, power,
 311         min_frequency=config['minimum-fit-frequency'],
 312         max_frequency=config['maximum-fit-frequency'],
 313         offset=config['model'] == _OffsetBreitWigner)
 314
 315     _LOG.debug('fit PSD(f) = C / ((A**2-f**2)**2 + (f*B)**2) with '
 316                'A = %g, B = %g, C = %g, D = %g' % (A, B, C, D))
 317
 318     if plot or _package_config['matplotlib']:
 319         _plot(deflection, freq_axis, power, A, B, C, D, config=config)
 320
 321     # Our A is in uV**2, so convert back to Volts**2
 322     fitted_variance = breit_wigner_area(A,B,C) * 1e-12
 323
 324     _LOG.debug('fitted deflection variance: %g V**2' % fitted_variance)
 325
 326     if _package_config['matplotlib']:
 327         plot(deflection, freq_axis, power, A, B, C, D,
 328                  config=config)
 329
 330     return min(fitted_variance, naive_variance)
 331
 332 def breit_wigner(f, A, B, C, D=0):
 333     """Breit-Wigner (sortof).
 334
 335     Inputs
 336       f  Frequency
 337       A  Resonant frequency
 338       B  Quality Q = A/B
 339       C  Scaling factor
 340       D  Optional white-noise offset
 341
 342     All parameters must be postive.
 343     """
 344     return abs(C) / ((A**2-f**2)**2 + (B*f)**2) + abs(D)
 345
 346 def fit_psd(freq_axis, psd_data, min_frequency=500, max_frequency=25000,
 347             offset=False):
 348     """Fit the FFTed vibration data to a Breit-Wigner.
 349
 350     Inputs
 351       freq_axis      array of frequencies in Hz
 352       psd_data       array of PSD amplitudes for the frequencies in freq_axis
 353       min_frequency  lower bound of Breit-Wigner fitting region
 354       max_frequency  upper bound of Breit-Wigner fitting region
 355       offset         add a white-noise offset to the Breit-Wigner model
 356     Output
 357       Breit-Wigner model fit parameters `A`, `B`, `C`, and `D`.
 358     """
 359     # cut out the relevent frequency range
 360     _LOG.debug('cut the frequency range %g to %g Hz'
 361                % (min_frequency, max_frequency))
 362     imin = 0
 363     while freq_axis[imin] < min_frequency : imin += 1
 364     imax = imin
 365     while freq_axis[imax] < max_frequency : imax += 1
 366     assert imax >= imin + 10 , 'less than 10 points in freq range (%g,%g)' % (
 367         min_frequency, max_frequency)
 368
 369     # generate guesses for Breit-Wigner parameters A, B, C, and D
 370     max_psd_index = _numpy.argmax(psd_data[imin:imax]) + imin
 371     max_psd = psd_data[max_psd_index]
 372     res_freq = freq_axis[max_psd_index]
 373
 374     # Breit-Wigner L(x) = C / ((A**2-x**2)**2 + (B*x)**2)
 375     # is expected power spectrum for
 376     # x'' + B x' + A^2 x'' = F_external(t)/m
 377     # (A = omega_0)
 378     # C = (2 k_B T B) / (pi m)
 379     #
 380     # A = resonant frequency
 381     # peak at  x_res = sqrt(A^2 - B^2/2)  (by differentiating)
 382     #  which could be complex if there isn't a peak (overdamped)
 383     # peak height    = C / (3 x_res^4 + A^4)
 384     #
 385     # Q = A/B
 386     #
 387     # Guessing Q = 1 is pretty safe.
 388
 389     Q_guess = 1
 390
 391     # so x_res^2 = B^2 Q^2 - B^2/2 = (Q^2-1/2)B^2
 392     #    B = x_res / sqrt(Q^2-1/2)
 393     B_guess = res_freq / _numpy.sqrt(Q_guess**2-0.5)
 394     A_guess = Q_guess*B_guess
 395     C_guess = max_psd * (-res_freq**4 + A_guess**4)
 396     if offset:
 397         D_guess = psd_data[max_psd_index]
 398         C_guess -= D_guess
 399     else:
 400         D_guess = 0
 401     _LOG.debug(('guessed params: resonant freq %g, max psd %g, Q %g, A %g, '
 402                 'B %g, C %g, D %g') % (
 403             res_freq, max_psd, Q_guess, A_guess, B_guess, C_guess, D_guess))
 404     # Half width w on lower side when L(a-w) = L(a)/2
 405     #  (a**2 - (a-w)**2)**2 + (b*(a-w))**2 = 2*(b*a)**2
 406     # Let W=(a-w)**2, A=a**2, and B=b**2
 407     #  (A - W)**2 + BW = 2*AB
 408     #  W**2 - 2AW + A**2 + BW = 2AB
 409     #  W**2 + (B-2A)W + (A**2-2AB) = 0
 410     #  W = (2A-B)/2 * [1 +/- sqrt(1 - 4(A**2-2AB)/(B-2A)**2]
 411     #    = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
 412     #  (a-w)**2 = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
 413     #  so w is a disaster ;)
 414     # For some values of A and B (non-underdamped), W is imaginary or negative.
 415     #
 416     # The mass m is given by m = k_B T / (pi A)
 417     # The spring constant k is given by k = m (omega_0)**2
 418     # The quality factor Q is given by Q = omega_0 m / gamma
 419
 420     # Fitting the PSD of V = photoSensitivity*x just rescales the parameters
 421
 422     # fit Breit-Wigner using scipy.optimize.leastsq
 423     def residual(p, y, x):
 424         return breit_wigner(x, *p) - y
 425     if offset:
 426         guess = _numpy.array((A_guess, B_guess, C_guess, D_guess))
 427     else:
 428         guess = _numpy.array((A_guess, B_guess, C_guess))
 429
 430     p,cov,info,mesg,ier = _leastsq(
 431         residual, guess,
 432         args=(psd_data[imin:imax], freq_axis[imin:imax]),
 433         full_output=True, maxfev=10000)
 434     _LOG.debug('fitted params: %s' % p)
 435     _LOG.debug('covariance matrix: %s' % cov)
 436     #_LOG.debug('info: %s' % info)
 437     _LOG.debug('message: %s' % mesg)
 438     if ier == 1:
 439         _LOG.debug('solution converged')
 440     else:
 441         _LOG.debug('solution did not converge')
 442     if offset:
 443         A,B,C,D = p
 444     else:
 445         A,B,C = p
 446         D = 0
 447     A=abs(A) # A and B only show up as squares in f(x)
 448     B=abs(B) # so ensure we get positive values.
 449     C=abs(C) # Only abs(C) is used in breit_wigner().
 450     return (A, B, C, D)
 451
 452 def breit_wigner_area(A, B, C):
 453     # Integrating the the power spectral density per unit time (power) over the
 454     # frequency axis [0, fN] returns the total signal power per unit time
 455     #  int_0^fN power(f)df = 1/T int_0^T |x(t)**2| dt
 456     #                      = <V(t)**2>, the variance for our AC signal.
 457     # The variance from our fitted Breit-Wigner is the area under the Breit-Wigner
 458     #  <V(t)**2> = (pi*C) / (2*B*A**2)
 459     return (_numpy.pi * C) / (2 * B * A**2)
 460
 461 def save(filename, group='/', raw=None, config=None,
 462          deflection_channel_config=None, processed=None):
 463     specs = [
 464         _SaveSpec(item=config, relpath='config/vibration',
 465                   config=_VibrationConfig),
 466         _SaveSpec(item=deflection_channel_config,
 467                   relpath='config/deflection',
 468                   config=_InputChannelConfig),
 469         _SaveSpec(item=raw, relpath='raw', units='bits'),
 470         _SaveSpec(item=processed, relpath='processed', units='V^2/Hz'),
 471         ]
 472     _save(filename=filename, group=group, specs=specs)
 473
 474 def load(filename=None, group='/'):
 475     specs = [
 476         _SaveSpec(key=('config', 'vibration'), relpath='config/vibration',
 477                   config=_VibrationConfig),
 478         _SaveSpec(key=('config', 'deflection'), relpath='config/deflection',
 479                   config=_InputChannelConfig),
 480         _SaveSpec(key=('raw',), relpath='raw', array=True, units='bits'),
 481         _SaveSpec(key=('processed',), relpath='processed', units='V^2/Hz'),
 482         ]
 483     return _load(filename=filename, group=group, specs=specs)
 484
 485 def plot(deflection=None, freq_axis=None, power=None, A=None, B=None,
 486              C=None, D=0, config=None, analyze=False):
 487     """Plot 3 subfigures displaying vibration data and analysis.
 488
 489      Time series (Vphoto vs sample index) (show any major drift events),
 490      A histogram of Vphoto, with a gaussian fit (demonstrate randomness), and
 491      FFTed Vphoto data (Vphoto vs Frequency) (show major noise components).
 492     """
 493     if not _matplotlib:
 494         raise _matplotlib_import_error
 495     figure = _matplotlib_pyplot.figure()
 496
 497     if power is None:
 498         assert deflection != None, (
 499             'must set at least one of `deflection` or `power`.')
 500         time_axes = figure.add_subplot(2, 1, 1)
 501         hist_axes = figure.add_subplot(2, 1, 2)
 502         freq_axes = None
 503     elif deflection is None:
 504         time_axes = hist_axes = None
 505         freq_axes = figure.add_subplot(1, 1, 1)
 506     else:
 507         time_axes = figure.add_subplot(3, 1, 1)
 508         hist_axes = figure.add_subplot(3, 1, 2)
 509         freq_axes = figure.add_subplot(3, 1, 3)
 510
 511     timestamp = _time.strftime('%H%M%S')
 512
 513     if deflection is not None:
 514         time_axes.plot(deflection, 'r.')
 515         time_axes.autoscale(tight=True)
 516         time_axes.set_title('free oscillation')
 517
 518         # plot histogram distribution and gaussian fit
 519         hist_axes.hold(True)
 520         n,bins,patches = hist_axes.hist(
 521             deflection, bins=30, normed=True, align='mid')
 522         gauss = _numpy.zeros((len(bins),), dtype=_numpy.float)
 523         mean = deflection.mean()
 524         std = deflection.std()
 525         pi = _numpy.pi
 526         exp = _numpy.exp
 527         gauss = _numpy.sqrt(2*pi)/std * exp(-0.5*((bins-mean)/std)**2)
 528         # Matplotlib's normed histogram uses bin heights of n/(len(x)*dbin)
 529         dbin = bins[1]-bins[0]
 530         hist_axes.plot(bins, gauss/dbin, 'r-')
 531         hist_axes.autoscale(tight=True)
 532     if power is not None:
 533         freq_axes.hold(True)
 534         freq_axes.set_yscale('log')
 535         freq_axes.plot(freq_axis, power, 'r.-')
 536         freq_axes.autoscale(tight=True)
 537         xmin,xmax = freq_axes.get_xbound()
 538         ymin,ymax = freq_axes.get_ybound()
 539
 540         # highlight the region we're fitting
 541         freq_axes.axvspan(
 542             config['minimum-fit-frequency'],
 543             config['maximum-fit-frequency'],
 544             facecolor='g', alpha=0.1, zorder=-2)
 545
 546         if A is not None:
 547             fitdata = breit_wigner(freq_axis, A, B, C, D)
 548             freq_axes.plot(freq_axis, fitdata, 'b-')
 549             noisefloor = D + 0*freq_axis;
 550             freq_axes.plot(freq_axis, noisefloor)
 551
 552             if B**2 < 2*A**2:
 553                 res_freq = _numpy.sqrt(A**2 - B**2/2)
 554                 freq_axes.axvline(res_freq, color='b', zorder=-1)
 555
 556         freq_axes.set_title('power spectral density %s' % timestamp)
 557         freq_axes.axis([xmin,xmax,ymin,ymax])
 558         freq_axes.set_xlabel('frequency (Hz)')
 559         freq_axes.set_ylabel('PSD')
 560     if hasattr(figure, 'show'):
 561         figure.show()
 562 _plot = plot  # alternative name for use inside analyze()