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

Question

问题

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

动机

下面我将提供问题的动机。如果您发现不需要回答上述问题来实现最终目标，我很高兴收到相应的 cmets。然后，我将创建一个单独的问题以保持排序。

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

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

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

在下面找到我的 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

“我正在寻找一种有效的方法来将 arr 的每个元素的 n 个最不重要的条目设置为 0 或 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 [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，但没有使用“官方”ufunc 的黑魔法。缺点：我很确定它会更慢，很可能是这样：

>>> 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])