logo标签列表

如何对齐两个不同大小的时间序列numpy数组？

pythonalignmenttime-seriesclosest-points
14
我有两个包含时间序列（Unix 时间戳）的 numpy 数组。我想要找到一对时间戳（每个数组中各一个），它们之间的差值在阈值范围内。
为了实现这个目标，我需要将两个时间序列数据对齐到两个数组中，使得每个索引都有最接近的一对。如果两个数组中有两个时间戳与另一个数组中的某个时间戳相距相等，那么我不介意选择其中任何一个，因为成对数比实际值更重要。
因此，对齐后的数据集将有两个相同大小的数组，以及一个填充了空数据的较小数组。
我考虑使用 timeseries 包和 align 函数。但是我不确定如何使用对齐函数来处理我的时间序列数据。
例如，考虑两个时间序列数组：
ts1=np.array([ 1311242821.0, 1311242882.0, 1311244025.0, 1311244145.0, 1311251330.0, 
               1311282555.0, 1311282614.0])
ts2=np.array([ 1311226761.0, 1311227001.0, 1311257033.0, 1311257094.0, 1311281265.0])

输出示例：

现在对于ts2[2] (1311257033.0)，它最接近的配对应该是ts1[4] (1311251330.0)，因为差值为5703.0，在threshold范围内，并且它是最小的。既然ts2[2]和ts1[4]已经配对，它们应该从其他计算中排除。

应找到这样的配对，因此输出数组可能比实际数组长。

abs(ts1[0]-ts2[0]) = 16060
abs(ts1[0]-ts2[1]) = 15820 //一对
abs(ts1[0]-ts2[2]) = 14212
abs(ts1[0]-ts2[3]) = 14273
abs(ts1[0]-ts2[4]) = 38444

abs(ts1[1]-ts2[0]) = 16121
abs(ts1[1]-ts2[1]) = 15881
abs(ts1[1]-ts2[2]) = 14151
abs(ts1[1]-ts2[3]) = 14212
abs(ts1[1]-ts2[4]) = 38383

abs(ts1[2]-ts2[0]) = 17264
abs(ts1[2]-ts2[1]) = 17024
abs(ts1[2]-ts2[2]) = 13008
abs(ts1[2]-ts2[3]) = 13069
abs(ts1[2]-ts2[4]) = 37240

abs(ts1[3]-ts2[0]) = 17384
abs(ts1[3]-ts2[1]) = 17144
abs(ts1[3]-ts2[2]) = 12888
abs(ts1[3]-ts2[3]) = 17144
abs(ts1[3]-ts2[4]) = 37120

abs(ts1[4]-ts2[0]) = 24569
abs(ts1[4]-ts2[1]) = 24329
abs(ts1[4]-ts2[2]) = 5703 //一对
abs(ts1[4]-ts2[3]) = 5764
abs(ts1[4]-ts2[4]) = 29935

abs(ts1[5]-ts2[0]) = 55794
abs(ts1[5]-ts2[1]) = 55554
abs(ts1[5]-ts2[2]) = 25522
abs(ts1[5]-ts2[3]) = 25461
abs(ts1[5]-ts2[4]) = 1290 //一对

abs(ts1[6]-ts2[0]) = 55853
abs(ts1[6]-ts2[1]) = 55613
abs(ts1[6]-ts2[2]) = 25581
abs(ts1[6]-ts2[3]) = 25520
abs(ts1[6]-ts2[4]) = 1349

所以这些对应关系是: (ts1[0], ts2[1]), (ts1[4], ts2[2]), (ts1[5], ts2[4])。其余的元素应该使用null作为它们的对应关系。最终的两个数组大小将为9。请告知我是否清楚理解了这个问题。
4
你能否发布几行代码，创建一个小的示例（或虚构）数据集和您期望的结果。例如：data_a = np.array([12345, 12846, 789789])等。这将有助于那些试图帮助您的人。 - YXD
@Neo 我知道这不是你的帖子，但你需要它来完成你正在处理的事情。在描述中并没有清楚地说明如何获取这些配对数据，因为在给定的示例中，所选配对的差异范围在1290-15820之间。这是一个相当大的范围，而且有比所选的更接近的配对数据。我理解一旦获得了这些配对数据就要将其删除，但除此之外，根据提供的信息无法回答这个问题。 - Matt
6个回答
2

使用numpy Mask数组输出对齐的时间序列(_ts1, _ts2)解决方案。
结果有3对，只有距离为1的对可以用于对齐时间序列，因此阈值为1。

def compute_diffs(threshold):
    dtype = [('diff', int), ('ts1', int), ('ts2', int), ('threshold', int)]
    diffs = np.empty((ts1.shape[0], ts2.shape[0]), dtype=dtype)
    pairs = np.ma.make_mask_none(diffs.shape)

    for i1, t1 in enumerate(ts1):
        for i2, t2 in enumerate(ts2):
            diffs[i1, i2] = (abs(t1 - t2), i1, i2, abs(i1-i2))

        d1 = diffs[i1][diffs[i1]['threshold'] == threshold]
        if d1.size == 1:
            (diff, y, x, t) = d1[0]
            pairs[y, x] = True
    return diffs, pairs

def align_timeseries(diffs):
    def _sync(ts, ts1, ts2, i1, i2, ii):
        while i1 < i2:
            ts1[ii] = ts[i1]; i1 +=1
            ts2[ii] = DTNULL
            ii += 1
        return ii, i1

    _ts1 = np.array([DTNULL]*9)
    _ts2 = np.copy(_ts1)
    ii = _i1 = _i2 = 0

    for n, (diff, i1, i2, t) in enumerate(np.sort(diffs, order='ts1')):
        ii, _i1 = _sync(ts1, _ts1, _ts2, _i1, i1, ii)
        ii, _i2 = _sync(ts2, _ts2, _ts1, _i2, i2, ii)

        if _i1 == i1:
            _ts1[ii] = ts1[i1]; _i1 += 1
            _ts2[ii] = ts2[i2]; _i2 += 1
            ii += 1

    ii, _i1 = _sync(ts1, _ts1, _ts2, _i1, ts1.size, ii)
    return _ts1, _ts2

主函数：

diffs, pairs = compute_diffs(threshold=1)
print('diffs[pairs]:{}'.format(diffs[pairs]))
_ts1, _ts2 = align_timeseries(diffs[pairs])
pprint(ts1, ts2, _ts1, _ts2)

Output:

diffs[pairs]:[(15820, 0, 1) ( 5703, 4, 2) ( 1290, 5, 4)]
           ts1                  ts2                    _ts1          diff          _ts2
0: 2011-07-21 12:07:01  2011-07-21 07:39:21     ---- -- -- -- -- --  ----   2011-07-21 07:39:21
1: 2011-07-21 12:08:02  2011-07-21 07:43:21     2011-07-21 12:07:01 15820   2011-07-21 07:43:21
2: 2011-07-21 12:27:05  2011-07-21 16:03:53     2011-07-21 12:08:02  ----   ---- -- -- -- -- --
3: 2011-07-21 12:29:05  2011-07-21 16:04:54     2011-07-21 12:27:05  ----   ---- -- -- -- -- --
4: 2011-07-21 14:28:50  2011-07-21 22:47:45     2011-07-21 12:29:05  ----   ---- -- -- -- -- --
5: 2011-07-21 23:09:15  ---- -- -- -- -- --     2011-07-21 14:28:50  5703   2011-07-21 16:03:53
6: 2011-07-21 23:10:14  ---- -- -- -- -- --     ---- -- -- -- -- --  ----   2011-07-21 16:04:54
7: ---- -- -- -- -- --  ---- -- -- -- -- --     2011-07-21 23:09:15  1290   2011-07-21 22:47:45
8: ---- -- -- -- -- --  ---- -- -- -- -- --     2011-07-21 23:10:14  ----   ---- -- -- -- -- --

使用Python进行测试：3.4.2

DTNULL是什么？变量没有在问题/回答的范围内介绍。 - undefined
1
除了问题中的小错误，我可以猜到问题的真正原因。
你所看到的是分配问题的一个经典例子。Scipy为你提供了匈牙利算法的实现，可以在这里查看文档here。它不必是时间戳，可以是任何数字（整数或浮点数）。
下面的代码片段将与两个不同大小的numpy数组一起使用阈值，为你提供成本数组（通过阈值进行过滤）或对应于两个numpy数组的索引对（再次，其成本由阈值过滤）。
注释将使用给定的时间戳数组作为示例来引导你。
import numpy as np
from scipy.optimize import linear_sum_assignment


def closest_pairs(inp1, inp2, threshold=np.inf):
    cost = np.zeros((inp1.shape[0], inp2.shape[0]), dtype=np.int64)

    for x in range(ts1.shape[0]):
        for y in range(ts2.shape[0]):
            cost[x][y] = abs(ts1[x] - ts2[y])

    print(cost)
    # cost for the above example:
    # [[16060 15820 14212 14273 38444]
    # [16121 15881 14151 14212 38383]
    # [17264 17024 13008 13069 37240]
    # [17384 17144 12888 12949 37120]
    # [24569 24329  5703  5764 29935]
    # [55794 55554 25522 25461  1290]
    # [55853 55613 25581 25520  1349]]

    # hungarian algorithm implementation provided by scipy
    row_ind, col_ind = linear_sum_assignment(cost)
    # row_ind = [0 1 3 4 5], col_ind = [1 0 3 2 4] 
    # where (ts1[5] - ts2[4]) = 1290

    # if you want the distances only
    out = [item 
           for item in cost[row_ind, col_ind] 
           if item < threshold]

    # if you want the pair of indices filtered by the threshold
    pairs = [(row, col) 
             for row, col in zip(row_ind, col_ind) 
             if cost[row, col] < threshold]

    return out, pairs


if __name__ == '__main__':
    out, pairs = closest_pairs(ts1, ts2, 6000)
    print(out, pairs)
    # out = [5703, 1290] 
    # pairs = [(4, 2), (5, 4)]

    out, pairs = closest_pairs(ts1, ts2)
    print(out, pairs)
    # out = [15820, 16121, 12949, 5703, 1290] 
    # pairs = [(0, 1), (1, 0), (3, 3), (4, 2), (5, 4)]
1
为了得到您的时间序列对，我建议您首先计算索引对（get_pairs），然后计算时间序列对（get_tspairs）。
在get_pairs中，我首先计算一个矩阵m，它表示两个时间序列之间每个点的差异。因此，矩阵的形状为(len(ts1), len(ts2))。然后我选择所有最小距离中的最小值。为了不选择多次相同的索引，我将选定索引的距离设置为np.inf。我继续这个过程，直到我们不能再选择元组索引为止。如果最小距离高于阈值，则停止该过程。
一旦我得到我的索引对，我调用get_tspairs以生成时间序列对。这里的第一步是使用所选的元组索引连接时间序列，然后添加未选择的索引并将其与None（Python中的NULL等效项）相关联。
这就是：
import numpy as np
import operator

ts1=np.array([ 1311242821.0, 1311242882.0, 1311244025.0, 1311244145.0, 1311251330.0, 
               1311282555.0, 1311282614.0])
ts2=np.array([ 1311226761.0, 1311227001.0, 1311257033.0, 1311257094.0, 1311281265.0])

def get_pairs(ts1, ts2, threshold=np.inf):
    m = np.abs(np.subtract.outer(ts1, ts2))
    indices = []
    while np.ma.masked_invalid(m).sum() != 'masked':
        ind = np.unravel_index(np.argmin(m), m.shape)
        if m[ind] < threshold:
            indices.append(ind)
            m[:,ind[1]] = np.inf
            m[ind[0],:] = np.inf
        else:
            m= np.inf
    return indices

def get_tspairs(pairs, ts1, ts2):
    ts_pairs = [(ts1[p[0]], ts2[p[1]]) for p in pairs]

    # We separate the selected indices from ts1 and ts2, then sort them
    ind_ts1 = sorted(map(operator.itemgetter(0), pairs))
    ind_ts2 = sorted(map(operator.itemgetter(1), pairs))

    # We only keep the non-selected indices
    l1 = np.delete(np.arange(len(ts1), dtype=np.int64), ind_ts1)
    l2 = np.delete(np.arange(len(ts2), dtype=np.int64), ind_ts1)

    ts_pairs.extend([(ts1[i], None) for i in l1])
    ts_pairs.extend([(ts2[i], None) for i in l2])

    return ts_pairs

if __name__ == '__main__':
    pairs = get_pairs(ts1, ts2)
    print(pairs)
    # [(5, 4), (4, 2), (3, 3), (0, 1), (1, 0)]

    ts_pairs = get_tspairs(pairs, ts1, ts2)
    print(ts_pairs)
    # [(1311282555.0, 1311281265.0), (1311251330.0, 1311257033.0), (1311244145.0, 1311257094.0), (1311242821.0, 1311227001.0), (1311242882.0, 1311226761.0), (1311244025.0, None), (1311282614.0, None), (1311257033.0, None)]
在这行代码的末尾有一个拼写错误："l2 = np.delete(np.arange(len(ts2), dtype=np.int64), ind_ts1)"；应该是ind_ts2而不是ind_ts1。 - undefined
已修复，谢谢 @jointfull - undefined
1

我不确定我是否正确理解了你的问题。如果是这样，并且假设你的数据已经排序，那么在使用迭代器时可以一次完成此操作。只需根据您的需求调整示例即可。

left = iter(range(15, 60, 3))
right = iter(range(0, 50, 5))

try:
    i = next(left)
    j = next(right)
    while True:
        if abs(i-j) < 1:
            print("pair", i, j)
            i = next(left)
            j = next(right)
        elif i <= j:
            print("left", i, None)
            i = next(left)
        else:
            print("right", None, j)
            j = next(right)
except StopIteration:
    pass
# one of the iterators may have leftover elements
for i in left:
    print("left", i, None)
for j in right:
    print("right", None, j)

打印

('right', None, 0)                                                                                                                                                                      
('right', None, 5)                                                                                                                                                                      
('right', None, 10)                                                                                                                                                                     
('pair', 15, 15)                                                                                                                                                                        
('left', 18, None)                                                                                                                                                                      
('right', None, 20)                                                                                                                                                                     
('left', 21, None)                                                                                                                                                                      
('left', 24, None)                                                                                                                                                                      
('right', None, 25)                                                                                                                                                                     
('left', 27, None)                                                                                                                                                                      
('pair', 30, 30)                                                                                                                                                                        
('left', 33, None)                                                                                                                                                                      
('right', None, 35)                                                                                                                                                                     
('left', 36, None)                                                                                                                                                                      
('left', 39, None)                                                                                                                                                                      
('right', None, 40)                                                                                                                                                                     
('left', 42, None)                                                                                                                                                                      
('pair', 45, 45)                                                                                                                                                                        
('left', 51, None)                                                                                                                                                                      
('left', 54, None)                                                                                                                                                                      
('left', 57, None)
你有一个bug..你左边缺了48个..（停止迭代总是错过了桶里的另一个） - undefined
1

我不知道你所说的时间戳对齐是什么意思。但是你可以使用time模块将时间戳表示为浮点数或整数。首先，您可以将任何用户格式转换为由time.struct_time定义的数组。然后，您可以将其转换为自纪元开始的秒数。然后，您就有了整数值来计算时间戳。

如何使用time.strptime()转换用户格式在文档中有很好的解释：

    >>> import time
    >>> t = time.strptime("30 Nov 00", "%d %b %y")
    >>> t
    time.struct_time(tm_year=2000, tm_mon=11, tm_mday=30, tm_hour=0, tm_min=0,
             tm_sec=0, tm_wday=3, tm_yday=335, tm_isdst=-1)
    >>> time.mktime(t)
    975538800.0
0

你有两个已排序的时间戳列表，需要将它们合并成一个列表，保持每个列表的元素彼此分开，并在切换或更改列表时计算差异。

我的第一个解决方案不使用numpy，包括以下步骤：1）为每个元素添加所属列表的ID，2）按时间戳排序，3）按列表ID分组，4）构建新列表并在需要时分隔每个元素并计算差异：

import numpy as np
from itertools import groupby
from operator import itemgetter

ts1 = np.array([1311242821.0, 1311242882.0, 1311244025.0, 1311244145.0, 1311251330.0, 1311282555.0, 1311282614.0])               
ts2 = np.array([1311226761.0, 1311227001.0, 1311257033.0, 1311257094.0, 1311281265.0])

def without_numpy():
    # 1) Add the list id to each element
    all_ts = [(_, 0) for _ in ts1] + [(_, 1) for _ in ts2] 
    merged_ts = [[], [], []]
    # 2) Sort by timestamp and 3) Group by list id
    groups = groupby(sorted(all_ts), key=itemgetter(1))
    # 4) Construct the new list
    diff = False
    for key, g in groups:
        group = list(g)
        ### See Note    
        for ts, k in group:
            if diff:
                merged_ts[key] = merged_ts[key][:-1]
                merged_ts[2][-1] = abs(end - ts)
                diff = False
            else:
                merged_ts[not key].append(None)
                merged_ts[2].append(None)
            merged_ts[key].append(ts)
        end = ts
        diff = True
    return merged_ts

使用numpy，过程有些不同，包括以下步骤：1）为每个元素添加所属列表的ID和一些辅助索引，2）按时间戳排序，3）标记每个开关或列表更改，4）对先前标记进行总和扫描，5）计算合并列表中每个元素的正确索引：
import numpy as np

ts1 = np.array([1311242821.0, 1311242882.0, 1311244025.0, 1311244145.0, 1311251330.0, 1311282555.0, 1311282614.0])
ts2 = np.array([1311226761.0, 1311227001.0, 1311257033.0, 1311257094.0, 1311281265.0])

    def with_numpy(): 
        dt = np.dtype([('ts', np.float), ('key', np.int), ('idx', np.int)])

        all_ts = np.sort(
                np.array(
                    [(_, 0, 1, 0) for _ in ts1] + [(_, 1, 1, 0) for _ in ts2], 
                    dtype=np.dtype([('ts', np.float), 
                                    ('key', np.int),    # list id
                                    ('index', np.int),  # index in result list
                                    ('single', np.int), # flag groups with only one element
                                   ])
                ), 
                order='ts'
            )


        #### See NOTE

        sh_dn = np.roll(all_ts, 1)

        all_ts['index'] = np.add.accumulate(all_ts['index']) - np.cumsum(
                            np.not_equal(all_ts['key'], sh_dn['key']))

        merged_ts = np.full(shape=(3, all_ts['index'][-1]+1), fill_value=np.nan)
        merged_ts[all_ts['key'], all_ts['index']] = all_ts['ts']
        merged_ts[2] = np.abs(merged_ts[0] - merged_ts[1])
        merged_ts = np.delete(merged_ts, -1, axis=1)
        merged_ts = np.transpose(merged_ts)

        return merged_ts

两个函数，使用或不使用numpy都会产生相同的结果。打印和格式化可以根据需要进行。哪个函数更好取决于手头的数据。

注意：如果在切换到另一个列表后，仅有一个值再切换回先前的列表，则上述函数将仅保留最后一个差异，可能会丢失较小的差异。在这种情况下，您可以在“####请参阅说明”处插入以下部分：

对于without_numpy函数，请插入：

if len(group) == 1:
    group.append(group[0])

对于with_numpy函数，插入以下内容：
sh_dn = np.roll(all_ts, 1)
sh_up = np.roll(all_ts, -1)
all_ts['single'] = np.logical_and( np.not_equal(all_ts['key'], sh_dn['key']), 
                                           np.equal(sh_dn['key'], sh_up['key']))
singles = np.where(all_ts['single']==1)[0]
all_ts = np.insert(all_ts, singles, all_ts[singles])
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接