45068


numpyで2つの三次元座標点列の引き当てをできるだけ高速に行いたいのですが、行き詰まっています。

＜前提＞  
`numpy`の`round()`と`builtins`の`round()`は挙動が異なります。  
（これ自体は質問の主旨ではありませんが、これについても解説いただけると大変ありがたいです。）

[code]

    >>> print(round(411105.886185, 5))
    411105.88619
    >>> print(np.round(411105.886185, 5))
    411105.88618
    
[/code]

＜本題＞

ある3次元座標点列 `V` を`builtins`の`round()`を使って小数点以下6桁で記録したa.tsvと  
`V` から一部を抽出して`numpy`の`round()`で小数点以下9桁で記録したb.tsvがあります。  
a.tsv、b.tsvともに元々の `V` の順番通りに記録されています。  
b.tsvの各行が、a.tsvでは何行目になるか、を特定したいです。

a.tsv  
977279.482707 734066.643064 662406.439074  
627635.945559 451974.042893 929737.099191  
1025463.393349 752819.302836 885502.793725  
971104.800369 916731.879454 475093.855238  
382780.576043 307121.604863 661611.845153  
...

b.tsv  
977279.482707500 734066.643064236 662406.439073611  
971104.800368500 916731.879453945 475093.855237700  
382780.576043500 307121.604863451 661611.845152960  
...

以下のコードで目的は達成できるのですが、計算量が `log(N**2)` となる点と、メモリ効率が悪い点が納得いきません。

[code]

    import numpy as np
    
    a = np.genfromtxt('a.tsv', delimiter='\t')
    b = np.genfromtxt('b.tsv', delimiter='\t')
    
    D = b-a[:, np.newaxis]
    S = np.sum(D**2, axis=2)
    I = np.argmin(S, axis=0)
    
    print(I)
    
[/code]

最初に `a` を辞書型に格納すれば `O(NlogN)` になると思ったのですが、先述の `round` の問題のため上手くいきませんでした。  
何卒ご教示お願いいたします。


CC BY-SA 4.0
2018-06-23T15:17:19.573
45022
2018-06-26T00:26:44.857
2018-06-23T15:59:23.060
28998
28998
question
python アルゴリズム numpy 浮動小数点数
numpyで2つの三次元座標点列の引き当て


a
を辞書型に格納して検索した場合に上手くいかないのは、bを四捨五入すると最初の四捨五入と両方で切り上げ又は切り捨てになってしまって１回の四捨五入であるaと同じ数字にならないためです。例えば、`0.0000014996`は、aでは`0.000001`になりますが、bでは`0.000001500`になりそれを四捨五入すると`0.000002`になります。

こういうケースが発生するのは、bで小数点7桁以下が`500`になる場合です。そのケースだけ、切り捨てになる場合と切り上げになる場合の両方で辞書を検索してやれば、どちらかで検索が一致するはずです。それで次のような式で計算できます。

[code]

    def func(a, b):
        delta = np.array([-0.0000000001, 0.0000000001])
    
        def func1(p):
            for d0 in delta:
                for d1 in delta:
                    for d2 in delta:
                        t = (np.round(p[0] + d0, 6), np.round(p[1] + d1, 6), np.round(p[2] + d2, 6))
                        if t in dic:
                            return list(t)
    
        dic = {(j[0], j[1], j[2]):i for i, j in enumerate(a)}
        result = np.zeros(b.shape)
        for i in range(len(b)):
            result[i] = func1(b[i])
        return result
    
[/code]

＜前提＞ の部分ですが、numpyのround()とbuiltinsのround()は、どちらも「偶数への丸め」が使われています。

[code]

    >>> print([round(0.5), round(1.5), round(2.5), round(3.5)])
    [0, 2, 2, 4]
    
[/code]

しかし、builtinsのround()を使って次の四捨五入をすると「偶数への丸め」とは違う結果になります。

[code]

    >>> print([round(0.05, 1), round(0.15, 1), round(0.25, 1), round(0.35, 1)])
    [0.1, 0.1, 0.2, 0.3]
    
[/code]

それは、'0.05'は2進数では有理数にならないため、次のような近似値になります。builtinsのround()をその近似値を使って「偶数への丸め」しているものと考えられます。一方numpyのround()は、「偶数への丸め」になるように調整しているので違いがでます。

[code]

    >>> print('{0:.19f}'.format(0.05))
    0.0500000000000000028
    
[/code]

詳しくは、Python チュートリアルの「[15\.
浮動小数点演算、その問題と制限](https://docs.python.jp/3/tutorial/floatingpoint.html)」のページをみてください。なお、「偶数への丸め」は、Wikipediaでは次のように説明されています。

> 偶数への丸め（round to
> even）は、端数が0.5より小さいなら切り捨て、端数が0.5より大きいならは切り上げ、端数がちょうど0.5なら切り捨てと切り上げのうち結果が偶数となる方へ丸める。JIS
> Z 8401で規則Aとして定められていて、規則B（四捨五入）より「望ましい」とされている。


CC BY-SA 4.0
2018-06-26T00:26:44.857
45068
2018-06-26T00:26:44.857
15171
45022
answer