Improve maps

[daanhb]

This is a work in progress to improve maps. It brings us closer to splitting off maps and reusing or merging with other packages - there is still a lot of overlap. The goal should probably remain to complete this code elsewhere, separate from DomainSets.

The main hierarchy is now: AbstractMap :> Map{T} :> TypedMap{T,U}. A Map{T} works like a Domain{T}, i.e., T plays a similar role. Most functionality assumes a Map{T}. We used to have AbstractMap{S,T}.

DomainSets is functionally equivalent to what it was before. Maps can now also take arbitrary vectors, and the functionality of jacobians was improved a tiny bit.

TODO: support nicer syntax, like 2 .* (0..1.0) to make a mapped domain.

[daanhb]

This touches on #21, #23 and #36.

The AbstractMap supertype is now close to the Transformation supertype in CoordinateTransformations.jl but Map{T} is different. I think Map{T} is a better match for Domain{T}, though a mapped domain typically stores the inverse map so the two T parameters involved may differ anyway.

[daanhb]

The new maps are not perfect, but they have certainly improved and they are simpler than before. Most maps have a proper jacobian - the jacobian of a map is again a map. There is some support even for jacobians of composite maps.

There is now also support for broadcast notation:

julia> using DomainSets

julia> 2*(0..1) .+ 1.5
1.5..3.5

julia> 2UnitDisk() .+ [1,2]
A mapped domain based on the 2-dimensional closed unit ball

The changes seem compatible with ApproxFun - at least the test suite passes on master.

[daanhb]

Should 2*d be allowed, or should it be 2 .* d?

[dlfivefifty]

If we follow the Set convention then....No, we should only support broadcast:

julia> 2Set([1,2,3])
ERROR: MethodError: no method matching *(::Int64, ::Set{Int64})
Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...) at operators.jl:529
  *(::T, ::T) where T<:Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8} at int.jl:54
  *(::Union{Int16, Int32, Int64, Int8}, ::BigInt) at gmp.jl:535
  ...
Stacktrace:
 [1] top-level scope at REPL[1]:1

julia> 2 .* Set([1,2,3])
3-element Array{Int64,1}:
 4
 6
 2

Though I'm actually surprised by this behaviour TBH, I'm wondering if it's just an oversight.

[daanhb]

Agree, so I implemented broadcast for multiplication as well. I had to implement that anyway. It's not really nice syntax though, since 2UnitDisk() now becomes 2 .* UnitDisk(). But it is consistent, and the broadcast framework is much more general than the various special cases we had before, so it is good to have. It would be easy to allow 2 UnitDisk() again if we want.

It is a bit of a pain to support older versions of Julia, since I had decided to use m(x) as syntax for applying a map to the argument x. A map is simply a function, after all, and using this syntax makes it easier to mix Map{T} with regular functions at some point later on. But the calling syntax can't be generically implemented on the abstract AbstractMap type. For now, I kept the syntax and added tests on the version of Julia in various places (as in VERSION > v"1.2").

[dlfivefifty]

I think it’s fine to drop earlier versions of Julia

[daanhb]

It was easy enough to add, but ugly - I'll take it out again later, I'm curious to see if the tests now pass first.

[dlfivefifty]

It would be easy to allow 2 UnitDisk() again if we want.

I think allow it since it’s convenient, and replicates Vector behaviour

[daanhb]

Okay, multiplication with a number is allowed again, and so is unary minus (that had become (-).(d) using broadcast instead of -d, yikes). Also, @vincentcp spotted a dangerous bug in the equality of mapped domains (the maps were not compared).

Types are promoted more carefully now:

julia> using DomainSets

julia> big(2) * (0..1)
0.0..2.0

julia> eltype(ans)
BigFloat

I don't think that was the case before.

No doubt many things will be wrong for types I did not think of and domains I did not test :-) But I think I'll stop here for now. Next steps (apart from testing) would be to take out the maps altogether, and possibly merge with other packages. The better maps will hopefully help further development of the new DomainIntegrals package, which I recently added to JuliaApproximation, for the (adaptive) computation of integrals on domains.

[daanhb]

If the tests pass, I'll merge this and tag 0.4. Tests pass with the ApproxFun packages: nothing really changed for any domain, it is only the maps. The last additions enable a jacobian for product maps. The affine mapping types are a bit numerous now, because I made some special cases in order to avoid having three type parameters even for a simple scalar affine map. Compared to some other packages, these maps can be rectangular, supporting a left and/or right inverse.

The inverse/jacobian/leftinverse/rightinverse of a map is again a map, but there is also a two-argument function leftinverse(mymap, x) which evaluates the left inverse at the point x without creating a new map in between (this can sometimes be implemented more efficiently, e.g. by solving a linear system rather than computing a matrix inverse). Of course the default is leftinverse(mymap)(x).

We can tune the construction and syntax later. For now, a linear map is made by LinearMap(rand(2,2)). This yields a VectorLinearMap, which stores the A in y=A*x by a dense array. An affine map of the form y=A*x+b is created using AffineMap(A, b). The type you get depends on the types of A and b. In the case of scalars, you'll get a ScalarAffineMap{T}. In the most generic case, you get a GenericAffineMap{T,AA,B} which can store any A and any b.