Add linear algebra routines

.github/workflows/documentation.yml

@@ -15,13 +15,17 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v2
+      - uses: julia-actions/setup-julia@latest
         with:
           version: '1.10'
       - name: Install dependencies
-        run: julia --project=docs/ -e 'using Pkg; Pkg.instantiate()'
+        run: |
+          import Pkg
+          Pkg.develop("DACE")
+          Pkg.instantiate()
+        shell: julia --color=yes --project=docs/ {0}
       - name: Build and deploy
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # If authenticating with GitHub Actions token
           DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # If authenticating with SSH deploy key
-        run: julia --project=docs/ docs/make.jl
+        run: julia --color=yes --project=docs/ docs/make.jl

.github/workflows/test.yml

@@ -9,7 +9,7 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        julia-version: ['1.9', '1.10']
+        julia-version: ['1.10', '1.11']
         julia-arch: [x64]
         os: [ubuntu-latest, macOS-latest, windows-latest]
         exclude:

Project.toml

@@ -7,19 +7,25 @@ version = "0.1.0"
 CxxWrap = "1f15a43c-97ca-5a2a-ae31-89f07a497df4"
 DACE_jll = "40de70a5-cf80-53d9-bda4-3aa67fea2f4f"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
-SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+libcxxwrap_julia_jll = "3eaa8342-bff7-56a5-9981-c04077f7cee7"
 
 [compat]
-julia = "1.6"
-CxxWrap = "~0.16"
-DACE_jll = "~0.5"
+CxxWrap = "^0.17"
+DACE_jll = "^0.7.1"
 DiffEqBase = "^6"
-SpecialFunctions = "^2"
 DocStringExtensions = "^0.9"
+SpecialFunctions = "^2"
+julia = ">=1.6"
 
 [extras]
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+DifferentiableEigen = "73a20539-4e65-4dcb-a56d-dc20f210a01b"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["DelimitedFiles", "DifferentiableEigen", "ForwardDiff", "GenericLinearAlgebra", "LinearAlgebra", "Test"]

docs/Project.toml

@@ -1,16 +1,15 @@
 [deps]
-Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 CxxWrap = "1f15a43c-97ca-5a2a-ae31-89f07a497df4"
+DACE = "f4439622-e757-4457-9ca0-f65a61d23f40"
 DACE_jll = "40de70a5-cf80-53d9-bda4-3aa67fea2f4f"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
-SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
-CxxWrap = "~0.16"
-DACE_jll = "~0.5"
 DiffEqBase = "^6"
-SpecialFunctions = "^2"
 DocStringExtensions = "^0.9"
+SpecialFunctions = "^2"

src/DACE.jl

@@ -16,7 +16,7 @@ module DACE
         m_ub::Float64
     end
 
-    # load the C++ interface and initialise it
+    # load the C++ interface and initialize it
     @wrapmodule(() -> libdace, :define_julia_module)
     function __init__()
         @initcxx
@@ -29,6 +29,9 @@ module DACE
     # extend DACE functionality #
     # ========================= #
 
+    # type alias for convenience
+    const DAOrF64 = Union{DA, Float64}
+
     # ------------------- #
     # custom constructors #
     # ------------------- #
@@ -40,36 +43,44 @@ module DACE
     DAAllocated() = DA(0.0)
     DAAllocated(x::Real) = DA(x)
 
-    # AlgebraicVector objects
-    # TODO is there a better way to do this?
-    AlgebraicVector(v::Vector{<:DA}) = AlgebraicVector{DA}(v)
-    AlgebraicVector(v::Vector{Float64}) = AlgebraicVector{Float64}(v)
-    AlgebraicVector{T}(v::Vector{<:T}) where T = begin
+    # AlgebraicVector and AlgebraicMatrix constructors from Julia vectors and matrices
+    AlgebraicVector(v::AbstractVector{<:DA}) = AlgebraicVector{DA}(v)
+    AlgebraicMatrix(m::AbstractMatrix{<:DA}) = AlgebraicMatrix{DA}(m)
+    AlgebraicVector(v::AbstractVector{Float64}) = AlgebraicVector{Float64}(v)
+    AlgebraicMatrix(m::AbstractMatrix{Float64}) = AlgebraicMatrix{Float64}(m)
+
+    AlgebraicVector{T}(v::AbstractVector{U}) where {T<:DAOrF64, U<:DAOrF64} = begin
         res = AlgebraicVector{T}(length(v))
-        for i in eachindex(v)
-            res[i] = v[i]
-        end
+        res .= v
+        return res
+    end
+
+    AlgebraicMatrix{T}(m::AbstractMatrix{U}) where {T<:DAOrF64, U<:DAOrF64} = begin
+        res = AlgebraicMatrix{T}(size(m)...)
+        res .= m
         return res
     end
 
     # -------------------------- #
     # overload functions in Base #
     # -------------------------- #
 
-    # addittive and multiplicative identities
+    # additive and multiplicative identities
     Base.zero(::Type{DA}) = DA(0.0)
     Base.one(::Type{DA}) = DA(1.0)
 
-    # similar function for AlgebraicVector
-    function Base.similar(foo::AlgebraicVector{DA})
-        sz = size(foo)[1]
-        bar = AlgebraicVector{DA}(sz)
-        return bar
-    end
+    # similar function for AlgebraicVector and AlgebraicMatrix
+    Base.similar(v::AlgebraicVector{<:DA}) = AlgebraicVector{DA}(v)
+    Base.similar(m::AlgebraicMatrix{<:DA}) = AlgebraicMatrix{DA}(m)
+    Base.similar(v::AlgebraicVector{Float64}) = AlgebraicVector{Float64}(v)
+    Base.similar(m::AlgebraicMatrix{Float64}) = AlgebraicMatrix{Float64}(m)
 
     # promotions and conversions to DA
     @inline Base.promote_rule(::Type{T}, ::Type{R}) where {T<:DA, R<:Real} = T
     @inline Base.promote_rule(::Type{R}, ::Type{T}) where {T<:DA, R<:Real} = T
+    # FIXME: how to handle these type promotions properly?
+    @inline Base.promote_rule(::Type{T}, ::Type{DAAllocated}) where T<:DA = DAAllocated
+    @inline Base.promote_rule(::Type{DAAllocated}, ::Type{T}) where T<:DA = DAAllocated
 
     for R in (AbstractFloat, AbstractIrrational, Integer, Rational)
         @eval Base.convert(::Type{<:DA}, x::$R) = DA(convert(Float64, x))
@@ -88,31 +99,17 @@ module DACE
         @eval Base.:^(a::$R, b::DA) = a^DACE.cons(b)
     end
 
-    # comparison operators (considering only the constant part)
-    for op = (:(==), :(!=), :<, :(<=), :>, :(>=))
-
-        # both arguments are DA objects
-        @eval Base.$op(a::DA, b::DA) = $op(DACE.cons(a), DACE.cons(b))
-
-        # one argument is a number
-        for R in (AbstractFloat, AbstractIrrational, Integer, Rational)
-            @eval Base.$op(a::DA, b::$R) = $op(DACE.cons(a), b)
-            @eval Base.$op(a::$R, b::DA) = $op(a, DACE.cons(b))
-        end
-    end
-
-    Base.isfinite(a::DA) = isfinite(DACE.cons(a))
-    Base.isinf(a::DA) = isinf(DACE.cons(a))
-    Base.isnan(a::DA) = isnan(DACE.cons(a))
+    Base.isfinite(a::DA) = !(isinf(a) || isnan(a))
 
     # assignment of vector and matrix elements
     Base.setindex!(v::AlgebraicVector{<:DA}, x::Real, i::Integer) = v[i] = convert(DA, x)
     Base.setindex!(m::AlgebraicMatrix{<:DA}, x::Real, i::Integer, j::Integer) = m[i,j] = convert(DA, x)
 
-    # custom show functions for DA and AlgebraicVector objects
+    # custom show functions for DA-related objects
+    Base.show(io::IO, m::Monomial) = print(io, toString(m))
     Base.show(io::IO, da::DA) = print(io, toString(da))
-    Base.show(io::IO, vec::AlgebraicVector) = print(io, toString(vec))
-    Base.show(io::IO, mat::AlgebraicMatrix) = print(io, toString(mat))
+    Base.show(io::IO, vec::AlgebraicVector{T}) where T<:DAOrF64 = print(io, toString(vec))
+    Base.show(io::IO, mat::AlgebraicMatrix{T}) where T<:DAOrF64 = print(io, toString(mat))
 
     # -------------------------------------- #
     # overload functions in SpecialFunctions #
@@ -133,26 +130,29 @@ module DACE
     # ---------------------------------------------- #
 
     # map inversion
-    invert(v::Vector{<:DA}) = Vector{DA}(invert(AlgebraicVector(v)))
+    invert(v::AbstractVector{<:DA}) = invert(AlgebraicVector(v))
 
     # linear part, Jacobian and Hessian
-    linear(v::Vector{<:DA}) = linear(AlgebraicVector(v))
-    jacobian(v::Vector{<:DA}) = jacobian(AlgebraicVector(v))
-    hessian(v::Vector{<:DA}) = stack(hessian(AlgebraicVector(v)), dims=3)
+    linear(v::AbstractVector{<:DA}) = linear(AlgebraicVector(v))
+    jacobian(v::AbstractVector{<:DA}) = jacobian(AlgebraicVector(v))
+    hessian(v::AbstractVector{<:DA}) = hessian(AlgebraicVector(v)) # vector of Hessian matrices
+    hess_stack(a::AbstractVector{<:DA}) = stack(hessian(a), dims=3)
 
     # compilation and evaluation of DA objects
-    compile(v::Vector{<:DA}) = compile(StdVector{DA}(v))
-    for R in (DA, Float64)
-        @eval eval(cda::compiledDA, v::Vector{<:$R}) = Vector{$R}(eval(cda, AlgebraicVector(v)))
-        @eval eval(da::DA, v::Vector{<:$R}) = eval(da, AlgebraicVector(v))
-        @eval eval(a::Vector{<:DA}, v::Vector{<:$R}) = Vector{$R}(eval(AlgebraicVector(a), AlgebraicVector(v)))
-        @eval eval(a::AlgebraicVector{DA}, v::Vector{<:$R}) = Vector{$R}(eval(a, AlgebraicVector(v)))
-    end
+    compile(v::AbstractVector{<:DA}) = compile(AlgebraicVector(v))
+
+    eval(cda::compiledDA, v::AbstractVector{<:DA}) = eval(cda, AlgebraicVector(v))
+    eval(a::AbstractVector{<:DA}, v::AbstractVector{<:DA}) = eval(AlgebraicVector(a), AlgebraicVector(v))
+    eval(a::AlgebraicVector{<:DA}, v::AbstractVector{<:DA}) = eval(a, AlgebraicVector(v))
+
+    eval(cda::compiledDA, v::AbstractVector{Float64}) = eval(cda, AlgebraicVector(v))
+    eval(a::AbstractVector{<:DA}, v::AbstractVector{Float64}) = eval(AlgebraicVector(a), AlgebraicVector(v))
+    eval(a::AlgebraicVector{<:DA}, v::AbstractVector{Float64}) = eval(a, AlgebraicVector(v))
 
     # ------- #
     # exports #
     # ------- #
 
-    export DA, AlgebraicVector, AlgebraicMatrix, compiledDA
+    export DA, AlgebraicVector, AlgebraicMatrix, compiledDA, Monomial
 
 end  # module DACE

test/comparison_operators.jl

@@ -13,7 +13,7 @@
 
     @testset "Test >=" begin
         DACE.init(10, 4)
-        
+
         x = 1 + DA(1, 1)
         y = DA(2)
         @test y >= x
@@ -22,4 +22,36 @@
 
     end
 
+    @testset "Test isnan" begin
+        DACE.init(2, 2)
+
+        x = 1.0 + DA(1, 1.0)
+        y = 1.0 + DA(2, Inf)
+        z = 1.0 + DA(2, NaN)
+        @test isnan(x) == false
+        @test isnan(y) == false
+        @test isnan(z) == true
+    end
+
+    @testset "Test isinf" begin
+        DACE.init(2, 2)
+
+        x = 1.0 + DA(1, 1.0)
+        y = 1.0 + DA(2, Inf)
+        z = 1.0 + DA(2, NaN)
+        @test isinf(x) == false
+        @test isinf(y) == true
+        @test isinf(z) == false
+    end
+
+    @testset "Test isfinite" begin
+        DACE.init(2, 2)
+
+        x = 1.0 + DA(1, 1.0)
+        y = 1.0 + DA(2, Inf)
+        z = 1.0 + DA(2, NaN)
+        @test isfinite(x) == true
+        @test isfinite(y) == false
+        @test isfinite(z) == false
+    end
 end