Python timeseries.TimeSeries类代码示例

def fftgram(timeseries, stride, pad=False):
    """
    calculates fourier-gram with automatic
    50% overlapping hann windowing.
    Parameters
    ----------
        timeseries : gwpy TimeSeries
            time series to take fftgram of
        stride : `int`
            number of seconds in single PSD
    Returns
    -------
        fftgram : gwpy spectrogram
            a fourier-gram
    """

    fftlength = stride
    dt = stride
    df = 1. / fftlength
    # number of values in a step
    stride *= timeseries.sample_rate.value
    # number of steps
    nsteps = 2 * int(timeseries.size // stride) - 1
    # only get positive frequencies
    if pad:
        nfreqs = int(fftlength * timeseries.sample_rate.value)
        df = df / 2
        f0 = df
    else:
        nfreqs = int(fftlength * timeseries.sample_rate.value) / 2
    dtype = np.complex
    # initialize the spectrogram
    out = Spectrogram(np.zeros((nsteps, nfreqs), dtype=dtype),
                      name=timeseries.name, epoch=timeseries.epoch,
                      f0=df, df=df, dt=dt, copy=True,
                      unit=1 / u.Hz**0.5, dtype=dtype)
    # stride through TimeSeries, recording FFTs as columns of Spectrogram
    for step in range(nsteps):
        # indexes for this step
        idx = (stride / 2) * step
        idx_end = idx + stride
        stepseries = timeseries[idx:idx_end]
        # zeropad, window, fft, shift zero to center, normalize
        # window
        stepseries = np.multiply(stepseries,
                                 np.hanning(stepseries.value.size))
        # take fft
        if pad:
            stepseries = TimeSeries(np.hstack((stepseries, np.zeros(stepseries.size))),
                                    name=stepseries.name, x0=stepseries.x0,
                                    dx=timeseries.dx)
            tempfft = stepseries.fft(stepseries.size)
        else:
            tempfft = stepseries.fft(stepseries.size)
        tempfft.override_unit(out.unit)

        out[step] = tempfft[1:]

    return out

def dump_calibrated_data(fname):
    data = numpy.load(fname)

    # Figure out the times covered by the file from the filename
    # I should start using HDF5 so I can store metadata
    temp = fname.split('.')[0]
    temp = temp.split('-')
    ifo = temp[0]
    st, dur = int(temp[-2]), int(temp[-1])
    et = st + dur

    maxidx = len(data)
    width = 45

    weights = 1. - ((numpy.arange(-width, width) / float(width))**2)

    # The VCO frequencies are integers so we could dither them
    # to avoid quantization error if we wanted to be fancy
    # but it seems to make no differece
    if False:
        from numpy.random import triangular
        data[:, 1] += triangular(-1., 0., 1., size=len(data))

    # Just fit the whole thing at once, to get a single coefficient
    a, b = numpy.polyfit(data[:, 0], data[:, 1], 1)
    print "%.1f %u" % (a, b)

    # Slide through the data fitting PSL to IMC for data around each sample
    coeffs = []
    for idx in xrange(maxidx):
        idx1 = max(0, idx - width)
        idx2 = min(idx + width, maxidx)
        coeffs.append(numpy.polyfit(data[idx1:idx2, 0], data[idx1:idx2, 1], 1,
                                    w=weights[idx1 - idx + width:idx2 - idx + width]))
    coeffs = numpy.array(coeffs)
    times = numpy.arange(len(coeffs)) + 0.5
    connection = datafind.GWDataFindHTTPConnection()
    cache = connection.find_frame_urls(
        ifo[0], '%s_R' % ifo, st, et, urltype='file')

    imc = TimeSeries.read(cache, "%s:IMC-F_OUT_DQ" % ifo, st, et)
    imc = imc[::16384 / 256]
    print imc
    samp_times = numpy.arange(len(imc)) / 256.

    coeffs0 = numpy.interp(samp_times, times, coeffs[:, 0])
    coeffs1 = numpy.interp(samp_times, times, coeffs[:, 1]) - 7.6e7

    vco_interp = coeffs0 * imc.data + coeffs1

    chan = "%s:IMC-VCO_PREDICTION" % (ifo,)
    vco_data = TimeSeries(vco_interp, epoch=st,
                          sample_rate=imc.sample_rate.value,
                          name=chan, channel=chan)
    vco_data.write("%s-vcoprediction-%u-%u.hdf" % (ifo, st, dur), format='hdf')

def noise_from_psd(length, sample_rate, psd, seed=0, name=None, unit=u.m):
    """ Create noise with a given psd.

    Return noise with a given psd. Note that if unique noise is desired
    a unique seed should be provided.
    Parameters
    ----------
    length : int
        The length of noise to generate in seconds.
    sample_rate : float
       the sample rate of the data
    stride : int
        Length of noise segments in seconds
    psd : FrequencySeries
        The noise weighting to color the noise.
    seed : {0, int}
        The seed to generate the noise.

    Returns
    --------
    noise : TimeSeries
        A TimeSeries containing gaussian noise colored by the given psd.
    """
    if name is None:
        name='noise'
    length *=sample_rate

    noise_ts = TimeSeries(np.zeros(length),
            sample_rate=sample_rate, name=name, unit=unit)

    randomness = lal.gsl_rng("ranlux", seed)

    N = int (sample_rate / psd.df.value)
    n = N/2+1
    stride = N/2

    if n > len(psd):
        raise ValueError("PSD not compatible with requested delta_t")
    psd = (psd[0:n]).to_lal()
    psd.data.data[n-1] = 0
    segment = TimeSeries(np.zeros(N), sample_rate=sample_rate).to_lal()
    length_generated = 0

    SimNoise(segment, 0, psd, randomness)
    while (length_generated < length):
        if (length_generated + stride) < length:
            noise_ts.data[length_generated:length_generated+stride] = segment.data.data[0:stride]
        else:
            noise_ts.data[length_generated:length] = segment.data.data[0:length-length_generated]

        length_generated += stride
        SimNoise(segment,stride, psd, randomness)
    return noise_ts

    def refresh(self):
        lines = [l for ax in self._fig.axes for l in ax.lines]
        axes = cycle(self._fig.get_axes(self.AXES_CLASS.name))
        params = self.params['draw']
        for i, channel in enumerate(self.channels):
            try:
                line = lines[i]
            except IndexError:
                # haven't plotted this channel before
                ax = next(axes)
                label = (hasattr(channel, 'label') and channel.label or
                         channel.texname)
                pparams = {}
                for key in params:
                    try:
                        if params[key][i]:
                            pparams[key] = params[key][i]
                    except IndexError:
                        pass

                ts = self.data[channel][0].copy()
                for t2 in self.data[channel][1:]:
                    ts.append(t2, pad=self.buffer.pad, gap=self.buffer.gap)
                l = ax.plot(ts, label=label, **pparams)
                ax.legend()
            else:
                # TODO: remove .copy() as soon as copy=True is fixed in gwpy
                ts = TimeSeries(line.get_ydata(), times=line.get_xdata(),
                                copy=True).copy()
                for t2 in self.buffer.get((ts.span[1], self.epoch), channel,
                                          fetch=False):
                    ts.append(t2, pad=self.buffer.pad, gap=self.buffer.gap)
                line.set_xdata(ts.times.value)
                line.set_ydata(ts.value)

        # format figure
        if 'ylim' not in self.params['refresh']:
            for ax in self._fig.get_axes(self.AXES_CLASS.name):
                ax.relim()
                ax.autoscale_view(scalex=False)
        self.logger.info('Figure data updated')
        for ax in self._fig.get_axes(self.AXES_CLASS.name):
            if float(self.epoch) > (self.gpsstart + self.duration):
                ax.set_xlim(float(self.epoch - self.duration),
                            float(self.epoch))
            else:
                ax.set_xlim(float(self.gpsstart),
                            float(self.gpsstart) + self.duration)
            ax.set_epoch(self.epoch)
        self.set_params('refresh')
        self._fig.refresh()
        self.logger.debug('Figure refreshed')

def get_data(channel,start_time,length):
    print 'Starting data transfer for channel: ' + str(channel)
    connection = datafind.GWDataFindHTTPConnection()
    cache = connection.find_frame_urls(ifo[0],ifo+'_R',start_time,start_time+length,urltype='file')
    data = TimeSeries.read(cache,channel)
    print "Got data for channel: " + str(channel)
    return data

 def _next(self):
     uchannels = self._unique_channel_names(self.channels)
     new = TimeSeriesDict()
     span = 0
     epoch = 0
     self.logger.debug('Waiting for next NDS2 packet...')
     while span < self.interval:
         try:
             buffers = next(self.iterator)
         except RuntimeError as e:
             self.logger.error('RuntimeError caught: %s' % str(e))
             self.restart()
             break
         for buff, c in zip(buffers, uchannels):
             ts = TimeSeries.from_nds2_buffer(buff)
             try:
                 new.append({c: ts}, gap=self.gap, pad=self.pad)
             except ValueError as e:
                 if 'discontiguous' in str(e):
                     e.args = ('NDS connection dropped data between %d and '
                               '%d' % (epoch, ts.span[0]),)
                 raise
             span = abs(new[c].span)
             epoch = new[c].span[-1]
             self.logger.debug('%ds data for %s received'
                               % (abs(ts.span), str(c)))
     out = type(new)()
     for chan in self.channels:
         out[chan] = new[self._channel_basename(chan)].copy()
     return out

 def _next(self):
     uchannels = self._unique_channel_names(self.channels)
     new = TimeSeriesDict()
     span = 0
     epoch = 0
     att = 0
     self.logger.debug('Waiting for next NDS2 packet...')
     while span < self.interval:
         try:
             buffers = next(self.iterator)
         except RuntimeError as e:
             self.logger.error('RuntimeError caught: %s' % str(e))
             if att < self.attempts:
                 att += 1
                 wait_time = att / 3 + 1
                 self.logger.warning(
                     'Attempting to reconnect to the nds server... %d/%d'
                     % (att, self.attempts))
                 self.logger.warning('Next attempt in minimum %d seconds' %
                                     wait_time)
                 self.restart()
                 sleep(wait_time - tconvert('now') % wait_time)
                 continue
             else:
                 self.logger.critical(
                     'Maximum number of attempts reached, exiting')
                 break
         att = 0
         for buff, c in zip(buffers, uchannels):
             ts = TimeSeries.from_nds2_buffer(buff)
             try:
                 new.append({c: ts}, gap=self.gap, pad=self.pad)
             except ValueError as e:
                 if 'discontiguous' in str(e):
                     e.message = (
                         'NDS connection dropped data between %d and '
                         '%d, restarting building the buffer from %d ') \
                         % (epoch, ts.span[0], ts.span[0])
                     self.logger.warning(str(e))
                     new = TimeSeriesDict()
                     new[c] = ts.copy()
                 elif ('starts before' in str(e)) or \
                         ('overlapping' in str(e)):
                     e.message = (
                         'Overlap between old data and new data in the '
                         'nds buffer, only the new data will be kept.')
                     self.logger.warning(str(e))
                     new = TimeSeriesDict()
                     new[c] = ts.copy()
                 else:
                     raise
             span = abs(new[c].span)
             epoch = new[c].span[-1]
             self.logger.debug('%ds data for %s received'
                               % (abs(ts.span), str(c)))
     out = type(new)()
     for chan in self.channels:
         out[chan] = new[self._channel_basename(chan)].copy()
     return out

def calibrate_imc_pslvco(ifo, start_time, dur, cache=None):
    st, et = start_time, start_time + dur
    if cache:
        pslvco = TimeSeries.read(cache, chan1_pat % ifo, start=st, end=et)
        imc = TimeSeries.read(cache, chan2_pat % ifo, start=st, end=et)
    else:
        imc = TimeSeries.fetch(chan2_pat % ifo, st, et)
        pslvco = TimeSeries.fetch(chan1_pat % ifo, st, et)

    arr_psl = pslvco[8::16]
    arr_imc = imc

    tmp1 = (arr_imc[8192::16384])[:-1]
    tmp2 = arr_psl[1:]
    a, b = numpy.polyfit(tmp1, tmp2, 1)

    return a, b

 def test_add_timeseries(self):
     a = TimeSeries([1, 2, 3, 4, 5], name='test name', epoch=0,
                    sample_rate=1)
     # test simple add using 'name'
     data.add_timeseries(a)
     self.assertIn('test name', globalv.DATA)
     self.assertEqual(globalv.DATA['test name'], [a])
     # test add using key kwarg
     data.add_timeseries(a, key='test key')
     self.assertIn('test key', globalv.DATA)
     self.assertEqual(globalv.DATA['test key'], [a])
     # test add to existing key with coalesce
     b = TimeSeries([6, 7, 8, 9, 10], name='test name 2', epoch=5,
                    sample_rate=1)
     data.add_timeseries(b, key='test key', coalesce=True)
     self.assertEqual(globalv.DATA['test key'],
                      [a.append(b, inplace=False)])

def psd():
    h5path = os.path.join(os.path.dirname(__file__), 'data',
                          'HLV-HW100916-968654552-1.hdf')
    try:
        data = TimeSeries.read(h5path, 'L1:LDAS-STRAIN', format='hdf5')
    except ImportError as e:
        pytest.skip(str(e))
    return data.psd(.4, overlap=.2, window=('kaiser', 24))

def _read_data(channel, st, et, frames=False):
    """
    get data, either from frames or from nds2
    """

    ifo = channel.split(':')[0]
    if frames:
        # read from frames
        connection = datafind.GWDataFindHTTPConnection()
	print ifo[0]
	if channel.split(':')[1] == 'GDS-CALIB_STRAIN':
	    cache = connection.find_frame_urls(ifo[0],ifo+'_HOFT_C00', st, et, urltype='file')
	else:
	    cache = connection.find_frame_urls(ifo[0], ifo + '_C', st, et,urltype='file')
        try:
            data = TimeSeries.read(cache, channel, st, et)
        except IndexError:
            cache = connection.find_frame_urls(ifo[0], ifo+'_R', st, et, urltype='file')
            data = TimeSeries.read(cache, channel, st, et)
    else:
        data = TimeSeries.fetch(channel, st, et)

    return data

def test_save_loudest_tile_features():
    # prepare input data
    channel = GW.channels[0]
    noise = TimeSeries(
        numpy.random.normal(loc=1, scale=.5, size=16384 * 68),
        sample_rate=16384, epoch=-34).zpk([], [0], 1)
    glitch = TimeSeries(
        signal.gausspulse(numpy.arange(-1, 1, 1./16384), bw=100),
        sample_rate=16384, epoch=-1) * 1e-4
    in_ = noise.inject(glitch)
    _, _, _, qgram, _, _, _ = core.scan(
        gps=0, channel=channel, xoft=in_, resample=4096, fftlength=8)

    # test loudest tiles
    channel.save_loudest_tile_features(qgram, correlate=glitch)
    assert channel.Q == numpy.around(qgram.plane.q, 1)
    assert channel.energy == numpy.around(qgram.peak['energy'], 1)
    assert channel.snr == numpy.around(qgram.peak['snr'], 1)
    assert channel.t == numpy.around(qgram.peak['time'], 3)
    assert channel.f == numpy.around(qgram.peak['frequency'], 1)
    assert channel.corr == numpy.around(glitch.max().value, 1)
    assert channel.delay == 0.0
    assert channel.stdev == glitch.std().value

def generate_fast_vco(ifo, segment, frames=False, fit=True):
    """
    Parameters:
    -----------
        ifo : start
            interferometer, e.g. 'L1'
        segment : array like
            time segment. first entry start second entry end
        frames : bool
            read from frames or nds2
        fit : bool
            fit from imc-f (default)
            or spline interpolation

    Returns:
    --------
        vco_data : saves file 'L1:IMC-VCO_PREDICTION-st-dur.hdf'
    """
    st = segment[0]
    et = segment[1]
    chan1_pat = '%s:SYS-TIMING_C_FO_A_PORT_11_SLAVE_CFC_FREQUENCY_5'
    chan2_pat = '%s:IMC-F_OUT_DQ'
    if frames:
        connection = datafind.GWDataFindHTTPConnection()
        cache = connection.find_frame_urls(
            ifo[0], '%s_R' % ifo, st, et + 1, urltype='file')
        if fit:
            imc = TimeSeries.read(cache, chan2_pat % ifo, st, et)
        else:
            imc = TimeSeries.read(cache, chan2_pat % ifo, st, st + 1)
        pslvco = TimeSeries.read(cache, chan1_pat % ifo, st, et + 1)
    else:
        if fit:
            imc = TimeSeries.fetch(chan2_pat % ifo, st, et)
        else:
	    print 'HI BEFORE LOADING IMC'
            imc = TimeSeries.fetch(chan2_pat % ifo, st, st + 1)
	print 'HI BEFORE LOADING PSL'
        pslvco = TimeSeries.fetch(chan1_pat % ifo, st, et + 1)
	print 'HI AFTER LOADING'

    pslvco = pslvco[16 + 8::16]

    if fit:
        imc_srate = int(imc.sample_rate.value)
        imc2 = imc[imc_srate / 2::imc_srate]
        data = np.array((imc2.value, pslvco.value)).T
        vco_interp = fit_with_imc(data, imc)
    else:
        vco_interp = interp_spline(pslvco)

    chan = "%s:IMC-VCO_PREDICTION" % (ifo,)
    vco_data = TimeSeries(vco_interp, epoch=st,
                          sample_rate=256,
                          name=chan, channel=chan)

    return vco_data

def _get_vco_data(vco_file, times):
    print('Finding mean values of %s' % channel)
    s = times.size
    amp = numpy.zeros(s)
    vco = TimeSeries.from_hdf5(vco_file)
    for i, t in enumerate(times):
        t = t.seconds + t.nanoseconds / 1e9
        idx1 = int(
            (t - vco.times.value[0] - 0.01) * vco.sample_rate.value)
        idx2 = int(
            (t - vco.times.value[0] + 0.01) * vco.sample_rate.value)
        temp = vco[idx1:idx2]
        amp[i] = (temp.mean().value - 3e6) / 1.e3
        print('    Processed trigger %d/%d' % (i + 1, s), end='\r')
    print('    Processed trigger %d/%d' % (s, s))
    return amp

 def test_fetch(self):
     try:
         nds_buffer = mockutils.mock_nds2_buffer(
             'X1:TEST', self.data, 1000000000, self.data.shape[0], 'm')
     except ImportError as e:
         self.skipTest(str(e))
     nds_connection = mockutils.mock_nds2_connection([nds_buffer])
     with mock.patch('nds2.connection') as mock_connection, \
          mock.patch('nds2.buffer', nds_buffer):
         mock_connection.return_value = nds_connection
         # use verbose=True to hit more lines
         ts = TimeSeries.fetch('X1:TEST', 1000000000, 1000000001,
                               verbose=True)
     nptest.assert_array_equal(ts.value, self.data)
     self.assertEqual(ts.sample_rate, self.data.shape[0] * units.Hz)
     self.assertTupleEqual(ts.span, (1000000000, 1000000001))
     self.assertEqual(ts.unit, units.meter)