Comparing approaches ‐ AlgebraicVector

Julia interface to AlgebraicVector

CxxWrap.jl approach

CxxWrap support template (parametric) types: https://github.com/JuliaInterop/CxxWrap.jl#template-parametric-types

    // adding AlgebraicVector
    mod.add_type<jlcxx::Parametric<jlcxx::TypeVar<1>>>("AlgebraicVector")
            .apply<AlgebraicVector<DA>, AlgebraicVector<double>>([](auto wrapped) {  // explicitly list all possible template types
        typedef typename decltype(wrapped)::type WrappedT;
        typedef typename WrappedT::value_type ScalarT;  // AlgebraicVector inherits from std::vector which sets value_type

        wrapped.template constructor<const size_t>();

        wrapped.method("sqr", [](const WrappedT& avec) { return avec.sqr(); });

        // add methods to Base
        wrapped.module().set_override_module(jl_base_module);
        // TODO: add show instead of print(ln)
        wrapped.module().method("print", [](const WrappedT& avec) { std::cout << avec; });
        wrapped.module().method("println", [](const WrappedT& avec) { std::cout << avec; });
        // implementing the indexing interface
        wrapped.module().method("getindex", [](const WrappedT& avec, const int i) { return avec[i-1]; });
        wrapped.module().method("setindex!", [](WrappedT& avec, const ScalarT& v, const int i) { avec[i-1] = v; });
        wrapped.module().method("firstindex", [](const WrappedT& avec) { return 1; });  // TODO: do we want to index from 1?
        wrapped.module().method("lastindex", [](const WrappedT& avec) { return avec.size(); });  // TODO: do we want to index from 1?
        // maths functions
        wrapped.module().method("sin", [](const WrappedT& avec) { return sin(avec); });
        wrapped.module().method("cos", [](const WrappedT& avec) { return cos(avec); });
        wrapped.module().method("tan", [](const WrappedT& avec) { return tan(avec); });
        // operators
        // TODO: ...
        // stop adding methods to base
        wrapped.module().unset_override_module();
    });

    // gradient of DA, returns AlgebraicVector
    mod.method("gradient", [](const DA& da)->AlgebraicVector<DA> { return da.gradient(); });

Notes:

• need to explicitly list all template types, here we specify DA and double
• implement indexing interface so in julia you can access elements of the vector using []
• could instead implement array interface so that AlgebraicVector behaves like a Julia array

Tutorial 1 Example 8

https://github.com/dacelib/dace/blob/master/Tutorials/Tutorial1/Example8.cpp

using DACE

DACE.init(1,2)

function somb(x::DACE.AlgebraicVector{DACE.DA})::DACE.DA
    return sin(sqrt(x[1]*x[1] + x[2]*x[2])) / sqrt(x[1]*x[1] + x[2]*x[2])
end

println("Initialize x as two-dim DA vector around (2,3)")

x = DACE.AlgebraicVector{DACE.DA}(2)
x[1] = 2.0 + DACE.DA(1, 1.0)
x[2] = 3.0 + DACE.DA(2, 1.0)
println("x")
println(x)

z = somb(x)
println("Sombrero function")
println(z)

grad_z = DACE.gradient(z)
println("Gradient of sombrero function")
println(grad_z)

ccall approach

We use a Julia Vector{DA} and override/add functions to implement the DACE AlgebraicVector functionality.

"""
    Base.sin(x::Vector{DA})::Vector{DA}

Compute the sine of a Vector{DA}. The result is
copied to a new Vector{DA}.
"""
function Base.sin(x::Vector{DA})::Vector{DA}
    temp = Vector{DA}(undef, length(x))
    for i in 1:length(x)
        temp[i] = sin(x[i])
    end

    return temp
end

"""
    deriv(da::DA, i::Integer)::DA

Compute the derivative of a DA object with respect to variable i.
The result is copied in a new DA object.
"""
function deriv(da::DA, i::Integer)::DA
    temp = DA()
    ccall((:daceDifferentiate, libdace), Cvoid, (Cuint, Ref{Variable}, Ref{Variable}), i, da.index, temp.index)
    exitondaceerror("Error: deriv of DA failed")

    return temp
end

"""
    gradient(x::DA)::Vector{DA}

Compute the gradient of the DA object. Returns a Vector{DA}
containing the derivatives of the DA object with respect to all
independent DA variables.
"""
function gradient(x::DA)::Vector{DA}
    nvar = getmaxvariables()
    temp = Vector{DA}(undef, nvar)
    for i in 1:nvar
        temp[i] = deriv(x, i)
    end

    return temp
end

"""
    Base.show(io::IO, vec::Vector{DA})

Print Vector of DA objects
"""
function Base.show(io::IO, vec::Vector{DA})
    #print(io, tostring(da))
    size = length(vec)
    print(io, "[[[ ")
    print(io, size)
    println(io, " vector")
    for i in 1:size
        println(io, vec[i])
    end
    println(io, "]]]")
end

Notes:

• any functionality on AlgebraicVector that we want to have in C++, we will need to implement for the Vector{DA} type in Julia
• using a native Julia vector

Tutorial 1 Example 8

https://github.com/dacelib/dace/blob/master/Tutorials/Tutorial1/Example8.cpp

using DACE

DACE.init(1, 2)

function somb(x::Vector{DACE.DA})::DACE.DA
    return sin(sqrt(x[1]*x[1] + x[2]*x[2])) / sqrt(x[1]*x[1] + x[2]*x[2])
end

println("Initialize x as two-dim DA vector around (2,3)")

x = Vector{DACE.DA}(undef, 2)
x[1] = 2.0 + DACE.DA(1, 1.0)
x[2] = 3.0 + DACE.DA(2, 1.0)
println("x")
println(x)

z = somb(x)
println("Sombrero function")
println(z)

grad_z = DACE.gradient(z)
println("Gradient of sombrero function")
println(grad_z)