如何对齐两个大小不等的时间序列numpy数组？
Posted
技术标签:
【中文标题】如何对齐两个大小不等的时间序列numpy数组？【英文标题】：How to align two unequal sized timeseries numpy array? 【发布时间】：2012-06-03 00:46:36 【问题描述】：我有两个包含时间序列（unix 时间戳）的 numpy 数组。 我想找到差异在阈值内的时间戳对（每个数组中的1个）。
为了实现这一点，我需要将两个时间序列数据对齐到两个数组中，这样每个索引都有最接近的对。 （如果数组中的两个时间戳与另一个数组中的另一个时间戳同样接近，我不介意选择其中一个，因为对的数量比实际值更重要。）
所以对齐的数据集将有两个大小相同的数组，外加一个用空数据填充的较小数组。
我正在考虑使用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)，它最接近的对应该是@987654325@，因为差异是5703.0，它在threshold 之内，并且是最小的。现在ts2[2] 和ts1[4] 已经配对，它们应该被排除在其他计算之外。
应该找到这样的对，因此输出数组可能比实际数组长
abs(ts1[0]-ts2[0]) = 16060abs(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]) = 24329abs(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]) = 25461abs(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。
如果这个问题清楚，请告诉我。
【问题讨论】：
您能否发布几行代码来创建一个小示例（或组成）数据集以及您期望的结果。即data_a = np.array([12345, 12846, 789789]) 等。会帮助那些试图帮助你的人。
@Neo 我知道这不是你的帖子，但你需要它来做你正在做的事情。在描述中并不清楚这些对是如何获得的，因为在给定示例中选择的对的差异范围在 1290-15820 之间。这是一个相当大的价差，而且有比选定的更接近的对。我知道一旦获得就删除对，但除此之外，这个问题不能用给出的信息来回答。
【参考方案1】：
使用numpy Mask arrays 输出对齐时间序列的解决方案（_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)
输出：
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
【讨论】：
【参考方案2】：我不知道您所说的对齐时间戳是什么意思。但是您可以使用 time 模块将时间戳表示为浮点数或整数。第一步，您可以将任何用户格式转换为由time.struct_time 定义的数组。在第二步中，您可以将其转换为纪元开始的秒数。然后你有整数值来计算时间戳。
如何使用time.strptime()转换用户格式在docs中有很好的解释：
>>> 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
【讨论】：
【参考方案3】：除了问题中的小错误，我可以猜出问题的真正含义。
您正在查看的是Assignment Problem 的经典示例。 Scipy 为您提供了Hungarian algorithm 的实现，请查看文档here。它不必是时间戳，可以是任何数字（整数或浮点数）。
下面的 sn-p 将使用 2 个不同大小的 numpy 数组以及一个阈值，为您提供成本数组（按阈值过滤）或对应于 2 个 numpy 数组的索引对（同样，成本对是按阈值过滤）。
cmets 将引导您完成，使用给定的时间戳数组作为示例。
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)]
【讨论】：
【参考方案4】：为了获得您的时间序列对，我建议您首先计算您的索引对 (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)]
【讨论】：
【参考方案5】：我不确定你的问题是否正确。如果是这样并假设您的数据已经排序，则可以在使用迭代器时一次性完成。只需根据您的需要调整示例即可。
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)
【讨论】：
【参考方案6】：你有两个排序的时间戳列表，你需要将它们合并为一个，保持每个列表的元素相互分离，计算列表切换或变化时的差异。
我的第一个解决方案是不使用 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，这两个函数都会产生相同的结果。可以根据需要进行打印和格式化。哪个功能更好取决于您手头的数据。
注意：如果切换到另一个列表，并且在只有一个值之后切换回上一个列表，则上面的函数将只保留最后一个差异，也许会失去更小的差异。在这种情况下，您可以在“#### See Note”所在的位置插入以下部分：
对于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])
【讨论】：
以上是关于如何对齐两个大小不等的时间序列numpy数组？的主要内容，如果未能解决你的问题，请参考以下文章