أسرع numpy الديكارت إلى تحويل الإحداثيات الكروية؟

Question

لدي مجموعة من 3 ملايين نقطة بيانات من مقياس Accellerometer 3-Axiz (XYZ) ، وأريد إضافة 3 أعمدة إلى الصفيف الذي يحتوي على الإحداثيات الكروية المكافئة (R ، Theta ، PHI). يعمل الرمز التالي ، ولكن يبدو بطيئًا جدًا. كيف يمكنني أن أفعل أفضل؟

import numpy as np
import math as m

def cart2sph(x,y,z):
    XsqPlusYsq = x**2 + y**2
    r = m.sqrt(XsqPlusYsq + z**2)               # r
    elev = m.atan2(z,m.sqrt(XsqPlusYsq))     # theta
    az = m.atan2(y,x)                           # phi
    return r, elev, az

def cart2sphA(pts):
    return np.array([cart2sph(x,y,z) for x,y,z in pts])

def appendSpherical(xyz):
    np.hstack((xyz, cart2sphA(xyz)))

Answer 1

هذا مشابه ل جاستن بيلإجابة ، ولكن باستخدام فقط numpy والاستفادة من التقييم المدمج:

import numpy as np

def appendSpherical_np(xyz):
    ptsnew = np.hstack((xyz, np.zeros(xyz.shape)))
    xy = xyz[:,0]**2 + xyz[:,1]**2
    ptsnew[:,3] = np.sqrt(xy + xyz[:,2]**2)
    ptsnew[:,4] = np.arctan2(np.sqrt(xy), xyz[:,2]) # for elevation angle defined from Z-axis down
    #ptsnew[:,4] = np.arctan2(xyz[:,2], np.sqrt(xy)) # for elevation angle defined from XY-plane up
    ptsnew[:,5] = np.arctan2(xyz[:,1], xyz[:,0])
    return ptsnew

لاحظ أنه ، كما هو مقترح في التعليقات ، أنا غيرت تعريف زاوية الارتفاع من وظيفتك الأصلية. على الجهاز الخاص بي ، والاختبار مع pts = np.random.rand(3000000, 3), ، ذهب الوقت من 76 ثانية إلى 3.3 ثانية. ليس لدي Cython ، لذا لم أتمكن من مقارنة التوقيت بهذا الحل.

Answer 2

إليك رمز Cython السريع الذي كتبته لهذا:

cdef extern from "math.h":
    long double sqrt(long double xx)
    long double atan2(long double a, double b)

import numpy as np
cimport numpy as np
cimport cython

ctypedef np.float64_t DTYPE_t

@cython.boundscheck(False)
@cython.wraparound(False)
def appendSpherical(np.ndarray[DTYPE_t,ndim=2] xyz):
    cdef np.ndarray[DTYPE_t,ndim=2] pts = np.empty((xyz.shape[0],6))
    cdef long double XsqPlusYsq
    for i in xrange(xyz.shape[0]):
        pts[i,0] = xyz[i,0]
        pts[i,1] = xyz[i,1]
        pts[i,2] = xyz[i,2]
        XsqPlusYsq = xyz[i,0]**2 + xyz[i,1]**2
        pts[i,3] = sqrt(XsqPlusYsq + xyz[i,2]**2)
        pts[i,4] = atan2(xyz[i,2],sqrt(XsqPlusYsq))
        pts[i,5] = atan2(xyz[i,1],xyz[i,0])
    return pts

استغرق الأمر الوقت لأسفل من 62.4 ثانية إلى 1.22 ثانية باستخدام 3،000،000 نقطة بالنسبة لي. هذا ليس رث جدا. أنا متأكد من أن هناك بعض التحسينات الأخرى التي يمكن إجراؤها.

Answer 3

! لا يزال هناك خطأ في جميع التعليمات البرمجية أعلاه .. وهذه نتيجة جوجل الأعلى .. TLDR: لقد اختبرت هذا باستخدام Vpython ، باستخدام ATAN2 لـ Theta (Elev) خاطئ ، استخدم ACOs! من الصحيح لـ PHI (AZIM). أوصي بوظيفة Sympy1.0 ACOS (لا تشكو حتى من ACOS (z/r) مع r = 0).

http://mathworld.wolfram.com/sphericalcoordinates.html

إذا قمنا بتحويل ذلك إلى نظام الفيزياء (R ، Theta ، Phi) = (R ، Elev ، Azimuth) لدينا:

r = sqrt(x*x + y*y + z*z)
phi = atan2(y,x)
theta = acos(z,r)

غير الأمثل ولكن صحيح رمز لنظام الفيزياء اليمنى:

from sympy import *
def asCartesian(rthetaphi):
    #takes list rthetaphi (single coord)
    r       = rthetaphi[0]
    theta   = rthetaphi[1]* pi/180 # to radian
    phi     = rthetaphi[2]* pi/180
    x = r * sin( theta ) * cos( phi )
    y = r * sin( theta ) * sin( phi )
    z = r * cos( theta )
    return [x,y,z]

def asSpherical(xyz):
    #takes list xyz (single coord)
    x       = xyz[0]
    y       = xyz[1]
    z       = xyz[2]
    r       =  sqrt(x*x + y*y + z*z)
    theta   =  acos(z/r)*180/ pi #to degrees
    phi     =  atan2(y,x)*180/ pi
    return [r,theta,phi]

يمكنك اختبارها بنفسك بوظيفة مثل:

test = asCartesian(asSpherical([-2.13091326,-0.0058279,0.83697319]))

بعض بيانات الاختبار الأخرى لبعض الأرباع:

[[ 0.          0.          0.        ]
 [-2.13091326 -0.0058279   0.83697319]
 [ 1.82172775  1.15959835  1.09232283]
 [ 1.47554111 -0.14483833 -1.80804324]
 [-1.13940573 -1.45129967 -1.30132008]
 [ 0.33530045 -1.47780466  1.6384716 ]
 [-0.51094007  1.80408573 -2.12652707]]

لقد استخدمت Vpython بالإضافة إلى تصور المتجهات:

test   = v.arrow(pos = (0,0,0), axis = vis_ori_ALA , shaftwidth=0.05, color=v.color.red)

Answer 4

لإكمال الإجابات السابقة ، إليك أ numExpr التنفيذ (مع احتمال احتياطي إلى Numpy) ،

import numpy as np
from numpy import arctan2, sqrt
import numexpr as ne

def cart2sph(x,y,z, ceval=ne.evaluate):
    """ x, y, z :  ndarray coordinates
        ceval: backend to use: 
              - eval :  pure Numpy
              - numexpr.evaluate:  Numexpr """
    azimuth = ceval('arctan2(y,x)')
    xy2 = ceval('x**2 + y**2')
    elevation = ceval('arctan2(z, sqrt(xy2))')
    r = eval('sqrt(xy2 + z**2)')
    return azimuth, elevation, r

بالنسبة لأحجام الصفيف الكبيرة ، يتيح ذلك عاملًا من السرعة 2 مقارنةً بتنفيذ نقي نقي ، وسيكون مماثلًا لسرعات C أو Cython. حل numpy الحالي (عند استخدامه مع ceval=eval الحجة) هو أيضا أسرع بنسبة 25 ٪ من appendSpherical_np وظيفة في mtrw إجابة لأحجام الصفيف الكبيرة ،

In [1]: xyz = np.random.rand(3000000,3)
   ...: x,y,z = xyz.T
In [2]: %timeit -n 1 appendSpherical_np(xyz)
1 loops, best of 3: 397 ms per loop
In [3]: %timeit -n 1 cart2sph(x,y,z, ceval=eval)
1 loops, best of 3: 280 ms per loop
In [4]: %timeit -n 1 cart2sph(x,y,z, ceval=ne.evaluate)
1 loops, best of 3: 145 ms per loop

على الرغم من الأحجام الصغيرة ، appendSpherical_np هو في الواقع أسرع ،

In [5]: xyz = np.random.rand(3000,3)
...: x,y,z = xyz.T
In [6]: %timeit -n 1 appendSpherical_np(xyz)
1 loops, best of 3: 206 µs per loop
In [7]: %timeit -n 1 cart2sph(x,y,z, ceval=eval)
1 loops, best of 3: 261 µs per loop
In [8]: %timeit -n 1 cart2sph(x,y,z, ceval=ne.evaluate)
1 loops, best of 3: 271 µs per loop