PDERoadmap
PDERoadmap copied to clipboard
SciML
Discussion on design of domain and grid interface
This issue isolates and continues discussions in https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/260 with @dlfivefifty and @ChrisRackauckas on a joint design of the interface for domains, etc. between the ApproxFun, DiffEqOperators, and ModelingToolkit libraries. We hope to implement it in the prototype: https://github.com/JuliaDiffEq/PDERoadmap/pull/1
Also see https://github.com/JuliaDiffEq/ModelingToolkit.jl/blob/master/test/domains.jl
For a starting point in one dimension,
d = Interval(1.0, 10.0) #Continuous domain
x = RegularGrid(d, N)
Note: I think it is better to give points than dx, but I am open. The reason is that if you give a dx then you never know how well things line up, or how many points it ends up creating.
When things move to an irregular grid,
d = Interval(1.0, 10.0) #Continuous domain
grid_points = 1.0:0.01:10.0 #This is regular, but you get the point.
x = IrregularGrid(d, grid_points)
For the product of a discrete domain with 2 values and continuous domain (which comes up a lot in economics applications),
d = DiscreteDomain(2) * Interval(0.0,5.0) #Or tensor product?
Or, when we later want to enable naming of types for the discrete domain (e.g. for a continuous-time markov chain of switching between L and H productivity in a model)
d1 = Interval(1.0, 2.0)
d2 = DiscreteDomain([:L, :H])
d = d1 ⊗ d2 #Is the tensor product the right one to use here?
Or, instead, does it make sense to use an enumeration? I really don't like how scoping works with enumerations in Julia, which limits the clarity of the code.
To get @MSeeker1340 up to speed, this is for being able to define the operators without the discretization. DiffEqOperators is all about how do you build the computational operator to compute, but there's a higher level question we should be asking as to, how can I just symbolically/mathematically specify some differential operator and have the correct computational operators be generated? At this level, it's not a finite difference, FEM, spectral, etc. discretization: just some operator on a domain which would make it easy to switch between forms (if we get it right).
ModelingToolkit.jl is an intermediate representation (IR) for models. You can think of it almost as a mini-CAS that DSLs can write into and compute on. The variables of this IR can hold domains (https://github.com/JuliaDiffEq/ModelingToolkit.jl/blob/master/src/variables.jl#L10) and there are differential operators. We can modify it at will to make the IR have what we need to generate operators.
However, I think we should get some operators working first before we talk about auto-generating. In some sense, this is the "FEniCS" part and we need dolphin first.
At this level, it's not a finite difference, FEM, spectral, etc. discretization: just some operator on a domain which would make it easy to switch between forms (if we get it right).
Well this is basically symbolic math right? Maybe we can get some inspiration from Mathematica? (I don't have much experience with it though).
The type system of lazy operator composition/combination already has a strong resemblance to a parse tree. So I wouldn't be surprised if we can write something that translates an abstract operator defined in ModelingToolkit to a concrete one that shares the same tree structure.
But yes I agree we should stick to simple finite-difference realizations of the operators at the moment, while maybe think a bit about incorporating ApproxFun.
Maybe we can get some inspiration from Mathematica?
We're building this with help from REDUCE.
But yes I agree we should stick to simple finite-difference realizations of the operators at the moment, while maybe think a bit about incorporating ApproxFun.
Yes, first things first.
The tree-vs-type hierarchy question is important to get right. In ApproxFun I overused templated variables, which puts a limit on how complicated the operator tree can be. However, building operators should be fast. It’s not clear to me what the right balance is.
The tree-vs-type hierarchy question is important to get right. In ApproxFun I overused templated variables, which puts a limit on how complicated the operator tree can be. However, building operators should be fast. It’s not clear to me what the right balance is.
We did it a little wrong once. https://github.com/JuliaDiffEq/ModelingToolkit.jl/issues/62 . I think the key is to separate "compile time" and "runtime" here. The resulting operator can be fully typed, just the computation to get it may not be inferred. That's fine because building the operator would either be at the top scope or above a function barrier, i.e. building it is a form of compilation so like how the Julia AST itself is untyped, the specification of the operator should be fully untyped but builds strongly typed code. We'll at least give it a try.
I understand the design concerns down the road when we might want to mess around with the IR, but does that change the basic types for writing the domains and grids?
In particular, isn't the design of https://github.com/JuliaDiffEq/PDERoadmap/issues/4#issue-318138562 style interfaces largely just going to help with dispatch and the user's API?
I understand the design concerns down the road when we might want to mess around with the IR, but does that change the basic types for writing the domains and grids?
Nope it shouldn't. It just carries those along but doesn't strictly use them. We can use them to do things like generate ODEs in local coordinates on known manifolds, but in terms of the IR its mostly for storing information.
@dlfivefifty @ChrisRackauckas What are the thoughts on using Interval(...) here, exactly as in ApproxFun? vs. 0.0:1.0 vs. 0.0 .. 1.0?
I was experimenting with instead using the return from the range built into Julia (i.e. 1.0:10.0) but
julia> 1.0:10.0 |> typeof
StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
julia> 1.0:10.0
1.0:1.0:10.0
As much as it would be great in notation to support defining intervals with 0.0:1.0, I do not think this usage scenario fits the Range concept, which is inherently iterable.
The better fit is probably https://github.com/JuliaIntervals/IntervalRootFinding.jl What about using the .. operator?:
julia> using IntervalArithmetic
julia> 1.0..10.0 |> typeof
IntervalArithmetic.Interval{Float64}
julia> 1.0..10.0
[1, 10]
Are there other precedents for using this notation for intervals?
a .. b is the agreed upon convention for intervals, see IntervalSets.jl
Abusing a UnitRange is a bad idea.
Perfect! a .. b as sugar to create the Interval it is. Though this tells me that ApproxFun and IntervalArithmetic are incompatible, and explains the error I saw. Should some sort of shared interval base library be established?
(Even in the longrun, with my ADL suggestion in https://discourse.julialang.org/t/function-name-conflict-adl-function-merging/10335 wouldn't help here, since the operator .. needs to work on two floats in both cases)
For the grids, are there any other libraries that would be nice to maintain compatibility with? For https://github.com/JuliaMath/Interpolations.jl it seems that they use the StepRangeLen type directly for the grid, but I wonder if we are better off maintaining a connection from a grid to a domain.
Also https://github.com/sisl/GridInterpolations.jl has its own grid objects
The whole point of IntervalSets.jl is to be a shared package for a .. b .
Done. What about grids? Any movement towards a coordinated library for those types?
Not on METADATA but we’ve been working on https://github.com/daanhb/Domains.jl
Unfortunately it was triggering a compile-time Julia bug when used in ApproxFun, which put a hold on development. We also decided “domain” should be an interface rather than an abstract type, and anything that supports in can be treated as a domain.
@dlfivefifty What are your thoughts on domains with discrete numbers of elements? Did you have that in mind with Domains.jl? In particular, having them "named"? I don't think using a enumeration is a good idea due to namespace clashes, but symbols work well.
@ChrisRackauckas Thoughts on getting uniform grids up and running for the prototype in DiffEqOperators.jl ? I think that IntervalSets.jl as a dependency is an easy addition? A temporary type to hold grids may be reasonable until things we are able to coordinate on a central grid type.
Making domain an interface has the benefit that, since 1 in 1 works, numbers can be interpreted as points and vectors can be interpreted as grids. So in 1D I’d just use linspace. For 2D one would want an array of SVector I think.
What is the connection between domains, intervals, etc.? It sounds like you are saying that the domain is the discretized grid that the function is defined on. Is that correct?
If so, then what do we need for the generic algorithms to work:
inis necessary, as you say- Some testing of equality is needed?
- The reason is ultimately to be able to ensure that operators we are composing are defined on the same domains, and potentially expand to the superset of the domain
- Don't we need iteration as well for generic algorithms to be implemented on the grids themselves? If so, then https://docs.julialang.org/en/stable/manual/interfaces/#man-interface-iteration-1 would also be a requirement of a domain?
Where I start to get confused in the interface, is how this would all work for algorithms with multiple dimensions. I see how the in works, but how do you write generic code? I think this becomes especially important with designing the operator interface, because I presume the iteration is used to go through the state with the update_coefficients! design.
In my usage domain is a possibly infinite set. So all of these would implement the domain interface:
a..b
a:b
1
[1,2,3]
Set([1,2,3])
The first is infinite (an interval) the rest are discrete.
Iteration is not possible with infinite sets but is for any of the discrete ones.
To generalise to higher dimensions, I interpret the “ambient space” of the domain as all elements of the type, and would use something like SVector{2,T} to represent 2D points.
@daanhb May have some thoughts.
OK, then for the "named" discrete sets, vectors of symbols would work perfectly.
[:L, :H]
Thinking about your domain, the generality makes complete sense for spectral methods. I guess the question with finite-difference methods is whether there is need to have separate generic types to for the continuous domain and the discretization of the continuous domain. That was the assumption I had with my suggested
N = 20
d = 1.0..10.0
x = RegularGrid(d, N)
But, is the d really needed to keep around, or is the x enough? For example, we could just go
x = 1.0:0.1:10.0 #For a regular grid
x = linspace(1.0, 10.0, 20) #Generates the same type
x_irregular = [1.0, 1.1, 1.5, 2.0, 5.0, 10.0] #Irregular grid
One worry in that approach would be whether you can have specialized code for regular vs. irregular grids. But I think it works completely fine, since the type of the range and the vector are different and it could dispatch on StepRangeLen. That is
julia> x = linspace(0.0, 1.0, 10)
0.0:0.1111111111111111:1.0
julia> typeof(x) <: StepRangeLen
true
julia> typeof(collect(x)) <: StepRangeLen
false
julia> step(x)
0.1111111111111111
Yes yes yes and yes.
Domains as an interface with in, iteration, array interface, and SVectors or scalars. SVector is important so that way (u,p,t,x) is universal for N-dimensional. We're only talking about grids now, continuous domains are for a higher level.
This isn't pertinent right now though. We can add grids to the operators quite easily once we have that setup. It will be important when defining the irregular grids, but it seems that interface works with arrays of scalars / SVector so it's fine.
We're only talking about grids now, continuous domains are for a higher level. OK, then no need to have the separation of the continuous domain and the discretization of it.
Some specific quetsion:
- So, lets say we are composing some operators. Would we then have then keep the "grid" with each of them (which could be used to ensure that the composed operators are defined on the same grid type). That is,
x = linspace(0.0, 1.0, 100) #100 node grid
sigma_2 = 1.1
mu = -0.2
L = mu * DriftOperator(x, :forwards) +sigma_2 * DiffusionOperator(x)
But what if there are the functions for the coefficients? Is the idea that the "functions" would conform to whatever the "grid"?
x = linspace(0.0, 1.0, 100) #100 node grid
sigma_2(u,p,t,x) = 1.1 * x^2
mu(u,p,t,x) = -0.2 * x
L = mu * DriftOperator(x, :forwards) +sigma_2 * DiffusionOperator(x)
- Alternatively, if we want the grid associated with every operator type, there could be a wrapper, as we had discussed,
x = linspace(0.0, 1.0, 100) #100 node grid
sigma_2_tilde(u,p,t,x) = 1.1 * x^2
mu_tilde(u,p,t,x) = -0.2 * x
sigma_2 = DiffEqFunction(sigma_2_tilde, x)
mu= DiffEqFunction(mu_tilde, x)
L = mu * DriftOperator(x, :forwards) +sigma_2 * DiffusionOperator(x)
SVector is important so that way (u,p,t,x) is universal for N-dimensional. I am not sure I understand you here. Are you saying that
xin the interface would be aSVector? I would think a tuple is more appropriate, especially when there are discrete states inx? (Yes, I knowSVectoris implemented as a tuple, but it is a tuple where everything is of the same type).
Although I would love it if for multiple dimensions we could use a named tuple. In particular, the canonical case in economics is a discrete i in {L,H} state which determines the drift of a continuous state, y. For example, if the drift is 0.1 * y if i = :L and 0.5 * y if i=:H, I would love to be able to write this as something like:
julia> x = @NT(y=0.1, i=:L)
(y = 0.1, i = :L)
julia> μ(u,p,t,x) = x[:i]==:L ? 0.1*x[:y] : 0.5*x[:y]
μ (generic function with 1 method)
julia> μ(0,0,0,x)
0.010000000000000002
julia> μ(0,0,0,@NT(y=0.5, i=:H))
0.25
(where the @NT is obviously dropped in v0.7
Although I would love it if for multiple dimensions we could use a named tuple. In particular, the canonical case in economics is a discrete i in {L,H} state which determines the drift of a continuous state, y. For example, if the drift is 0.1 * y if i = :L and 0.5 * y if i=:H, I would love to be able to write this as something like:
NamedTuples satisfy the interface. You can do that if you want.
But what if there are the functions for the coefficients? Is the idea that the "functions" would conform to whatever the "grid"?
Yes, it would just be on the user side the x is of the right dimension. It would be like an FEM node listing.
NamedTuples satisfy the interface. You can do that if you want.
Is it that simple? Presumably the interface of what is passed into the x of the function has something to do with the types of the dimensions of the grid? Forgetting about the named tuples for a second, it seems like in the construction of the tuple itself the underlying code in the library is going to have to take the Cartesian product of all of the elements of the types?
It seems like if I want to have an operator which works on the product of a continuous and discrete grid, I would have to put it in a tuple, or something like that:
domain = ( linspace(0.0, 1.0, 10), [:L, :H] ) #A tuple of the domain
...or add syntactic sugar for the cartesian product. Then generic code calling the functions in the library could construct a tuple to call the underlying functions of type Tuple{eltype(domain[1]), eltype(domain[2])}
And in order to have names, wouldn't we need to have those named associated with the underlying domain. That is, something like AxisArrays?
i_domain = Axis{:i}([:L, :H])
y_domain= Axis{:y}(1.0:0.1:2.0)
That way, a named tuple could be created to call the user functions by getting the names of each axis with AxisArrays.axisname(i_domain) etc. and then the type with eltype(i_domain), etc.
Non-lazy Cartesian products of grids would be SVector.(x, y’). Lazy could be done with BroadcastArray currently in InfiniteArrays.jl, or via a special type.
OK, but
julia> x = collect(linspace(0.0, 1.0, 5))
5-element Array{Float64,1}:
0.0
0.25
0.5
0.75
1.0
julia> y = [:L, :H]
2-element Array{Symbol,1}:
:L
:H
julia> z = [1, 2]
2-element Array{Int64,1}:
1
2
julia> SVector.(x, z')
5×2 Array{StaticArrays.SArray{Tuple{2},Float64,1,2},2}:
[0.0, 1.0] [0.0, 2.0]
[0.25, 1.0] [0.25, 2.0]
[0.5, 1.0] [0.5, 2.0]
[0.75, 1.0] [0.75, 2.0]
[1.0, 1.0] [1.0, 2.0]
julia> SVector.(x, y')
ERROR: MethodError: no method matching transpose(::Symbol)
Basically, this approach of the cartesian products seems to promote discrete dimensions, and doesn't allow symbols?
If I add in a transpose,
julia> Base.transpose(x::Symbol) = x
julia> SVector.(x, y')
5×2 Array{StaticArrays.SArray{Tuple{2},Any,1,2},2}:
Any[0.0, :L] Any[0.0, :H]
Any[0.25, :L] Any[0.25, :H]
Any[0.5, :L] Any[0.5, :H]
Any[0.75, :L] Any[0.75, :H]
Any[1.0, :L] Any[1.0, :H]
It looks the "promotion" turns these into "Any", which is not what we want.
Isn't the basic problem that SArray are designed to hold a single type?
OK. So
julia> using IterTools
julia> grids = (1.0:0.1:2.0, [:L, :H])
(1.0:0.1:2.0, Symbol[:L, :H])
julia> [i for i in product(grids...)]
22-element Array{Tuple{Float64,Symbol},1}:
(1.0, :L)
(1.1, :L)
(1.2, :L)
(1.3, :L)
...
So at least IterTools could be used to get an enumeration on the cartesian problem? I suppose something similar with a Axis generating a NamedTuple would also be possible?
Yes, any of these AbstractVector with scalars or AbstractVectors are a good way to specify the points of a domain. It doesn't matter which you choose, though being able to determine dx will matter for operators.