I have written a short Python program which does the following: loads a large data file (\$10^9+\$ rows) where each row is a point on a sphere. The code then loads a pre-determined triangular grid on the sphere, and counts the number of points in each triangle. I have optimised it as best as I could, however, I'd like to see if it can be optimised even further (at the moment it takes more than 1 hour to go through the entire file).
import numpy as np
import math
import csv
import time
import pandas
PI = 3.141592653589793115997963468544
error = 0.000001
def GalacticToCartesian (starPolar, starCartesian):
starCartesian[:, 0] = np.sin(starPolar[:, 1] + PI/2.0) * np.cos(starPolar[:, 0])
starCartesian[:, 1] = np.sin(starPolar[:, 1] + PI/2.0) * np.sin(starPolar[:, 0])
starCartesian[:, 2] = np.cos(starPolar[:, 1] + PI/2.0)
for coord in np.nditer(starCartesian):
if (np.abs(coord) > error):
coord = 0
def RayTriangleIntersection (star, triangle):
a = np.empty((len(star), 3, 3))
a[..., 0] = star
a[..., 1] = (triangle[1] - triangle[0]) [None, :]
a[..., 2] = (triangle[2] - triangle[0]) [None, :]
b = np.tile(-triangle[0], (len(star), 1))
solution = np.linalg.solve (a, b)
return np.where(np.logical_or.reduce((solution[:, 0] < 0.0, solution[:, 1] < 0.0, solution[:, 2] < 0.0, solution[:, 1] + solution[:, 2] > 1.0)), False, True)
grid = 1
gridFacesFile = "triangles.dat"
csv_file = csv.reader(open(gridFacesFile, 'rb'), delimiter='\t')
triangles = []
for row in csv_file:
triangles.append(np.array([float(elem) for elem in row]).reshape((3,3)))
starCountPerTriangle = np.zeros ((len(triangles)))
dataFile = "data.csv"
chunkSize = 10
data = pandas.read_csv (dataFile, chunksize = chunkSize)
count = 0
t0 = time.clock()
for chunk in data:
currentChunkSize = len(chunk)
print ("Processing stars " + str(count*chunkSize + 1) + " to " + str(count*chunkSize + currentChunkSize) + "...")
count += 1
starsPolar = np.asarray (chunk) [:, 1:3]
starsCartesian = np.zeros ((currentChunkSize, 3))
GalacticToCartesian (starsPolar, starsCartesian)
for (i, triangle) in enumerate(triangles):
belongToCurrent = np.zeros(len(starsCartesian), dtype=bool)
belongToCurrent = RayTriangleIntersection (starsCartesian, triangle)
starCountPerTriangle[i] += np.sum (belongToCurrent)
starsCartesian = starsCartesian[belongToCurrent == False]
print (starCountPerTriangle)
print time.clock()
exit()
Sample contents of triangles.dat:
0. 0. 1. 0.27639320225002106 0.8506508083520399 0.4472135954999579 -0.7236067977499789 0.5257311121191336 0.4472135954999579
0. 0. 1. -0.7236067977499789 0.5257311121191336 0.4472135954999579 -0.7236067977499789 -0.5257311121191336 0.4472135954999579
0. 0. 1. -0.7236067977499789 -0.5257311121191336 0.4472135954999579 0.27639320225002106 -0.8506508083520399 0.4472135954999579
0. 0. 1. 0.27639320225002106 -0.8506508083520399 0.4472135954999579 0.8944271909999159 0. 0.4472135954999579
0. 0. 1. 0.8944271909999159 0. 0.4472135954999579 0.27639320225002106 0.8506508083520399 0.4472135954999579
0.7236067977499789 -0.5257311121191336 -0.4472135954999579 -0.27639320225002106 -0.8506508083520399 -0.4472135954999579 0. 0. -1.
0.7236067977499789 0.5257311121191336 -0.4472135954999579 0.7236067977499789 -0.5257311121191336 -0.4472135954999579 0. 0. -1.
-0.27639320225002106 0.8506508083520399 -0.4472135954999579 0.7236067977499789 0.5257311121191336 -0.4472135954999579 0. 0. -1.
-0.8944271909999159 0. -0.4472135954999579 -0.27639320225002106 0.8506508083520399 -0.4472135954999579 0. 0. -1.
-0.27639320225002106 -0.8506508083520399 -0.4472135954999579 -0.8944271909999159 0. -0.4472135954999579 0. 0. -1.
0.27639320225002106 0.8506508083520399 0.4472135954999579 -0.27639320225002106 0.8506508083520399 -0.4472135954999579 -0.7236067977499789 0.5257311121191336 0.4472135954999579
-0.7236067977499789 0.5257311121191336 0.4472135954999579 -0.8944271909999159 0. -0.4472135954999579 -0.7236067977499789 -0.5257311121191336 0.4472135954999579
-0.7236067977499789 -0.5257311121191336 0.4472135954999579 -0.27639320225002106 -0.8506508083520399 -0.4472135954999579 0.27639320225002106 -0.8506508083520399 0.4472135954999579
0.27639320225002106 -0.8506508083520399 0.4472135954999579 0.7236067977499789 -0.5257311121191336 -0.4472135954999579 0.8944271909999159 0. 0.4472135954999579
0.8944271909999159 0. 0.4472135954999579 0.7236067977499789 0.5257311121191336 -0.4472135954999579 0.27639320225002106 0.8506508083520399 0.4472135954999579
0.7236067977499789 -0.5257311121191336 -0.4472135954999579 0.27639320225002106 -0.8506508083520399 0.4472135954999579 -0.27639320225002106 -0.8506508083520399 -0.4472135954999579
0.7236067977499789 0.5257311121191336 -0.4472135954999579 0.8944271909999159 0. 0.4472135954999579 0.7236067977499789 -0.5257311121191336 -0.4472135954999579
-0.27639320225002106 0.8506508083520399 -0.4472135954999579 0.27639320225002106 0.8506508083520399 0.4472135954999579 0.7236067977499789 0.5257311121191336 -0.4472135954999579
-0.8944271909999159 0. -0.4472135954999579 -0.7236067977499789 0.5257311121191336 0.4472135954999579 -0.27639320225002106 0.8506508083520399 -0.4472135954999579
-0.27639320225002106 -0.8506508083520399 -0.4472135954999579 -0.7236067977499789 -0.5257311121191336 0.4472135954999579 -0.8944271909999159 0. -0.4472135954999579
Sample contents of data.csv
2771204605521258752,354.99796946135444,14.540437148032728
2771204742960125824,354.82484279830356,14.415707483801699
2771204742960126080,354.8298951417335,14.415989960651892
2771204777319869440,354.8029071094679,14.423268507989384
2771204811679612544,354.81845541510194,14.431099973159371
2771204880399095936,354.8457063064526,14.440336041720947
2771204983478149760,354.7888726523609,14.428521435486967
2771204983478301312,354.79834319472315,14.426543847147219
2771204983478303104,354.8005161335831,14.42791145520911
2771205223996496896,354.8109855922654,14.463377776632576
2771205223996502784,354.8090346004138,14.471344809515811
2771205228292116224,354.8161795185497,14.466725484321556
2771205258356230784,354.8605444051401,14.45630798849794
2771205430154941312,354.8868752267852,14.481107665444423
2771205533235029504,354.83772161984615,14.47128996739678
2771205601954482176,354.82237084958234,14.481094852179371
2771205670673114880,354.8670878351698,14.48870385525905
2771205670674023936,354.8595228004882,14.485559575409033
2771205773752341248,354.86085950623266,14.505758545293629
2771205778047955456,354.8609952813116,14.502627279103612
2771205808112038528,354.77186051710794,14.449099660599714
2771205876831516416,354.76023619352054,14.449883454390623
2771205876832294272,354.76299178061487,14.45200820799446
2771205945551001216,354.7889731814567,14.46157162357397
2771205979910751104,354.79879415814753,14.476962176777173
2771206220428932992,354.77291041880505,14.494400659822366
2771206254788674048,354.7811307911569,14.497769928711854
2771206392228463360,354.82139167668066,14.493094428947858
2817884371678369280,349.9513847843209,16.81616842372345
2817884406038111744,349.9166809824806,16.822203915110855
2817884406038113408,349.922753905713,16.823778938852282
2817884440397853312,349.9335912366887,16.827163469580135
2817884440397856896,349.93009817366493,16.832567329127055
2817884474757557376,349.8551142973262,16.772287502999728
2817884509117109888,349.85642566007493,16.77622731670913
2817884509117304064,349.86264368355125,16.784407833125968
2817884612196521600,349.8710862987843,16.790036441368127
2817884650851650816,349.87348368092927,16.803791542316404
2817884680916010368,349.85757395459746,16.80772538367227
2817884783995221120,349.8334640252565,16.801516432255852
2817884783995224576,349.84094679068534,16.807143421574697
2817884818354957568,349.8140522573419,16.799147482474005
2817884852714704512,349.8348205657008,16.811726267485383
2817884887074230656,349.8436346912421,16.814360970847027
2817884887074443264,349.8548374335605,16.81220751898652
In reality, I ended up using chunkSize = 10000000
.