Rethink maps about domainsets.jl HOT 9 CLOSED

Maps as an interface is definitely the right approach. We would probably still need a special type for affine maps however as x -> (A*x) has to recompile for each A.

I think using composition notion ∘ is preferable to * as maps are not naturally an "algebra", and also, we want "right-associativity" (see https://discourse.julialang.org/t/why-is-multiplication-a-b-c-left-associative-foldl-not-right-associative-foldr/17552)

This question goes hand-in-hand with the question of broadcast and domains. That is, can we write exp(im*(0..1)) or should we require exp.(im .* (0..1)). Broadcast notation should be the "default" as its the easiest to implement, and support for other operations would be implemented by lowering to their broadcast variants, as in base (e.g., +(a::AbstractArray, b::AbstractArray) calls broadcast(+, a, b).)

I believe there is a package out there that has something like TypedFunction, but I couldn't find it.

from domainsets.jl.

Okay, good points! Let's use ∘ then for composition. I'm also very much in favour of using broadcast notation to easily make mapped domains. Perhaps that could even lead to fusion, but I have little experience with the broadcast framework so far.

About affine maps, I was thinking of the following:

struct CompositeMap{A} <: AbstractMap
    maps  :: A
end

where maps would be, say, a tuple of an array and an int. If we can not get this to work without having to define x -> A*x then I agree it is not worth it.

Perhaps it depends on what we define the operation of applying a map to be. If the syntax is y = f(x), it would be awkward to define the calling operation on abstract arrays in this package (even impossible because the type has to be concrete?). The syntax * has disadvantages too, because then we would like to define cos * 0.5 to be cos(0.5).

On the other hand, if we continue to use y = applymap(A, x), then we could define applymap for arrays, numbers, etcetera. Or we could have two fallbacks applymap(A, x) = A*x and applymap(f::Function, x) = f(x). The interface of what we consider to be a map would include an implementation of applymap.

As an interface and from the point of view that maps are functions, y = f(x) is of course nicer syntax.

I'm talking in thin air now, perhaps some exploring is called for.

from domainsets.jl.

My current feelings are:

1. * should be reserved only for algebras.
2. Operators and functions are not naturally an algebra. Operators are "functions of functions".
3. So ApproxFun's notation D*f is wrong: it should be D(f). and operator composition should be D ∘ D, not D*D. Note that this is the same change Chebfun did from Linop to Chebop.
4. Operators and functions are equivalent to (Quasi-)Matrices and (Quasi-)vectors, that do form an algebra and should support *. But "functions as functions" and "functions as vectors" should be two different types.

How does this fit into maps? Something like a rotation is both a function and a matrix, depending on what it's acting on. I guess I'm arguing these should be different types...

from domainsets.jl.

I completely agree here.

from domainsets.jl.

So, this would lead to using the y=f(x) syntax for calling a map, which is the most natural one, and wrapping matrices and numbers into a map type, is that right?

For Domains, each domain is free to implement in. However, in is also implemented at the general level of Domain{T} with some checks and conversions, calling indomain for concrete domains who can then assume that the argument is of type T and need not worry about conversions. We can't do this for the calling operation on an abstract Map type. (Though perhaps we don't want to.)

from domainsets.jl.

While we are rethinking maps, we may want to rethink MappedDomain as well. The only sensible definition of a mapped domain in DomainSets.jl is where you think of the mapped domain as the image of the original domain under a map. In particular, it should hold regardless of any distribution of the points within the domain or of any metric in the space involved. Thus, mapping 0..1 using y=x^2 should equal 0..1.

The equality of a MappedDomain with a concrete domain that represents the same image will of course be difficult to check in general. However, in all my current use cases of mapped domains, I am actually interested in retaining the map, since its jacobian is of interest. That is, I'm implicitly assuming the map alters some corresponding metric space. The latter clearly has no place in DomainSets.jl. (Or does it?)

Currently, there is no specific problem, but the two use-cases are incompatible and it may be good to keep the distinction in mind.

from domainsets.jl.

I think this needs to go hand in hand with broadcasting, and essentially the storing the map means a lazy broadcast. I think we could mirror the relationship between broadcast and BroadcastArray in LazyArrays.jl. That is, broadcast(+, 1:3, 1) returns 2:4 (and forgets the map) while BroadcastArray(+, 1:3, 1) is left “unmaterialised” and thereby remembers the map.

So I guess I’d rename MappedDomain to BroadcastDomain . Broadcasting is better than mapping since it’s more type stable, consider the arc:

BroadcastDomain(+, BroadcastDomain(exp, im*(0..1)), 2+im)

One issue is to support in we also need a maps inverse. Probably we could do this cleanly via a broadcastinv machinery.

from domainsets.jl.

Sounds sensible, it is a type of lazy domain. Most importantly, a user who is interested in retaining the map needs a safe way to store the map and the original domain. It is still somewhat fishy to do that in DomainSets itself, since it remains the case that as far as in is concerned 0..1 and BroadcastDomain( x->x^2, 0..1) are the same. If there is ever a package that has metric spaces, it should perhaps be in there.

from domainsets.jl.

I think this remains the largest open issue with DomainSets: to remove maps from the package itself, enable the use of various other packages that define maps, and to come up with a convenient interface for making mapped domains, for example using broadcast and lazy domains as outlined above.

The definition of a MappedDomain is already largely independent of what a map exactly is. Specific requirements for mapped domains are that ideally we know both the map and its inverse, and we want to support maps that are not fully invertible such as embeddings (i.e. a map from R to R2 to represent a line segment in the plane).

from domainsets.jl.