type promotion for composite maps

Author: daanhb

Project.toml

@@ -1,7 +1,7 @@
 name = "FunctionMaps"
 uuid = "a85aefff-f8ca-4649-a888-c8e5398bc76c"
 authors = ["Daan Huybrechs <[email protected]> and contributors"]
-version = "0.1.1"
+version = "0.1.2"
 
 [deps]
 CompositeTypes = "b152e2b5-7a66-4b01-a709-34e65c35f657"

src/FunctionMaps.jl

@@ -12,6 +12,9 @@ import CompositeTypes: component, components
 
 # Exhaustive list of exports:
 
+# from CompositeTypes
+export component, components, ncomponents
+
 # from util/common.jl
 export prectype, numtype
 

src/concrete/coordinates.jl

@@ -1,4 +1,3 @@
-
 # This file contains a few specific maps, included mainly to test
 # functionality on non-trivial maps
 
@@ -12,6 +11,9 @@ end
 
 CartToPolarMap() = CartToPolarMap{Float64}()
 
+similarmap(::CartToPolarMap, ::Type{SVector{2,T}}) where {T<:Number} = CartToPolarMap{T}()
+isequalmap(m1::CartToPolarMap, m2::CartToPolarMap) = true
+
 mapsize(m::CartToPolarMap) = (2,2)
 
 applymap(map::CartToPolarMap{T}, x) where {T} =
@@ -28,10 +30,6 @@ inverse(m::CartToPolarMap, x) = inverse(m)(x)
 
 isrealmap(m::CartToPolarMap) = true
 
-convert(::Type{Map{SVector{2,T}}}, ::CartToPolarMap) where {T} = CartToPolarMap{T}()
-
-isequalmap(m1::CartToPolarMap, m2::CartToPolarMap) = true
-
 
 """
 A Polar to Cartesian map. The angle is mapped to the second dimension,
@@ -43,6 +41,9 @@ end
 
 PolarToCartMap() = PolarToCartMap{Float64}()
 
+similarmap(::PolarToCartMap, ::Type{SVector{2,T}}) where {T<:Number} = PolarToCartMap{T}()
+isequalmap(m1::PolarToCartMap, m2::PolarToCartMap) = true
+
 mapsize(m::PolarToCartMap) = (2,2)
 
 applymap(map::PolarToCartMap{T}, x) where {T} = SVector{2,T}((x[1]+1)/2*cos(pi*x[2]), (x[1]+1)/2*sin(pi*x[2]))
@@ -55,10 +56,6 @@ inverse(m::PolarToCartMap, x) = inverse(m)(x)
 
 isrealmap(m::PolarToCartMap) = true
 
-convert(::Type{Map{SVector{2,T}}}, ::PolarToCartMap) where {T} = PolarToCartMap{T}()
-
-isequalmap(m1::PolarToCartMap, m2::PolarToCartMap) = true
-
 
 """
 The map `[cos(2πt), sin(2πt)]` from `[0,1]` to the unit circle in `ℝ^2`.
@@ -67,6 +64,9 @@ struct UnitCircleMap{T} <: Map{T} end
 
 UnitCircleMap() = UnitCircleMap{Float64}()
 
+similarmap(::UnitCircleMap, ::Type{T}) where {T<:Number} = UnitCircleMap{T}()
+isequalmap(::UnitCircleMap, ::UnitCircleMap) = true
+
 mapsize(m::UnitCircleMap) = (2,)
 
 applymap(m::UnitCircleMap{T}, t) where {T} = SVector(cos(2*T(pi)*t), sin(2*T(pi)*t))
@@ -89,6 +89,9 @@ end
 
 AngleMap() = AngleMap{Float64}()
 
+similarmap(::AngleMap, ::Type{SVector{2,T}}) where {T<:Number} = AngleMap{T}()
+isequalmap(::AngleMap, ::AngleMap) = true
+
 function applymap(m::AngleMap{T}, x) where {T}
     twopi = 2*convert(T, pi)
     θ = atan(x[2],x[1])
@@ -124,6 +127,9 @@ struct UnitDiskMap{T} <: Map{SVector{2,T}} end
 
 UnitDiskMap() = UnitDiskMap{Float64}()
 
+similarmap(::UnitDiskMap, ::Type{SVector{2,T}}) where {T<:Number} = UnitDiskMap{T}()
+isequalmap(::UnitDiskMap, ::UnitDiskMap) = true
+
 mapsize(m::UnitDiskMap) = (2,2)
 
 applymap(m::UnitDiskMap{T}, x) where {T} =

src/generic/composite.jl

@@ -1,16 +1,41 @@
+# ensure that the domaintype of each map in the sequence agrees with the codomain type of the
+# preceding map
+function match_domain_codomain_types(::Type{T}, map, maps...) where {T}
+    Y = promote_type(T, domaintype(map))
+    m1 = convert_domaintype(Y, map)
+    (m1, match_domain_codomain_types(codomaintype(m1), maps...)...)
+end
+
+match_domain_codomain_types(::Type{T}, map) where {T} = (convert_domaintype(T, map),)
+
+"""
+    ComposedMap{T,MAPS}
 
-"The composition of several maps."
+The composition of several maps.
+
+The `components` of a `ComposedMap` are the maps in the order that they are applied
+to the input.
+"""
 struct ComposedMap{T,MAPS} <: CompositeLazyMap{T}
     maps    ::  MAPS
 end
 
-ComposedMap(maps...) = ComposedMap{domaintype(maps[1])}(maps...)
-ComposedMap{T}(maps...) where {T} = ComposedMap{T,typeof(maps)}(maps)
+function ComposedMap(map1, maps...)
+    P = prectype(map1, maps...)
+    if P == Any
+        # don't try to promote types
+        _ComposedMap(domaintype(map1), map1, maps...)
+    else
+        T = to_prectype(P, domaintype(map1))
+        _ComposedMap(T, match_domain_codomain_types(T, map1, maps...)...)
+    end
+end
+ComposedMap{T}(map1, maps...) where {T} = _ComposedMap(T, match_domain_codomain_types(T, map1, maps...)...)
+_ComposedMap(::Type{T}, maps...) where {T} = ComposedMap{T,typeof(maps)}(maps)
 
-# TODO: make proper conversion
-similarmap(m::ComposedMap, ::Type{T}) where {T} = ComposedMap{T}(m.maps...)
+similarmap(m::ComposedMap, ::Type{T}) where {T} = ComposedMap{T}(components(m)...)
 
-codomaintype(m::ComposedMap) = codomaintype(m.maps[end])
+codomaintype(m::ComposedMap) = codomaintype(last(m.maps))
 
 # Maps are applied in the order that they appear in m.maps
 applymap(m::ComposedMap, x) = applymap_rec(x, m.maps...)
@@ -19,7 +44,7 @@ applymap_rec(x, map1, maps...) = applymap_rec(applymap(map1, x), maps...)
 
 # The size of a composite map depends on the first and the last map to be applied
 # We check whether they are scalar_to_vector, vector_to_vector, etcetera
-mapsize(m::ComposedMap) = _composed_mapsize(m, m.maps[end], m.maps[1], mapsize(m.maps[end]), mapsize(m.maps[1]))
+mapsize(m::ComposedMap) = _composed_mapsize(m, last(m.maps), first(m.maps), mapsize(last(m.maps)), mapsize(first(m.maps)))
 _composed_mapsize(m, m_end, m1, S_end::Tuple{Int,Int}, S1::Tuple{Int,Int}) = (S_end[1],S1[2])
 _composed_mapsize(m, m_end, m1, S_end::Tuple{Int,Int}, S1::Tuple{Int}) =
     is_vector_to_scalar(m_end) ? () : (S_end[1],)
@@ -109,9 +134,10 @@ struct MulMap{T,MAPS} <: CompositeLazyMap{T}
     maps    ::  MAPS
 end
 
+MulMap(map1::Map, maps::Map...) = MulMap(promote_maps(map1, maps...)...)
 MulMap(map1::Map{T}, maps::Map{T}...) where {T} = MulMap{T}(map1, maps...)
 MulMap{T}(maps::Map{T}...) where {T} = MulMap{T,typeof(maps)}(maps)
-MulMap{T}(maps...) where {T} = _mulmap(T, convert.(Map{T}, maps)...)
+MulMap{T}(maps...) where {T} = _mulmap(T, convert_domaintype.(Ref(T), maps)...)
 _mulmap(::Type{T}, maps...) where {T} = MulMap{T,typeof(maps)}(maps)
 
 similarmap(m::MulMap, ::Type{T}) where {T} = MulMap{T}(m.maps...)
@@ -155,9 +181,10 @@ struct SumMap{T,MAPS} <: CompositeLazyMap{T}
     maps    ::  MAPS
 end
 
+SumMap(map1::Map, maps::Map...) = SumMap(promote_maps(map1, maps...)...)
 SumMap(map1::Map{T}, maps::Map{T}...) where {T} = SumMap{T}(map1, maps...)
 SumMap{T}(maps::Map{T}...) where {T} = SumMap{T,typeof(maps)}(maps)
-SumMap{T}(maps...) where {T} = _summap(T, convert.(Map{T}, maps)...)
+SumMap{T}(maps...) where {T} = _summap(T, convert_domaintype.(Ref(T), maps)...)
 _summap(::Type{T}, maps...) where {T} = SumMap{T,typeof(maps)}(maps)
 
 similarmap(m::SumMap, ::Type{T}) where {T} = SumMap{T}(m.maps...)

src/generic/isomorphism.jl

@@ -34,6 +34,11 @@ See also: [`VectorToNumber`](@ref).
 """
 NumberToVector() = NumberToVector{Float64}()
 
+similarmap(::NumberToVector, ::Type{T}) where {T<:Number} = NumberToVector{T}()
+similarmap(::VectorToNumber, ::Type{SVector{1,T}}) where {T<:Number} = VectorToNumber{T}()
+isequalmap(::NumberToVector, ::NumberToVector) = true
+isequalmap(::VectorToNumber, ::VectorToNumber) = true
+
 mapsize(::VectorToNumber) = (1,1)
 mapsize(::NumberToVector) = (1,)
 
@@ -77,6 +82,11 @@ See also: [`VectorToComplex`](@ref).
 """
 ComplexToVector() = ComplexToVector{Float64}()
 
+similarmap(::ComplexToVector, ::Type{Complex{T}}) where {T<:Number} = ComplexToVector{T}()
+similarmap(::VectorToComplex, ::Type{SVector{2,T}}) where {T<:Number} = VectorToComplex{T}()
+isequalmap(::ComplexToVector, ::ComplexToVector) = true
+isequalmap(::VectorToComplex, ::VectorToComplex) = true
+
 mapsize(::VectorToComplex) = (1,2)
 mapsize(::ComplexToVector) = (2,)
 
@@ -115,6 +125,12 @@ end
 struct FlatToNested{N,T,U,DIM} <: Isomorphism{SVector{N,T},U}
 end
 
+# TODO: not all type combinations make sense
+similarmap(::NestedToFlat{N,T,U,DIM}, ::Type{V}) where {N,T,U,DIM,V} = NestedToFlat{N,T,V,DIM}()
+similarmap(::FlatToNested{N,T,U,DIM}, ::Type{SVector{N,V}}) where {N,T,U,DIM,V<:Number} = FlatToNested{N,V,U,DIM}()
+isequalmap(::NestedToFlat{N,T,U,DIM}, ::NestedToFlat{N,S,V,DIM}) where {N,T,U,DIM,S,V} = true
+isequalmap(::FlatToNested{N,T,U,DIM}, ::FlatToNested{N,S,V,DIM}) where {N,T,U,DIM,S,V} = true
+
 inverse(::NestedToFlat{N,T,U,DIM}) where {N,T,U,DIM} = FlatToNested{N,T,U,DIM}()
 inverse(::FlatToNested{N,T,U,DIM}) where {N,T,U,DIM} = NestedToFlat{N,T,U,DIM}()
 inverse(m::NestedToFlat, x) = inverse(m)(x)

src/generic/jacobian.jl

@@ -1,10 +1,12 @@
-
 "A lazy Jacobian `J` stores a map `m` and returns `J(x) = jacobian(m, x)`."
 struct LazyJacobian{T,M} <: SimpleLazyMap{T}
 	map	::	M
 end
 
-LazyJacobian(m) = LazyJacobian{domaintype(m),typeof(m)}(m)
+LazyJacobian(m) = LazyJacobian{domaintype(m)}(m)
+LazyJacobian{T}(m) where T = _LazyJacobian(T, convert_domaintype(T, m))
+_LazyJacobian(::Type{T}, m) where T = LazyJacobian{T,typeof(m)}(m)
+
 applymap(m::LazyJacobian, x) = jacobian(supermap(m), x)
 
 mapsize(m::LazyJacobian) = mapsize(supermap(m))
@@ -29,9 +31,8 @@ struct DeterminantMap{T,M} <: SimpleLazyMap{T}
 end
 
 DeterminantMap(m) = DeterminantMap{domaintype(m)}(m)
-DeterminantMap{T}(m) where {T} = DeterminantMap{T,typeof(m)}(m)
-DeterminantMap{T}(m::Map{T}) where {T} = DeterminantMap{T,typeof(m)}(m)
-DeterminantMap{T}(m::Map{S}) where {S,T} = DeterminantMap{T}(convert(Map{T}, m))
+DeterminantMap{T}(m) where T = _DeterminantMap(T, convert_domaintype(T, m))
+_DeterminantMap(::Type{T}, m) where T = DeterminantMap{T,typeof(m)}(m)
 
 determinantmap(m) = DeterminantMap(m)
 
@@ -56,9 +57,8 @@ struct AbsMap{T,M} <: SimpleLazyMap{T}
 end
 
 AbsMap(m) = AbsMap{domaintype(m)}(m)
-AbsMap{T}(m) where {T} = AbsMap{T,typeof(m)}(m)
-AbsMap{T}(m::Map{T}) where {T} = AbsMap{T,typeof(m)}(m)
-AbsMap{T}(m::Map{S}) where {S,T} = AbsMap{T}(convert(Map{T}, m))
+AbsMap{T}(m) where T = _AbsMap(T, convert_domaintype(T, m))
+_AbsMap(::Type{T}, m) where T = AbsMap{T,typeof(m)}(m)
 
 absmap(m) = AbsMap(m)
 
@@ -73,9 +73,8 @@ struct LazyDiffVolume{T,M} <: SimpleLazyMap{T}
 end
 
 LazyDiffVolume(m) = LazyDiffVolume{domaintype(m)}(m)
-LazyDiffVolume{T}(m) where {T} = LazyDiffVolume{T,typeof(m)}(m)
-LazyDiffVolume{T}(m::Map{T}) where {T} = LazyDiffVolume{T,typeof(m)}(m)
-LazyDiffVolume{T}(m::Map{S}) where {S,T} = LazyDiffVolume{T}(convert(Map{T}, m))
+LazyDiffVolume{T}(m) where T = _LazyDiffVolume(T, convert_domaintype(T, m))
+_LazyDiffVolume(::Type{T}, m) where T = LazyDiffVolume{T,typeof(m)}(m)
 
 applymap(m::LazyDiffVolume, x) = diffvolume(supermap(m), x)
 

src/generic/map.jl

@@ -159,6 +159,28 @@ isvectorvalued_type(::Type{T}) where {T} = false
 isvectorvalued(m) =
     isvectorvalued_type(domaintype(m)) && isvectorvalued_type(codomaintype(m))
 
+"""
+    mapsize(m[, i])
+
+The size of a vectorvalued map is the size of its jacobian matrix.
+
+A map with size (3,5) maps vectors of length 5 to vectors of length 3. Its jacobian matrix
+is a 3-by-5 matrix.
+
+Optionally, similar to the `size` function, a second argument `i` can be provided. Choose `i`
+to be `1` or `2` to return only the corresponding element of the size. The dimension of the
+domain type is `mapsize(m, 2)`, while that of the codomain type is `mapsize(m, 1)`.
+
+If a map is scalar valued, its `mapsize` equals `()`. This reflects the convention in Julia that
+the size of a number is an empty tuple: `size(5) == ()`. If a map is scalar-to-vector, its `mapsize`
+is the tuple `(m,)` of length one. In this case the jacobian is a vector with size `(m,)`. If a map
+is vector-to-scalar, its `mapsize` is `(1,n)`.
+
+Thus, in all cases of scalarvalued and vectorvalued maps, the `mapsize` of a map agrees with the
+size of the Jacobian matrix.
+"""
+function mapsize end
+
 # mapsize should be defined for vector valued maps
 # The size of a map equals the size of its jacobian
 # The jacobian can be a number, a vector, an adjoint vector, or a matrix

src/generic/product.jl

@@ -114,8 +114,6 @@ function VcatMap{T}(maps...) where T
 	VcatMap{T,M,N}(maps...)
 end
 
-mapdim(map) = mapsize(map,2)
-
 VcatMap{T,M,N}(maps::Union{Tuple,Vector}) where {T,M,N} = VcatMap{T,M,N}(maps...)
 function VcatMap{T,M,N}(maps...) where {T,M,N}
 	DIM1 = map(t->mapsize(t,1), maps)

test/runtests.jl

@@ -15,5 +15,7 @@ include("test_basic.jl")
 include("test_affine.jl")
 include("test_canonical.jl")
 include("test_product.jl")
+include("test_composite.jl")
+
 include("test_maps.jl")
 include("test_arithmetics.jl")

test/test_composite.jl

@@ -0,0 +1,77 @@
+function test_composite_maps()
+    T = Float64
+    a = T(0)
+    b = T(1)
+    c = T(2)
+    d = T(3)
+    ma = IdentityMap{T}()
+    mb = interval_map(a, b, c, d)
+
+    r = suitable_point_to_map(ma)
+    m1 = ComposedMap(mb, ma)
+    test_generic_map(m1)
+    @test m1(r) ≈ ma(mb(r))
+    m2 = ComposedMap(mb, m1)
+    test_generic_map(m2)
+
+    @test ncomponents(composedmap(m1,m1)) == 4
+
+    m5 = ComposedMap(LinearMap(rand(T,2,2)), AffineMap(rand(T,2,2),rand(T,2)))
+    test_generic_map(m5)
+    @test jacobian(m5) isa ConstantMap
+    @test m5[Component(1)] isa LinearMap
+    @test m5[Component(2)] isa AffineMap
+    @test ComposedMap(m5[Component(1:2)]...) == m5
+    @test_throws BoundsError m5[Component(3)]
+
+    m6 = multiply_map(ma,ma)
+    @test m6(one(T)/2) == one(T)/4
+    @test jacobian(m6) isa SumMap
+    @test jacobian(m6)(one(T)) == 2
+    @test jacobian(m6, one(T)) == 2
+
+    m7 = sum_map(ma,ma)
+    @test m7(one(T)) == 2
+    @test jacobian(m7) isa ConstantMap
+    @test jacobian(m7, one(T)) == 2
+    @test sum_map() == ()
+    @test sum_map(ma) == ma
+    @test sum_map(ma, mb, ma) isa SumMap{T,Tuple{X,Y,Z}} where X where Y where Z
+    @test sum_map(SumMap(ma, mb), ma) isa SumMap{T,Tuple{X,Y,Z}} where X where Y where Z
+    @test sum_map(ma, SumMap(ma, mb)) isa SumMap{T,Tuple{X,Y,Z}} where X where Y where Z
+    @test sum_map(SumMap(ma, mb), SumMap(ma,mb)) isa SumMap{T,Tuple{W,X,Y,Z}} where W where X where Y where Z
+    @test FunctionMaps.SumMap{Float64}(LinearMap(2), LinearMap(2.0)) isa FunctionMaps.SumMap{Float64}
+    @test allequal(map(domaintype,components(FunctionMaps.SumMap{Float64}(LinearMap(2), LinearMap(2.0)))))
+
+    @test mapsize(ComposedMap(LinearMap(2),LinearMap(rand(T,2)),LinearMap(rand(T,2,2)))) == (2,)
+    @test mapsize(ComposedMap(LinearMap(2),LinearMap(rand(T,2)))) == (2,)
+    @test mapsize(ComposedMap(LinearMap(rand(T,2)),LinearMap(rand(T,2)'),LinearMap(rand(T,2)))) == (2,)
+    @test mapsize(ComposedMap(LinearMap(rand(T,2,2)),LinearMap(rand(T,2)'),LinearMap(2))) == (1,2)
+    @test mapsize(ComposedMap(LinearMap(one(T)),LinearMap(one(T)))) == ()
+
+    @test composedmap() == ()
+    @test composedmap(ma) == ma
+    @test composedmap(ma,ma) == ma
+    @test composedmap(ma,ma,ma) == ma
+
+    @test composite_jacobian(ma) == jacobian(ma)
+
+    @test multiply_map() == ()
+    @test multiply_map(ma) == ma
+    @test multiply_map(ma,ma)(2*one(T)) == 4
+    @test multiply_map(ma,ma,ma)(2*one(T)) == 8
+    @test FunctionMaps.mul_jacobian() == ()
+    @test FunctionMaps.mul_jacobian(ma) == jacobian(ma)
+    @test jacobian(multiply_map(ma, ma, ma)) isa FunctionMaps.SumMap
+    
+    @test ncomponents(multiply_map(multiply_map(ma,ma),multiply_map(ma,ma))) == 4
+    @test multiply_map(LinearMap(2), LinearMap(2.0)) isa FunctionMaps.MulMap{Float64}
+    @test FunctionMaps.MulMap{Float64}(LinearMap(2), LinearMap(2.0)) isa FunctionMaps.MulMap{Float64}
+    @test allequal(map(domaintype,components(FunctionMaps.MulMap{Float64}(LinearMap(2), LinearMap(2.0)))))
+
+    @test sum_jacobian() == ()
+    @test sum_jacobian(ma) == jacobian(ma)
+
+    @test convert_domaintype(Float64, ComposedMap(LinearMap(2), LinearMap(2))) isa Map{Float64}
+    @test allequal(map(domaintype,components(convert_domaintype(Float64, ComposedMap(LinearMap(2), LinearMap(2))))))
+end