Svector (#11)

* replace 0.5 by 1//2

* support static arrays

* fix appveyor

Author: jw3126

.travis.yml

@@ -4,7 +4,7 @@ os:
   - linux
   - osx
 julia:
-  - 0.5
+  - 0.6
   - nightly
 notifications:
   email: false

REQUIRE

@@ -1,2 +1 @@
-julia 0.5
-Compat 0.17.0
+julia 0.6

appveyor.yml

@@ -1,7 +1,7 @@
 environment:
   matrix:
-  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.5/julia-0.5-latest-win32.exe"
-  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.5/julia-0.5-latest-win64.exe"
+  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.6/julia-0.6-latest-win32.exe"
+  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.6/julia-0.6-latest-win64.exe"
   - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x86/julia-latest-win32.exe"
   - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x64/julia-latest-win64.exe"
 

src/NumericalIntegration.jl

@@ -2,56 +2,59 @@ __precompile__()
 
 module NumericalIntegration
 
-using Compat
-
 export integrate
 export Trapezoidal, TrapezoidalEven, TrapezoidalFast, TrapezoidalEvenFast
 export SimpsonEven, SimpsonEvenFast
 export IntegrationMethod
 
-@compat abstract type IntegrationMethod end
+abstract type IntegrationMethod end
 
-immutable Trapezoidal         <: IntegrationMethod end
-immutable TrapezoidalEven     <: IntegrationMethod end
-immutable TrapezoidalFast     <: IntegrationMethod end
-immutable TrapezoidalEvenFast <: IntegrationMethod end
-immutable SimpsonEven         <: IntegrationMethod end # https://en.wikipedia.org/wiki/Simpson%27s_rule#Alternative_extended_Simpson.27s_rule
-immutable SimpsonEvenFast     <: IntegrationMethod end
+struct Trapezoidal         <: IntegrationMethod end
+struct TrapezoidalEven     <: IntegrationMethod end
+struct TrapezoidalFast     <: IntegrationMethod end
+struct TrapezoidalEvenFast <: IntegrationMethod end
+struct SimpsonEven         <: IntegrationMethod end # https://en.wikipedia.org/wiki/Simpson%27s_rule#Alternative_extended_Simpson.27s_rule
+struct SimpsonEvenFast     <: IntegrationMethod end
 
 const HALF = 1//2
 
-function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector{Y}, ::Trapezoidal)
-  @assert length(x) == length(y) "x and y vectors must be of the same length!"
-  retval = zero(promote(x[1], y[1])[1])
+function _zero(x,y)
+    ret = zero(eltype(x)) + zero(eltype(y))
+    ret / 2
+end
+
+function integrate(x::AbstractVector, y::AbstractVector, ::Trapezoidal)
+    @assert length(x) == length(y) "x and y vectors must be of the same length!"
+    retval = _zero(x,y)
     for i in 1 : length(y)-1
-      retval += (x[i+1] - x[i]) * (y[i] + y[i+1])
+        retval += (x[i+1] - x[i]) * (y[i] + y[i+1])
     end
-  return HALF * retval
+    return HALF * retval
 end
 
-function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector{Y}, ::TrapezoidalEven)
+function integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEven)
     @assert length(x) == length(y) "x and y vectors must be of the same length!"
     return (x[2] - x[1]) * (HALF * (y[1] + y[end]) + sum(y[2:end-1]))
 end
 
-function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector{Y}, ::TrapezoidalFast)
-    retval = zero(promote(x[1], y[1])[1])
+function integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalFast)
+    retval = _zero(x,y)
     @fastmath @simd for i in 1 : length(y)-1
         @inbounds retval += (x[i+1] - x[i]) * (y[i] + y[i+1])
     end
     return HALF * retval
 end
 
-function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector{Y}, ::TrapezoidalEvenFast)
-    retval = zero(promote(x[1], y[1])[1])
+function integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEvenFast)
+    retval = _zero(x,y)
     N = length(y) - 1
     @fastmath @simd for i in 2 : N
         @inbounds retval += y[i]
     end
     @inbounds return (x[2] - x[1]) * (retval + HALF*y[1] + HALF*y[end])
 end
 
-function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector{Y}, ::SimpsonEven)
+function integrate(x::AbstractVector, y::AbstractVector, ::SimpsonEven)
     @assert length(x) == length(y) "x and y vectors must be of the same length!"
     retval = (17*y[1] + 59*y[2] + 43*y[3] + 49*y[4] + 49*y[end-3] + 43*y[end-2] + 59*y[end-1] + 17*y[end]) / 48
     for i in 5 : length(y) - 1
@@ -60,7 +63,7 @@ function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector
     return (x[2] - x[1]) * retval
 end
 
-function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector{Y}, ::SimpsonEvenFast)
+function integrate(x::AbstractVector, y::AbstractVector, ::SimpsonEvenFast)
     @inbounds retval = 17*y[1] + 59*y[2] + 43*y[3] + 49*y[4]
     @inbounds retval += 49*y[end-3] + 43*y[end-2] + 59*y[end-1] + 17*y[end]
     retval /= 48
@@ -70,6 +73,6 @@ function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector
     @inbounds return (x[2] - x[1]) * retval
 end
 
-integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector{Y}) = integrate(x, y, TrapezoidalFast())
+integrate(x::AbstractVector, y::AbstractVector) = integrate(x, y, TrapezoidalFast())
 
 end

test/REQUIRE

@@ -0,0 +1 @@
+StaticArrays

test/runtests.jl

@@ -1,20 +1,36 @@
 using NumericalIntegration
 using Base.Test
+using StaticArrays
 
-x = collect(-π : π/1000 : π)
-y = sin.(x)
-p = collect(0 : π/1000 : π)
-q = sin.(p)
-for M in subtypes(IntegrationMethod)
-    println("Testing method: ", M)
-    for T in [Float32, Float64, BigFloat]
-        for (xs,ys,val,atol) in [
-                                 (x,y,0,1e-4),
-                                 (p,q,2,1e-4),
-                              ]
-            result = @inferred integrate(T.(xs), T.(ys),M())
-            @test isapprox(result, val, atol=atol)
-            @test typeof(result) == T
+@testset "compare with analytic result" begin
+    x = collect(-π : π/1000 : π)
+    y = sin.(x)
+    p = collect(0 : π/1000 : π)
+    q = sin.(p)
+    for M in subtypes(IntegrationMethod)
+        for T in [Float32, Float64, BigFloat]
+            for (xs,ys,val,atol) in [
+                                     (x,y,0,1e-4),
+                                     (p,q,2,1e-4),
+                                  ]
+                result = @inferred integrate(T.(xs), T.(ys),M())
+                @test isapprox(result, val, atol=atol)
+                @test typeof(result) == T
+            end
         end
     end
 end
+
+@testset "SVector" begin
+    xs  = linspace(0,1,10)
+    ys1 = randn(10)
+    ys2 = randn(10)
+    ys = map(SVector, ys1, ys2)
+    for M in subtypes(IntegrationMethod)
+        m = M()
+        res1 = integrate(xs,ys1,m)
+        res2 = integrate(xs,ys2,m)
+        res = @inferred integrate(xs,ys,m)
+        @test res ≈ @SVector [res1, res2]
+    end
+end