n-dimension integration (#17)

* n-dimensional integration

* added Project and Manifest

* removed Manifest and Project

* fixing versioninfo not defined appveyor issue

* added HCubature to REQUIRE (for testing n-dimensional integration)

* moved HCubature to test/REQUIRE, and making range 0.7 friendly

* added using Pkg to work with nightly

Author: m-wells

appveyor.yml

@@ -28,8 +28,8 @@ install:
 build_script:
 # Need to convert from shallow to complete for Pkg.clone to work
   - IF EXIST .git\shallow (git fetch --unshallow)
-  - C:\projects\julia\bin\julia -e "versioninfo();
+  - C:\projects\julia\bin\julia -e "using InteractiveUtils; versioninfo(); using Pkg;
       Pkg.clone(pwd(), \"NumericalIntegration\"); Pkg.build(\"NumericalIntegration\")"
 
 test_script:
-  - C:\projects\julia\bin\julia -e "Pkg.test(\"NumericalIntegration\")"
+  - C:\projects\julia\bin\julia -e "using Pkg; Pkg.test(\"NumericalIntegration\")"

src/NumericalIntegration.jl

@@ -176,6 +176,22 @@ function integrate(x::AbstractVector, y::AbstractVector, m::RombergEven)
     @inbounds return rombaux[maxsteps, prevrow]
 end
 
+"""
+    integrate(Y::AbstractArray{T,N}, Varargs{AbstractVector,N}, method)
+
+Given an n-dimensional grid of values, compute the total integral along each dim
+"""
+function integrate(X::NTuple{N,AbstractVector}, Y::AbstractArray{T,N}, M::IntegrationMethod) :: T where {T,N}
+    dims = size(Y)
+    d,n = length(dims), dims[end]
+    if d == 1
+        return integrate(X[1], Y, M)
+    else
+        return integrate(X[end], [integrate(X[1:d-1], selectdim(Y,d,i), M) for i in 1:n], M)
+    end
+end
+
+
 
 # cumulative integrals
 
@@ -242,12 +258,13 @@ function cumul_integrate(x::AbstractVector, y::AbstractMatrix, M::IntegrationMet
     return hcat([cumul_integrate(x,selectdim(y,dims,j),M) for j=1:size(y,dims)]...)
 end
 
-
 #default behaviour
 integrate(x::AbstractVector, y::AbstractVector) = integrate(x, y, Trapezoidal())
 
 integrate(x::AbstractVector, y::AbstractMatrix; dims=2) = integrate(x, y, Trapezoidal(); dims=dims)
 
+integrate(X::NTuple, Y::AbstractArray) = integrate(X, Y, Trapezoidal())
+
 cumul_integrate(x::AbstractVector, y::AbstractVector) = cumul_integrate(x, y, Trapezoidal())
 
 cumul_integrate(x::AbstractVector, y::AbstractMatrix; dims=2) = cumul_integrate(x, y, Trapezoidal(); dims=dims)

test/REQUIRE

@@ -1 +1,2 @@
 StaticArrays
+HCubature

test/runtests.jl

@@ -2,6 +2,7 @@ using NumericalIntegration
 using Test
 using InteractiveUtils # for subtypes
 using StaticArrays
+using HCubature # for testing n-dimensional integration
 
 @testset "compare with analytic result" begin
     x = collect(-π : 2*π/2048 : π)
@@ -22,6 +23,25 @@ using StaticArrays
     end
 end
 
+@testset "n-dimensional integration testing" begin
+    X = range(0,stop=2π,length=10)
+    Y = range(-π,stop=π,length=10)
+    Z = range(0,stop=2,length=10)
+
+    A = Array{Float64}(undef,length(X),length(Y),length(Z))
+    f(x) = sin(x[1]) + cos(x[2]) + 2x[3]
+
+    for (k,z) in enumerate(Z)
+        for (j,y) in enumerate(Y)
+            for (i,x) in enumerate(X)
+                A[i,j,k] = f([x,y,z])
+            end
+        end
+    end
+
+    @test isapprox(integrate((X,Y,Z),A), hcubature(f,[0,-π,0],[2π,π,2])[1], atol=1e-4)
+end
+
 @testset "SVector" begin
     xs = range(0, stop=1, length=9)
     ys1 = randn(9)