Rethink maps

[daanhb]

We are converging to Domain as an interface: it is any object type that supports eltype and in. We make methods with which one can combine domains, where the individual domains don't have to be Domain types. This is clearly very flexible. Perhaps we can do the same for maps and make the framework more generic/applicable?

A map is really just a function: y = f(x). Thus, perhaps any type that supports the calling operator can be seen as a map. Or perhaps any type that supports * can be a map. We could treat a matrix A as a map, the map then being the matrix-vector product y = A*x. Even a number could be a map: 2 represents multiplication by 2. As long as * is associative (which I think Julia already assumes in many cases), the interpretation of * as applying a map is consistent with composition: (A*B)*x = A*(B*x). With this approach, I think we would not even have to implement affine maps. We can just use the existing linear algebra functions and vectors/matrices.

Julia functions have a known output type if you give it a (combination of) input type(s). The current definition of AbstractMap{S,T} really is like a typed function. It corresponds to a function for which you restrict the input type to be S. The output type T is perhaps superfluous. From this point of view, a better name might be TypedMap{S} or TypedFunction{S}. It could have an implementation like:

(*)(f::TypedFunction{S}, x) where {S} = *(f, convert(S, x))

Any concrete subtype would implement multiplication for arguments of type S only.

In order to generically support maps like 2, cos and rand(2,2), it is somewhat restrictive to assume that the input type is known. We can in general assume that, given the map and the input type, the output type is known. (Worst case it is Any.) Thus, we could have TypedMap{S} <: AbstractMap where AbstractMap has no type parameters at all. What the map does, depends on the type of the input and is expressed by implementing (*)(m::MyMap, x).

Ideally, eventually, computing Jacobians etcetera is done generically via automatic differentiation. Composition of maps ideally would use fusion (like broadcast), rather than layered CompositeMap's or CompositeMap's with many elements.

As with domains, as long as you are manipulating variables of type AbstractMap we can use standard function names. If you intend to manipulate any object that implements the map interface, whatever we decide it to be, the intention has to be implicit in the name of the function that is being called (like convert_map and promote_map) or the operation has to be consistent with the generic meaning of the function name (as with *).

I'm sure there'll be some problems I'm not seeing at the moment too.

[dlfivefifty]

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.

[daanhb]

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.

[dlfivefifty]

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...

[daanhb]

I completely agree here.

[daanhb]

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.)

[daanhb]

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.

[dlfivefifty]

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.

[daanhb]

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.

[daanhb]

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).