switch to immutable arrays for speed

JuliaAstro · Jan 6, 2015 · 9e55067 · 9e55067
1 parent eae6587
commit 9e55067
Show file tree

Hide file tree

Showing 4 changed files with 89 additions and 46 deletions.
diff --git a/README.md b/README.md
@@ -49,40 +49,42 @@ ICRS(0.0,-8.673617379884037e-19)  # pretty close to round-tripping
 
 ## Speed
 
-For small numbers of points, this can be much faster than
-`astropy.coordinates`. Here is an example transforming 100 points
-from ICRS to Galactic:
+For small numbers of points, this package can be much faster than
+`astropy.coordinates` in Python. The following is an example of
+transforming points from ICRS to Galactic. The speed difference is
+factor of approximately 10,000 for 10 coordinates, 80 for 1000
+coordinates and 4 for 100,000 coordinates.
 
 **astropy.coordinates**
 
 ```python
-In [1]: from astropy.coordinates import ICRS
+In [1]: from numpy import pi
 
-In [2]: import numpy as np
+In [2]: from numpy.random import rand
 
-In [3]: ra = 2.* np.pi * np.random.rand(100)
+In [3]: from astropy.coordinates import SkyCoord
 
-In [4]: dec = np.pi * np.random.rand(100) - np.pi/2.
-
-In [5]: c1 = ICRS(ra, dec, unit=('rad', 'rad'))
-
-In [6]: %timeit c1.galactic
-100 loops, best of 3: 6.68 ms per loop
+In [4]: for n in [10, 1000, 100000]:
+   ...:     c = SkyCoord(2.*pi*rand(n), pi*rand(n)-pi/2, unit=('rad', 'rad'))
+   ...:     %timeit c.galactic
+   ...: 
+100 loops, best of 3: 12.5 ms per loop
+100 loops, best of 3: 13 ms per loop
+10 loops, best of 3: 66.7 ms per loop
 ```
 
-**SkyCoords**
+**SkyCoords.jl**
 
-```julia
+```jlcon
 julia> using SkyCoords
 
 julia> using TimeIt
 
-julia> ra = 2pi*rand(100)
-
-julia> dec = pi*rand(100) - pi/2.
-
-julia> c1 = [ICRS(ra[i], dec[i]) for i=1:100]
-
-julia> @timeit to_galactic(c1)
-10000 loops, best of 3: 61.16 µs per loop
-```
+julia> for n in [10, 1000, 100000]
+           c = [ICRS(2pi*rand(), pi*(rand() - 0.5)) for i=1:n]
+           @timeit to_galactic(c)
+       end
+100000 loops, best of 3: 1.33 µs per loop
+1000 loops, best of 3: 163.94 µs per loop
+10 loops, best of 3: 16.23 ms per loop
+```
diff --git a/bench.jl b/bench.jl
@@ -1,6 +1,9 @@
+#!/usr/bin/env julia
 using SkyCoords
 using TimeIt
 
-c1 = [ICRS(2pi*rand(), pi*(rand() - 0.5)) for i=1:100]
-
-@timeit to_galactic(c1)
+for n in [10, 1000, 100000]
+    print("$n coordinates: ")
+    c = [ICRS(2pi*rand(), pi*(rand() - 0.5)) for i=1:n]
+    @timeit to_galactic(c)
+end
diff --git a/src/SkyCoords.jl b/src/SkyCoords.jl
@@ -11,36 +11,77 @@ export ICRS,
 
 # Helper functions -----------------------------------------------------------
 
+# Use an immutable array to avoid memory allocations all over the place.
+immutable Matrix33
+    a11::Float64
+    a12::Float64
+    a13::Float64
+    a21::Float64
+    a22::Float64
+    a23::Float64
+    a31::Float64
+    a32::Float64
+    a33::Float64
+end
+
+immutable Vector3
+    x::Float64
+    y::Float64
+    z::Float64
+end
+
+function *(x::Matrix33, y::Matrix33)
+    Matrix33(x.a11 * y.a11 + x.a12 * y.a21 + x.a13 * y.a31,
+             x.a11 * y.a12 + x.a12 * y.a22 + x.a13 * y.a32,
+             x.a11 * y.a13 + x.a12 * y.a23 + x.a13 * y.a33,
+             x.a21 * y.a11 + x.a22 * y.a21 + x.a23 * y.a31,
+             x.a21 * y.a12 + x.a22 * y.a22 + x.a23 * y.a32,
+             x.a21 * y.a13 + x.a22 * y.a23 + x.a23 * y.a33,
+             x.a31 * y.a11 + x.a32 * y.a21 + x.a33 * y.a31,
+             x.a31 * y.a12 + x.a32 * y.a22 + x.a33 * y.a32,
+             x.a31 * y.a13 + x.a32 * y.a23 + x.a33 * y.a33)
+end
+
+transpose(m::Matrix33) = Matrix33(m.a11, m.a21, m.a31,
+                                  m.a12, m.a22, m.a32,
+                                  m.a13, m.a23, m.a33)
+
+function *(m::Matrix33, v::Vector3)
+    Vector3(m.a11 * v.x + m.a12 * v.y + m.a13 * v.z,
+            m.a21 * v.x + m.a22 * v.y + m.a23 * v.z,
+            m.a31 * v.x + m.a32 * v.y + m.a33 * v.z)
+end
+
 # Create rotation matrix about a given axis (x, y, z)
 function xrotmat(angle)
     s = sin(angle)
     c = cos(angle)
-    [1.  0.  0.;
-     0.  c   s ;
-     0. -s   c ]
+    Matrix33(1., 0., 0.,
+             0.,  c,  s,
+             0., -s,  c)
 end
 
 function yrotmat(angle)
     s = sin(angle)
     c = cos(angle)
-    [c   0. -s ;
-     0.  1.  0.;
-     s   0.  c ]
+    Matrix33(c,  0., -s,
+             0., 1., 0.,
+             s,  0.,  c)
 end
 
 function zrotmat(angle)
     s = sin(angle)
     c = cos(angle)
-    [ c   s   0.;
-     -s   c   0.;
-      0.  0.  1.]
+    Matrix33(c,   s,  0.,
+             -s,  c,  0.,
+             0., 0.,  1.)
 end
 
 # (lon, lat) -> [x, y, z] unit vector
-coords2cart(lon, lat) = [cos(lat)*cos(lon); cos(lat)*sin(lon); sin(lat)]
+coords2cart(lon, lat) = Vector3(cos(lat)*cos(lon), cos(lat)*sin(lon), sin(lat))
 
 # [x, y, z] unit vector -> (lon, lat)
-cart2coords(r) = atan2(r[2], r[1]), atan2(r[3], sqrt(r[1]*r[1] + r[2]*r[2]))
+cart2coords(v) = atan2(v.y, v.x), atan2(v.z, sqrt(v.x*v.x + v.y*v.y))
 
 # Computes the precession matrix from J2000 to the given Julian equinox.
 # Expression from from Capitaine et al. 2003 as expressed in the USNO
@@ -72,10 +113,10 @@ end
 # Constants --------------------------------------------------------------
 
 # Delete once 0.2 is no longer supported:
-if !isdefined(:rad2deg)
-  const rad2deg = radians2degrees
-  const deg2rad = degrees2radians
-end
+#if !isdefined(:rad2deg)
+#  const rad2deg = radians2degrees
+#  const deg2rad = degrees2radians
+#end
 
 # ICRS --> FK5 at J2000 (See USNO Circular 179, section 3.5)
 eta0 = deg2rad(-19.9 / 3600000.)

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -1,10 +1,7 @@
-include("SkyCoords.jl")
+#!/usr/bin/env julia
 using SkyCoords
 using Base.Test
 
-rad2deg = radians2degrees
-deg2rad = degrees2radians
-
 # Angular separation between two points (angles in radians)
 #
 #   The angular separation is calculated using the Vincenty formula [1]_,
@@ -35,7 +32,7 @@ tol = deg2rad(0.03 / 3600.)
 # The testdata file was generated with the ref_icrs_fk5.py script in
 # astropy.coordinates.test.accuracy. The reference values were computed
 # using AST.
-data, hdr = readcsv("testdata/icrs_fk5.csv", has_header=true)
+data, hdr = readcsv("testdata/icrs_fk5.csv", header=true)
 
 for i=1:size(data, 1)
     equinox = float(data[i,1][2:end])