Merge pull request #88 from SciML/linear_algebra

support linear algebra on ArrayPartition

Author: ChrisRackauckas

.travis.yml

@@ -3,11 +3,11 @@ language: julia
 os:
   - linux
 julia:
-  - 1.0
-  - nightly
-matrix:
-  allow_failures:
-    - julia: nightly
+  - 1
+#  - nightly
+#matrix:
+#  allow_failures:
+#    - julia: nightly
 notifications:
   email: false
 # uncomment the following lines to override the default test script

Project.toml

@@ -5,19 +5,22 @@ version = "2.2.0"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 ZygoteRules = "700de1a5-db45-46bc-99cf-38207098b444"
 
 [compat]
-ArrayInterface = "1.2, 2.0"
+ArrayInterface = "2.7"
 RecipesBase = "0.7, 0.8, 1.0"
 Requires = "0.5, 1.0"
 StaticArrays = "0.10, 0.11, 0.12"
 ZygoteRules = "0.2"
-julia = "1"
+julia = "1.3"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/RecursiveArrayTools.jl

@@ -3,7 +3,7 @@ __precompile__()
 module RecursiveArrayTools
 
   using Requires, RecipesBase, StaticArrays, Statistics,
-        ArrayInterface, ZygoteRules
+        ArrayInterface, ZygoteRules, LinearAlgebra
 
   abstract type AbstractVectorOfArray{T, N, A} <: AbstractArray{T, N} end
   abstract type AbstractDiffEqArray{T, N, A} <: AbstractVectorOfArray{T, N, A} end

src/array_partition.jl

@@ -233,6 +233,8 @@ ArrayPartitionStyle(::Val{N}) where N = ArrayPartitionStyle{Broadcast.DefaultArr
 function Broadcast.BroadcastStyle(::ArrayPartitionStyle{AStyle}, ::ArrayPartitionStyle{BStyle}) where {AStyle, BStyle}
     ArrayPartitionStyle(Broadcast.BroadcastStyle(AStyle(), BStyle()))
 end
+Broadcast.BroadcastStyle(::ArrayPartitionStyle, ::Broadcast.DefaultArrayStyle{0}) = Broadcast.DefaultArrayStyle{1}()
+Broadcast.BroadcastStyle(::ArrayPartitionStyle, ::Broadcast.DefaultArrayStyle{N}) where N = Broadcast.DefaultArrayStyle{N}()
 
 combine_styles(args::Tuple{})         = Broadcast.DefaultArrayStyle{0}()
 combine_styles(args::Tuple{Any})      = Broadcast.result_style(Broadcast.BroadcastStyle(args[1]))
@@ -252,7 +254,7 @@ end
     ArrayPartition(f, N)
 end
 
-@inline function Base.copyto!(dest::ArrayPartition, bc::Broadcast.Broadcasted)
+@inline function Base.copyto!(dest::ArrayPartition, bc::Broadcast.Broadcasted{ArrayPartitionStyle{Style}}) where Style
     N = npartitions(dest, bc)
     for i in 1:N
         copyto!(dest.x[i], unpack(bc, i))
@@ -293,3 +295,10 @@ common_number(a, b) =
     (b == 0 ? a :
      (a == b ? a :
       throw(DimensionMismatch("number of partitions must be equal"))))
+
+## Linear Algebra
+
+ArrayInterface.zeromatrix(A::ArrayPartition) = ArrayInterface.zeromatrix(reduce(vcat,vec.(A.x)))
+LinearAlgebra.ldiv!(A::LinearAlgebra.LU,b::ArrayPartition) = ldiv!(A,Array(b))
+LinearAlgebra.ldiv!(A::LinearAlgebra.QR,b::ArrayPartition) = ldiv!(A,Array(b))
+LinearAlgebra.ldiv!(A::LinearAlgebra.SVD,b::ArrayPartition) = ldiv!(A,Array(b))

test/partitions_test.jl

@@ -26,7 +26,21 @@ a = 5
 @. p = p*p2
 K = p.*p2
 
-@test_broken p.*rand(10)
+x = rand(10)
+y = p.*x
+@test y[1:5] == p.x[1] .* x[1:5]
+@test y[6:10] == p.x[2] .* x[6:10]
+y = p.*x'
+for i in 1:10
+  @test y[1:5,i] == p.x[1] .* x[i]
+  @test y[6:10,i] == p.x[2] .* x[i]
+end
+y = p .* p'
+@test y[1:5,1:5] == p.x[1] .* p.x[1]'
+@test y[6:10,6:10] == p.x[2] .* p.x[2]'
+@test y[1:5,6:10] == p.x[1] .* p.x[2]'
+@test y[6:10,1:5] == p.x[2] .* p.x[1]'
+
 b = rand(10)
 c = rand(10)
 copyto!(b,p)
@@ -94,13 +108,11 @@ S = [
     ]
 
 for sizes in S
-  x = ArrayPartition( randn.(sizes[1]) ) 
+  x = ArrayPartition( randn.(sizes[1]) )
   y = ArrayPartition( zeros.(sizes[2]) )
   y_array = zeros(length(x))
   copyto!(y,x)           #testing Base.copyto!(dest::ArrayPartition,A::ArrayPartition)
   copyto!(y_array,x)     #testing Base.copyto!(dest::Array,A::ArrayPartition)
   @test all([x[i] == y[i] for i in eachindex(x)])
   @test all([x[i] == y_array[i] for i in eachindex(x)])
 end
-  
-

test/runtests.jl

@@ -7,4 +7,5 @@ using Test
   @time @testset "VecOfArr Indexing Tests" begin include("basic_indexing.jl") end
   @time @testset "VecOfArr Interface Tests" begin include("interface_tests.jl") end
   @time @testset "StaticArrays Tests" begin include("copy_static_array_test.jl") end
+  @time @testset "Upstream Tests" begin include("upstream.jl") end
 end

test/upstream.jl

@@ -0,0 +1,30 @@
+using OrdinaryDiffEq, NLsolve, RecursiveArrayTools, Test, ArrayInterface
+function lorenz(du,u,p,t)
+    du[1] = 10.0*(u[2]-u[1])
+    du[2] = u[1]*(28.0-u[3]) - u[2]
+    du[3] = u[1]*u[2] - (8/3)*u[3]
+end
+u0 = ArrayPartition([1.0,0.0],[0.0])
+@test ArrayInterface.zeromatrix(u0) isa Matrix
+tspan = (0.0,100.0)
+prob = ODEProblem(lorenz,u0,tspan)
+sol = solve(prob,Tsit5())
+sol = solve(prob,AutoTsit5(Rosenbrock23(autodiff=false)))
+sol = solve(prob,AutoTsit5(Rosenbrock23()))
+
+function f!(F, vars)
+    x = vars.x[1]
+    F.x[1][1] = (x[1]+3)*(x[2]^3-7)+18
+    F.x[1][2] = sin(x[2]*exp(x[1])-1)
+    y=vars.x[2]
+    F.x[2][1] = (y[1]+3)*(y[2]^3-7)+18
+    F.x[2][2] = sin(y[2]*exp(y[1])-1)
+end
+
+# To show that the function works
+F = ArrayPartition([0.0 0.0],[0.0, 0.0])
+u0= ArrayPartition([0.1; 1.2], [0.1; 1.2])
+result = f!(F, u0)
+
+# To show the NLsolve error that results with ArrayPartitions:
+nlsolve(f!, ArrayPartition([0.1; 1.2], [0.1; 1.2]))