Basic operator composition

[MSeeker1340]

The aim is to implement the most basic subset of #60 so that basic operator compositions like addition and multiplication can be constructed, which would then enable the implementation of lazy jacobian/W operator from SciML/DifferentialEquations.jl#220.

The main improvement I'm considering is how to implement DiffEqOperatorCombination. The prototype in the notebook is still influenced by LinearMaps.jl in that both the component operators and their coefficients are stored. While this is perfectly fine for LinearMaps where the coefficients are constant, it would become a bit messy if we allow the coefficients to also be time-dependent. For example, in the notebook I need to also include an update_func for the coefficients in addition to the update functions of the component operators.

My solution to this is to add an additional type DiffEqScaledOperator, which is constructed by multiplying a DiffEqScalar and another operator. The coefficients in DiffEqOperatorCombination will be dropped. For example, a composite operator constructed as L = a * A + B for will be constructed as

DiffEqOperatorCombination(
  DiffEqScaledOperator(DiffEqScalar(a), A),
  B
)

In this way, calling update_coefficients! to a composite operator, however complex it is, boils down to doing a tree transversal and calling the embedded update functions only on the leaf nodes (in the above case: DiffEqScalar(a), A and B).

I have also added a setval! method to the scalar and array operators as a way to bypass the update_coefficients! interface and directly mutate their value. This should be used sparsely and only in situation where outside information other than (u,p,t) is needed to update the operator (for example, dt for the W operator).

Finally, since full is deprecated in v0.7 I have switched to using Matrix(L). For the as_array function in the notebook, maybe convert(::AbstractArray, L)?

[ChrisRackauckas]

Why make it an overloaded singleton instead of a function?

[ChrisRackauckas]

For the as_array function in the notebook, maybe convert(::AbstractArray, L)?

Good choice.

[MSeeker1340]

Some notes before moving on:

1. I'm not sure why it's necessary to keep both update_coefficients and update_coefficients!. Despite the name difference both are mutating, and the only difference is that the former returns the updated operator and the latter returns nothing. The return value isn't really used anyways. I think having just one update_coefficients! method should be enough.

2. While this is not necessary, I should probably write a pprint method for the composite operators. The default outputs can be extremely hard to decipher even for simple compositions. I'm thinking about doing something like

   alpha = DiffEqScalar(2.0)
   A = DiffEqArrayOperator(rand(2,2))
   update_func = ...
   B = DiffEqArrayOperator(rand(2,2); update_func=update_func)
   L = alpha * A + B
   
   linear combination of
     scaled operator with
       coefficient: constant scalar operator with value 2.0
       operator: constant array operator with value [...]
     array operator with value [...] and update function ...

   This may take a while but should be possible using the recursive structure of the composite operators.

[MSeeker1340]

Ah and also I'm keeping the "AbstractDiffEqCompositeOperator" name for now. The "Abstract" should be enough to distinguish this from the concrete operator composition type.

[ChrisRackauckas]

All good ideas.

[jlperla]

should probably write a pprint method for the composite operators

A thought: if you associate with each type a "description" and a "default symbol" then you could recursively walk the tree and display as a typical linear non-commutative structure. For example,

• DiffEqArrayOperator -> C and description differential equation array
• Number -> c and scalar constant
• DifferentialOperator (or whatever for the derivative operators) -> D and differential operator

alpha = DiffEqScalar(2.0)
A = DiffEqArrayOperator(rand(2,2))
A2 = DiffEqArrayOperator(rand(2,2))
update_func = ...
B = DiffEqArrayOperator(rand(2,2); update_func=update_func)
D = DifferentialOperator(1, :forward) #or whatever those are called/
L = alpha * (A + A2) + B + D

Then you could traverse the tree (numbering them with subscripts arbitrarily) to display this as

Linear combination with L = c_1 (C_1 + C_2) + C_3 + D_1
where
c_i: scalar constant
C_i: differential equation array
D_1: differential operator

[MSeeker1340]

A thought: if you associate with each type a "description" and a "default symbol" then you could recursively walk the tree and display as a typical linear non-commutative structure.

I like this idea! I think a further improvement could be to use multi-index the deeper the tree goes, so in your example C_1 -> C_11 and C_2 -> C_12.

Definitely want to try my hands on this. Seems to be a good change of style from what I've been doing the past few weeks ;)

[MSeeker1340]

I still need to add some documentation and tests but this should cover all the basic compositions. In particular all arithmetics operations on the operators (including left division, mul! and ldiv! via FactorizedDiffEqArrayOperator) is now supported. For composition I also overloaded the ∘ operator because it's equivalent to * in the linear operator setting.

Like I said about the decision to remove the decision in DiffEqLinearCombination and instead use DiffEqScaledOperator, none of the composite operators types now include data of its own other than the component operators. This allows me to create a psuedo-iterator interface in the form of getops, which can be used by update_coefficients! and similar functions (e.g. the pprint method).

(By the way, in v0.7 update_coefficients!(L,u,p,t) will be more versatile because we can now pass p as a named tuple. If a complex operator is composed of several components, all of which have their own parameter sets, we can simply pass in all the parameters and let the components access their parameters by name. In the past this is not possible because the components need to know about the exact position of the parameters in p, or use a custom struct for the purpose of emulating a named tuple).

[ChrisRackauckas]

(By the way, in v0.7 update_coefficients!(L,u,p,t) will be more versatile because we can now pass p as a named tuple. If a complex operator is composed of several components, all of which have their own parameter sets, we can simply pass in all the parameters and let the components access their parameters by name. In the past this is not possible because the components need to know about the exact position of the parameters in p, or use a custom struct for the purpose of emulating a named tuple).

Very good point. That will likely make a lot of things easier.

[jlperla]

For named tuples: If you are creating a lot of variations on parameters, where having defaults would be useful, then note that Parameters.jl has upcoming support: https://mauro3.github.io/Parameters.jl/latest/manual.html#Named-Tuple-Support-1

It is pretty slick and works with 0.6 and 0.7. I suspect that Mauro will release something soon if you need it. I also think that it is worth being dependent on the package just to use that and the @unpack feature for named tuples, which is a great combination.

[ChrisRackauckas]

Will we need this in DiffEqBase to be the default mass matrix?

[MSeeker1340]

I think so.

[MSeeker1340]

Tests are now finished. @ChrisRackauckas do you think it's ok to close this PR for now and create a new release for DiffEqOperators? I know there are still a ton of depwarns for the derivative operators part but I want to focus on jacobians in OrdinaryDiffEq for the moment, the basis of which should be covered by this PR.

[ChrisRackauckas]

Yup don't worry about the depwarns. I'll get femtocleaner to work through them.