代码之家  ›  专栏  ›  技术社区  ›  Brendan

从向量的numpy数组中删除重复项(在给定的公差范围内)

  •  5
  • Brendan  · 技术社区  · 16 年前

    我有一个nx5数组,包含n个形式为'id'、'x'、'y'、'z'和'energy'的向量。我需要删除重复的点(即x,y,z都匹配的地方)在0.1的公差范围内。理想情况下,我可以创建一个函数,在其中传递数组、需要匹配的列和匹配的公差。

    跟随 this thread on Scipy-user ,我可以使用记录数组移除基于完整数组的重复项,但只需要匹配数组的一部分。而且,这在一定的公差范围内是不匹配的。

    我可以费力地用一个 for 在python中循环,但是有更好的numponic方法吗?

    3 回复  |  直到 16 年前
        1
  •  2
  •   denis    16 年前

    你可以看看 scipy.spatial.KDTree . n有多大?
    加上:哎呀, tree.query_pairs 不在scipy 0.7.1中。

    当有疑问时,使用暴力:把空间(这边^3)分成几个小单元, 每个单元格一个点:

    """ scatter points to little cells, 1 per cell """
    from __future__ import division         
    import sys                              
    import numpy as np                      
    
    side = 100                              
    npercell = 1  # 1: ~ 1/e empty          
    exec "\n".join( sys.argv[1:] )  # side= ...
    N = side**3 * npercell                  
    print "side: %d  npercell: %d  N: %d" % (side, npercell, N)
    np.random.seed( 1 )                     
    points = np.random.uniform( 0, side, size=(N,3) )
    
    cells = np.zeros( (side,side,side), dtype=np.uint )
    id = 1
    for p in points.astype(int):
        cells[tuple(p)] = id                
        id += 1                             
    
    cells = cells.flatten()
        # A C, an E-flat, and a G walk into a bar. 
        # The bartender says, "Sorry, but we don't serve minors."
    nz = np.nonzero(cells)[0]               
    print "%d cells have points" % len(nz)
    print "first few ids:", cells[nz][:10]
    
        2
  •  0
  •   Vincent    16 年前

    还没有测试过这个,但是如果你把数组按x和y排序,那么z应该会得到重复的列表。然后你需要选择保留哪一个。

    def find_dup_xyz(anarray, x, y, z): #for example in an data = array([id,x,y,z,energy]) x=1 y=2 z=3
        dup_xyz=[]
        for i, row in enumerated(sortedArray):
            nx=1
            while (abs(row[x] - sortedArray[i+nx[x])<0.1) and (abs(row[z] and sortedArray[i+nx[y])<0.1) and (abs(row[z] - sortedArray[i+nx[z])<0.1):
                  nx=+1
                  dup_xyz.append(row)
    return dup_xyz
    

    还发现了这个 http://mail.scipy.org/pipermail/scipy-user/2008-April/016504.html

        3
  •  0
  •   Brendan    16 年前

    我终于有了一个解决方案,我很高兴,这是一个稍微清理剪贴从我自己的代码。可能还有一些虫子。

    注意:它仍然使用“for”循环。我可以利用丹尼斯上面的kdtree思想和舍入来得到完整的解决方案。

    import numpy as np
    
    def remove_duplicates(data, dp_tol=None, cols=None, sort_by=None):
        '''
        Removes duplicate vectors from a list of data points
        Parameters:
            data        An MxN array of N vectors of dimension M 
            cols        An iterable of the columns that must match 
                        in order to constitute a duplicate 
                        (default: [1,2,3] for typical Klist data array) 
            dp_tol      An iterable of three tolerances or a single 
                        tolerance for all dimensions. Uses this to round 
                        the values to specified number of decimal places 
                        before performing the removal. 
                        (default: None)
            sort_by     An iterable of columns to sort by (default: [0])
    
        Returns:
            MxI Array   An array of I vectors (minus the 
                        duplicates)
    
        EXAMPLES:
    
        Remove a duplicate
    
        >>> import wien2k.utils
        >>> import numpy as np
        >>> vecs1 = np.array([[1, 0, 0, 0],
        ...     [2, 0, 0, 0],
        ...     [3, 0, 0, 1]])
        >>> remove_duplicates(vecs1)
        array([[1, 0, 0, 0],
               [3, 0, 0, 1]])
    
        Remove duplicates with a tolerance
    
        >>> vecs2 = np.array([[1, 0, 0, 0  ],
        ...     [2, 0, 0, 0.001 ],
        ...     [3, 0, 0, 0.02  ],
        ...     [4, 0, 0, 1     ]])
        >>> remove_duplicates(vecs2, dp_tol=2)
        array([[ 1.  ,  0.  ,  0.  ,  0.  ],
               [ 3.  ,  0.  ,  0.  ,  0.02],
               [ 4.  ,  0.  ,  0.  ,  1.  ]])
    
        Remove duplicates and sort by k values
    
        >>> vecs3 = np.array([[1, 0, 0, 0],
        ...     [2, 0, 0, 2],
        ...     [3, 0, 0, 0],
        ...     [4, 0, 0, 1]])
        >>> remove_duplicates(vecs3, sort_by=[3])
        array([[1, 0, 0, 0],
               [4, 0, 0, 1],
               [2, 0, 0, 2]])
    
        Change the columns that constitute a duplicate
    
        >>> vecs4 = np.array([[1, 0, 0, 0],
        ...     [2, 0, 0, 2],
        ...     [1, 0, 0, 0],
        ...     [4, 0, 0, 1]])
        >>> remove_duplicates(vecs4, cols=[0])
        array([[1, 0, 0, 0],
               [2, 0, 0, 2],
               [4, 0, 0, 1]])
    
        '''
        # Deal with the parameters
        if sort_by is None:
            sort_by = [0]
        if cols is None:
            cols = [1,2,3]
        if dp_tol is not None:
            # test to see if already an iterable
            try:
                null = iter(dp_tol)
                tols = np.array(dp_tol)
            except TypeError:
                tols = np.ones_like(cols) * dp_tol
            # Convert to numbers of decimal places
            # Find the 'order' of the axes
        else:
            tols = None
    
        rnd_data = data.copy()
        # set the tolerances
        if tols is not None:
            for col,tol in zip(cols, tols):
                rnd_data[:,col] = np.around(rnd_data[:,col], decimals=tol)
    
        # TODO: For now, use a slow Python 'for' loop, try to find a more
        # numponic way later - see: http://stackoverflow.com/questions/2433882/
        sorted_indexes = np.lexsort(tuple([rnd_data[:,col] for col in cols]))
        rnd_data = rnd_data[sorted_indexes]
        unique_kpts = []
        for i in xrange(len(rnd_data)):
            if i == 0:
                unique_kpts.append(i)    
            else:
                if (rnd_data[i, cols] == rnd_data[i-1, cols]).all():
                    continue
                else:
                    unique_kpts.append(i)    
    
        rnd_data =  rnd_data[unique_kpts]
        # Now sort
        sorted_indexes = np.lexsort(tuple([rnd_data[:,col] for col in sort_by]))
        rnd_data = rnd_data[sorted_indexes]
        return rnd_data
    
    
    
    if __name__ == '__main__':
        import doctest
        doctest.testmod()