Add function to return gradient of continuous system

[ueliwechsler]

As mentioned in this post by @schillic in #96 (comment)
I suggest adding a functionality corresponding to the discrete successor method for continuous systems. Possible names are:

• derivative
• tangent
• gradient
• vectorfield

[mforets]

What would be the signature of these methods, i mean in terms of the inputs/outputs?

About the name, I think that vector_field, tangent_field or direction_field sounds good to me. Perhaps I would reserve derivative/tangent/gradient for the point-wise functions.

[ueliwechsler]

What would be the signature of these methods, i mean in terms of the inputs/outputs?

I am not sure, if I understand.
I was thinking about not specifying the types and use the default operators for the types, e.g.

function vector_field(sys::LinearContinuousSystem, x, u)
return state_matrix(sys)*x + input_matrix(sys)*u
end

[mforets]

To clarify my previous comment, what i had in mind was the distinction between the field and the function-like behavior (ie. the field applied to a given point / tuple). As it turns out in applications having the field is needed (eg. for what you asked here, we can just wrap the field and pass it to a solver).

Here is a proposal:

using Revise, MathematicalSystems, LazySets

MathematicalSystems.system(sys::AbstractSystem) = sys
MathematicalSystems.system(sys::InitialValueProblem) = sys.s

struct VectorField{T}
    field::T
end

# function-like evaluation
@inline function (V::VectorField)(args...)
    evaluate(V, args...)
end

function evaluate(V::VectorField, args...)
    return V.field(args...)
end

function VectorField(sys::AbstractSystem)
    sys = system(sys)
    if islinear(sys)
        if inputdim(sys) == 0
            field = (x) -> state_matrix(sys) * x
        else
            field = (x, u) -> state_matrix(sys) * x + input_matrix(sys) * u
        end
    elseif isaffine(sys)
        if inputdim(sys) == 0
            field = (x) -> state_matrix(sys) * x + affine_term(sys)
        else
            field = (x, u) -> state_matrix(sys) * x + affine_term(sys) + input_matrix(sys) * u
        end
    else
        error("the vector field for a system of type $sys is not implemented yet")
    end
    
    return VectorField(field)
end

function vector_field(sys::AbstractSystem, args...)
    return evaluate(VectorField(sys), args...)
end

To be used as:

julia> V = VectorField(@system(x' = x));

julia> V([0.0]), V([1.0])
([0.0], [1.0])

julia> V = VectorField(@system(x' = 2x + 1));
​
julia> V([0.0]), V([1.0])
([1.0], [3.0])

julia> vector_field(@system(x' = x), [0.0])
1-element Array{Float64,1}:
 0.0

[ueliwechsler]

I implemented another way to get the same result.
Is there a disadvantage to the structure proposed in the following?

MathematicalSystems.jl/src/successor.jl

Lines 215 to 237 in 0bbe3df

	struct VectorField{T}
	field::T
	end
	# function-like evaluation
	@inline function (V::VectorField)(args...)
	evaluate(V, args...)
	end
	function evaluate(V::VectorField, args...)
	return V.field(args...)
	end
	function VectorField(sys::AbstractContinuousSystem)
	if inputdim(sys) == 0 && noisedim(sys) == 0
	field = (x) -> vector_field(sys, x)
	elseif inputdim(sys) == 0 || noisedim(sys) == 0
	field = (x, u) -> vector_field(sys, x, u)
	else
	field = (x, u, w) -> vector_field(sys, x, u, w)
	end
	return VectorField(field)
	end

As a result, the vector_field / successor method could be concisely defined in MathematicalSystems and the VectorField type e.g. in a down-stream package.

[ueliwechsler]

In the same branch, I also propose a more generic way to generate the vector_field and successor method than in the current implementation. It is a tradeoff between lines of code and performance. (if you like it, maybe #142 could be done in a similar way).

[mforets]

ok, i see, and sorry for my late reply.. i was too busy the last couple of days :/

i checked the branch you propose, thanks! do you mind making to follow up with a pull request?

and the VectorField type e.g. in a down-stream package.

i find it useful to define VectorField in this package.

It is a tradeoff between lines of code and performance.

not sure that i follow, why do you mention performance? can you elaborate?

Is there a disadvantage to the structure proposed in the following?

i don't know, but i think that apply and vector_field have different meanings: apply is to apply a map; because they are functions in -> out := map(in), while vector_field evaluates the right-hand side at a given point. at least, i don't see why should we mix these terms..

but apart from some comments,eg. do not assume that what is neither linear or affine has a field f etc, the code looks good to me, we can continue the discussion in your PR.

[ueliwechsler]

ok, i see, and sorry for my late reply.. i was too busy the last couple of days :/

No worries, same here :)

not sure that i follow, why do you mention performance? can you elaborate?

I assumed that adding the additional logic for distinguishing between the cases will increase the runtime compared to just defining successor for every type individually. But benchmarks show that the effect is neglectable.

i don't know, but i think that apply and vector_field have different meanings: apply is to apply a map; because they are functions in -> out := map(in), while vector_field evaluates the right-hand side at a given point. at least, i don't see why should we mix these terms..

I used apply because it can be used for both vector_field and successor since its the same computation.
I am not sure if I understand this distinction. Isn't apply a more general term and vector_field and successor ar a subset of it?

do not assume that what is neither linear or affine has a field f

What about defining a getter for the functional fields, e.g. f, p like mapping similar to state_matrix, etc? An maybe also add an hasmapping method for checking if a mapping exists.