Numpy：在双数组中设置尾数的最后 n 个元素

Question

问题

假设我们有一个 numpy 数组 arr 双精度数和一个小的正整数 n。我正在寻找一种有效的方法来将 arr 的每个元素的 n 最低有效条目设置为 0 或 1。有 ufunc 吗？如果没有，是否有合适的 C 函数可以应用于 Cython 的元素？

动机

下面我将提供问题的动机。如果您发现上述问题的答案不需要实现最终目标，我很乐意收到相应的评论。然后，我将创建一个单独的问题以使事情井然有序。

这个问题的动机是实现一个接受相对公差参数的 np.unique(arr, True) 版本。因此，np.unique 的第二个参数很重要：我需要知道原始数组中唯一元素（第一次出现！）的索引。因此，元素排序并不重要。

我知道 questions and solutions on np.unique with tolerance。但是，我还没有找到一个解决方案，它也 returns 原始数组中第一次出现的唯一元素的索引。此外，我看到的解决方案是基于排序的，运行s in O(arr.size log(arr.size))。但是，使用哈希映射可以实现恒定时间的解决方案。

想法是将arr中的每个元素上下舍入，并将这些元素放入哈希映射中。如果任一值已在哈希映射中，则忽略条目。否则，该元素包含在结果中。由于hash map的插入和查找运行平均时间恒定，理论上这种方法应该比基于排序的方法更快。

在下面找到我的 Cython 实现：

import numpy as np
cimport numpy as np
import cython
from libcpp.unordered_map cimport unordered_map

@cython.boundscheck(False)
@cython.wraparound(False)
def unique_tol(np.ndarray[DOUBLE_t, ndim=1] lower,
               np.ndarray[DOUBLE_t, ndim=1] higher):
    cdef long i, count
    cdef long endIndex = lower.size
    cdef unordered_map[double, short] vals = unordered_map[double, short]()
    cdef np.ndarray[DOUBLE_t, ndim=1] result_vals = np.empty_like(lower)
    cdef np.ndarray[INT_t, ndim=1] result_indices = np.empty_like(lower, 
                                                                  dtype=int)

    count = 0
    for i in range(endIndex): 
        if not vals.count(lower[i]) and not vals.count(higher[i]):

            # insert in result
            result_vals[count] = lower[i]
            result_indices[count] = i

            # put lowerVal and higherVal in the hashMap
            vals[lower[i]]
            vals[higher[i]]

            # update the index in the result
            count += 1

    return result_vals[:count], result_indices[:count]

通过适当舍入调用的这个方法可以完成这项工作。例如，如果小于 10^-6 的差异应该被忽略，我们会写

unique_tol(np.round(a, 6), np.round(a+1e-6, 6))

现在我想用基于尾数操作的相对舍入程序替换 np.round。我知道 alternative ways of relative rounding，但我认为直接操作尾数应该更高效和优雅。（诚然，我认为性能提升并不显着。但我会对解决方案感兴趣。）

编辑

Warren Weckesser 的解决方案非常有效。然而，结果并不像我希望的那样适用，因为两个相差很小的数字可能有不同的指数。统一尾数将不会导致相似的数字。我想我必须坚持现有的相对舍入解决方案。

Answer 1

"I am looking for an efficient way to set the n least significant entries of each element of arr to 0 or to 1."

您可以创建数据类型为 numpy.uint64 的数组视图，然后根据需要操作该视图中的位。

比如我把这个数组的尾数最低21位设为0

In [46]: np.set_printoptions(precision=15)                                                            

In [47]: x = np.array([0.0, -1/3, 1/5, -1/7, np.pi, 6.02214076e23])                                   

In [48]: x                                                                                            
Out[48]: 
array([ 0.000000000000000e+00, -3.333333333333333e-01,
        2.000000000000000e-01, -1.428571428571428e-01,
        3.141592653589793e+00,  6.022140760000000e+23])

创建 x 中数据类型为 numpy.uint64 的数据视图：

In [49]: u = x.view(np.uint64)

看一下值的二进制表示。

In [50]: [np.binary_repr(t, width=64) for t in u]                                                     
Out[50]: 
['0000000000000000000000000000000000000000000000000000000000000000',
 '1011111111010101010101010101010101010101010101010101010101010101',
 '0011111111001001100110011001100110011001100110011001100110011010',
 '1011111111000010010010010010010010010010010010010010010010010010',
 '0100000000001001001000011111101101010100010001000010110100011000',
 '0100010011011111111000011000010111001010010101111100010100010111']

把低n位设为0，再看一下

In [51]: n = 21                                                                                       

In [52]: u &= ~np.uint64(2**n-1)                                                              

In [53]: [np.binary_repr(t, width=64) for t in u]                                                     
Out[53]: 
['0000000000000000000000000000000000000000000000000000000000000000',
 '1011111111010101010101010101010101010101010000000000000000000000',
 '0011111111001001100110011001100110011001100000000000000000000000',
 '1011111111000010010010010010010010010010010000000000000000000000',
 '0100000000001001001000011111101101010100010000000000000000000000',
 '0100010011011111111000011000010111001010010000000000000000000000']

因为u是与x中相同数据的视图，所以x也被修改了in-place。

In [54]: x                                                                      
Out[54]: 
array([ 0.000000000000000e+00, -3.333333332557231e-01,
        1.999999999534339e-01, -1.428571428405121e-01,
        3.141592653468251e+00,  6.022140758954589e+23])

Answer 2

与@WarrenWeckesser 类似，但没有黑魔法，而是使用 "official" ufuncs。缺点：我很确定它比较慢，很可能是这样：

>>> a = np.random.normal(size=10)**5
>>> a
array([ 9.87664561e-12, -1.79654870e-03,  4.36740261e-01,  7.49256141e+00,
       -8.76894617e-01,  2.93850753e+00, -1.44149959e-02, -1.03026094e-03,
        3.18390143e-03,  3.05521581e-03])
>>> 
>>> mant,expn = np.frexp(a)
>>> mant
array([ 0.67871792, -0.91983293,  0.87348052,  0.93657018, -0.87689462,
        0.73462688, -0.92255974, -0.5274936 ,  0.81507877,  0.78213525])
>>> expn
array([-36,  -9,  -1,   3,   0,   2,  -6,  -9,  -8,  -8], dtype=int32)
>>> a_binned = np.ldexp(np.round(mant,5),expn)
>>> a_binned
array([ 9.87667590e-12, -1.79654297e-03,  4.36740000e-01,  7.49256000e+00,
       -8.76890000e-01,  2.93852000e+00, -1.44150000e-02, -1.03025391e-03,
        3.18390625e-03,  3.05523437e-03])

Numpy：在双数组中设置尾数的最后 n 个元素

Numpy: Setting n last elements of mantissas in double array

python

floating-point

numpy

rounding