Merge pull request #6 from sisl/TW_julia07

Julia 0.7 Compatability

.travis.yml

@@ -4,7 +4,7 @@ os:
   - linux
   # - osx
 julia:
-  - 0.6
+  - 0.7
   - nightly
 matrix:
   allow_failures:
@@ -13,6 +13,6 @@ notifications:
   email: false
 script:
   - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
-  - julia --check-bounds=yes -e 'Pkg.clone(pwd()); Pkg.test("Vec"; coverage=true)'
+  - julia --check-bounds=yes -e 'using Pkg; Pkg.add(pwd()); Pkg.test("Vec"; coverage=true)'
 after_success:
-  - julia -e 'cd(Pkg.dir("Vec")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'
+  - julia -e 'using Pkg; cd(Pkg.dir("Vec")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'

REQUIRE

@@ -1,2 +1,2 @@
-julia 0.6
-StaticArrays
+julia 0.7
+StaticArrays

src/Vec.jl

@@ -3,6 +3,9 @@ __precompile__()
 module Vec
 
 using StaticArrays
+import LinearAlgebra
+import Printf: @printf
+import LinearAlgebra: ⋅, ×
 
 export
     AbstractVec,

src/coordinate_transforms.jl

@@ -131,13 +131,12 @@ end
 Base.convert(::Type{VecE3}, p::ECEF) = VecE3(p.x, p.y, p.z)
 Base.convert(::Type{ECEF}, p::VecE3) = ECEF(p.x, p.y, p.z)
 
-"""
-Universal Transverse Mercator coordinate system
-"""
-
 const UTM_LATITUDE_LIMIT_NORTH = deg2rad(84)
 const UTM_LATITUDE_LIMIT_SOUTH = deg2rad(-80)
 
+"""
+Universal Transverse Mercator coordinate system
+"""
 struct UTM <: AbstractCoordinate
     e::Float64 # [m]
     n::Float64 # [m]
@@ -264,25 +263,25 @@ function Base.convert(::Type{UTM}, lla::LatLonAlt, datum::GeodeticDatum=WGS_84,
         error("latitude $(rad2deg(lat)) is out of limits!")
     end
 
-    const FN = 0.0      # false northing, zero in the northern hemisphere
-    const FE = 500000.0 # false easting
-    const ko = 0.9996   # central scale factor
+    FN = 0.0      # false northing, zero in the northern hemisphere
+    FE = 500000.0 # false easting
+    ko = 0.9996   # central scale factor
 
-    zone_centers = -177.0*pi/180 + 6.0*pi/180*collect(0:59) # longitudes of the zone centers
+    zone_centers = collect(0:59).*(6.0*π/180) .- (177.0*π/180) # longitudes of the zone centers
     if zone == -1
-        zone = indmin(map(x->abs(lon - x), zone_centers)) # index of min zone center
+        zone = argmin(map(x->abs(lon - x), zone_centers)) # index of min zone center
     end
     long0 = zone_centers[zone]
 
-    s  = sin(lat)
-    c  = cos(lat)
-    t  = tan(lat)
+    s = sin(lat)
+    c = cos(lat)
+    t = tan(lat)
 
     a = datum.a
     b = datum.b
 
     n  = (a-b)/(a+b)
-    e₁  = sqrt((a^2-b^2)/a^2) # first eccentricity
+    e₁ = sqrt((a^2-b^2)/a^2) # first eccentricity
     ep = sqrt((a^2-b^2)/b^2) # second eccentricity
     nu = a/(1-e₁^2*s^2)^0.5   # radius of curvature in the prime vertical
     dh = lon - long0         # longitudinal distance from central meridian

src/geom/hyperplanes.jl

@@ -22,39 +22,39 @@ struct Plane3
     Constructs a plane from its normal n and distance to the origin d
     such that the algebraic equation of the plane is n⋅x + d = 0.
     """
-    Plane3(normal::VecE3, offset::Real) = new(normalize(normal), convert(Float64, offset))
+    Plane3(normal::VecE3, offset::Real) = new(LinearAlgebra.normalize(normal), convert(Float64, offset))
 
 
     """
     Construct a plane from its normal and a point on the plane.
     """
     function Plane3(normal::VecE3, P::VecE3)
-        n = normalize(normal)
+        n = LinearAlgebra.normalize(normal)
         offset = -n⋅P
         return new(n, offset)
     end
 end
 
 """
-Constructs a plane passing through the three points. 
+Constructs a plane passing through the three points.
 """
 function Plane3(p0::VecE3, p1::VecE3, p2::VecE3)
     v0 = p2 - p0
     v1 = p1 - p0
 
     normal = v0×v1
-    nnorm = norm(normal)
-    if nnorm <= norm(v0)*norm(v1)*eps()
+    nnorm = LinearAlgebra.norm(normal)
+    if nnorm <= LinearAlgebra.norm(v0)*LinearAlgebra.norm(v1)*eps()
         M = hcat(convert(Vector{Float64}, v0),
                  convert(Vector{Float64}, v1))'
-        s = svdfact(M, thin=false)
+        s = LinearAlgebra.svdfact(M, thin=false)
         normal = convert(VecE3, s[:V][:,2])
     else
-        normal = normalize(normal)
+        normal = LinearAlgebra.normalize(normal)
     end
-    
+
     offset = -p0⋅normal
-    
+
     return Plane3(normal, offset)
 end
 

src/geom/line_segments.jl

@@ -26,14 +26,14 @@ function get_distance(seg::LineSegment, P::VecE2)
         return 0.0
     end
 
-    r = dot(ab, pb)/denom
+    r = (ab⋅pb)/denom
 
     if r ≤ 0.0
-        norm(P - seg.A)
+        LinearAlgebra.norm(P - seg.A)
     elseif r ≥ 1.0
-        norm(P - seg.B)
+        LinearAlgebra.norm(P - seg.B)
     else
-        norm(P - (seg.A + r*ab))
+        LinearAlgebra.norm(P - (seg.A + r*ab))
     end
 end
 
@@ -49,11 +49,11 @@ The angular distance between the two line segments
 function angledist(segA::LineSegment, segB::LineSegment)
     u = segA.B - segA.A
     v = segB.B - segB.A
-    sqdenom = dot(u,u)*dot(v,v)
+    sqdenom = (u⋅u)*(v⋅v)
     if isapprox(sqdenom, 0.0, atol=1e-10)
         return NaN
     end
-    return acos(dot(u,v) / sqrt(sqdenom))
+    return acos((u⋅v) / sqrt(sqdenom))
 end
 
 """

src/geom/lines.jl

@@ -32,8 +32,8 @@ function get_distance(line::Line, P::VecE2)
         return 0.0
     end
 
-    r = dot(ab, pb)/denom
-    return norm(P - (line.C + r*ab))
+    r = (ab⋅pb)/denom
+    return LinearAlgebra.norm(P - (line.C + r*ab))
 end
 
 

src/geom/projectiles.jl

@@ -26,13 +26,13 @@ function closest_time_of_approach(A::Projectile, B::Projectile)
     if aΔ ≈ 0.0
         return 0.0
     else
-        return -dot(W, Δ) / aΔ
+        return -(W⋅Δ) / aΔ
     end
 end
 function closest_approach_distance(A::Projectile, B::Projectile, t_CPA::Float64=closest_time_of_approach(A, B))
     Aₜ = convert(VecE2, propagate(A, t_CPA).pos)
     Bₜ = convert(VecE2, propagate(B, t_CPA).pos)
-    return norm(Aₜ - Bₜ)
+    return LinearAlgebra.norm(Aₜ - Bₜ)
 end
 function closest_time_of_approach_and_distance(A::Projectile, B::Projectile)
     t_CPA = closest_time_of_approach(A, B)
@@ -47,11 +47,11 @@ function get_intersection_time(A::Projectile, seg::LineSegment)
     v₂ = seg.B - seg.A
     v₃ = polar(1.0, A.pos.θ + π/2)
 
-    denom = dot(v₂, v₃)
+    denom = (v₂⋅v₃)
 
     if !isapprox(denom, 0.0, atol=1e-10)
-        d₁ = cross(v₂, v₁) / denom # time for projectile (0 ≤ t₁)
-        t₂ = dot(v₁, v₃) / denom # time for segment (0 ≤ t₂ ≤ 1)
+        d₁ = (v₂×v₁) / denom # time for projectile (0 ≤ t₁)
+        t₂ = (v₁⋅v₃) / denom # time for segment (0 ≤ t₂ ≤ 1)
         if 0.0 ≤ d₁ && 0.0 ≤ t₂ ≤ 1.0
             return d₁/A.v
         end
@@ -74,12 +74,12 @@ function closest_time_of_approach_and_distance(P::Projectile, seg::LineSegment,
     v₂ = seg.B - seg.A
     v₃ = polar(1.0, P.pos.θ + π/2)
 
-    denom = dot(v₂, v₃)
+    denom = (v₂⋅v₃)
 
     if !isapprox(denom, 0.0, atol=1e-10)
 
-        d₁ = cross(v₂, v₁) / denom # time for projectile (0 ≤ d₁)
-        t₂ = dot(v₁, v₃) / denom # time for segment (0 ≤ t₂ ≤ 1)
+        d₁ = (v₂×v₁) / denom # time for projectile (0 ≤ d₁)
+        t₂ = (v₁⋅v₃) / denom # time for segment (0 ≤ t₂ ≤ 1)
 
         if 0.0 ≤ d₁ && 0.0 ≤ t₂ ≤ 1.0
             t = d₁/P.v
@@ -90,8 +90,8 @@ function closest_time_of_approach_and_distance(P::Projectile, seg::LineSegment,
                 r = polar(1.0, P.pos.θ)
                 projA = proj(seg.A - o, r, VecE2)
                 projB = proj(seg.B - o, r, VecE2)
-                tA = max(norm(projA)/P.v * sign(dot(r, projA)), 0.0)
-                tB = max(norm(projB)/P.v * sign(dot(r, projB)), 0.0)
+                tA = max(LinearAlgebra.norm(projA)/P.v * sign(r⋅projA), 0.0)
+                tB = max(LinearAlgebra.norm(projB)/P.v * sign(r⋅projB), 0.0)
                 pA = VecE2(propagate(P, tA).pos)
                 pB = VecE2(propagate(P, tB).pos)
                 distA = normsquared(seg.A - pA)
@@ -121,8 +121,8 @@ function closest_time_of_approach_and_distance(P::Projectile, seg::LineSegment,
                 r = polar(1.0, P.pos.θ)
                 projA = proj(seg.A - o, r, VecE2)
                 projB = proj(seg.B - o, r, VecE2)
-                tA = max(norm(projA)/P.v * sign(dot(r, projA)), 0.0)
-                tB = max(norm(projB)/P.v * sign(dot(r, projB)), 0.0)
+                tA = max(LinearAlgebra.norm(projA)/P.v * sign(r⋅projA), 0.0)
+                tB = max(LinearAlgebra.norm(projB)/P.v * sign(r⋅projB), 0.0)
                 pA = VecE2(propagate(P, tA).pos)
                 pB = VecE2(propagate(P, tB).pos)
                 distA = normsquared(seg.A - pA)

src/geom/rays.jl

@@ -23,18 +23,17 @@ function intersects(A::Ray, B::Ray)
     return u > 0.0 && v > 0.0
 end
 function intersects(ray::Ray, line::Line;
-    ε::Float64 = 1e-10;
+    ε::Float64 = 1e-10,
     )
 
     v₁ = VecE2(ray) - line.C
     v₂ = polar(1.0, line.θ)
     v₃ = polar(1.0, ray.θ + π/2)
 
-    denom = dot(v₂, v₃)
+    denom = v₂⋅v₃
 
     if !(denom ≈ 0.0)
-        t₁ = cross(v₂, v₁) / denom # time for ray (0 ≤ t₁)
-        # t₂ = dot(v₁, v₃) / denom # time for line can be anything
+        t₁ = (v₂×v₁) / denom # time for ray (0 ≤ t₁)
         return -ε ≤ t₁
     else
         # denom is zero if the line and the ray are parallel
@@ -48,19 +47,19 @@ function intersects(ray::Ray, seg::LineSegment)
     v₂ = seg.B - seg.A
     v₃ = polar(1.0, ray.θ + π/2)
 
-    denom = dot(v₂, v₃)
+    denom = v₂⋅v₃
 
     if !isapprox(denom, 0.0, atol=1e-10)
-        t₁ = cross(v₂, v₁) / denom # time for ray (0 ≤ t₁)
-        t₂ = dot(v₁, v₃) / denom # time for segment (0 ≤ t₂ ≤ 1)
+        t₁ = (v₂×v₁) / denom # time for ray (0 ≤ t₁)
+        t₂ = (v₁⋅v₃) / denom # time for segment (0 ≤ t₂ ≤ 1)
         return 0 ≤ t₁ && 0 ≤ t₂ ≤ 1
     else
         # denom is zero if the segment and the ray are parallel
         # only collide if they are perfectly aligned
         # must ensure that at least one point is in the positive ray direction
         r = polar(1.0, ray.θ)
         return are_collinear(R, seg.A, seg.B) &&
-               (dot(r, seg.A - R) ≥ 0 || dot(r, seg.B - R) ≥ 0)
+               (r⋅(seg.A - R) ≥ 0 || r⋅(seg.B - R) ≥ 0)
     end
 end
 
@@ -95,11 +94,11 @@ function Base.intersect(ray::Ray, seg::LineSegment)
     v₂ = seg.B - seg.A
     v₃ = polar(1.0, ray.θ + π/2)
 
-    denom = dot(v₂, v₃)
+    denom = v₂⋅v₃
 
     if !isapprox(denom, 0.0, atol=1e-10)
-        t₁ = cross(v₂, v₁) / denom # time for ray (0 ≤ t₁)
-        t₂ = dot(v₁, v₃) / denom # time for segment (0 ≤ t₂ ≤ 1)
+        t₁ = (v₂×v₁) / denom # time for ray (0 ≤ t₁)
+        t₂ = (v₁⋅v₃) / denom # time for segment (0 ≤ t₂ ≤ 1)
         if 0 ≤ t₁ && 0 ≤ t₂ ≤ 1
             return R + polar(t₁, ray.θ)
         end
@@ -109,7 +108,7 @@ function Base.intersect(ray::Ray, seg::LineSegment)
         # must ensure that at least one point is in the positive ray direction
         r = polar(1.0, ray.θ)
         if are_collinear(R, seg.A, seg.B) &&
-               (dot(r, seg.A - R) ≥ 0 || dot(r, seg.B - R) ≥ 0)
+               (r⋅(seg.A - R) ≥ 0 || r⋅(seg.B - R) ≥ 0)
             return R
         end
     end

src/quat.jl

@@ -73,16 +73,16 @@ The angle (in radians) between two rotations
 """
 function angledist(a::Quat, b::Quat)
     d = a * conj(b)
-    return 2*atan2(norm(imag(d)), abs(d.w))
+    return 2*atan2(LinearAlgebra.norm(imag(d)), abs(d.w))
 end
 
 """
 The quaternion that transforms a into b through a rotation
 Based on Eigen's implementation for setFromTwoVectors
 """
 function quat_for_a2b(a::VecE3, b::VecE3)
-    v0 = normalize(a)
-    v1 = normalize(b)
+    v0 = LinearAlgebra.normalize(a)
+    v1 = LinearAlgebra.normalize(b)
     c = v1⋅v0
 
     # if dot == -1, vectors are nearly opposites
@@ -97,7 +97,7 @@ function quat_for_a2b(a::VecE3, b::VecE3)
         c = max(c, -1.0)
         M = hcat(convert(Vector{Float64}, v0),
                  convert(Vector{Float64}, v1))'
-        s = svdfact(M, thin=false)
+        s = LinearAlgebra.svdfact(M, thin=false)
         axis = convert(VecE3, s[:V][:,2])
         w2 = (1.0 + c)*0.5
         return Quat(axis * sqrt(1.0 - w2), sqrt(w2))
@@ -116,7 +116,7 @@ a and b for t ∈ [0,1].
 Based on the Eigen implementation for slerp.
 """
 function lerp(a::Quat, b::Quat, t::Float64)
-  
+
     d = a⋅b
     absD = abs(d)
 
@@ -182,7 +182,7 @@ struct RPY <: FieldVector{3, Float64}
 end
 function Base.convert(::Type{RPY}, q::Quat)
 
-    q2 = normalize(q)
+    q2 = LinearAlgebra.normalize(q)
     x = q2[1]
     y = q2[2]
     z = q2[3]
@@ -203,7 +203,7 @@ function Base.convert(::Type{RPY}, q::Quat)
 
     # yaw (z-axis rotation)
     siny = 2(w * z + x * y)
-    cosy = 1.0 - 2(y * y + z * z) 
+    cosy = 1.0 - 2(y * y + z * z)
     yaw = atan2(siny, cosy)
 
     RPY(roll, pitch, yaw)

src/vecSE2.jl

@@ -48,9 +48,9 @@ Base.:%(a::VecSE2, b::Real) = error(op_overload_error)
 scale_euclidean(a::VecSE2, b::Real) = VecSE2(b*a.x, b*a.y, a.θ)
 clamp_euclidean(a::VecSE2, lo::Real, hi::Real) = VecSE2(clamp(a.x, lo, hi), clamp(a.y, lo, hi), a.θ)
 
-Base.norm(a::VecSE2, p::Real=2) = error("norm is not defined for VecSE2 - use norm(VecE2(v)) to get the norm of the Euclidean part")
+norm(a::VecSE2, p::Real=2) = error("norm is not defined for VecSE2 - use norm(VecE2(v)) to get the norm of the Euclidean part")
 normsquared(a::VecSE2) = norm(a)^2 # this will correctly throw an error
-Base.normalize(a::VecSE2, p::Real=2) = error("normalize is not defined for VecSE2. Use normalize_euclidean(v) to normalize the euclidean part of the vector")
+normalize(a::VecSE2, p::Real=2) = error("normalize is not defined for VecSE2. Use normalize_euclidean(v) to normalize the euclidean part of the vector")
 
 function normalize_euclidian(a::VecSE2, p::Real=2)
     n = norm(VecE2(a))

test/runtests.jl

@@ -1,5 +1,7 @@
 using Vec
-using Base.Test
+using Test
+import LinearAlgebra: norm, normalize, dot
+import Random: srand
 
 @test invlerp(0.0,1.0,0.5) ≈ 0.5
 @test invlerp(0.0,1.0,0.4) ≈ 0.4