think of interface / dastructures to handle linear-PILS

[lucaferranti]

Let us first focus on symmetric and linear PILS. In a linear PILS we have

A(p) = A0 + A1p1 + A2p2 + .... + An*pn

A few alternatives off the top of my head

Alternative 1

struct IntervalAffineArray{T, MT<: AbstractVecOrMat{T}}
  coeffs::Vector{M}
  params::Vector{Interval{T}
end

Alternative 2

using Symbolics.jl

struct IntervalParametricArray{T, MT <: AbstractVecOrMat{Expr}}
  A::MT
  params::Vector{Interval{T}}
end

I guess this would generalize to more complex expressions, but it would require to analyze the expressions in the matrix to figure out the case, i.e. for the linear case we would need a function

is_affine_multivariate_polynomial(ex::Expr)::Bool

e.g.

is_affine_multivariate_polynomial(x^2+ y) == false
is_affine_multivariate_polynomial(x + y + z + 1) == true
is_affine_multivariate_polynomial(xy + x + y) == false

maybe if we are planning to focus on linear PILS atm, it could be good to go for option 1 (unless that is worse than 2 anyway)

in both cases, the solve interface could be

AbstractParametricIntervalLinearSolver <: AbstractIntervalLinearSolver

solve(Ap, bp, method::AbstractParametricIntervalLinearSolver) = ....

[lucaferranti]

even when going with alternative 1, there could still be a constructor to construct the matrix from the symbolic expression

[schillic]

If you do not need to rely on the concrete type, you can keep the eltype of A parametric.

struct IntervalParametricArray{T, MT<:AbstractVecOrMat}
    A::MT
    params::Vector{Interval{T}}
end

But I agree to go for Option 1 because it can always be changed later. A short note in the documentation that this is work in progress and the interface may change would be good, though.

[lucaferranti]

Looking at skalna2006 paper, an alternative structure could be

Alternative 3

struct IntervalParametricArray where {T}
  Ω::Matrix{Vector{T}}
  params::Vector{Interval{T}}
end

so that A(i, j) = Ω(i, j)ᵀp (p being the vector of symbolic parameters).

that's basically the same of alternative 1, except that instead of having a vector of matrices you have a matrix of vectors.

Maybe not the most optimal solution, but would be at least straightforward to implement with DynamicPolynomials

julia> @polyvar x y z
(x, y, z)

julia> A = [x+2y+z+1 2x+z+2;y/2 z+x+y+1]
2×2 Matrix{Polynomial{true, Float64}}:
 x + 2.0y + z + 1.0  2.0x + z + 2.0
 0.5y                x + y + z + 1.0

julia> map(a -> coefficients(a, [x, y, z, 1]), A)
2×2 Matrix{Vector{Float64}}:
 [1.0, 2.0, 1.0, 1.0]  [2.0, 0.0, 1.0, 2.0]
 [0.0, 0.5, 0.0, 0.0]  [1.0, 1.0, 1.0, 1.0]

[lucaferranti]

With a quick googling, it seems Symbolics.jl doesn't have (yet) the needed functionalities (extract coefficients wrt given variable, check if expression is linear etc). Alternatives could be 1) use DynamicPolynomials (maybe an overkill) 2) if for the moment we restrict to linear expressions, could just create a simple type LinearExpression with addition and subtraction overloaded.

good to notice that the interfaces are just for human-readable output/-writable input, at least for the algorithm in Skalna06 symbolic manipulation is not needed

[mforets]

If possible let's use the AbstractIntervalMatrix interface (conceptually, a parametric interval matrix resp system is a special case of a general interval matrix, resp system). Moreover, I would rather leave the params vector type generic (algebraic ops and using eg. using statically sized arrays is a reasonable thing to do).

These observations lead to:

import IntervalMatrices: AbstractIntervalMatrix

struct ParametricIntervalMatrix{N, MN<: AbstractMatrix{N}, VMN<:AbstractVector{MN}, T, IT, VI<: AbstractVector{IT}} <: AbstractIntervalMatrix{IT}
  coeffs::VMN
  params::VI
end

We may decide that the "A0" matrix is another field, or that the length of coeffs is 1 plus the length of params and checked in a constructor.

[mforets]

good to notice that the interfaces are just for human-readable output/-writable input, at least for the algorithm in Skalna06 symbolic manipulation is not needed

Indeed, the discussion of using symbolic expressions to initialize such objects is nice but also somehow orthogonal, so we may split it to another issue.

it seems Symbolics.jl doesn't have (yet) the needed functionalities (extract coefficients wrt given variable, check if expression is linear etc)

But it shouldn't be hard to do, at least naively (now that I think about it, there is some code in LazySets to parse hyperplanes and half-spaces, both linear expressions).

The long-term benefits of using an established and active project such as Symbolics.jl (or MultivariatePolynomials, or both) exceed the short-term benefits of cooking custom symbolic structs as in https://github.com/JuliaIntervals/IntervalLinearAlgebra.jl/pull/104 -- which I find it is a fantastic demo!

[lucaferranti]

Yes parametrizing on the constructor type would indeed be good, I was just lazy when writing the comment 😄

Now that I think a little better about it, maybe it makes sense to keep the list of parameters separated from the matrix and make ParametricMatrix callable, that is

struct ParametricMatrix{T, VT<:AbstractVector{T}, MT<:AbstractMatrix{VT}} <: AbstractMatrix{VT}
  coeffs::MT
end

# define array interface here, getindex, setindex, size etc

function (A::ParametricIntervalMatrix)(p::AbstractVector)
   return [dot(a[1:end-1], p) + a[end] for a in A]
end

or something like that, so that one could simply do

A = ParametricMatrix([1+p1 p2;p2 p1-1])
p1 = [1..2, 3..4]
p2 = [1, 2]
Ap = A(p)
Amid = A(p1)

This would make it easier to use the matrix over different parameters domains or getting an instance of the system at a specific parameter value.

[mforets]

ok :+1: something like that but also allowing to store the coeffs in a struct seems like a nice approach.

i mean, after it's been created, it has non-zero cost to "extract" the coeffs from Ap, even though it'll possibly be required to use the params p separately from the coeffs matrix in downstream computations.

[mforets]

i mean the "callable" approach also works with the definition in https://github.com/JuliaIntervals/IntervalLinearAlgebra.jl/issues/100#issuecomment-933592992. and provided that the container type VI is mutable, one can always update such coeffs after creation:

@variables p1 p2

A = ParametricIntervalMatrix([1+p1 p2;p2 p1-1]) # params field is initialized with interval(1)
p = [1..2, 3..4]

Ap = A(p) # instance with old A and new p

setparams!(A, p) # in-place update of A's params field 

In this way, ParametricIntervalMatrix works like a "template".