Computation on large array in python

user_7909061 0
Viewed 483

i have a question concerning computation of large array in python( >=3.6.8 ) on Ubuntu SO. I have a large array A, about 20 000 000 (or more) of points (each point is a 3d coordinate x,y,z). I would like to know which of these points are close to a given plane. I use sympy python library for representing point, plane and to compute distance. I have another array T of planes (as sympy object), tipically an array of 7 or 10 elements. For example, for each point p of my A array, and for each t plane, i do something like:

        for p in A:
            for t in T:
                if  t.distance(p) <= Delta:
                    do something

Obviously, if i loop this check for all points and for all plane, i achieve the result after several hours. I done this loop also using processes splitting my array of points in chunks, but also in this way, computation requires a lot of time (hours). I read numpy has where method, but i don't know how to set condition like:

        for t in T:
            B = A[np.where(t.distance(Point3D(A)) <= 0.1)]

where: t is a sympy plane; Point3D is sympy class for casting generic x,y,z point into sympy point; B is array of A indexes which satisfy condition.

where condition as above reported causes exception :"Nonzero coordinates cannot be removed". Is there some fast approach to compute the loop into 10 or 15 minutes (max)? Thank you in advance. Bye.

sympy numpy arrays python
2 Answers

Using map seems to be very fast.

Here a code test example with One Million points in 4 seconds .

in this case map needs a function and two iterables , so we need iterable points (one million) and iterable planes (one million , all the same in this case).

import math
import numpy as np
import time

# Function to find point plane distance
def ppdist(p, plane ):
    x1,y1,z1 = p
    a,b,c,d = plane
    d = abs((a * x1 + b * y1 + c * z1 + d))
    e = (math.sqrt(a * a + b * b + c * c))
    return d/e


# num points
numpts = 1000000
#creating random 3d points
pts = np.random.randint(0,1234,size=(numpts,3))

#plane 
a = 20
b = -20
c = 50
d = 80

plane=[a,b,c,d]
#creating a list of identical planes    
planes = [ plane for n in range( numpts ) ]
# Function call

start_time = time.time()
res_map = list(map( ppdist , pts , planes))
end_time = time.time()
elapsed = end_time-start_time
print(len(res_map), " distances in ", elapsed, "seconds")

output for One Million Point Plane distance computations :

1000000  distances in  4.155193328857422 seconds

I think 20M distance computations needs 4.15*20 -> 83 seconds ~
Try it and let me know.

Here an extended example to abtain a new list with the points that are at a distance less than a value , ** processing 1M points in less than 5 seconds **

import math
import numpy as np
import time

# Function to find point plane distance
def ppdist(p, plane ):
    x1,y1,z1 = p
    a,b,c,d = plane
    d = abs((a * x1 + b * y1 + c * z1 + d))
    e = (math.sqrt(a * a + b * b + c * c))
    return d/e

# Function returns point if distance is less than value
def pdistmin(point , dist , val):
    return point if dist < val else None

# num points
numpts = 1000000
#creating random 3d points
pts = np.random.randint(0,1234,size=(numpts,3))


#plane 
a = 20
b = -20
c = 50
d = 80

plane=[a,b,c,d]
#creating a list of identical planes    
planes = [ plane for n in range( numpts ) ]
# Function call

start_time = time.time()
res_dist = list(map( ppdist , pts , planes))
end_time = time.time()
elapsed = end_time-start_time
print(len(res_dist), " distances in ", elapsed, "seconds")


start_time = time.time()
less = [ 10 for n in range( numpts ) ]
xxx = list(map( pdistmin, pts, res_dist , less ))
end_time = time.time()
elapsed = end_time-start_time
print(len(xxx), " points or None computation in ", elapsed, "seconds")


start_time = time.time()
needed_points = [ p for p in xxx if p is not None ]
end_time = time.time()
elapsed = end_time-start_time
print(len(needed_points), " points in ", elapsed, "seconds")

output:

1000000  distances in  4.1579577922821045 seconds
1000000  points or None computation in  0.42284178733825684 seconds
   9372  points found in  0.014847517013549805 seconds

So less than 5 seconds for 1M points . Probably 100 seconds for 20M Points

user_106386520

SymPy is a symbolic library. The idea is that you would use it to symbolically derive a formula that you then compute using more efficient numerical libraries e.g.:

In [20]: x, y, z = symbols('x, y, z')

In [21]: arbitrary_point = Point3D(x, y, z)

In [22]: delta = 2

In [23]: plane = Plane(Point3D(1, 1, 1), normal_vector=(1,4,7))

In [24]: condition = plane.distance(arbitrary_point) < delta

In [25]: condition
Out[25]: 
│√66⋅(x - 1)   2⋅√66⋅(y - 1)   7⋅√66⋅(z - 1)│    
│─────────── + ───────────── + ─────────────│ < 2
│     66             33              66     │    

In [26]: f = lambdify((x, y, z), condition)

In [27]: f(2, 3, 4)
Out[27]: False

In [28]: f(1, 1, 2)
Out[28]: True

You can see the code used to make the lambdified function like this:

In [29]: import inspect

In [31]: print(inspect.getsource(f))
def _lambdifygenerated(x, y, z):
    return (less(abs((1/66)*sqrt(66)*(x - 1) + (2/33)*sqrt(66)*(y - 1) + (7/66)*sqrt(66)*(z - 1)), 2))

You can now use this function with numpy arrays for the x, y and z values:

In [33]: xvals = np.arange(0, 10)

In [34]: yvals = np.arange(0, 10)

In [35]: zvals = np.arange(0, 10)

In [36]: f(xvals, yvals, zvals)
Out[36]: 
array([ True,  True,  True, False, False, False, False, False, False,
       False])

This will handle large arrays of inputs more efficiently e.g. 20 million points can be processed in 1 second:

In [38]: xvals = yvals = zvals = np.arange(0, 20*10**6)

In [39]: %time ok = f(xvals, yvals, zvals)
CPU times: user 698 ms, sys: 553 ms, total: 1.25 s
Wall time: 1.26 s
user_94509910
Related