diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 0000000..700707c --- /dev/null +++ b/.github/dependabot.yml @@ -0,0 +1,7 @@ +# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates +version: 2 +updates: + - package-ecosystem: "github-actions" + directory: "/" # Location of package manifests + schedule: + interval: "weekly" diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml new file mode 100644 index 0000000..7a5bc0e --- /dev/null +++ b/.github/workflows/CI.yml @@ -0,0 +1,47 @@ +name: CI +on: + push: + branches: + - main + tags: ['*'] + pull_request: + workflow_dispatch: +concurrency: + # Skip intermediate builds: always. + # Cancel intermediate builds: only if it is a pull request build. + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} +jobs: + test: + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} + runs-on: ${{ matrix.os }} + timeout-minutes: 60 + permissions: # needed to allow julia-actions/cache to proactively delete old caches that it has created + actions: write + contents: read + strategy: + fail-fast: false + matrix: + version: + - '1.10' + - '1.11' + - 'pre' + os: + - ubuntu-latest + arch: + - x64 + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: julia-actions/cache@v2 + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-runtest@v1 + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v4 + with: + files: lcov.info + token: ${{ secrets.CODECOV_TOKEN }} + fail_ci_if_error: false diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml new file mode 100644 index 0000000..cba9134 --- /dev/null +++ b/.github/workflows/CompatHelper.yml @@ -0,0 +1,16 @@ +name: CompatHelper +on: + schedule: + - cron: 0 0 * * * + workflow_dispatch: +jobs: + CompatHelper: + runs-on: ubuntu-latest + steps: + - name: Pkg.add("CompatHelper") + run: julia -e 'using Pkg; Pkg.add("CompatHelper")' + - name: CompatHelper.main() + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + COMPATHELPER_PRIV: ${{ secrets.DOCUMENTER_KEY }} + run: julia -e 'using CompatHelper; CompatHelper.main()' diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml new file mode 100644 index 0000000..0cd3114 --- /dev/null +++ b/.github/workflows/TagBot.yml @@ -0,0 +1,31 @@ +name: TagBot +on: + issue_comment: + types: + - created + workflow_dispatch: + inputs: + lookback: + default: "3" +permissions: + actions: read + checks: read + contents: write + deployments: read + issues: read + discussions: read + packages: read + pages: read + pull-requests: read + repository-projects: read + security-events: read + statuses: read +jobs: + TagBot: + if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot' + runs-on: ubuntu-latest + steps: + - uses: JuliaRegistries/TagBot@v1 + with: + token: ${{ secrets.GITHUB_TOKEN }} + ssh: ${{ secrets.DOCUMENTER_KEY }} diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..95731a5 --- /dev/null +++ b/.gitignore @@ -0,0 +1,6 @@ +*.jl.*.cov +*.jl.cov +*.jl.mem +/Manifest.toml +/docs/Manifest.toml +/docs/build/ diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..ad3142c --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2025 Daan Huybrechs and contributors + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..3f93729 --- /dev/null +++ b/Project.toml @@ -0,0 +1,24 @@ +name = "FunctionMaps" +uuid = "a85aefff-f8ca-4649-a888-c8e5398bc76c" +authors = ["Daan Huybrechs and contributors"] +version = "0.1.0" + +[deps] +CompositeTypes = "b152e2b5-7a66-4b01-a709-34e65c35f657" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + +[compat] +Aqua = "0.8" +CompositeTypes = "0.1.3" +LinearAlgebra = "1" +StaticArrays = "1" +Test = "1" +julia = "1.10" + +[extras] +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test", "Aqua"] diff --git a/README.md b/README.md new file mode 100644 index 0000000..c24472d --- /dev/null +++ b/README.md @@ -0,0 +1,5 @@ +# FunctionMaps + +[![Build Status](https://github.com/daanhb/FunctionMaps.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/daanhb/FunctionMaps.jl/actions/workflows/CI.yml?query=branch%3Amain) +[![Coverage](https://codecov.io/gh/daanhb/FunctionMaps.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/daanhb/FunctionMaps.jl) +[![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 0000000..be2a388 --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,3 @@ +[deps] +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +FunctionMaps = "a85aefff-f8ca-4649-a888-c8e5398bc76c" diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 0000000..9c4f36f --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,17 @@ +using FunctionMaps +using Documenter + +DocMeta.setdocmeta!(FunctionMaps, :DocTestSetup, :(using FunctionMaps); recursive=true) + +makedocs(; + modules=[FunctionMaps], + authors="Daan Huybrechs and contributors", + sitename="FunctionMaps.jl", + format=Documenter.HTML(; + edit_link="main", + assets=String[], + ), + pages=[ + "Home" => "index.md", + ], +) diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 0000000..8901584 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,14 @@ +```@meta +CurrentModule = FunctionMaps +``` + +# FunctionMaps + +Documentation for [FunctionMaps](https://github.com/daanhb/FunctionMaps.jl). + +```@index +``` + +```@autodocs +Modules = [FunctionMaps] +``` diff --git a/src/FunctionMaps.jl b/src/FunctionMaps.jl new file mode 100644 index 0000000..931d0e0 --- /dev/null +++ b/src/FunctionMaps.jl @@ -0,0 +1,58 @@ +module FunctionMaps + +using CompositeTypes, CompositeTypes.Display, CompositeTypes.Indexing +using LinearAlgebra +using StaticArrays + +import Base: + convert, show, + ==, hash # for the equality of maps + +import CompositeTypes: component, components + +# Exhaustive list of exports: + +# from util/common.jl +export prectype, numtype + +# from generic/map.jl +export Map, MapRef, + applymap, isequalmap, + domaintype, codomaintype, + inverse, leftinverse, rightinverse, + mapsize, jacobian, jacdet, diffvolume, + isrealmap +# from generic/canonical.jl +export canonicalmap +# from generic/composite.jl +export composedmap +# from generic/product.jl +export productmap + +# from types/basic.jl +export IdentityMap, + StaticIdentityMap, VectorIdentityMap, + ZeroMap, UnityMap, ConstantMap, + isconstantmap, mapconstant +# from types/affine.jl +export AffineMap, Translation, LinearMap, + affinematrix, affinevector, + islinearmap, isaffinemap + +include("util/common.jl") + +include("generic/map.jl") +include("generic/interface.jl") +include("generic/canonical.jl") +include("generic/lazy.jl") +include("generic/inverse.jl") +include("generic/jacobian.jl") +include("generic/composite.jl") +include("generic/product.jl") +include("generic/isomorphism.jl") +include("types/basic.jl") +include("types/affine.jl") +include("types/coordinates.jl") +include("types/arithmetics.jl") + +end diff --git a/src/generic/canonical.jl b/src/generic/canonical.jl new file mode 100644 index 0000000..b951982 --- /dev/null +++ b/src/generic/canonical.jl @@ -0,0 +1,94 @@ +""" + canonicalmap([ctype::CanonicalType, ]map) + +Return an associated canonical map, if any, of the given map. + +Optionally, a canonical type argument may specify an alternative canonical map. +Canonical maps help with converting between equal maps of different types. +""" +canonicalmap(m) = m + +"Does the map have a canonical map?" +hascanonicalmap(m) = !(m === canonicalmap(m)) + + +"Supertype of different kinds of canonical objects." +abstract type CanonicalType end + +canonicalmap(ctype::CanonicalType, m) = m +hascanonicalmap(ctype::CanonicalType, m) = !(m === canonicalmap(ctype, m)) + + +""" + Equal <: CanonicalType + +A canonical object that is equal but simpler. +""" +struct Equal <: CanonicalType end + +canonicalmap(::Equal, m) = m +equalmap(m) = canonicalmap(Equal(), m) +hasequalmap(m) = hascanonicalmap(Equal(), m) + +""" + simplify(m) + +Simplify the given map to an equal map. +""" +simplify(m) = equalmap(m) +""" + simplifies(m) + +Does the map simplify? +""" +simplifies(m) = hasequalmap(m) + +"Convert the given map to a map defined in FunctionMaps.jl." +tofunctionmap(m::Map) = m +tofunctionmap(m::MapRef) = tofunctionmap(functionmap(m)) + + +""" + Equivalent <: CanonicalType + +A canonical object that is equivalent but may have different type. +""" +struct Equivalent <: CanonicalType end + +canonicalmap(::Equivalent, m) = canonicalmap(Equal(), m) + +canonicalmap(::Equivalent, m::Map{<:StaticVector{1,T}}) where {T} = convert(Map{T}, m) +canonicalmap(::Equivalent, m::Map{NTuple{N,T}}) where {N,T} = convert(Map{SVector{N,T}}, m) + +equivalentmap(m) = canonicalmap(Equivalent(), m) +hasequivalentmap(m) = hascanonicalmap(Equivalent(), m) + +""" + CanonicalExtensionType <: CanonicalType + +Canonical types used to translate between packages. +""" +abstract type CanonicalExtensionType <: CanonicalType +end + +"Return the extension type associated with the given object." +canonicalextensiontype(m) = canonicalextensiontype(typeof(m)) + + +==(m1::Map, m2::Map) = isequalmap(m1, m2) # Method from Base + +""" + isequalmap(map1, map2) + +Are the two given maps equal? +""" +isequalmap(m1, m2) = isequalmap1(m1, m2) +isequalmap1(m1, m2) = simplifies(m1) ? isequalmap(simplify(m1), m2) : isequalmap2(m1, m2) +isequalmap2(m1, m2) = simplifies(m2) ? isequalmap(m1, simplify(m2)) : default_isequalmap(m1, m2) +default_isequalmap(m1, m2) = m1===m2 + +# Associated with == we have to define the hashes of maps +Base.hash(m::Map, h::UInt) = map_hash(m, h) + +map_hash(m) = map_hash(m, zero(UInt)) +map_hash(m, h::UInt) = invoke(hash, Tuple{Any,UInt}, m, h) diff --git a/src/generic/composite.jl b/src/generic/composite.jl new file mode 100644 index 0000000..d6923ac --- /dev/null +++ b/src/generic/composite.jl @@ -0,0 +1,228 @@ + +"The composition of several maps." +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) + +# TODO: make proper conversion +similarmap(m::ComposedMap, ::Type{T}) where {T} = ComposedMap{T}(m.maps...) + +codomaintype(m::ComposedMap) = codomaintype(m.maps[end]) + +# Maps are applied in the order that they appear in m.maps +applymap(m::ComposedMap, x) = applymap_rec(x, m.maps...) +applymap_rec(x) = x +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])) +_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],) +_composed_mapsize(m, m_end, m1, S_end::Tuple{Int,Int}, S1::Tuple{}) = + is_vector_to_scalar(m_end) ? () : (S_end[1],) +_composed_mapsize(m, m_end, m1, S_end::Tuple{Int}, S1::Tuple{Int,Int}) = (S_end[1],S1[2]) +_composed_mapsize(m, m_end, m1, S_end::Tuple{Int}, S1::Tuple{Int}) = (S_end[1],) +_composed_mapsize(m, m_end, m1, S_end::Tuple{Int}, S1::Tuple{}) = (S_end[1],) +_composed_mapsize(m, m_end, m1, S_end::Tuple{}, S1::Tuple{Int,Int}) = (1,S1[2]) +_composed_mapsize(m, m_end, m1, S_end::Tuple{}, S1::Tuple{Int}) = () +_composed_mapsize(m, m_end, m1, S_end::Tuple{}, S1::Tuple{}) = () + +function jacobian(m::ComposedMap, x) + f, fd = backpropagate(x, reverse(components(m))...) + fd +end +backpropagate(x, m1) = (applymap(m1, x), jacobian(m1, x)) +function backpropagate(x, m2, ms...) + f, fd = backpropagate(x, ms...) + applymap(m2, f), jacobian(m2, f) * fd +end + +for op in (:inverse, :leftinverse, :rightinverse) + @eval $op(cmap::ComposedMap) = ComposedMap(reverse(map($op, components(cmap)))...) +end + +inverse(m::ComposedMap, x) = inverse_rec(x, reverse(components(m))...) +inverse_rec(x) = x +inverse_rec(x, map1, maps...) = inverse_rec(inverse(map1, x), maps...) + +leftinverse(m::ComposedMap, x) = leftinverse_rec(x, reverse(components(m))...) +leftinverse_rec(x) = x +leftinverse_rec(x, map1, maps...) = leftinverse_rec(leftinverse(map1, x), maps...) + +rightinverse(m::ComposedMap, x) = rightinverse_rec(x, reverse(components(m))...) +rightinverse_rec(x) = x +rightinverse_rec(x, map1, maps...) = rightinverse_rec(rightinverse(map1, x), maps...) + + +promote_composing_maps(m1, m2) = + _promote_composing_maps(codomaintype(m1), domaintype(m2), m1, m2) +function _promote_composing_maps(::Type{S}, ::Type{T}, m1, m2) where {S,T} + U = promote_type(S, T) + convert_codomaintype(U, m1), convert_domaintype(U, m2) +end + +composedmap() = () +composedmap(m) = m +composedmap(m1, m2) = composedmap1(promote_composing_maps(m1, m2)...) +composedmap1(m1, m2) = composedmap2(m1, m2) +composedmap2(m1, m2) = default_composedmap(m1, m2) +default_composedmap(m1, m2) = ComposedMap(m1, m2) + +composedmap(m1, m2, maps...) = composedmap(composedmap(m1, m2), maps...) + +composedmap(m1::ComposedMap, m2::ComposedMap) = + ComposedMap(components(m1)..., components(m2)...) +function composedmap1(m1::ComposedMap, m2) + T = promote_type(codomaintype(m1), domaintype(m2)) + ComposedMap(components(m1)..., convert_domaintype(T, m2)) +end +function composedmap2(m1, m2::ComposedMap) + T = promote_type(codomaintype(m1), domaintype(m2)) + ComposedMap(convert_codomaintype(T, m1), components(m2)...) +end + +# Arguments to ∘ should be reversed before passing on to mapcompose +Base.:∘(map1::AbstractMap, map2::AbstractMap) = composedmap(map2, map1) + + +isequalmap(m1::ComposedMap, m2::ComposedMap) = + ncomponents(m1) == ncomponents(m2) && all(map(isequalmap, components(m1), components(m2))) +map_hash(m::ComposedMap, h::UInt) = hashrec("ComposedMap", collect(components(m)), h) + +Display.combinationsymbol(m::ComposedMap) = Display.Symbol('∘') +Display.displaystencil(m::ComposedMap) = + composite_displaystencil(m; reversecomponents=true) +map_object_parentheses(m::ComposedMap) = true +map_stencil_parentheses(m::ComposedMap) = true +show(io::IO, mime::MIME"text/plain", m::ComposedMap) = composite_show(io, mime, m) +show(io::IO, m::ComposedMap) = composite_show_compact(io, m) + +## Lazy multiplication + +"The lazy multiplication of one or more maps." +struct MulMap{T,MAPS} <: CompositeLazyMap{T} + maps :: MAPS +end + +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(::Type{T}, maps...) where {T} = MulMap{T,typeof(maps)}(maps) + +similarmap(m::MulMap, ::Type{T}) where {T} = MulMap{T}(m.maps...) + +applymap(m::MulMap, x) = reduce(*, applymap.(components(m), Ref(x))) + +multiply_map() = () +multiply_map(m) = m +multiply_map(m1, m2) = multiply_map1(m1, m2) +multiply_map1(m1, m2) = multiply_map2(m1, m2) +multiply_map2(m1, m2) = default_multiply_map(m1, m2) +default_multiply_map(m1, m2) = MulMap(m1, m2) + +multiply_map(m1, m2, maps...) = multiply_map(multiply_map(m1, m2), maps...) + +multiply_map(m1::MulMap, m2::MulMap) = + MulMap(components(m1)..., components(m2)...) +multiply_map1(m1::MulMap, m2) = MulMap(components(m1)..., m2) +multiply_map2(m1, m2::MulMap) = MulMap(m1, components(m2)...) + +function Display.displaystencil(m::MulMap) + A = Any[] + list = components(m) + push!(A, "x -> ") + push!(A, Display.SymbolObject(list[1])) + push!(A, "(x)") + for i in 2:length(list) + push!(A, " * ") + push!(A, Display.SymbolObject(list[i])) + push!(A, "(x)") + end + A +end +show(io::IO, mime::MIME"text/plain", m::MulMap) = composite_show(io, mime, m) + + +## Lazy sum + +"The lazy sum of one or more maps." +struct SumMap{T,MAPS} <: CompositeLazyMap{T} + maps :: MAPS +end + +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(::Type{T}, maps...) where {T} = SumMap{T,typeof(maps)}(maps) + +similarmap(m::SumMap, ::Type{T}) where {T} = SumMap{T}(m.maps...) + +applymap(m::SumMap, x) = reduce(+, applymap.(components(m), Ref(x))) + +sum_map() = () +sum_map(m) = m +sum_map(m1, m2) = sum_map1(m1, m2) +sum_map1(m1, m2) = sum_map2(m1, m2) +sum_map2(m1, m2) = default_sum_map(m1, m2) +default_sum_map(m1, m2) = SumMap(m1, m2) + +sum_map(m1, m2, maps...) = sum_map(sum_map(m1, m2), maps...) + +sum_map(m1::SumMap, m2::SumMap) = + SumMap(components(m1)..., components(m2)...) +sum_map1(m1::SumMap, m2) = SumMap(components(m1)..., m2) +sum_map2(m1, m2::SumMap) = SumMap(m1, components(m2)...) + +function Display.displaystencil(m::SumMap) + A = Any[] + list = components(m) + push!(A, "x -> ") + push!(A, Display.SymbolObject(list[1])) + push!(A, "(x)") + for i in 2:length(list) + push!(A, " + ") + push!(A, Display.SymbolObject(list[i])) + push!(A, "(x)") + end + A +end +show(io::IO, mime::MIME"text/plain", m::SumMap) = composite_show(io, mime, m) + + +# Define the jacobian of a composite map +jacobian(m::ComposedMap) = composite_jacobian(reverse(components(m))...) +composite_jacobian(map1) = jacobian(map1) +composite_jacobian(map1, map2) = multiply_map(composedmap(map2, jacobian(map1)), jacobian(map2)) +function composite_jacobian(map1, map2, maps...) + rest = ComposedMap(reverse(maps)..., map2) + f1 = composedmap(rest, jacobian(map1)) + f2 = composite_jacobian(map2, maps...) + multiply_map(f1, f2) +end + +jacobian(m::MulMap) = mul_jacobian(components(m)...) +mul_jacobian() = () +mul_jacobian(map1) = jacobian(map1) +mul_jacobian(map1, map2) = sum_map(multiply_map(jacobian(map1), map2), multiply_map(map1, jacobian(map2))) +function mul_jacobian(map1, map2, maps...) + rest = multiply_map(map2, maps...) + mul_jacobian(map1, rest) +end +function jacobian(m::MulMap, x) + z = map(t -> applymap(t,x), components(m)) + zd = map(t -> jacobian(t, x), components(m)) + sum(prod(z[1:i-1]) * zd[i] * prod(z[i+1:end]) for i in 1:ncomponents(m)) +end + + +jacobian(m::SumMap) = sum_jacobian(components(m)...) +sum_jacobian() = () +sum_jacobian(map1) = jacobian(map1) +sum_jacobian(maps...) = sum_map(map(jacobian, maps)...) + +jacobian(m::SumMap, x) = sum(jacobian(mc, x) for mc in components(m)) diff --git a/src/generic/interface.jl b/src/generic/interface.jl new file mode 100644 index 0000000..24ac3cc --- /dev/null +++ b/src/generic/interface.jl @@ -0,0 +1,73 @@ + +""" + MapStyle(m) + +Trait to indicate whether or not `m` implements the map interface. +""" +abstract type MapStyle end + +""" + IsMap() + +indicates an object implements the map interface. +""" +struct IsMap <: MapStyle end + +""" + NotMap() + +indicates an object does not implement the Map interface. +""" +struct NotMap <: MapStyle end + + +MapStyle(m) = MapStyle(typeof(m)) +# the default is no map +MapStyle(::Type) = NotMap() +# subtypes of Map are maps +MapStyle(::Type{<:Map}) = IsMap() + +""" + functionmap(m) + +Return a map associated with the object `m`. +""" +functionmap(m::Map) = m + +""" + MapRef(m) + +A reference to a map. + +In a function call, `MapRef(x)` can be used to indicate that `x` should be +treated as a map, e.g., `foo(x, MapRef(m))`. +""" +struct MapRef{M} + map :: M +end + +functionmap(m::MapRef) = m.map +domaintype(m::MapRef) = domaintype(functionmap(m)) + + +""" +`AnyMap` is the union of `Map` and `MapRef`. + +In both cases `map(m::AnyMap)` returns the map itself. +""" +const AnyMap = Union{Map,MapRef} + +""" + checkmap(m) + +Checks that `m` is a map or refers to a map and if so returns that map, +throws an error otherwise. +""" +checkmap(m::Map) = m +# we trust the explicit intention of a user providing a map reference +checkmap(m::MapRef) = functionmap(m) +# for other objects we check MapStyle +checkmap(m) = _checkmap(m, MapStyle(m)) +_checkmap(m, ::IsMap) = m +_checkmap(m, ::NotMap) = + throw(ArgumentError("Map does not implement map interface as indicated by MapStyle.")) diff --git a/src/generic/inverse.jl b/src/generic/inverse.jl new file mode 100644 index 0000000..3af081e --- /dev/null +++ b/src/generic/inverse.jl @@ -0,0 +1,56 @@ + + +"A lazy inverse stores a map `m` and returns `inverse(m, x)`." +struct LazyInverse{T,M} <: SimpleLazyMap{T} + map :: M +end + +LazyInverse(m) = LazyInverse{codomaintype(m)}(m) +LazyInverse{T}(m) where {T} = LazyInverse{T,typeof(m)}(m) + +applymap(m::LazyInverse, x) = inverse(supermap(m), x) + +Display.displaystencil(m::LazyInverse) = ["LazyInverse(", supermap(m), ")"] +show(io::IO, mime::MIME"text/plain", m::LazyInverse) = composite_show(io, mime, m) + + +""" + inverse(m[, x]) + +Return the inverse of `m`. The two-argument function evaluates the inverse +at the point `x`. +""" +inverse(m) = LazyInverse(m) +inverse(m::LazyInverse) = supermap(m) +# Concrete maps should implement inverse(m, x) + +Base.:\(m::AbstractMap, x) = inverse(m, x) + +implements_inverse(m) = !(inverse(m) isa LazyInverse) + +""" + leftinverse(m[, x]) + +Return a left inverse of the given map. This left inverse `mli` is not unique, +but in any case it is such that `(mli ∘ m) * x = x` for each `x` in the domain +of `m`. + +The two-argument function applies the left inverse to the point `x`. +""" +leftinverse(m) = inverse(m) +leftinverse(m, x) = inverse(m, x) + +""" + rightinverse(m[, x]) + +Return a right inverse of the given map. This right inverse `mri` is not unique, +but in any case it is such that `(m ∘ mri) * y = y` for each `y` in the range +of `m`. + +The two-argument function applies the right inverse to the point `x`. +""" +rightinverse(m) = inverse(m) +rightinverse(m, x) = inverse(m, x) + +isequalmap(m1::LazyInverse, m2::LazyInverse) = isequalmap(supermap(m1), supermap(m2)) +map_hash(m::LazyInverse, h::UInt) = hashrec("LazyInverse", supermap(m), h) diff --git a/src/generic/isomorphism.jl b/src/generic/isomorphism.jl new file mode 100644 index 0000000..5155f03 --- /dev/null +++ b/src/generic/isomorphism.jl @@ -0,0 +1,124 @@ + +""" + Isomorphism{T,U} <: TypedMap{T,U} + +An isomorphism is a bijection between types that preserves norms. +""" +abstract type Isomorphism{T,U} <: TypedMap{T,U} end + +show(io::IO, m::Isomorphism{T,U}) where {T,U} = print(io, "x : $(T) -> x : $(U)") +map_object_parentheses(m::Isomorphism) = true + +struct VectorToNumber{T} <: Isomorphism{SVector{1,T},T} +end +struct NumberToVector{T} <: Isomorphism{T,SVector{1,T}} +end + +""" + VectorToNumber() + VectorToNumber{T}() + +Map a length 1 vector `x` to the number `x[1]`. + +See also: [`NumberToVector`](@ref). +""" +VectorToNumber() = VectorToNumber{Float64}() + +""" + NumberToVector() + NumberToVector{T}() + +Map a number `x` to the length 1 vector `[x]`. + +See also: [`VectorToNumber`](@ref). +""" +NumberToVector() = NumberToVector{Float64}() + +mapsize(::VectorToNumber) = (1,1) +mapsize(::NumberToVector) = (1,) + +inverse(::VectorToNumber{T}) where {T} = NumberToVector{T}() +inverse(::NumberToVector{T}) where {T} = VectorToNumber{T}() +inverse(m::VectorToNumber, x) = inverse(m)(x) +inverse(m::NumberToVector, x) = inverse(m)(x) + +applymap(::VectorToNumber, x) = x[1] +applymap(::NumberToVector, x) = SVector(x) + +jacobian(::VectorToNumber{T}, x) where {T} = transpose(SVector(one(T))) +jacobian(::VectorToNumber{T}) where {T} = ConstantMap{SVector{1,T}}(transpose(SVector(one(T)))) + +jacobian(::NumberToVector{T}, x) where {T} = SVector(one(T)) +jacobian(::NumberToVector{T}) where {T} = ConstantMap{T}(SVector(one(T))) + + +struct VectorToComplex{T} <: Isomorphism{SVector{2,T},Complex{T}} +end +struct ComplexToVector{T} <: Isomorphism{Complex{T},SVector{2,T}} +end + +""" + VectorToComplex() + VectorToComplex{T}() + +Map a length 2 vector ``[a;b]`` to the complex number ``a+bi``. + +See also: [`ComplexToVector`](@ref). +""" +VectorToComplex() = VectorToComplex{Float64}() + +""" + ComplexToVector() + ComplexToVector{T}() + +Map a complex number ``a+bi`` to the length 2 vector ``[a; b]``. + +See also: [`VectorToComplex`](@ref). +""" +ComplexToVector() = ComplexToVector{Float64}() + +mapsize(::VectorToComplex) = (1,2) +mapsize(::ComplexToVector) = (2,) + +applymap(::VectorToComplex, x) = x[1] + im*x[2] +applymap(::ComplexToVector, x) = SVector(real(x), imag(x)) + +inverse(::VectorToComplex{T}) where {T} = ComplexToVector{T}() +inverse(::ComplexToVector{T}) where {T} = VectorToComplex{T}() +inverse(m::VectorToComplex, x) = inverse(m)(x) +inverse(m::ComplexToVector, x) = inverse(m)(x) + +jacobian(::VectorToComplex{T}, x) where {T} = transpose(SVector(one(T),one(T)*im)) +jacobian(::VectorToComplex{T}) where {T} = ConstantMap{SVector{2,T}}(transpose(SVector(one(T),one(T)*im))) + + +"Map a static vector to a tuple." +struct VectorToTuple{N,T} <: Isomorphism{SVector{N,T},NTuple{N,T}} +end +"Map a tuple to a static vector." +struct TupleToVector{N,T} <: Isomorphism{NTuple{N,T},SVector{N,T}} +end + +inverse(::VectorToTuple{N,T}) where {N,T} = TupleToVector{N,T}() +inverse(::TupleToVector{N,T}) where {N,T} = VectorToTuple{N,T}() +inverse(m::VectorToTuple, x) = inverse(m)(x) +inverse(m::TupleToVector, x) = inverse(m)(x) + +applymap(::VectorToTuple, x) = tuple(x...) +applymap(::TupleToVector, x) = SVector(x) + + +"Map a nested vector or tuple to a flat vector." +struct NestedToFlat{N,T,U,DIM} <: Isomorphism{U,SVector{N,T}} +end +"Map a flattened vector to a nested one." +struct FlatToNested{N,T,U,DIM} <: Isomorphism{SVector{N,T},U} +end + +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) +inverse(m::FlatToNested, x) = inverse(m)(x) + +applymap(::NestedToFlat{N,T,U,DIM}, x) where {N,T,U,DIM} = convert_tocartesian(x, Val{DIM}()) +applymap(::FlatToNested{N,T,U,DIM}, x) where {N,T,U,DIM} = convert_fromcartesian(x, Val{DIM}()) diff --git a/src/generic/jacobian.jl b/src/generic/jacobian.jl new file mode 100644 index 0000000..a6b30dd --- /dev/null +++ b/src/generic/jacobian.jl @@ -0,0 +1,179 @@ + +"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) +applymap(m::LazyJacobian, x) = jacobian(supermap(m), x) + +mapsize(m::LazyJacobian) = mapsize(supermap(m)) + +Display.displaystencil(m::LazyJacobian) = ["LazyJacobian(", supermap(m), ")"] +show(io::IO, mime::MIME"text/plain", m::LazyJacobian) = composite_show(io, mime, m) + + +""" + jacobian(m[, x]) + +Return the jacobian map. The two-argument version evaluates the jacobian +at a point `x`. +""" +jacobian(m) = LazyJacobian(m) +jacobian!(y, m, x) = y .= jacobian(m, x) + + +"A `DeterminantMap` returns the determinant of the result of a given map." +struct DeterminantMap{T,M} <: SimpleLazyMap{T} + map :: M +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(m) = DeterminantMap(m) + +applymap(m::DeterminantMap, x) = det(supermap(m)(x)) + + +""" + jacdet(m[, x]) + +Return the determinant of the jacobian as a map. The two-argument version +evaluates the jacobian determinant at a point `x`. +""" +jacdet(m) = determinantmap(jacobian(m)) +jacdet(m, x) = det(jacobian(m, x)) + + +"An `AbsMap` returns the absolute value of the result of a given map." +struct AbsMap{T,M} <: SimpleLazyMap{T} + map :: M +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(m) = AbsMap(m) + +applymap(m::AbsMap, x) = abs(supermap(m)(x)) + + +"A lazy volume element evaluates to `diffvolume(m, x)` on the fly." +struct LazyDiffVolume{T,M} <: SimpleLazyMap{T} + map :: M +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)) + +applymap(m::LazyDiffVolume, x) = diffvolume(supermap(m), x) + +Display.displaystencil(m::LazyDiffVolume) = ["LazyDiffVolume(", supermap(m), ")"] +show(io::IO, mime::MIME"text/plain", m::LazyDiffVolume) = composite_show(io, mime, m) + +""" + diffvolume(m[, x]) + +Compute the differential volume (at a point `x`). If `J` is the Jacobian matrix, +possibly rectangular, then the differential volume is `sqrt(det(J'*J))`. + +If the map is square, then the differential volume is the absolute value of the +Jacobian determinant. +""" +function diffvolume(m) + if issquaremap(m) + if mapsize(m,1) > 1 + absmap(jacdet(m)) + else + jacdet(m) + end + else + LazyDiffVolume(m) + end +end + +function diffvolume(m, x) + if issquaremap(m) + if mapsize(m,1) > 1 + abs(jacdet(m, x)) + else + jacdet(m, x) + end + else + jac = jacobian(m, x) + sqrt(det(adjoint(jac)*jac)) + end +end + + + +const NumberLike = Union{Number,UniformScaling} + +""" + to_matrix(::Type{T}, A[, b]) + +Convert the `A` in the affine map `A*x` or `A*x+b` with domaintype `T` to a matrix. +""" +to_matrix(::Type{T}, A) where {T} = A +to_matrix(::Type{T}, A::AbstractMatrix) where {T} = A +to_matrix(::Type{T}, A::NumberLike) where {T<:Number} = A +to_matrix(::Type{<:StaticVector{N,T}}, A::Number) where {N,T} = A * one(SMatrix{N,N,T}) +to_matrix(::Type{<:StaticVector{N,T}}, A::UniformScaling) where {N,T} = A.λ * one(SMatrix{N,N,T}) +to_matrix(::Type{T}, A::Number) where {T<:AbstractVector} = A * I +to_matrix(::Type{T}, A::UniformScaling) where {T<:Number} = one(T) +to_matrix(::Type{T}, A::UniformScaling) where {T<:AbstractVector} = A + +to_matrix(::Type{T}, A, b) where {T} = A +to_matrix(::Type{T}, A::AbstractMatrix, b) where {T} = A +to_matrix(::Type{T}, A::Number, b::Number) where {T<:Number} = A +to_matrix(::Type{T}, A::UniformScaling{S}, b::Number) where {S,T<:Number} = + convert(promote_type(S,T,typeof(b)), A.λ) +to_matrix(::Type{<:StaticVector{N,T}}, A::NumberLike, b::StaticVector{N,T}) where {N,T} = A * one(SMatrix{N,N,T}) +to_matrix(::Type{T}, A::NumberLike, b::AbstractVector) where {S,T<:AbstractVector{S}} = + A * Array{S,2}(I, length(b), length(b)) + +""" + to_vector(::Type{T}, A[, b]) + +Convert the `b` in the affine map `A*x` or `A*x+b` with domaintype `T` to a vector. +""" +to_vector(::Type{T}, A) where {T} = zero(T) +to_vector(::Type{T}, A::StaticVector{M,S}) where {T,M,S} = zero(SVector{M,S}) +to_vector(::Type{T}, A::StaticMatrix{M,N,S}) where {T,M,N,S} = zero(SVector{M,S}) +to_vector(::Type{T}, A::AbstractArray) where {T} = zeros(eltype(T),size(A,1)) +to_vector(::Type{T}, A, b) where {T} = b + + +"Return a zero matrix of the same size as the map." +zeromatrix(m) = zeromatrix(m, domaintype(m), codomaintype(m)) +zeromatrix(m, ::Type{T}, ::Type{U}) where {T,U} = zeros(numtype(T),mapsize(m)) +zeromatrix(m, ::Type{T}, ::Type{U}) where {T<:Number,U<:Number} = zero(promote_type(T,U)) +zeromatrix(m, ::Type{T}, ::Type{U}) where {N,T<:StaticVector{N},U<:Number} = + transpose(zero(SVector{N,eltype(T)})) +zeromatrix(m, ::Type{T}, ::Type{U}) where {N,M,T<:StaticVector{N},U<:StaticVector{M}} = + zero(SMatrix{M,N,promote_type(eltype(T),eltype(U))}) +zeromatrix(m, ::Type{T}, ::Type{U}) where {T<:Number,M,U<:StaticVector{M}} = + zero(SVector{M,promote_type(T,eltype(U))}) +zeromatrix(m, ::Type{T}, ::Type{U}) where {T<:AbstractVector,U<:Number} = + transpose(zeros(promote_type(eltype(T),U), mapsize(m,1))) + + +"Return a zero vector of the same size as the codomain of the map." +zerovector(m) = zerovector(m, codomaintype(m)) +zerovector(m, ::Type{U}) where {U} = zero(U) +zerovector(m, ::Type{<:StaticVector{M,T}}) where {M,T} = zero(SVector{M,T}) +# If the output type is a vector, the map itself should store the size information. +zerovector(m, ::Type{<:AbstractVector{T}}) where {T} = zeros(T, mapsize(m,1)) + +"Return an identity matrix with the same size as the map." +identitymatrix(m) = identitymatrix(m, codomaintype(m)) +identitymatrix(m, ::Type{T}) where {T} = one(T) +identitymatrix(m, ::Type{<:StaticVector{N,T}}) where {N,T} = one(SMatrix{N,N,T}) +identitymatrix(m, ::Type{<:AbstractVector{T}}) where {T} = Diagonal{T}(ones(mapsize(m,1))) diff --git a/src/generic/lazy.jl b/src/generic/lazy.jl new file mode 100644 index 0000000..65fb950 --- /dev/null +++ b/src/generic/lazy.jl @@ -0,0 +1,53 @@ + +""" +A lazy map has an action that is defined in terms of other maps. Those maps are +stored internally, and the action of the lazy map is computed on-the-fly and only +when invoked. +""" +abstract type LazyMap{T} <: Map{T} end + + +"A composite lazy map is defined in terms of several other maps." +abstract type CompositeLazyMap{T} <: LazyMap{T} end +"A simple lazy map derives from a single other map." +abstract type SimpleLazyMap{T} <: LazyMap{T} end + +supermap(m::SimpleLazyMap) = m.map +components(m::CompositeLazyMap) = m.maps + +Base.getindex(m::CompositeLazyMap, I::ComponentIndex...) = component(m, map(Indexing.to_index, I)...) + +isrealmap(m::SimpleLazyMap) = isrealmap(supermap(m)) +isrealmap(m::CompositeLazyMap) = all(map(isrealmap, components(m))) + +"A `DerivedMap` inherits all of its properties from another map, but has its own type." +abstract type DerivedMap{T} <: SimpleLazyMap{T} end + +applymap(m::DerivedMap, x) = supermap(m)(x) +appymap!(y, m::DerivedMap, x) = applymap!(y, supermap(m), x) + +jacobian(m::DerivedMap) = jacobian(supermap(m)) + + +"A `WrappedMap{T}` takes any object and turns it into a `Map{T}`." +struct WrappedMap{T,M} <: DerivedMap{T} + map :: M +end +WrappedMap{T}(map) where {T} = WrappedMap{T,typeof(map)}(map) +WrappedMap(map) = WrappedMap{Float64}(map) + +similarmap(m::WrappedMap, ::Type{T}) where {T} = WrappedMap{T}(m) + +convert(::Type{Map}, m::Map) = m +convert(::Type{Map}, m) = WrappedMap(m) +convert(::Type{Map{T}}, m) where {T} = WrappedMap{T}(m) +convert(::Type{Map{T}}, m::WrappedMap) where {T} = WrappedMap{T}(m.map) + +isequalmap(m1::WrappedMap, m2::Function) = m1.map == m2 +isequalmap(m1::Function, m2::WrappedMap) = m1 == m2.map + +map_hash(m::WrappedMap, h::UInt) = hashrec("WrappedMap", supermap(m), h) + +Display.displaystencil(m::WrappedMap{T}) where {T} = + ["WrappedMap{$T}(", supermap(m), ")"] +show(io::IO, mime::MIME"text/plain", m::WrappedMap) = composite_show(io, mime, m) diff --git a/src/generic/map.jl b/src/generic/map.jl new file mode 100644 index 0000000..00c39f6 --- /dev/null +++ b/src/generic/map.jl @@ -0,0 +1,175 @@ + +"An `AbstractMap` represents a function `y=f(x)` of a single variable." +abstract type AbstractMap end + +"A `Map{T}` is a map of a single variable of type `T`." +abstract type Map{T} <: AbstractMap end + +Map(m) = convert(Map, m) +Map{T}(m) where {T} = convert(Map{T}, m) + +"A `TypedMap{T,U}` maps a variable of type `T` to a variable of type `U`." +abstract type TypedMap{T,U} <: Map{T} end + +const EuclideanMap{N,T} = Map{<:StaticVector{N,T}} +const VectorMap{T} = Map{Vector{T}} + +CompositeTypes.Display.displaysymbol(m::Map) = 'F' + +""" + domaintype(m) + +What is the expected type of a point in the domain of the function map `m`? +""" +domaintype(m) = domaintype(typeof(m)) +domaintype(::Type{M}) where {M} = Any +domaintype(::Type{<:Map{T}}) where {T} = T + +""" + codomaintype(m[, T]) + +What is the codomain type of the function map `m`, given that `T` is its domain type? +""" +codomaintype(m) = codomaintype(m, domaintype(m)) +codomaintype(m, ::Type{T}) where {T} = _codomaintype(typeof(m), T) +_codomaintype(::Type{M}, ::Type{T}) where {M,T} = Base.promote_op(applymap, M, T) +_codomaintype(::Type{M}, ::Type{Any}) where {M} = Any +_codomaintype(M::Type{<:TypedMap{T,U}}, ::Type{T}) where {T,U} = U +_codomaintype(M::Type{<:TypedMap{T,U}}, ::Type{Any}) where {T,U} = U + +prectype(::Type{<:Map{T}}) where T = prectype(T) +numtype(::Type{<:Map{T}}) where T = numtype(T) + +isrealmap(m) = isrealtype(domaintype(m)) && isrealtype(codomaintype(m)) +isrealmap(::UniformScaling{T}) where {T} = isrealtype(T) + +convert(::Type{AbstractMap}, m::AbstractMap) = m +convert(::Type{Map{T}}, m::Map{T}) where {T} = m +convert(::Type{Map{T}}, m::Map{S}) where {S,T} = similarmap(m, T) +convert(::Type{TypedMap{T,U}}, m::TypedMap{T,U}) where {T,U} = m +convert(::Type{TypedMap{T,U}}, m::TypedMap) where {T,U} = similarmap(m, T, U) + +convert_domaintype(::Type{T}, map::Map{T}) where {T} = map +convert_domaintype(::Type{U}, map::Map{T}) where {T,U} = convert(Map{U}, map) +convert_domaintype(::Type{Any}, map) = map +convert_domaintype(::Type{Any}, map::Map{T}) where T = map +convert_domaintype(::Type{T}, map) where T = _convert_domaintype(T, map, domaintype(map)) +_convert_domaintype(::Type{T}, map, ::Type{T}) where T = map +_convert_domaintype(::Type{U}, map, ::Type{T}) where {T,U} = Map{U}(map) + +convert_numtype(::Type{U}, map::Map{T}) where {T,U} = convert(Map{to_numtype(U,T)}, map) +convert_prectype(::Type{U}, map::Map{T}) where {T,U} = convert(Map{to_prectype(U,T)}, map) + +convert_codomaintype(::Type{U}, map) where U = + _convert_codomaintype(U, map, domaintype(map), codomaintype(map)) +# types match: we can return map +_convert_codomaintype(::Type{U}, map, ::Type{T}, ::Type{U}) where {T,U} = map +# ...or, outputs are vectors with the same eltype: that's good enough, further +# conversion might be unnecessary and expensive +_convert_codomaintype(::Type{<:AbstractVector{S}}, map, ::Type{T}, ::Type{<:AbstractVector{S}}) where {T,S} = map +# no match yet: we now try to convert the numtype +_convert_codomaintype(::Type{U}, map, ::Type{T}, ::Type{V}) where {T,U,V} = + _convert_codomaintype2(U, convert_numtype(numtype(U), map)) +_convert_codomaintype2(::Type{U}, map::Map) where U = + _convert_codomaintype2(U, map, domaintype(map), codomaintype(map)) +# better match this time around: we can return map +_convert_codomaintype2(::Type{U}, map, ::Type{T}, ::Type{U}) where {T,U} = map +_convert_codomaintype2(::Type{<:AbstractVector{S}}, map, ::Type{T}, ::Type{<:AbstractVector{S}}) where {T,S} = map +# still no match: we can't do this automatically +_convert_codomaintype2(::Type{U}, map, ::Type{T}, ::Type{V}) where {T,U,V} = + throw(ArgumentError("Don't know how to convert the codomain type of $(map) to $(U).")) + +convert_codomaintype(::Type{Any}, map) = map + + +# Users may call a map, concrete subtypes specialize the `applymap` function +(m::AbstractMap)(x) = applymap(m, x) + +# For Map{T}, we allow invocation with multiple arguments by conversion to T +(m::Map{T})(x) where {T} = promote_and_apply(m, x) +(m::Map{T})(x...) where {T} = promote_and_apply(m, convert(T, x)) + +"Promote map and point to compatible types." +promote_map_point_pair(m, x) = (m, x) +promote_map_point_pair(m::Map, x) = _promote_map_point_pair(m, x, promote_type(domaintype(m), typeof(x))) +# This is the line where we promote both the map and the point: +_promote_map_point_pair(m, x, ::Type{T}) where {T} = + convert(Map{T}, m), convert(T, x) +_promote_map_point_pair(m, x, ::Type{Any}) = m, x + +# Some exceptions: +# - types match, do nothing +promote_map_point_pair(m::Map{T}, x::T) where {T} = m, x +# - an Any map, do nothing +promote_map_point_pair(m::Map{Any}, x) = m, x +# - tuples: these are typically composite maps, promotion may happen later +promote_map_point_pair(m::Map{<:Tuple}, x::Tuple) = m, x +# - abstract vectors: promotion may be expensive +promote_map_point_pair(m::Map{<:AbstractVector}, x::AbstractVector) = m, x +# - SVector: promotion is likely cheap +promote_map_point_pair(m::EuclideanMap{N,T}, x::AbstractVector{S}) where {N,S,T} = + _promote_map_point_pair(m, x, SVector{N,promote_type(S,T)}) + +# For maps of type Map{T}, we call promote_map_point_pair and then applymap +promote_and_apply(m::Map, x) = applymap(promote_map_point_pair(m,x)...) + +applymap!(y, m, x) = y .= m(x) + +# Fallback for functions that are not of type AbstractMap +applymap(m, x) = m(x) +applymap(m::AbstractMap, x) = error("Please implement applymap for map $(m)") + + +"Can the maps be promoted to a common domain type without throwing an error?" +promotable_maps(maps...) = promotable_eltypes(map(domaintype, maps)...) + +"Promote the given maps to have a common domain type." +promote_maps() = () +promote_maps(m1) = m1 +function promote_maps(m1, m2) + T = promote_type(domaintype(m1), domaintype(m2)) + convert_domaintype(T, m1), convert_domaintype(T, m2) +end +function promote_maps(m1, m2, m3, maps...) + T = mapreduce(domaintype, promote_type, (m1,m2,m3,maps...)) + convert_domaintype.(T, (m1,m2,m3,maps...)) +end + +isvectorvalued_type(::Type{T}) where {T<:Number} = true +isvectorvalued_type(::Type{T}) where {T<:AbstractVector} = true +isvectorvalued_type(::Type{T}) where {T} = false + + +"Is the map a vector-valued function, i.e., a function from Rn to Rm?" +isvectorvalued(m) = + isvectorvalued_type(domaintype(m)) && isvectorvalued_type(codomaintype(m)) + +# 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 +mapsize(m, i) = _mapsize(m, i, mapsize(m)) +_mapsize(m, i, size::Tuple{Int,Int}) = i <= 2 ? size[i] : 1 +_mapsize(m, i, size::Tuple{Int}) = i <= 1 ? size[i] : 1 +_mapsize(m, i, size::Tuple{}) = 1 + +"Is the given map a square map?" +issquaremap(m) = isvectorvalued(m) && (mapsize(m,1) == mapsize(m,2)) + +isoverdetermined(m) = mapsize(m,1) >= mapsize(m,2) +isunderdetermined(m) = mapsize(m,1) <= mapsize(m,2) + +is_scalar_to_vector(m) = mapsize(m) isa Tuple{Int} +is_vector_to_scalar(m) = mapsize(m) isa Tuple{Int,Int} && codomaintype(m)<:Number +is_scalar_to_scalar(m) = mapsize(m) == () +is_vector_to_vector(m) = mapsize(m) isa Tuple{Int,Int} && !is_vector_to_scalar(m) + +# Display routines +map_stencil(m, x) = [Display.SymbolObject(m), '(', x, ')'] +map_stencil_broadcast(m, x) = [Display.SymbolObject(m), ".(", x, ')'] +map_stencil_broadcast(m::Function, x) = [repr(m), ".(", x, ')'] + +Display.object_parentheses(m::Map) = map_object_parentheses(m) +Display.stencil_parentheses(m::Map) = map_stencil_parentheses(m) + +map_object_parentheses(m) = false +map_stencil_parentheses(m) = false diff --git a/src/generic/product.jl b/src/generic/product.jl new file mode 100644 index 0000000..8cc6e03 --- /dev/null +++ b/src/generic/product.jl @@ -0,0 +1,201 @@ + +""" +A product map is diagonal and acts on each of the components of x separately: +`y = f(x)` becomes `y_i = f_i(x_i)`. +""" +abstract type ProductMap{T} <: CompositeLazyMap{T} end + +components(m::ProductMap) = m.maps +factors(d::ProductMap) = components(d) + +VcatMapElement = Union{Map{<:StaticVector},Map{<:Number}} + +ProductMap(maps::Tuple) = ProductMap(maps...) +ProductMap(maps::StaticVector) = ProductMap(maps...) +ProductMap(maps...) = TupleProductMap(maps...) +ProductMap(maps::VcatMapElement...) = VcatMap(maps...) +ProductMap(maps::AbstractVector) = VectorProductMap(maps) + +ProductMap{T}(maps...) where {T} = _TypedProductMap(T, maps...) +_TypedProductMap(::Type{T}, maps...) where {T<:Tuple} = TupleProductMap(maps...) +_TypedProductMap(::Type{SVector{N,T}}, maps...) where {N,T} = VcatMap{T}(maps...) +_TypedProductMap(::Type{T}, maps...) where {T<:AbstractVector} = VectorProductMap{T}(maps...) + +compatibleproductdims(d1::ProductMap, d2::ProductMap) = + mapsize(d1) == mapsize(d2) && + all(map(==, map(mapsize, components(d1)), map(mapsize, components(d2)))) + +isconstantmap(m::ProductMap) = mapreduce(isconstantmap, &, components(m)) +islinearmap(m::ProductMap) = mapreduce(islinearmap, &, components(m)) +isaffinemap(m::ProductMap) = mapreduce(isaffinemap, &, components(m)) + +affinematrix(m::ProductMap) = toexternalmatrix(m, map(affinematrix, components(m))) +affinevector(m::ProductMap) = toexternalpoint(m, map(affinevector, components(m))) +mapconstant(m::ProductMap) = toexternalpoint(m, map(mapconstant, components(m))) + +jacobian(m::ProductMap, x) = + toexternalmatrix(m, map(jacobian, components(m), tointernalpoint(m, x))) +function jacobian(m::ProductMap{T}) where {T} + if isaffinemap(m) + ConstantMap{T}(affinematrix(m)) + else + LazyJacobian(m) + end +end + +# diffvolume(m::ProductMap) = multiply_map(map(diffvolume, factors(m))...) + +similarmap(m::ProductMap, ::Type{T}) where {T} = ProductMap{T}(components(m)) + +tointernalpoint(m::ProductMap, x) = x +toexternalpoint(m::ProductMap, y) = y + +applymap(m::ProductMap, x) = + toexternalpoint(m, map(applymap, components(m), tointernalpoint(m, x))) + +productmap(m1, m2) = productmap1(m1, m2) +productmap1(m1, m2) = productmap2(m1, m2) +productmap2(m1, m2) = default_productmap(m1, m2) +default_productmap(m1, m2) = ProductMap(m1, m2) + +productmap(m1::ProductMap, m2::ProductMap) = + ProductMap(components(m1)..., components(m2)...) +productmap1(m1::ProductMap, m2) = ProductMap(components(m1)..., m2) +productmap2(m1, m2::ProductMap) = ProductMap(m1, components(m2)...) + +for op in (:inverse, :leftinverse, :rightinverse) + @eval $op(m::ProductMap) = ProductMap(map($op, components(m))) + @eval $op(m::ProductMap, x) = toexternalpoint(m, map($op, components(m), tointernalpoint(m, x))) +end + +function composedmap(m1::ProductMap, m2::ProductMap) + if compatibleproductdims(m1, m2) + ProductMap(map(composedmap, components(m1), components(m2))) + else + ComposedMap(m1,m2) + end +end + +mapsize(m::ProductMap) = (sum(t->mapsize(t,1), components(m)), sum(t->mapsize(t,2), components(m))) + +isequalmap(m1::ProductMap, m2::ProductMap) = all(map(isequalmap, components(m1), components(m2))) +map_hash(m::ProductMap, h::UInt) = hashrec("ProductMap", collect(components(m)), h) + +canonicalmap(m::ProductMap) = any(map(hascanonicalmap, factors(m))) ? + ProductMap(map(canonicalmap, factors(m))) : m +canonicalmap(::Equal, m::ProductMap) = any(map(hasequalmap, factors(m))) ? + ProductMap(map(equalmap, factors(m))) : m +canonicalmap(::Equivalent, m::ProductMap) = any(map(hasequivalentmap, factors(m))) ? + ProductMap(map(equivalentmap, factors(m))) : m + +Display.combinationsymbol(m::ProductMap) = Display.Symbol('⊗') +Display.displaystencil(m::ProductMap) = composite_displaystencil(m) +map_object_parentheses(m::ProductMap) = true +map_stencil_parentheses(m::ProductMap) = true +show(io::IO, mime::MIME"text/plain", m::ProductMap) = composite_show(io, mime, m) +show(io::IO, m::ProductMap) = composite_show_compact(io, m) + +""" +A `VcatMap` is a product map with domain and codomain vectors +concatenated (`vcat`) into a single vector. +""" +struct VcatMap{T,M,N,DIM1,DIM2,MAPS} <: ProductMap{SVector{N,T}} + maps :: MAPS +end + +VcatMap(maps::Union{Tuple,Vector}) = VcatMap(maps...) +VcatMap(maps...) = VcatMap{numtype(maps...)}(maps...) + +VcatMap{T}(maps::Union{Tuple,Vector}) where T = VcatMap{T}(maps...) +function VcatMap{T}(maps...) where T + M = sum(t->mapsize(t,1), maps) + N = sum(t->mapsize(t,2), maps) + 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) + DIM2 = map(t->mapsize(t,2), maps) + VcatMap{T,M,N,DIM1,DIM2}(convert_numtype.(Ref(T), maps)...) +end + +VcatMap{T,M,N,DIM1,DIM2}(maps...) where {T,M,N,DIM1,DIM2} = + VcatMap{T,M,N,DIM1,DIM2,typeof(maps)}(maps) + +mapsize(m::VcatMap{T,M,N}) where {T,M,N} = (M,N) + +tointernalpoint(m::VcatMap{T,M,N,DIM1,DIM2}, x) where {T,M,N,DIM1,DIM2} = + convert_fromcartesian(x, Val{DIM2}()) +toexternalpoint(m::VcatMap{T,M,N,DIM1,DIM2}, y) where {T,M,N,DIM1,DIM2} = + convert_tocartesian(y, Val{DIM1}()) + +size_as_matrix(A::AbstractArray) = size(A) +size_as_matrix(A::Number) = (1,1) + +# The Jacobian is block-diagonal +function toexternalmatrix(m::VcatMap{T,M,N}, matrices) where {T,M,N} + A = zeros(T, M, N) + k = 0 + l = 0 + for el in matrices + m,n = size_as_matrix(el) + A[k+1:k+m,l+1:l+n] .= el + k += m + l += n + end + SMatrix{M,N}(A) +end + + +""" +A `VectorProductMap` is a product map where all components are univariate maps, +with inputs and outputs collected into a `Vector`. +""" +struct VectorProductMap{T<:AbstractVector,M} <: ProductMap{T} + maps :: Vector{M} +end + +VectorProductMap(maps...) = VectorProductMap(maps) +VectorProductMap(maps) = VectorProductMap(collect(maps)) +function VectorProductMap(maps::Vector) + T = mapreduce(numtype, promote_type, maps) + VectorProductMap{Vector{T}}(maps) +end + +VectorProductMap{T}(maps...) where {T} = VectorProductMap{T}(maps) +VectorProductMap{T}(maps) where {T} = VectorProductMap{T}(collect(maps)) +function VectorProductMap{T}(maps::Vector) where {T} + Tmaps = convert.(Map{eltype(T)}, maps) + VectorProductMap{T,eltype(Tmaps)}(Tmaps) +end + +# the Jacobian is a diagonal matrix +toexternalmatrix(m::VectorProductMap, matrices) = Diagonal(matrices) + +mapsize(m::VectorProductMap) = (length(m.maps), length(m.maps)) + +""" +A `TupleProductMap` is a product map with all components collected in a tuple. +There is no vector-valued function associated with this map. +""" +struct TupleProductMap{T<:Tuple,MM} <: ProductMap{T} + maps :: MM +end + +TupleProductMap(maps::Vector) = TupleProductMap(maps...) +TupleProductMap(maps...) = TupleProductMap(maps) +function TupleProductMap(maps::Tuple) + T = Tuple{map(domaintype, maps)...} + TupleProductMap{T}(maps) +end + +TupleProductMap{T}(maps::Vector) where {T<:Tuple} = TupleProductMap{T}(maps...) +TupleProductMap{T}(maps...) where {T<:Tuple} = TupleProductMap{T}(maps) +function TupleProductMap{T}(maps::NTuple{N,<:AbstractMap}) where {N,T <: Tuple} + Tmaps = map((t,d) -> convert(Map{t},d), tuple(T.parameters...), maps) + TupleProductMap{T,typeof(Tmaps)}(Tmaps) +end +TupleProductMap{T}(maps) where {T <: Tuple} = TupleProductMap{T,typeof(maps)}(maps) diff --git a/src/types/affine.jl b/src/types/affine.jl new file mode 100644 index 0000000..a7b97d9 --- /dev/null +++ b/src/types/affine.jl @@ -0,0 +1,589 @@ + +""" +An affine map has the general form `y = A*x + b`. + +We use `affinematrix(m)` and `affinevector(m)` to denote `A` and `b` respectively. Concrete +subtypes include linear maps of the form `y = A*x` and translations of the +form `y = x + b`. +""" +abstract type AbstractAffineMap{T} <: Map{T} end + +unsafe_matrix(m::AbstractAffineMap) = m.A +unsafe_vector(m::AbstractAffineMap) = m.b + +"Return the matrix `A` in the affine map `Ax+b`." +affinematrix(m::AbstractAffineMap) = to_matrix(domaintype(m), unsafe_matrix(m), unsafe_vector(m)) + +"Return the vector `b` in the affine map `Ax+b`." +affinevector(m::AbstractAffineMap) = to_vector(domaintype(m), unsafe_matrix(m), unsafe_vector(m)) + +applymap(m::AbstractAffineMap, x) = _affine_applymap(m, x, unsafe_matrix(m), unsafe_vector(m)) +_affine_applymap(m, x, A, b) = A*x + b + +applymap!(y, m::AbstractAffineMap, x) = _affine_applymap!(y, m, x, unsafe_matrix(m), unsafe_vector(m)) +function _affine_applymap!(y, m, x, A, b) + mul!(y, A, x) + y .+= b + y +end + +_affine_A_isreal(A) = isreal(A) +_affine_A_isreal(::UniformScaling{T}) where T = isrealtype(T) + +isrealmap(m::AbstractAffineMap) = _affine_isrealmap(m, unsafe_matrix(m), unsafe_vector(m)) +_affine_isrealmap(m, A, b) = _affine_A_isreal(A) && isreal(b) + +jacobian(m::AbstractAffineMap{T}) where {T} = ConstantMap{T}(affinematrix(m)) +jacobian(m::AbstractAffineMap, x) = affinematrix(m) + +jacdet(m::AbstractAffineMap, x) = _affine_jacdet(m, x, unsafe_matrix(m)) +_affine_jacdet(m, x, A) = det(A) +_affine_jacdet(m, x::Number, A::UniformScaling) = A.λ +_affine_jacdet(m, x::AbstractVector, A::Number) = A^length(x) +_affine_jacdet(m, x::AbstractVector, A::UniformScaling) = A.λ^length(x) + +function diffvolume(m::AbstractAffineMap{T}) where T + J = jacobian(m) + c = sqrt(det(affinevector(J)'*affinevector(J))) + ConstantMap{T}(c) +end + +islinearmap(m::AbstractMap) = false +islinearmap(m::AbstractAffineMap) = _affine_islinearmap(m, unsafe_vector(m)) +_affine_islinearmap(m, b) = all(b .== 0) + +isaffinemap(m) = false +isaffinemap(m::Map) = islinearmap(m) || isconstantmap(m) +isaffinemap(m::AbstractAffineMap) = true + +isequalmap(m1::AbstractAffineMap, m2::AbstractAffineMap) = + affinematrix(m1) == affinematrix(m2) && affinevector(m1) == affinevector(m2) + +isequalmap(m1::AbstractAffineMap, m2::IdentityMap) = + islinearmap(m1) && affinematrix(m1) == affinematrix(m2) +isequalmap(m1::IdentityMap, m2::AbstractAffineMap) = isequalmap(m2, m1) + +map_hash(m::AbstractAffineMap, h::UInt) = hashrec("AbstractAffineMap", affinematrix(m), affinevector(m), h) + +mapsize(m::AbstractAffineMap) = _affine_mapsize(m, domaintype(m), unsafe_matrix(m), unsafe_vector(m)) +_affine_mapsize(m, T, A::AbstractArray, b) = size(A) +_affine_mapsize(m, T, A::AbstractVector, b::AbstractVector) = (length(A),) +_affine_mapsize(m, T, A::Number, b::Number) = () +_affine_mapsize(m, T, A::Number, b::AbstractVector) = (length(b),length(b)) +_affine_mapsize(m, T, A::UniformScaling, b::Number) = () +_affine_mapsize(m, T, A::UniformScaling, b) = (length(b),length(b)) + + +Display.displaystencil(m::AbstractAffineMap) = vcat(["x -> "], map_stencil(m, 'x')) +show(io::IO, mime::MIME"text/plain", m::AbstractAffineMap) = composite_show(io, mime, m) + +map_stencil(m::AbstractAffineMap, x) = _affine_map_stencil(m, x, unsafe_matrix(m), unsafe_vector(m)) +_affine_map_stencil(m, x, A, b) = [A, " * ", x, " + ", b] +_affine_map_stencil(m, x, A, b::Real) = + b >= 0 ? [A, " * ", x, " + ", b] : [A, " * ", x, " - ", abs(b)] + +map_stencil_broadcast(m::AbstractAffineMap, x) = _affine_map_stencil_broadcast(m, x, unsafe_matrix(m), unsafe_vector(m)) +_affine_map_stencil_broadcast(m, x, A, b) = [A, " .* ", x, " .+ ", b] +_affine_map_stencil_broadcast(m, x, A::Number, b) = [A, " * ", x, " .+ ", b] + +map_object_parentheses(m::AbstractAffineMap) = true +map_stencil_parentheses(m::AbstractAffineMap) = true + +######################## +# Linear maps: y = A*x +######################## + +""" +The supertype of all linear maps `y = A*x`. +Concrete subtypes may differ in how `A` is represented. +""" +abstract type LinearMap{T} <: AbstractAffineMap{T} end + +mapsize(m::LinearMap) = _linearmap_size(m, domaintype(m), unsafe_matrix(m)) +_linearmap_size(m, T, A) = size(A) +_linearmap_size(m, ::Type{T}, A::Number) where {T<:Number} = () +_linearmap_size(m, ::Type{T}, A::AbstractVector) where {T<:Number} = (length(A),) +_linearmap_size(m, ::Type{T}, A::Number) where {N,T<:StaticVector{N}} = (N,N) + +affinematrix(m::LinearMap) = to_matrix(domaintype(m), unsafe_matrix(m)) +affinevector(m::LinearMap) = to_vector(domaintype(m), unsafe_matrix(m)) + +LinearMap(A::Number) = ScalarLinearMap(A) +LinearMap(A::SMatrix) = StaticLinearMap(A) +LinearMap(A::MMatrix) = StaticLinearMap(A) +LinearMap(A::Matrix) = VectorLinearMap(A) +LinearMap(A) = GenericLinearMap(A) + +LinearMap{T}(A::Number) where {T<:Number} = _LinearMap(A, T, promote_type(T,typeof(A))) +_LinearMap(A::Number, ::Type{T}, ::Type{T}) where {T<:Number} = ScalarLinearMap{T}(A) +_LinearMap(A::Number, ::Type{T}, ::Type{S}) where {T<:Number,S} = GenericLinearMap{T}(A) +LinearMap{T}(A::SMatrix{M,N}) where {M,N,S,T <: SVector{N,S}} = StaticLinearMap{S}(A) +LinearMap{T}(A::MMatrix{M,N}) where {M,N,S,T <: SVector{N,S}} = StaticLinearMap{S}(A) +LinearMap{T}(A::Matrix) where {S,T <: Vector{S}} = VectorLinearMap{S}(A) +LinearMap{T}(A) where {T} = GenericLinearMap{T}(A) + +# convenience functions +LinearMap(a::Number...) = LinearMap(promote(a...)) +LinearMap(a::NTuple{N,T}) where {N,T <: Number} = LinearMap{SVector{N,T}}(Diagonal(SVector{N,T}(a))) + +applymap(m::LinearMap, x) = _linear_applymap(m, x, unsafe_matrix(m)) +_linear_applymap(m, x, A) = A*x + +applymap!(y, m::LinearMap, x) = _linear_applymap!(y, m, x, unsafe_matrix(m)) +_linear_applymap!(y, m, x, A) = mul!(y, A, x) + +islinearmap(m::LinearMap) = true + +isrealmap(m::LinearMap) = _linear_isrealmap(m, unsafe_matrix(m)) +_linear_isrealmap(m, A) = isreal(A) + +isequalmap(m1::LinearMap, m2::LinearMap) = affinematrix(m1) == affinematrix(m2) + +# inverse should be called only on square maps, otherwise use +# leftinverse or rightinverse in order to use pinv instead of inv +inverse(m::LinearMap) = (@assert issquaremap(m); LinearMap(inv(m.A))) +inverse(m::LinearMap, x) = (@assert issquaremap(m); m.A \ x) + +function leftinverse(m::LinearMap) + @assert isoverdetermined(m) + LinearMap(matrix_pinv(m.A)) +end +function rightinverse(m::LinearMap) + @assert isunderdetermined(m) + LinearMap(matrix_pinv(m.A)) +end +function leftinverse(m::LinearMap, x) + @assert isoverdetermined(m) + m.A \ x +end +function rightinverse(m::LinearMap, x) + @assert isunderdetermined(m) + m.A \ x +end + +similarmap(m::LinearMap, ::Type{T}) where {T} = LinearMap{T}(m.A) + +convert(::Type{Map}, a::Number) = LinearMap(a) +convert(::Type{Map{T}}, a::Number) where {T} = LinearMap{T}(a) + +convert(::Type{LinearMap}, ::StaticIdentityMap{T}) where {T} = LinearMap{T}(1) +convert(::Type{LinearMap{T}}, ::StaticIdentityMap) where {T} = LinearMap{T}(1) +convert(::Type{AbstractAffineMap}, m::StaticIdentityMap) = convert(LinearMap, m) +convert(::Type{AbstractAffineMap{T}}, m::StaticIdentityMap) where {T} = convert(LinearMap, m) + + +map_stencil(m::LinearMap, x) = [unsafe_matrix(m), " * ", x] +map_stencil_broadcast(m::LinearMap, x) = _linear_map_stencil_broadcast(m, x, unsafe_matrix(m)) +_linear_map_stencil_broadcast(m, x, A) = [A, " .* ", x] +_linear_map_stencil_broadcast(m, x, A::Number) = [A, " * ", x] + + +"A `GenericLinearMap` is a linear map `y = A*x` for any type of `A`." +struct GenericLinearMap{T,AA} <: LinearMap{T} + A :: AA +end + +""" +What is the suggested domaintype for a generic linear map `A*x` with +the given argument 'A'? +""" +glm_domaintype(A) = Any +glm_domaintype(A::Number) = typeof(A) +glm_domaintype(A::AbstractMatrix{T}) where T = Vector{T} +glm_domaintype(A::StaticMatrix{M,N,T}) where {M,N,T} = SVector{N,T} +glm_domaintype(A::AbstractVector{T}) where {T} = T +glm_domaintype(A::Diagonal{T,<:StaticVector{N,T}}) where {N,T} = SVector{N,T} + +GenericLinearMap(A) = GenericLinearMap{glm_domaintype(A)}(A) + +# Allow any A +GenericLinearMap{T}(A) where {T} = GenericLinearMap{T,typeof(A)}(A) + +# Promote some eltypes if applicable +GenericLinearMap{T}(A::Number) where {T <: Number} = + _GenericLinearMap(A, T, promote_type(T,typeof(A))) +_GenericLinearMap(A::Number, ::Type{T}, ::Type{T}) where {T <: Number} = + GenericLinearMap{T,T}(A) +_GenericLinearMap(A::Number, ::Type{T}, ::Type{S}) where {S,T<:Number} = + GenericLinearMap{T,typeof(A)}(A) + +GenericLinearMap{T}(A::Number) where {T <: AbstractVector} = + _GenericLinearMap(A, T, promote_type(eltype(T),typeof(A))) +_GenericLinearMap(A::Number, ::Type{T}, ::Type{S}) where {S,T<:AbstractVector{S}} = + GenericLinearMap{T,S}(A) +_GenericLinearMap(A::Number, ::Type{T}, ::Type{S}) where {S,T<:AbstractVector} = + GenericLinearMap{T,typeof(A)}(A) + +GenericLinearMap{T}(A::AbstractMatrix{S}) where {S,T <: AbstractVector{S}} = + GenericLinearMap{T,typeof(A)}(A) +GenericLinearMap{T}(A::AbstractMatrix{U}) where {S,T <: AbstractVector{S},U} = + GenericLinearMap{T}(convert(AbstractMatrix{S}, A)) +GenericLinearMap{T}(A::AbstractVector{T}) where {T} = + GenericLinearMap{T,typeof(A)}(A) +GenericLinearMap{T}(A::AbstractVector{S}) where {S,T} = + GenericLinearMap{T}(convert(AbstractVector{T}, A)) +GenericLinearMap{T}(A::AbstractMatrix) where {T<:Number} = + throw(ArgumentError("Linear map with matrix A can not have scalar domaintype.")) + +# Preserve the action on vectors with a number type +inverse(m::GenericLinearMap{T,AA}) where {T<:AbstractVector,AA<:Number} = + LinearMap{T}(inv(m.A)) +leftinverse(m::GenericLinearMap{T,AA}) where {T<:AbstractVector,AA<:Number} = + LinearMap{T}(inv(m.A)) +rightinverse(m::GenericLinearMap{T,AA}) where {T<:AbstractVector,AA<:Number} = + LinearMap{T}(inv(m.A)) + +convert(::Type{Map}, A::UniformScaling) = GenericLinearMap{Vector{Any}}(A) +convert(::Type{Map{T}}, A::UniformScaling) where {T} = GenericLinearMap{T}(A) + +"A `ScalarLinearMap` is a linear map `y = A*x` for scalars." +struct ScalarLinearMap{T} <: LinearMap{T} + A :: T +end + +ScalarLinearMap(A::Number) = ScalarLinearMap{typeof(A)}(A) + +isrealmap(m::ScalarLinearMap{T}) where {T} = isrealtype(T) + +show(io::IO, m::ScalarLinearMap) = show_scalar_linear_map(io, m.A) +show_scalar_linear_map(io, A::Real) = print(io, "x -> $(A) * x") +show_scalar_linear_map(io, A::Complex) = print(io, "x -> ($(A)) * x") + + +"A `VectorLinearMap` is a linear map `y = A*x` using vectors and matrices." +struct VectorLinearMap{T} <: LinearMap{Vector{T}} + A :: Matrix{T} +end + +VectorLinearMap(A::AbstractMatrix{T}) where {T} = + VectorLinearMap{T}(A) + + +"A `StaticLinearMap` is a linear map `y = A*x` using static arrays." +struct StaticLinearMap{T,N,M,L} <: LinearMap{SVector{N,T}} + A :: SMatrix{M,N,T,L} +end + +StaticLinearMap(A::AbstractMatrix{T}) where {T} = + StaticLinearMap{T}(A) + +StaticLinearMap{T}(A::StaticMatrix{M,N,S}) where {M,N,T,S} = + StaticLinearMap{T}(convert(AbstractMatrix{T}, A)) +StaticLinearMap{T}(A::StaticMatrix{M,N,T}) where {M,N,T} = + StaticLinearMap{T,M,N}(A) +StaticLinearMap{T,N,M}(A::AbstractMatrix) where {T,N,M} = + StaticLinearMap{T,N,M,M*N}(A) + +convert(::Type{Map{SVector{N,T}}}, m::VectorLinearMap) where {N,T} = StaticLinearMap{T,N,N}(m.A) + + +# Implement the interface for abstract arrays, +# representing the linear map x->A*x +MapStyle(A::AbstractArray) = IsMap() + +Map(A::AbstractArray) = LinearMap(A) +Map{T}(A::AbstractArray) where T = LinearMap{T}(A) + +domaintype(A::AbstractArray) = domaintype(Map(A)) + +applymap(A::AbstractArray, x) = A*x +mapsize(A::AbstractArray) = size(A) + +islinearmap(A::AbstractArray) = true +isaffinemap(A::AbstractArray) = true +affinematrix(A::AbstractArray) = A +affinevector(A::AbstractArray) = zerovector(A) + +inverse(A::AbstractMatrix) = inv(A) +inverse(A::AbstractMatrix, x) = A \ x + +jacobian(A::AbstractMatrix) = ConstantMap{glm_domaintype(A)}(A) +jacobian(A::AbstractMatrix, x) = A + +canonicalmap(A::AbstractArray) = Map(A) +canonicalmap(::Equal, A::AbstractArray) = Map(A) + + + + +########################## +# Translations: y = x + b +########################## + + +"A `Translation` represents the map `y = x + b`." +abstract type Translation{T} <: AbstractAffineMap{T} end + +# unsafe_matrix(m::Translation) = I +affinematrix(m::Translation) = to_matrix(domaintype(m), LinearAlgebra.I, unsafe_vector(m)) +affinevector(m::Translation) = unsafe_vector(m) + +mapsize(m::Translation) = _translation_mapsize(m, domaintype(m), unsafe_vector(m)) +_translation_mapsize(m, ::Type{T}, b::Number) where {T} = () +_translation_mapsize(m, ::Type{T}, b) where {T} = (length(b),length(b)) + +map_stencil(m::Translation, x) = _translation_map_stencil(m, x, unsafe_vector(m)) +_translation_map_stencil(m, x, b) = [x, " + ", b] +_translation_map_stencil(m, x, b::Real) = + b >= 0 ? [x, " + ", b] : [x, " - ", abs(b)] +map_stencil_broadcast(m::Translation, x) = _translation_map_stencil_broadcast(m, x, unsafe_vector(m)) +_translation_map_stencil_broadcast(m, x, b) = [x, " .+ ", b] + +"Translation by a scalar value." +struct ScalarTranslation{T} <: Translation{T} + b :: T +end + +isrealmap(m::ScalarTranslation{T}) where {T} = isrealtype(T) + +show(io::IO, m::ScalarTranslation) = show_scalar_translation(io, m.b) +show_scalar_translation(io, b::Real) = print(io, "x -> x", b < 0 ? " - " : " + ", abs(b)) +show_scalar_translation(io, b) = print(io, "x -> x + ", b) + + +"Translation by a static vector." +struct StaticTranslation{T,N} <: Translation{SVector{N,T}} + b :: SVector{N,T} +end + +"Translation by a vector." +struct VectorTranslation{T} <: Translation{Vector{T}} + b :: Vector{T} +end + +"Translation by a generic vectorlike object." +struct GenericTranslation{T,B} <: Translation{T} + b :: B +end + +Translation(b::Number) = ScalarTranslation(b) +Translation(b::StaticVector) = StaticTranslation(b) +Translation(b::Vector) = VectorTranslation(b) +Translation(b) = GenericTranslation(b) + +Translation{T}(b::Number) where {T<:Number} = ScalarTranslation{T}(b) +Translation{T}(b::AbstractVector) where {N,S,T<:StaticVector{N,S}} = StaticTranslation{S,N}(b) +Translation{T}(b::Vector) where {S,T<:Vector{S}} = VectorTranslation{S}(b) +Translation{T}(b) where {T} = GenericTranslation{T}(b) + +jacdet(m::Translation{T}) where {T} = UnityMap{T,eltype(T)}() +jacdet(m::Translation{T}, x) where {T} = one(eltype(T)) + +isrealmap(m::Translation) = isreal(unsafe_vector(m)) + +isequalmap(m1::Translation, m2::Translation) = unsafe_vector(m1)==unsafe_vector(m2) + +similarmap(m::Translation, ::Type{T}) where {T} = Translation{T}(m.b) + +applymap(m::Translation, x) = _translation_applymap(m, x, unsafe_vector(m)) +_translation_applymap(m, x, b) = x + b +applymap!(y, m::Translation, x) = _translation_applymap!(y, m, x, unsafe_vector(m)) +_translation_applymap!(y, m, x, b) = y .= x .+ m.b + +inverse(m::Translation{T}) where {T} = Translation{T}(-m.b) +inverse(m::Translation, x) = x - m.b + +ScalarTranslation(b::Number) = ScalarTranslation{typeof(b)}(b) + +StaticTranslation(b::AbstractVector{T}) where {T} = StaticTranslation{T}(b) + +StaticTranslation{T}(b::StaticVector{N}) where {N,T} = + StaticTranslation{T,N}(b) + +VectorTranslation(b::AbstractVector{T}) where {T} = VectorTranslation{T}(b) + +GenericTranslation(b) = GenericTranslation{typeof(b)}(b) +GenericTranslation(b::AbstractVector{T}) where {T} = + GenericTranslation{Vector{T}}(b) +GenericTranslation(b::StaticVector{N,T}) where {N,T} = + GenericTranslation{SVector{N,T}}(b) + +GenericTranslation{T}(b) where {T} = GenericTranslation{T,typeof(b)}(b) +GenericTranslation{T}(b::Number) where {T<:Number} = + GenericTranslation{T,T}(b) + + +############################ +# Affine maps: y = A*x + b +############################ + + +""" +The supertype of all affine maps that store `A` and `b`. +Concrete subtypes differ in how `A` and `b` are represented. +""" +abstract type AffineMap{T} <: AbstractAffineMap{T} end + +AffineMap(A::Number, b::Number) = ScalarAffineMap(A, b) +AffineMap(A::StaticMatrix, b::StaticVector) = StaticAffineMap(A, b) +AffineMap(A::Matrix, b::Vector) = VectorAffineMap(A, b) +AffineMap(A::UniformScaling{Bool}, b::Number) = ScalarAffineMap(one(b), b) +AffineMap(A, b) = GenericAffineMap(A, b) + +AffineMap{T}(A::Number, b::Number) where {T<:Number} = ScalarAffineMap{T}(A, b) +AffineMap{T}(A::AbstractMatrix, b::AbstractVector) where {N,S,T<:SVector{N,S}} = StaticAffineMap{S,N}(A, b) +AffineMap{T}(A::Matrix, b::Vector) where {S,T<:Vector{S}} = VectorAffineMap{S}(A, b) +AffineMap{T}(A::UniformScaling{Bool}, b::Number) where {T} = ScalarAffineMap{T}(one(T), b) +AffineMap{T}(A, b) where {T} = GenericAffineMap{T}(A, b) + +similarmap(m::AffineMap, ::Type{T}) where {T} = AffineMap{T}(m.A, m.b) + +convert(::Type{AffineMap}, m) = (@assert isaffinemap(m); AffineMap(affinematrix(m), affinevector(m))) +convert(::Type{AffineMap{T}}, m) where {T} = (@assert isaffinemap(m); AffineMap{T}(affinematrix(m), affinevector(m))) +# avoid ambiguity errors with convert(::Type{T}, x::T) in Base: +convert(::Type{AffineMap}, m::AffineMap) = m +convert(::Type{AffineMap{T}}, m::AffineMap{T}) where T = m + +# If y = A*x+b, then x = inv(A)*(y-b) = inv(A)*y - inv(A)*b +inverse(m::AffineMap) = (@assert issquaremap(m); AffineMap(inv(m.A), -inv(m.A)*m.b)) +inverse(m::AffineMap, x) = (@assert issquaremap(m); m.A \ (x-m.b)) + +function leftinverse(m::AffineMap) + @assert isoverdetermined(m) + pA = matrix_pinv(m.A) + AffineMap(pA, -pA*m.b) +end +function rightinverse(m::AffineMap) + @assert isunderdetermined(m) + pA = matrix_pinv(m.A) + AffineMap(pA, -pA*m.b) +end +function leftinverse(m::AffineMap, x) + @assert isoverdetermined(m) + m.A \ (x-m.b) +end +function rightinverse(m::AffineMap, x) + @assert isunderdetermined(m) + m.A \ (x-m.b) +end + + +"An affine map for any combination of types of `A` and `b`." +struct GenericAffineMap{T,AA,B} <: AffineMap{T} + A :: AA + b :: B +end + +GenericAffineMap(A, b) = GenericAffineMap{typeof(b)}(A, b) +GenericAffineMap(A::AbstractVector{S}, b::AbstractVector{T}) where {S,T} = + GenericAffineMap{promote_type(S,T)}(A, b) +GenericAffineMap(A::AbstractArray{S}, b::AbstractVector{T}) where {S,T} = + GenericAffineMap{Vector{promote_type(S,T)}}(A, b) +GenericAffineMap(A::StaticMatrix{M,N,S}, b::StaticVector{M,T}) where {M,N,S,T} = + GenericAffineMap{SVector{N,promote_type(S,T)}}(A, b) +GenericAffineMap(A::StaticMatrix{M,N,S}, b::AbstractVector{T}) where {M,N,S,T} = + GenericAffineMap{SVector{N,promote_type(S,T)}}(A, b) +GenericAffineMap(A::S, b::AbstractVector{T}) where {S<:Number,T} = + GenericAffineMap{Vector{promote_type(S,T)}}(A, b) +GenericAffineMap(A::S, b::StaticVector{N,T}) where {S<:Number,N,T} = + GenericAffineMap{SVector{N,promote_type(S,T)}}(A, b) +GenericAffineMap(A::UniformScaling{Bool}, b) = + GenericAffineMap(UniformScaling{eltype(b)}(1), b) + + +# Fallback routine for generic A and b, special cases follow +GenericAffineMap{T}(A, b) where {T} = GenericAffineMap{T,typeof(A),typeof(b)}(A, b) + +GenericAffineMap{T}(A::AbstractVector{S}, b::AbstractVector{U}) where {T<:Number,S,U} = + GenericAffineMap{T}(convert(AbstractVector{T}, A), convert(AbstractVector{T}, b)) +GenericAffineMap{T}(A::AbstractVector{T}, b::AbstractVector{T}) where {T<:Number} = + GenericAffineMap{T,typeof(A),typeof(b)}(A, b) +GenericAffineMap{T}(A::Number, b) where {T} = GenericAffineMap{T,eltype(T),typeof(b)}(A, b) +GenericAffineMap{T}(A::Number, b::AbstractVector) where {N,S,T <: StaticVector{N,S}} = + GenericAffineMap{T,S,SVector{N,S}}(A, b) +# Promote element types of abstract arrays +GenericAffineMap{T}(A::AbstractArray, b::AbstractVector) where {S,T<:AbstractVector{S}} = + GenericAffineMap{T}(convert(AbstractArray{eltype(T)},A), convert(AbstractVector{eltype(T)}, b)) +GenericAffineMap{T}(A::AbstractArray{S}, b::AbstractVector{S}) where {S,T<:AbstractVector{S}} = + GenericAffineMap{T,typeof(A),typeof(b)}(A, b) +GenericAffineMap{T}(A::UniformScaling{Bool}, b::AbstractVector) where {S,T<:AbstractVector{S}} = + GenericAffineMap{T}(A*one(S), convert(AbstractVector{S}, b)) +GenericAffineMap{T}(A::UniformScaling{S}, b::AbstractVector{S}) where {S,T<:AbstractVector{S}} = + GenericAffineMap{T,typeof(A),typeof(b)}(A, b) + + +similarmap(m::GenericAffineMap, ::Type{T}) where {T} = AffineMap{T}(m.A, m.b) + +convert(::Type{GenericAffineMap{T}}, m::GenericAffineMap) where {T} = + GenericAffineMap{T}(m.A, m.b) + + + +"An affine map with scalar representation." +struct ScalarAffineMap{T} <: AffineMap{T} + A :: T + b :: T +end + +ScalarAffineMap(A, b) = ScalarAffineMap(promote(A, b)...) + +isrealmap(m::ScalarAffineMap{T}) where {T} = isrealtype(T) + +show(io::IO, m::ScalarAffineMap) = show_scalar_affine_map(io, m.A, m.b) +show_scalar_affine_map(io, A::Real, b::Real) = print(io, "x -> $(A) * x", b < 0 ? " - " : " + ", abs(b)) +show_scalar_affine_map(io, A::Complex, b::Complex) = print(io, "x -> ($(A)) * x + ", b) +show_scalar_affine_map(io, A, b) = print(io, "x -> ($(A)) * x + $(b)") + + +convert(::Type{ScalarAffineMap{T}}, m::ScalarAffineMap) where {T} = + ScalarAffineMap{T}(m.A, m.b) + +"An affine map with array and vector representation." +struct VectorAffineMap{T} <: AffineMap{Vector{T}} + A :: Matrix{T} + b :: Vector{T} +end + +VectorAffineMap(A::AbstractArray{T}, b::AbstractVector{T}) where {T} = + VectorAffineMap{T}(A, b) +function VectorAffineMap(A::AbstractArray{S}, b::AbstractVector{T}) where {S,T} + U = promote_type(S,T) + VectorAffineMap(convert(AbstractArray{U}, A), convert(AbstractVector{U}, b)) +end + +convert(::Type{VectorAffineMap{T}}, m::VectorAffineMap) where {T} = + VectorAffineMap{T}(m.A, m.b) + + + +"An affine map with representation using static arrays." +struct StaticAffineMap{T,N,M,L} <: AffineMap{SVector{N,T}} + A :: SMatrix{M,N,T,L} + b :: SVector{M,T} +end + +# Constructors: +# - first, we deduce T +StaticAffineMap(A::AbstractMatrix{T}, b::AbstractVector{T}) where {T} = + StaticAffineMap{T}(A, b) +function StaticAffineMap(A::AbstractMatrix{S}, b::AbstractVector{T}) where {S,T} + U = promote_type(S,T) + StaticAffineMap(convert(AbstractMatrix{U}, A), convert(AbstractVector{U}, b)) +end + +StaticAffineMap{T}(A::AbstractMatrix, b::AbstractVector) where {T} = + StaticAffineMap{T}(convert(AbstractMatrix{T}, A), convert(AbstractVector{T}, b)) + +# - then, we determine N and/or M, from the arguments +function StaticAffineMap{T}(A::AbstractMatrix{T}, b::StaticVector{M,T}) where {T,M} + @assert size(A) == (M,M) + StaticAffineMap{T,M,M}(A, b) +end +StaticAffineMap{T}(A::StaticMatrix{M,N,T}, b::AbstractVector) where {T,N,M} = + StaticAffineMap{T,N,M}(A, b) +StaticAffineMap{T}(A::StaticMatrix{M,N,T}, b::StaticVector{M,T}) where {T,N,M} = + StaticAffineMap{T,N,M}(A, b) +# line below catches ambiguity error +StaticAffineMap{T}(A::StaticMatrix{M1,N,T}, b::StaticVector{M2,T}) where {T,N,M1,M2} = + throw(ArgumentError("Non-matching dimensions")) +StaticAffineMap{T,N}(A::AbstractMatrix, b::AbstractVector) where {T,N} = + StaticAffineMap{T,N,N}(A, b) +StaticAffineMap{T,N}(A::StaticMatrix{M,N}, b::AbstractVector) where {T,N,M} = + StaticAffineMap{T,N,M}(A, b) + +# - finally invoke the constructor (and implicitly convert the data if necessary) +StaticAffineMap{T,N,M}(A::AbstractMatrix, b::AbstractVector) where {T,N,M} = + StaticAffineMap{T,N,M,M*N}(A, b) + +convert(::Type{Map{SVector{N,T}}}, m::VectorAffineMap) where {N,T} = + StaticAffineMap{T,N}(m.A, m.b) diff --git a/src/types/arithmetics.jl b/src/types/arithmetics.jl new file mode 100644 index 0000000..d768b12 --- /dev/null +++ b/src/types/arithmetics.jl @@ -0,0 +1,141 @@ + +# Simplifications go here + +# Basic maps + +composedmap1(m1::IdentityMap, m2) = m2 +composedmap1(m1::IdentityMap{T}, m2::Map{T}) where {T} = m2 +composedmap1(m1::IdentityMap{T}, m2::Map{S}) where {S,T} = convert(Map{T}, m2) + +composedmap2(m1, m2::IdentityMap) = m1 +composedmap2(m1::Map{T}, m2::IdentityMap{T}) where {T} = m1 +composedmap2(m1::Map{S}, m2::IdentityMap{T}) where {S,T} = convert(Map{T}, m1) + +composedmap2(m1, m2::ConstantMap) = m2 +composedmap2(m1::Map{T}, m2::ConstantMap) where {T} = ConstantMap{T}(mapconstant(m2)) +composedmap1(m1::ConstantMap{T}, m2) where {T} = ConstantMap{T}(m2(mapconstant(m1))) + +composedmap2(m1, m2::ZeroMap) = m2 +composedmap2(m1::Map{T}, m2::ZeroMap{S,U}) where {S,T,U} = ZeroMap{T,U}() +composedmap1(m1::ConstantMap{T}, m2::ZeroMap{S,U}) where {S,T,U} = ZeroMap{T,U}() + +multiply_map1(m1::ZeroMap, m2) = m1 +multiply_map2(m1, m2::ZeroMap) = m2 +multiply_map2(m1, m2::ConstantMap) = _constant_multiply_map(m1, m2) +_constant_multiply_map(m1::ConstantMap{T}, m2::ConstantMap{S}) where {S,T} = + ConstantMap{promote_type(S,T)}(mapconstant(m1)*mapconstant(m2)) +_constant_multiply_map(m1, m2::ConstantMap) = default_multiply_map(m1, m2) + +multiply_map(m1::Function, m2::Function) = t -> m1(t)*m2(t) + +sum_map1(m1::ZeroMap, m2) = m2 +sum_map2(m1, m2::ZeroMap) = m1 +sum_map2(m1, m2::ConstantMap) = _constant_sum_map(m1, m2) +_constant_sum_map(m1::ConstantMap{T}, m2::ConstantMap{S}) where {S,T} = + ConstantMap{promote_type(S,T)}(mapconstant(m1)+mapconstant(m2)) +_constant_sum_map(m1, m2::ConstantMap) = default_sum_map(m1, m2) + +## Affine maps + +composedmap(m1::AbstractAffineMap, m2::AbstractAffineMap) = affine_composition(m1, m2) + +""" +Compute the affine map that represents map2 after map1, that is: +`y = a2*(a1*x+b1)+b2 = a2*a1*x + a2*b1 + b2`. +""" +affine_composition(map1::AbstractAffineMap, map2::AbstractAffineMap) = + AffineMap(affinematrix(map2) * affinematrix(map1), + affinematrix(map2)*affinevector(map1) + affinevector(map2)) + +affine_composition(map1::AffineMap, map2::AffineMap) = + AffineMap(unsafe_matrix(map2) * unsafe_matrix(map1), + unsafe_matrix(map2)*unsafe_vector(map1) + unsafe_vector(map2)) + +affine_composition(map1::LinearMap, map2::LinearMap) = + LinearMap(unsafe_matrix(map2) * unsafe_matrix(map1)) + +affine_composition(map1::LinearMap, map2::AffineMap) = + AffineMap(unsafe_matrix(map2) * unsafe_matrix(map1), unsafe_vector(map2)) + +affine_composition(map1::AffineMap, map2::LinearMap) = + AffineMap(unsafe_matrix(map2) * unsafe_matrix(map1), unsafe_matrix(map2)*unsafe_vector(map1)) + +affine_composition(map1::Translation, map2::Translation) = + Translation(unsafe_vector(map2) + unsafe_vector(map1)) + +affine_composition(map1::Translation, map2::LinearMap) = + AffineMap(unsafe_matrix(map2), unsafe_matrix(map2)*unsafe_vector(map1)) + +affine_composition(map1::LinearMap, map2::Translation) = + AffineMap(unsafe_matrix(map1), unsafe_vector(map2)) + +# The sum of two affine maps is again an affine map +sum_map(map1::AbstractAffineMap, map2::AbstractAffineMap) = + AffineMap(affinematrix(map1)+affinematrix(map2), affinevector(map1)+affinevector(map2)) + + +isequalmap(m1::ProductMap, m2::IdentityMap) = all(map(isidentitymap, components(m1))) +isequalmap(m1::IdentityMap, m2::ProductMap) = isequalmap(m2, m1) + + +""" +Map the interval `[a,b]` to the interval `[c,d]`. + +This function deals with infinite intervals, and the type of the +map returned may depend on the value (finiteness) of the given endpoints. +""" +interval_map(a, b, c, d) = interval_map(promote(a,b,c,d)...) + +function interval_map(a::T, b::T, c::T, d::T) where {T} + FT = float(T) + if isfinite(a) && isfinite(b) && isfinite(c) && isfinite(d) + bounded_interval_map(a, b, c, d) + elseif isfinite(a) && !isfinite(b) && isfinite(c) && !isfinite(d) + # (a,Inf) to (c,Inf) + AffineMap(one(FT), c-a) + elseif isfinite(a) && !isfinite(b) && !isfinite(c) && isfinite(d) + # (a,Inf) to (Inf,d) + AffineMap(-one(FT), d+a) + elseif !isfinite(a) && isfinite(b) && isfinite(c) && !isfinite(d) + # (Inf,b) to (c,Inf) + AffineMap(-one(FT), c+b) + elseif !isfinite(a) && isfinite(b) && !isfinite(c) && isfinite(d) + # (Inf,b) to (Inf,d) + AffineMap(one(FT), d-b) + elseif !isfinite(a) && !isfinite(b) && !isfinite(c) && !isfinite(d) + if (a < 0) && (b > 0) && (c < 0) && (d > 0) + # (-Inf,Inf) to (-Inf,Inf) + StaticIdentityMap{FT}() + elseif (a < 0) && (b > 0) && (c > 0) && (d < 0) + # (-Inf,Inf) to (Inf,-Inf) + LinearMap(-one(FT)) + elseif (a > 0) && (b < 0) && (c < 0) && (d > 0) + # (Inf,-Inf) to (-Inf,Inf) + LinearMap(-one(FT)) + elseif (a > 0) && (b < 0) && (c > 0) && (d < 0) + # (Inf,-Inf) to (Inf,-Inf) + StaticIdentityMap{FT}() + elseif (a > 0) && (b > 0) && (c > 0) && (d > 0) + # (Inf,Inf) to (Inf,Inf) + StaticIdentityMap{FT}() + elseif (a < 0) && (b < 0) && (c < 0) && (d < 0) + # (-Inf,-Inf) to (-Inf,-Inf) + StaticIdentityMap{FT}() + elseif (a < 0) && (b < 0) && (c > 0) && (d > 0) + # (-Inf,-Inf) to (Inf,Inf) + LinearMap(-one(FT)) + elseif (a > 0) && (b > 0) && (c < 0) && (d < 0) + # (Inf,Inf) to (-Inf,-Inf) + LinearMap(-one(FT)) + else + throw(ArgumentError("Requested affine map is unbounded")) + end + else + throw(ArgumentError("Requested affine map is unbounded")) + end +end + +"Like interval_map, but guaranteed to return a scalar affine map." +bounded_interval_map(a, b, c, d) = bounded_interval_map(promote(a,b,c,d)...) +bounded_interval_map(a::T, b::T, c::T, d::T) where {T} = + AffineMap((d-c)/(b-a), c - a*(d-c)/(b-a)) diff --git a/src/types/basic.jl b/src/types/basic.jl new file mode 100644 index 0000000..57c34fc --- /dev/null +++ b/src/types/basic.jl @@ -0,0 +1,160 @@ + +"Supertype of identity maps." +abstract type IdentityMap{T} <: Map{T} end + +IdentityMap(n::Int) = DynamicIdentityMap(n) +IdentityMap() = StaticIdentityMap() +IdentityMap(::Val{N}) where {N} = StaticIdentityMap(Val(N)) + +IdentityMap{T}(n::Int) where {T} = DynamicIdentityMap{T}(n) +IdentityMap{T}(n::Int) where {T<:StaticTypes} = StaticIdentityMap{T}() +IdentityMap{T}(::Val{N}) where {N,T} = StaticIdentityMap{T}(Val(N)) +IdentityMap{T}() where {T} = StaticIdentityMap{T}() + +applymap(map::IdentityMap, x) = x +applymap!(y, map::IdentityMap, x) = y .= x + +inverse(m::IdentityMap) = m +inverse(m::IdentityMap, x) = x + +islinearmap(::IdentityMap) = true +isrealmap(::IdentityMap{T}) where {T} = isrealtype(T) + +isidentitymap(::IdentityMap) = true +isidentitymap(m::Map{T}) where {T} = m == StaticIdentityMap{T}() + +mapsize(m::IdentityMap{T}) where {T<:Number} = () +mapsize(m::IdentityMap{T}) where {T} = (euclideandimension(T),euclideandimension(T)) + +affinematrix(m::IdentityMap) = identitymatrix(m) +affinevector(m::IdentityMap) = zerovector(m) + +jacobian(m::IdentityMap) = ConstantMap(affinematrix(m)) +jacobian(m::IdentityMap, x) = affinematrix(m) + +jacdet(m::IdentityMap, x) = 1 + +determinantmap(m::IdentityMap{T}) where {T} = UnityMap{T,prectype(T)}() + +mapcompose(m1::IdentityMap) = m1 +mapcompose(m1::IdentityMap, maps...) = mapcompose(maps...) +mapcompose2(m1, m2::IdentityMap, maps...) = mapcompose(m1, maps...) + +show(io::IO, m::IdentityMap{T}) where {T} = print(io, "x -> x") +map_object_parentheses(m::IdentityMap) = true + +"The identity map for variables of type `T`." +struct StaticIdentityMap{T} <: IdentityMap{T} +end + +StaticIdentityMap() = StaticIdentityMap{Float64}() +StaticIdentityMap(::Val{N}) where {N} = StaticIdentityMap{SVector{N,Float64}}() + +StaticIdentityMap{T}(n::Int) where {T} = + (@assert n == euclideandimension(T); StaticIdentityMap{T}()) +StaticIdentityMap{T}(::Val{N}) where {N,T} = + (@assert N == euclideandimension(T); StaticIdentityMap{T}()) + +similarmap(m::StaticIdentityMap, ::Type{T}) where {T<:StaticTypes} = StaticIdentityMap{T}() +similarmap(m::StaticIdentityMap, ::Type{T}) where {T} = + DynamicIdentityMap{T}(euclideandimension(T)) + +convert(::Type{StaticIdentityMap{T}}, ::StaticIdentityMap) where {T} = StaticIdentityMap{T}() + +isequalmap(m1::StaticIdentityMap, m2::StaticIdentityMap) = true +map_hash(m::StaticIdentityMap, h::UInt) = hash("StaticIdentityMap", h) + +"Identity map with dynamic size determined by a dimension field." +struct DynamicIdentityMap{T} <: IdentityMap{T} + dimension :: Int +end + +const EuclideanIdentityMap{N,T} = StaticIdentityMap{SVector{N,T}} +const VectorIdentityMap{T} = DynamicIdentityMap{Vector{T}} + +DynamicIdentityMap(dimension::Int) = VectorIdentityMap(dimension) +VectorIdentityMap(dimension::Int) = VectorIdentityMap{Float64}(dimension) + +mapsize(m::DynamicIdentityMap) = (m.dimension, m.dimension) + +similarmap(m::DynamicIdentityMap, ::Type{T}) where {T} = + DynamicIdentityMap{T}(m.dimension) +similarmap(m::DynamicIdentityMap, ::Type{T}) where {T<:StaticTypes} = + StaticIdentityMap{T}() + +isequalmap(m1::DynamicIdentityMap, m2::DynamicIdentityMap) = m1.dimension == m2.dimension +map_hash(m::DynamicIdentityMap, h::UInt) = hashrec("DynamicIdentityMap", m.dimension, h) + + +"The supertype of constant maps from `T` to `U`." +abstract type ConstantMap{T,U} <: TypedMap{T,U} end + +applymap(m::ConstantMap, x) = mapconstant(m) + +isconstantmap(m::AbstractMap) = false +isconstantmap(m::ConstantMap) = true + +isrealmap(m::ConstantMap{T,U}) where {T,U} = + isrealtype(T) && isrealtype(U) && isreal(mapconstant(m)) + +mapsize(m::ConstantMap) = _constant_mapsize(m, mapconstant(m)) +_constant_mapsize(m::ConstantMap{T,U}, c) where {T<:Number,U<:Number} = () +_constant_mapsize(m::ConstantMap{T,U}, c) where {T<:Number,U} = (length(c),) +_constant_mapsize(m::ConstantMap{T,U}, c) where {T,U<:Number} = (1,euclideandimension(T)) +_constant_mapsize(m::ConstantMap{T,U}, c) where {T,U} = (length(c), euclideandimension(T)) + +affinematrix(m::ConstantMap) = zeromatrix(m) +affinevector(m::ConstantMap) = mapconstant(m) + +jacobian(m::ConstantMap{T}) where {T} = ConstantMap{T}(affinematrix(m)) +jacobian(m::ConstantMap, x) = affinematrix(m) + +jacdet(::ConstantMap, x) = 0 + +determinantmap(m::ConstantMap{T}) where {T} = ConstantMap{T}(det(mapconstant(m))) +absmap(m::ConstantMap{T}) where {T} = ConstantMap{T}(abs(mapconstant(m))) + +diffvolume(m::ConstantMap{T,U}) where {T,U} = ZeroMap{T,U}() + +isequalmap(m1::ConstantMap, m2::ConstantMap) = mapconstant(m1)==mapconstant(m2) +map_hash(m::ConstantMap, h::UInt) = hashrec("ConstantMap", mapconstant(m), h) + +similarmap(m::ConstantMap, ::Type{T}) where {T} = ConstantMap{T}(mapconstant(m)) +similarmap(m::ConstantMap, ::Type{T}, ::Type{U}) where {T,U} = ConstantMap{T,U}(m.c) + +ConstantMap() = ConstantMap{Float64}() +ConstantMap(c) = FixedConstantMap(c) +ConstantMap{T}() where {T} = UnityMap{T}() +ConstantMap{T}(c) where {T} = FixedConstantMap{T}(c) +ConstantMap{T,U}() where {T,U} = UnityMap{T,U}() +ConstantMap{T,U}(c) where {T,U} = FixedConstantMap{T,U}(c) + +show(io::IO, m::ConstantMap{T}) where {T} = print(io, "x -> $(mapconstant(m))") +map_object_parentheses(m::ConstantMap) = true + + +"The zero map `f(x) = 0`." +struct ZeroMap{T,U} <: ConstantMap{T,U} +end +ZeroMap{T}() where {T} = ZeroMap{T,T}() +mapconstant(m::ZeroMap{T,U}) where {T,U} = zero(U) +similarmap(m::ZeroMap{S,U}, ::Type{T}) where {T,S,U} = ZeroMap{T,U}() +similarmap(m::ZeroMap, ::Type{T}, ::Type{U}) where {T,U} = ZeroMap{T,U}() + + +"The unity map `f(x) = 1`." +struct UnityMap{T,U} <: ConstantMap{T,U} +end +UnityMap{T}() where {T} = UnityMap{T,real(numtype(T))}() +mapconstant(m::UnityMap{T,U}) where {T,U} = one(U) +similarmap(m::UnityMap{S,U}, ::Type{T}) where {T,S,U} = UnityMap{T,U}() +similarmap(m::UnityMap, ::Type{T}, ::Type{U}) where {T,U} = UnityMap{T,U}() + + +"The constant map `f(x) = c`." +struct FixedConstantMap{T,U} <: ConstantMap{T,U} + c :: U +end +FixedConstantMap{T}(c::U) where {T,U} = FixedConstantMap{T,U}(c) +FixedConstantMap(c::T) where {T} = FixedConstantMap{T}(c) +mapconstant(m::FixedConstantMap) = m.c diff --git a/src/types/coordinates.jl b/src/types/coordinates.jl new file mode 100644 index 0000000..1e4d93e --- /dev/null +++ b/src/types/coordinates.jl @@ -0,0 +1,138 @@ + +# This file contains a few specific maps, included mainly to test +# functionality on non-trivial maps + +""" +A Cartesion to Polar map. First dimension is interpreted as radial distance, +second as an angle. +The unit circle is mapped to the square `[-1,1]x[-1,1]`. +""" +struct CartToPolarMap{T} <: Map{SVector{2,T}} +end + +CartToPolarMap() = CartToPolarMap{Float64}() + +mapsize(m::CartToPolarMap) = (2,2) + +applymap(map::CartToPolarMap{T}, x) where {T} = + SVector{2,T}(sqrt(x[1]^2+x[2]^2)*2-1, atan(x[2],x[1])/pi) + +function jacobian(m::CartToPolarMap{T}, x) where {T} + d = sqrt(x[1]^2+x[2]^2) + r = 1/(1+(x[2]/x[1])^2) + SMatrix{2,2,T}(2*x[1]/d, 1/pi*r*(-x[2]/x[1]^2), 2*x[2]/d, 1/pi*r*1/x[1]) +end + +inverse(m::CartToPolarMap{T}) where {T} = PolarToCartMap{T}() +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, +radius to the first. +The square `[-1,1]x[-1,1]` is mapped to the unit circle. +""" +struct PolarToCartMap{T} <: Map{SVector{2,T}} +end + +PolarToCartMap() = PolarToCartMap{Float64}() + +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])) + +jacobian(m::PolarToCartMap{T}, x) where {T} = + SMatrix{2,2,T}(cos(pi*x[2])/2, sin(pi*x[2])/2, -pi*(x[1]+1)/2*sin(pi*x[2]), pi*(x[1]+1)/2*cos(pi*x[2])) + +inverse(m::PolarToCartMap{T}) where {T} = CartToPolarMap{T}() +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`. +""" +struct UnitCircleMap{T} <: Map{T} end + +UnitCircleMap() = UnitCircleMap{Float64}() + +mapsize(m::UnitCircleMap) = (2,) + +applymap(m::UnitCircleMap{T}, t) where {T} = SVector(cos(2*T(pi)*t), sin(2*T(pi)*t)) + +function jacobian(m::UnitCircleMap{T}, t) where {T} + a = 2*T(pi) + SVector(-a*sin(a*t), a*cos(a*t)) +end + +# we know that the differential volume is 2*pi +diffvolume(m::UnitCircleMap{T}) where {T} = ConstantMap{T}(2*T(pi)) + +""" +`AngleMap` is a left inverse of `UnitCircleMap`. A 2D vector `x` is projected onto +the intersection point with the unit circle of the line connecting `x` to the +origin. The angle of this point, scaled to the interval `[0,1)`, is the result. +""" +struct AngleMap{T} <: Map{SVector{2,T}} +end + +AngleMap() = AngleMap{Float64}() + +function applymap(m::AngleMap{T}, x) where {T} + twopi = 2*convert(T, pi) + θ = atan(x[2],x[1]) + if θ < 0 + # atan2 returns an angle in (-π,π], convert to [0,2π) using periodicity. + θ += twopi + end + # And divide by 2π to scale to [0,1) + θ / twopi +end + +mapsize(m::AngleMap) = (1,2) + +function jacobian(m::AngleMap{T}, t) where {T} + x = t[1]; y = t[2] + twopi = 2*convert(T, pi) + v = SVector(one(T) / (1+(y/x)^2) * (-y/x^2) * 1/twopi, + one(T) / (1+(y/x)^2) * one(T)/x * 1/twopi) + transpose(v) +end + + +leftinverse(m::UnitCircleMap{T}) where {T} = AngleMap{T}() +leftinverse(m::UnitCircleMap, x) = leftinverse(m)(x) +rightinverse(m::AngleMap{T}) where {T} = UnitCircleMap{T}() +rightinverse(m::AngleMap, x) = rightinverse(m)(x) + + +""" +The map `r*[cos(2πt), sin(2πt)]` from `[0,1]^2` to the unit disk in `ℝ^2`. +""" +struct UnitDiskMap{T} <: Map{SVector{2,T}} end + +UnitDiskMap() = UnitDiskMap{Float64}() + +mapsize(m::UnitDiskMap) = (2,2) + +applymap(m::UnitDiskMap{T}, x) where {T} = + SVector(x[1]*cos(2*T(pi)*x[2]), x[1]*sin(2*T(pi)*x[2])) + +function jacobian(m::UnitDiskMap{T}, x) where {T} + a = 2*T(pi) + SMatrix{2,2}(cos(a*x[2]), sin(a*x[2]), -a*x[1]*sin(a*x[2]), a*x[1]*cos(a*x[2])) +end + +# we know that the differential volume is 2*r*pi +diffvolume(m::UnitDiskMap, x) = 2*x[1]*pi diff --git a/src/util/common.jl b/src/util/common.jl new file mode 100644 index 0000000..0c962d6 --- /dev/null +++ b/src/util/common.jl @@ -0,0 +1,195 @@ + +isrealtype(::Type{<:Real}) = true +isrealtype(::Type{<:Complex}) = false +isrealtype(::Type{T}) where {T} = isrealtype(eltype(T)) + +const StaticTypes = Union{Number,<:StaticVector{N} where N,<:NTuple{N,Any} where N} + +"Apply the `hash` function recursively to the given arguments." +hashrec(x) = hash(x) +hashrec(x, args...) = hash(x, hashrec(args...)) + +# Workaround for #88, manually compute the hash of an array using all its elements +hashrec() = zero(UInt) +function hashrec(A::AbstractArray, args...) + h = hash(size(A)) + for x in A + h = hash(x, h) + end + hash(h, hashrec(args...)) +end + +"What is the euclidean dimension of the given type (if applicable)?" +euclideandimension(::Type{T}) where {T <: Number} = 1 +euclideandimension(::Type{T}) where {N,T <: StaticVector{N}} = N +euclideandimension(::Type{T}) where {N,T <: NTuple{N,Any}} = N +euclideandimension(::Type{T}) where {T} = + throw(ArgumentError("Don't know the euclidean dimension of $(T).")) + +################# +# Element type +################# + +""" + convert_eltype(T, x) + +Convert `x` such that its `eltype` equals `T`. +""" +convert_eltype(::Type{T}, d::AbstractArray) where {T} = convert(AbstractArray{T}, d) +convert_eltype(::Type{T}, d::AbstractRange) where {T} = map(T, d) +convert_eltype(::Type{T}, d::Set) where {T} = convert(Set{T}, d) +convert_eltype(::Type{T}, d::Number) where {T} = convert(T, d) + +"Are the given element types promotable to a concrete supertype?" +promotable_eltypes(types...) = isconcretetype(promote_type(types...)) +promotable_eltypes(::Type{S}, ::Type{T}) where {S<:AbstractVector,T<:AbstractVector} = + promotable_eltypes(eltype(S), eltype(T)) + +################# +# Precision type +################# + +""" + prectype(x[, ...]) + +The floating point precision type associated with the argument(s). +""" +prectype(x) = prectype(typeof(x)) +prectype(::Type{T}) where T = Any # fallback +prectype(::Type{T}) where {T<:AbstractFloat} = T +prectype(::Type{T}) where {T<:Number} = prectype(float(T)) +prectype(::Type{<:AbstractArray{T}}) where T = prectype(T) +prectype(::Type{<:Complex{T}}) where {T} = prectype(T) +prectype(::Type{NTuple{N,T}}) where {N,T} = prectype(T) +prectype(::Type{Tuple{A}}) where {A} = prectype(A) +prectype(::Type{Tuple{A,B}}) where {A,B} = prectype(A,B) +prectype(::Type{Tuple{A,B,C}}) where {A,B,C} = prectype(A,B,C) +@generated function prectype(T::Type{<:Tuple}) + quote $(promote_type(map(prectype, T.parameters[1].parameters)...)) end +end + +prectype(a...) = promote_type(map(prectype, a)...) + +""" + convert_prectype(T, x) + +Convert `x` such that its `prectype` equals `T`. +""" +convert_prectype(::Type{T}, x) where {T} = + convert_eltype(to_prectype(T, eltype(x)), x) + +""" + to_prectype(U, T) + +Return the type to which `T` can be converted, such that the `prectype` becomes `U`. +""" +to_prectype(::Type{U}, ::Type{T}) where {T,U} = throw(ArgumentError("Don't know how to convert the prectype of $(T) to $(U).")) +to_prectype(::Type{U}, ::Type{T}) where {T <: Real,U <: Real} = U +to_prectype(::Type{U}, ::Type{Complex{T}}) where {T <: Real,U <: Real} = Complex{U} +to_prectype(::Type{U}, ::Type{Vector{T}}) where {T,U} = Vector{to_prectype(U,T)} +to_prectype(::Type{U}, ::Type{<:StaticVector{N,T}}) where {N,T,U} = SVector{N,to_prectype(U,T)} +to_prectype(::Type{U}, ::Type{MVector{N,T}}) where {N,T,U} = MVector{N,to_prectype(U,T)} +to_prectype(::Type{U}, ::Type{<:StaticMatrix{M,N,T}}) where {M,N,T,U} = SMatrix{M,N,to_prectype(U,T)} +to_prectype(::Type{U}, ::Type{MMatrix{M,N,T}}) where {M,N,T,U} = MMatrix{M,N,to_prectype(U,T)} + +""" + promote_prectype(a, b[, ...]) + +Promote the precision types of the arguments to a joined supertype. +""" +promote_prectype(a) = a +promote_prectype(a, b) = _promote_prectype(prectype(a,b), a, b) +promote_prectype(a, b, c...) = _promote_prectype(prectype(a,b,c...), a, b, c...) +_promote_prectype(U, a) = convert_prectype(U, a) +_promote_prectype(U, a, b) = convert_prectype(U, a), convert_prectype(U, b) +_promote_prectype(U, a, b, c...) = + (convert_prectype(U, a), convert_prectype(U, b), _promote_prectype(U, c...)...) + +################# +# Numeric type +################# + +"The numeric element type of `x` in a Euclidean space." +numtype(x) = numtype(typeof(x)) +numtype(::Type{T}) where T = Any +numtype(::Type{T}) where {T<:Number} = T +numtype(::Type{<:AbstractArray{T}}) where T = T +numtype(::Type{NTuple{N,T}}) where {N,T} = T +numtype(::Type{Tuple{A,B}}) where {A,B} = promote_type(numtype(A), numtype(B)) +numtype(::Type{Tuple{A,B,C}}) where {A,B,C} = promote_type(numtype(A), numtype(B), numtype(C)) +numtype(::Type{Tuple{A,B,C,D}}) where {A,B,C,D} = promote_type(numtype(A), numtype(B), numtype(C), numtype(D)) +@generated function numtype(T::Type{<:Tuple}) + quote $(promote_type(T.parameters[1].parameters...)) end +end + +numtype(a...) = promote_type(map(numtype, a)...) + +""" + convert_numtype(T, x) + +Convert `x` such that its `numtype` equals `T`. +""" +convert_numtype(::Type{T}, x) where {T} = + convert_eltype(to_numtype(T, eltype(x)), x) + +""" + to_numtype(U, T) + +Return the type to which `T` can be converted, such that the `numtype` becomes `U`. +""" +to_numtype(::Type{U}, ::Type{T}) where {T,U} = throw(ArgumentError("Don't know how to convert the numtype of $(T) to $(U).")) +to_numtype(::Type{U}, ::Type{T}) where {T <: Number,U <: Number} = U +to_numtype(::Type{U}, ::Type{Vector{T}}) where {T,U} = Vector{to_numtype(U,T)} +to_numtype(::Type{U}, ::Type{<:StaticVector{N,T}}) where {N,T,U} = SVector{N,to_numtype(U,T)} +to_numtype(::Type{U}, ::Type{MVector{N,T}}) where {N,T,U} = MVector{N,to_numtype(U,T)} +to_numtype(::Type{U}, ::Type{<:StaticMatrix{M,N,T}}) where {M,N,T,U} = SMatrix{M,N,to_numtype(U,T)} +to_numtype(::Type{U}, ::Type{MMatrix{M,N,T}}) where {M,N,T,U} = MMatrix{M,N,to_numtype(U,T)} + +""" + promote_numtype(a, b[, ...]) + +Promote the numeric types of the arguments to a joined supertype. +""" +promote_numtype(a) = a +promote_numtype(a, b) = _promote_numtype(numtype(a,b), a, b) +promote_numtype(a, b, c...) = _promote_numtype(numtype(a,b,c...), a, b, c...) +_promote_numtype(U, a) = convert_numtype(U, a) +_promote_numtype(U, a, b) = convert_numtype(U, a), convert_numtype(U, b) +_promote_numtype(U, a, b, c...) = + (convert_numtype(U, a), convert_numtype(U, b), _promote_numtype(U, c...)...) + +## Conversion from nested vectors to flat vectors and back + +""" +Convert a vector from a cartesian format to a nested tuple according to the +given dimensions. + +For example: +`convert_fromcartesian([1,2,3,4,5], Val{(2,2,1)}()) -> ([1,2],[3,4],5)` +""" +@generated function convert_fromcartesian(x::AbstractVector, ::Val{DIM}) where {DIM} + dimsum = [0; cumsum([d for d in DIM])] + E = Expr(:tuple, [ (dimsum[i+1]-dimsum[i] > 1 ? Expr(:call, :SVector, [:(x[$j]) for j = dimsum[i]+1:dimsum[i+1]]...) : :(x[$(dimsum[i+1])])) for i in 1:length(DIM)]...) + return quote $(E) end +end + +"The inverse function of `convert_fromcartesian`." +@generated function convert_tocartesian(x, ::Val{DIM}) where {DIM} + dimsum = [0; cumsum([d for d in DIM])] + E = vcat([[:(x[$i][$j]) for j in 1:DIM[i]] for i in 1:length(DIM)]...) + quote SVector($(E...)) end +end + +# we use matrix_pinv rather than pinv to preserve static matrices +matrix_pinv(A) = pinv(A) +# matrix_pinv(A::SMatrix{M,N}) where {M,N} = SMatrix{N,M}(pinv(A)) +matrix_pinv(A::StaticVector{N}) where {N} = convert(Transpose{Float64, SVector{N,Float64}}, pinv(A)) + +## Composite objects + +"Factors of a product-like composite object (equivalent to `components(d)`)." +function factors end +"The number of factors of a product-like composite object." +nfactors(d) = length(factors(d)) +"Factor `I...` of a product-like composite object." +factor(d, I...) = getindex(factors(d), I...) diff --git a/test/aqua.jl b/test/aqua.jl new file mode 100644 index 0000000..76d32d6 --- /dev/null +++ b/test/aqua.jl @@ -0,0 +1,14 @@ +using Aqua + +# Aqua.test_all(FunctionMaps) + +Aqua.test_ambiguities(FunctionMaps) +# The unbound tests annoyingly flag several valid uses of NTuple{N,T} +# Aqua.test_unbound_args(FunctionMaps) +Aqua.test_undefined_exports(FunctionMaps) +Aqua.test_project_extras(FunctionMaps) +Aqua.test_stale_deps(FunctionMaps) +Aqua.test_deps_compat(FunctionMaps) +# Piracies to be removed in next breaking release +Aqua.test_piracies(FunctionMaps) +Aqua.test_persistent_tasks(FunctionMaps) diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..5253792 --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,20 @@ +module FunctionMapsTests + +using Test, LinearAlgebra, StaticArrays +using CompositeTypes, CompositeTypes.Indexing + +include("using_fmaps.jl") + +include("aqua.jl") + +include("test_common.jl") +include("test_interface.jl") + +include("test_generic.jl") +include("test_basic.jl") +include("test_affine.jl") +include("test_product.jl") +include("test_maps.jl") +include("test_arithmetics.jl") + +end diff --git a/test/test_affine.jl b/test/test_affine.jl new file mode 100644 index 0000000..10a9d4d --- /dev/null +++ b/test/test_affine.jl @@ -0,0 +1,237 @@ +function test_affine_maps(T) + A = rand(T,2,2) + @test FunctionMaps.to_matrix(Vector{T}, A) == A + @test FunctionMaps.to_matrix(T, 2) == 2 + @test FunctionMaps.to_matrix(SVector{2,T}, 2) == SMatrix{2,2}(2,0,0,2) + @test FunctionMaps.to_matrix(SVector{2,T}, LinearAlgebra.I) == SMatrix{2,2}(1,0,0,1) + @test FunctionMaps.to_matrix(Vector{T}, 2) == UniformScaling(2) + @test FunctionMaps.to_matrix(Vector{T}, LinearAlgebra.I) == LinearAlgebra.I + # test fallback with nonsensical call + @test FunctionMaps.to_matrix(Tuple{Int}, 2) == 2 + + @test FunctionMaps.to_matrix(T, A, 2) == A + @test FunctionMaps.to_matrix(T, 2, 3) == 2 + @test FunctionMaps.to_matrix(T, UniformScaling(2), 3) == 2 + @test FunctionMaps.to_matrix(T, LinearAlgebra.I, zero(T)) isa T + @test FunctionMaps.to_matrix(SVector{2,T}, 2, SVector(1,1)) == SMatrix{2,2}(2,0,0,2) + @test FunctionMaps.to_matrix(Vector{T}, 2, [1,2]) == [2 0 ; 0 2] + + @test FunctionMaps.to_vector(T, 2) == 0 + @test FunctionMaps.to_vector(SVector{2,T}, 2) == SVector(0,0) + @test FunctionMaps.to_vector(Vector{T}, A) == [0,0] + @test FunctionMaps.to_vector(T, 2, 3) == 3 + + if T != BigFloat # BigFloat's make pinv fail for StaticArrays + @test FunctionMaps.matrix_pinv(SMatrix{2,2}(rand(T),rand(T),rand(T),rand(T))) isa SMatrix{2,2} + @test FunctionMaps.matrix_pinv(SVector(rand(T),rand(T))) isa Transpose{T,SVector{2,T}} + end + + test_linearmap(T) + test_translation(T) + test_affinemap(T) + test_abstractarraymap(T) +end + + + +function test_linearmap(T) + m1 = LinearMap(2one(T)) + @test m1 isa LinearMap{T} + @test domaintype(m1) == T + @test islinearmap(m1) + @test isaffinemap(m1) + @test m1(3) == 6 + @test m1(3one(T)) == 6 + @test m1(3one(widen(T))) isa widen(T) + @test m1(SVector(1,2)) == SVector(2,4) + @test m1([1.0,2.0]) == [2.0,4.0] + mw1 = convert(Map{widen(T)}, m1) + @test mw1 isa Map{widen(T)} + @test mw1.A isa widen(T) + @test jacobian(m1) isa ConstantMap{T} + @test jacobian(m1, 1) == 2 + @test LinearMap(one(T)) == StaticIdentityMap{T}() + + m2 = LinearMap(2) + @test domaintype(m2) == Int + @test m2 === ScalarLinearMap(2) + @test m2(one(T)) isa T + @test m2 == convert(Map, 2) + @test jacobian(m2, 1) == 2 + @test jacobian(m2) isa ConstantMap{Int} + @test jacobian(m2, 1) == 2 + @test LinearMap(1) == StaticIdentityMap{T}() + + m3 = LinearMap(SMatrix{2,2}(one(T), 2one(T), 3one(T), 4one(T))) + @test m3 isa LinearMap{SVector{2,T}} + @test m3 === StaticLinearMap(m3.A) + @test m3(SVector(1,2)) == SVector(7, 10) + @test m3(SVector{2,T}(1,2)) == SVector{2,T}(7, 10) + @test m3 ∘ m3 isa LinearMap + @test LinearMap(SA[1 0; 0 1]) == StaticIdentityMap{domaintype(m3)}() + + A = rand(T,2,2) + m4 = LinearMap(A) + @test m4 isa LinearMap{Vector{T}} + @test m4([1,2]) == A * [1,2] + y = zeros(T,2) + @test (FunctionMaps.applymap!(y, m4, [1,2]); y == A * [1,2]) + @test jacobian(m4, [1,2]) == A + + m5 = LinearMap{Vector{T}}(UniformScaling(2*one(T))) + @test m5 isa LinearMap{Vector{T}} + @test m5([1,2]) == 2 * [1,2] + @test jacobian(m5, [1,2]) == UniformScaling(2) + + m6 = LinearMap{SVector{2,T}}(2) + @test m6 isa GenericLinearMap{SVector{2,T},T} + @test m6(SA[one(T),one(T)]) == [2,2] + @test leftinverse(m6) isa LinearMap{SVector{2,T}} + @test leftinverse(m6)([2,2]) == [1,1] + + # Test construction and conversion + @test LinearMap{T}(1) isa ScalarLinearMap{T} + @test LinearMap{Vector{T}}(rand(Int,5,5)) isa VectorLinearMap{T} + @test LinearMap{SVector{2,T}}(SMatrix{2,2}(1, 2, 3, 4)) isa StaticLinearMap{T,2,2,4} + + @test GenericLinearMap(one(T)) isa GenericLinearMap{T,T} + @test GenericLinearMap(SMatrix{2,2,T}(1,1,1,1)) isa GenericLinearMap{SVector{2,T},SMatrix{2,2,T,4}} + + @test VectorLinearMap(rand(T,3,3)) isa VectorLinearMap{T} + @test StaticLinearMap(SMatrix{2,2,T}(1,2,3,4)) isa StaticLinearMap{T,2,2,4} + + @test convert(Map{SVector{2,T}}, LinearMap(rand(T,2,2))) isa StaticLinearMap{T,2,2,4} + @test convert(Map{T}, 1) isa LinearMap{T} +end + + + +function test_translation(T) + v = randvec(T,3) + m = Translation(v) + @test !islinearmap(m) + @test isaffinemap(m) + @test inverse(inverse(m)) == m + @test jacobian(m) isa ConstantMap + @test affinevector(jacobian(m)) == [1 0 0; 0 1 0; 0 0 1] + @test jacdet(m) isa UnityMap{SVector{3,T},T} + + @test mapsize(Translation(one(T))) == () + + # Test construction and conversion + @test Translation(one(T)) isa ScalarTranslation{T} + @test Translation{T}(1) isa ScalarTranslation{T} + @test Translation(SVector{2,T}(1,2)) isa StaticTranslation{T,2} + @test Translation{SVector{2,T}}(rand(T,2)) isa StaticTranslation{T,2} + @test Translation(rand(Int,2)) isa VectorTranslation{Int} + @test Translation{Vector{T}}(rand(Int,2)) isa VectorTranslation{T} + @test Translation(1:10) isa GenericTranslation{Vector{Int}} + @test Translation{Vector{T}}(1:10) isa GenericTranslation{Vector{T}} + + @test StaticTranslation(MVector(one(T),2)) isa StaticTranslation{T,2} + @test StaticTranslation{T}(SVector(1,2)) isa StaticTranslation{T,2} + + @test VectorTranslation(1:10) isa VectorTranslation{Int} + @test VectorTranslation{T}(1:10) isa VectorTranslation{T} + + @test GenericTranslation(1:10) isa GenericTranslation{Vector{Int}} + @test GenericTranslation{Vector{T}}(1:10) isa GenericTranslation{Vector{T}} + + @test GenericTranslation(one(T)) isa GenericTranslation{T,T} + @test GenericTranslation{T}(1) isa GenericTranslation{T,T} + @test_throws InexactError GenericTranslation{Int}(one(T)/2) + + @test convert(Map{SVector{2,T}}, Translation(rand(T,2))) isa StaticTranslation{T,2} +end + + +function test_affinemap(T) + m1 = AffineMap(T(2), T(3)) + @test m1 isa ScalarAffineMap{T} + @test !islinearmap(m1) + @test isaffinemap(m1) + @test m1(2) == 7 + + @test convert(ScalarAffineMap{BigFloat}, m1) == m1 + + @test m1 ∘ m1 isa AffineMap + @test (m1 ∘ m1)(2) == 2*(2*2+3)+3 + + m2 = AffineMap(T(2), SVector{2,T}(1,2)) + @test m2 isa GenericAffineMap{SVector{2,T}} + @test m2(SVector(1,2)) == SVector(3,6) + @test mapsize(m2) == (2,2) + + m3 = AffineMap(UniformScaling(2*one(T)), [one(T),2*one(T)]) + @test m3 isa AffineMap{Vector{T}} + @test mapsize(m3) == (2,2) + @test m3([1,2]) == 2 * [1,2] + [1,2] + y = zeros(T,2) + @test (FunctionMaps.applymap!(y, m3, [1,2]); y == m3([1,2])) + @test jacobian(m3, [1,2]) == [2 0; 0 2] + @test jacdet(m3, [1,2]) == 4 + + @test AffineMap(LinearAlgebra.I, 2.0) isa ScalarAffineMap + m4 = GenericAffineMap(LinearAlgebra.I, 2.0) + @test m4.A isa UniformScaling{Float64} + @test jacdet(m4, 3.0) == 1.0 + + @test mapsize(AffineMap(rand(3),rand(3))) == (3,) + @test mapsize(AffineMap(LinearAlgebra.I,2)) == () + + # Test construction and conversion + @test AffineMap(1, 2*one(T)) isa ScalarAffineMap{T} + @test AffineMap{BigFloat}(1, 2.0) isa ScalarAffineMap{BigFloat} + @test AffineMap(rand(T,2),rand(T,2)) isa GenericAffineMap{T} + @test AffineMap(rand(T,2,2),rand(T,2)) isa VectorAffineMap{T} + @test AffineMap{Vector{T}}(rand(Int,2,2),rand(Int,2)) isa VectorAffineMap{T} + @test AffineMap{SVector{2,T}}(rand(T,2,2),rand(T,2)) isa StaticAffineMap{T,2,2,4} + @test AffineMap{SVector{2,T}}(rand(Int,2,2),rand(T,2)) isa StaticAffineMap{T,2,2,4} + @test AffineMap{SVector{2,T}}(rand(Int,2,2),rand(Int,2)) isa StaticAffineMap{T,2,2,4} + @test AffineMap{SVector{2,T}}(1,rand(Int,2)) isa GenericAffineMap{SVector{2,T},T,SVector{2,T}} + @test AffineMap{SVector{2,T}}(SMatrix{3,2}(3,2,1,4,5,6),SVector(1,2,3)) isa StaticAffineMap{T,2,3,6} + + @test GenericAffineMap(SMatrix{2,2,T}(1,2,3,4), [1,2]) isa GenericAffineMap{SVector{2,T},SMatrix{2,2,T,4}} + @test GenericAffineMap(SMatrix{1,1,Int}(1),SVector{1,T}(1)) isa GenericAffineMap{SVector{1,T}} + @test GenericAffineMap(rand(T,2,2),rand(Int,2)) isa GenericAffineMap{Vector{T},Matrix{T},Vector{T}} + @test GenericAffineMap{Vector{T}}(1, [1,2]) isa GenericAffineMap{Vector{T},T,Vector{Int}} + @test GenericAffineMap{Vector{T}}(rand(Int,2,2), rand(Int,2)) isa GenericAffineMap{Vector{T},Matrix{T},Vector{T}} + @test convert(GenericAffineMap{Vector{T}}, GenericAffineMap(rand(Int,2,2), rand(Int,2))) isa GenericAffineMap{Vector{T}} + @test GenericAffineMap(rand(T,2),rand(T,2)) isa GenericAffineMap{T} + @test GenericAffineMap{T}(rand(Int,2),rand(Int,2)) isa GenericAffineMap{T} + + @test VectorAffineMap(rand(T,2,2),rand(T,2)) isa VectorAffineMap{T} + + @test StaticAffineMap(SMatrix{1,1,T}(1),SVector{1,T}(1)) isa StaticAffineMap{T} + @test StaticAffineMap(SMatrix{1,1,T}(1),SVector{1,Int}(1)) isa StaticAffineMap{T} + + @test StaticAffineMap(rand(T,2,2),SVector{2,T}(1,1)) isa StaticAffineMap{T} + @test_throws AssertionError StaticAffineMap(rand(T,3,2),SVector{2,T}(1,1)) + @test StaticAffineMap(SMatrix{2,2,T}(1,2,3,4),[1,1]) isa StaticAffineMap{T,2} + @test StaticAffineMap(SMatrix{3,2,T}(1,2,3,4,5,6),SVector{3,T}(1,2,3)) isa StaticAffineMap{T,2,3,6} + + @test convert(Map{SVector{2,T}}, AffineMap(rand(2,2),rand(2))) isa StaticAffineMap{T,2,2,4} +end + +function test_abstractarraymap(T) + A = rand(T, 3, 3) + @test FunctionMaps.MapStyle(A) isa FunctionMaps.IsMap + @test FunctionMaps.checkmap(A) == A + @test Map(A) isa LinearMap + @test Map(A) isa VectorLinearMap + @test Map(SA[1 2; 3 4]) isa StaticLinearMap + @test Map(Diagonal(rand(3))) isa GenericLinearMap + @test mapsize(A) == size(A) + @test applymap(A, 1:3) == A*(1:3) + @test islinearmap(A) + @test isaffinemap(A) + @test affinematrix(A) == A + @test affinevector(A) == [0,0,0] + @test jacobian(A) isa ConstantMap + @test jacobian(A, 1:3) == A + @test inverse(A) ≈ inv(A) + @test inverse(A, 1:3) ≈ A \ (1:3) + @test canonicalmap(A) isa LinearMap + @test FunctionMaps.equalmap(A) isa LinearMap + @test isequalmap(A, LinearMap(A)) +end diff --git a/test/test_arithmetics.jl b/test/test_arithmetics.jl new file mode 100644 index 0000000..d1503ac --- /dev/null +++ b/test/test_arithmetics.jl @@ -0,0 +1,22 @@ +@testset "map interval" begin + @test interval_map(1.0, Inf, 2.0, Inf) == AffineMap(1.0, 1.0) + @test interval_map(1.0, Inf, Inf, 2.0) == AffineMap(-1.0, 3.0) + @test interval_map(Inf, 1.0, 2.0, Inf) == AffineMap(-1.0, 3.0) + @test interval_map(Inf, 1.0, Inf, 2.0) == AffineMap(1.0, 1.0) + @test interval_map(-Inf, Inf, -Inf, Inf) == IdentityMap() + @test interval_map(-Inf, Inf, Inf, -Inf) == LinearMap(-1) + @test interval_map(Inf, -Inf, -Inf, Inf) == LinearMap(-1) + @test interval_map(Inf, -Inf, Inf, -Inf) == IdentityMap() + @test interval_map(Inf, Inf, Inf, Inf) == IdentityMap() + @test interval_map(-Inf, -Inf, -Inf, -Inf) == IdentityMap() + @test interval_map(Inf, Inf, -Inf, -Inf) == LinearMap(-1) + @test interval_map(-Inf, -Inf, Inf, Inf) == LinearMap(-1) + @test_throws ArgumentError interval_map(-Inf, Inf, -Inf, -Inf) + @test_throws ArgumentError interval_map(-1, 1, -1, Inf) + @test_throws ArgumentError interval_map(-1, 1, -1, Inf) + + a = 1.0; b = 2.0; c = 3.0; d = 5; + @test interval_map(a, b, c, d)(a) == c + @test interval_map(a, b, c, d)(b) == d + @test interval_map(a, b, c, d) == bounded_interval_map(a, b, c, d) +end diff --git a/test/test_basic.jl b/test/test_basic.jl new file mode 100644 index 0000000..a24c35f --- /dev/null +++ b/test/test_basic.jl @@ -0,0 +1,114 @@ +function test_isomorphisms(T) + m1 = FunctionMaps.VectorToNumber{T}() + @test m1(SA[1.0]) == 1.0 + m1b = FunctionMaps.NumberToVector{T}() + @test inverse(m1) == m1b + @test inverse(m1, 0.4) == m1b(0.4) + @test inverse(m1b) == m1 + @test inverse(m1b, SA[0.4]) == m1(SA[0.4]) + @test m1b(1.0) == SA[1.0] + @test mapsize(m1) == (1,1) + @test mapsize(m1b) == (1,) + @test mapsize(m1 ∘ inverse(m1)) == () + @test mapsize(inverse(m1) ∘ m1) == (1,1) + @test jacobian(m1 ∘ inverse(m1), 1) == one(T) + @test jacobian(inverse(m1) ∘ m1, [1]) == ones(T,1,1) + @test jacobian(m1) isa ConstantMap + @test jacobian(m1b) isa ConstantMap + + m2 = FunctionMaps.VectorToComplex{T}() + @test m2(SA[one(T), one(T)]) == 1 + im + m2b = FunctionMaps.ComplexToVector{T}() + @test inverse(m2) == m2b + @test inverse(m2, zero(T)+0im) == m2b(zero(T)+0im) + @test inverse(m2b) == m2 + @test inverse(m2b, SA[one(T),one(T)]) == one(T)+one(T)*im + @test m2b(one(T)+one(T)*im) == SA[one(T),one(T)] + @test mapsize(m2) == (1,2) + @test mapsize(m2b) == (2,) + @test mapsize(m2 ∘ m2b) == () + @test mapsize(m2b ∘ m2) == (2,2) + + m3 = FunctionMaps.VectorToTuple{2,T}() + m3b = FunctionMaps.TupleToVector{2,T}() + @test m3(SA[one(T), one(T)]) == (one(T),one(T)) + @test inverse(m3b, SA[one(T), one(T)]) == (one(T),one(T)) + @test inverse(m3) == m3b + @test inverse(m3b) == m3 + @test m3b( (one(T),one(T)) ) == SA[one(T),one(T)] + @test inverse(m3, (one(T),one(T)) ) == SA[one(T),one(T)] + + m4 = FunctionMaps.NestedToFlat{3,T,Tuple{Tuple{T,T},T},(2,1)}() + m4b = FunctionMaps.FlatToNested{3,T,Tuple{Tuple{T,T},T},(2,1)}() + x4 = ([T(1),T(2)],T(3)) + y4 = T[1,2,3] + @test m4(x4) == y4 + @test m4b(y4) == x4 + @test inverse(m4) == m4b + @test inverse(m4b) == m4 + @test inverse(m4, y4) == x4 + @test inverse(m4b, x4) == y4 +end + +function test_scaling_maps(T) + test_generic_map(LinearMap(T(2))) + test_generic_map(LinearMap(T(2), T(3))) + test_generic_map(LinearMap(T(2), T(3), T(4))) + test_generic_map(LinearMap(T(2), T(3), T(4), T(5))) +end + +function test_identity_map(T) + i1 = StaticIdentityMap{T}() + i2 = StaticIdentityMap{SVector{2,T}}() + test_generic_map(i1) + test_generic_map(i2) + @test i1 == i2 + @test hash(i1) == hash(i2) + @test islinearmap(i1) + @test isaffinemap(i1) + @test convert(StaticIdentityMap{SVector{2,T}}, i1) === i2 + @test jacobian(i1) isa ConstantMap + @test jacobian(i1, 1) == 1 + @test jacdet(i1, 1) == 1 + m1 = convert(FunctionMaps.AbstractAffineMap{T}, i1) + @test m1 isa LinearMap{T} + @test jacdet(m1, 1) == 1 + @test convert(FunctionMaps.AbstractAffineMap, i1) isa LinearMap{T} + m2 = convert(FunctionMaps.LinearMap{T}, i1) + @test m2 isa LinearMap{T} + @test jacdet(m2, 1) == 1 + + i3 = VectorIdentityMap{T}(10) + test_generic_map(i3) + r = rand(T, 10) + @test i3(r) ≈ r + @test hash(i3) == hash(VectorIdentityMap{Int}(10)) + + @test IdentityMap() ∘ LinearMap(2) == LinearMap(2.0) +end + +function test_basic_maps(T) + @test UnityMap{SVector{2,Float64}}() == UnityMap{SVector{2,Float64},Float64}() + @test hash(ConstantMap(2)) == hash(ConstantMap(2.0)) + @test FunctionMaps.absmap(ConstantMap(-2)) == ConstantMap(2) +end + +function test_wrapped_maps(T) + m1 = WrappedMap{T}(cos) + m2 = WrappedMap{T}(sin) + @test m1(one(T)) ≈ cos(one(T)) + @test m2(one(T)) ≈ sin(one(T)) + @test isequalmap(m1, cos) + @test isequalmap(sin, m2) + m3 = m1 ∘ m2 + @test m3(one(T)) ≈ cos(sin(one(T))) + + @test WrappedMap(cos) isa WrappedMap{Float64} + @test convert(Map, cos) isa WrappedMap + @test convert(Map{BigFloat}, m1) isa WrappedMap{BigFloat} + + @test convert(Map{T}, cos) isa FunctionMaps.WrappedMap{T,typeof(cos)} + + @test convert(Map, LinearAlgebra.I) isa GenericLinearMap{Vector{Any}} + @test convert(Map{Vector{T}}, LinearAlgebra.I) isa GenericLinearMap{Vector{T}} +end diff --git a/test/test_common.jl b/test/test_common.jl new file mode 100644 index 0000000..b573fee --- /dev/null +++ b/test/test_common.jl @@ -0,0 +1,97 @@ +function test_dimension() + @test euclideandimension(Int) == 1 + @test euclideandimension(Float64) == 1 + @test euclideandimension(ComplexF64) == 1 + @test euclideandimension(SVector{2,Float64}) == 2 + @test euclideandimension(MVector{2,Float64}) == 2 + @test euclideandimension(Tuple{Int,Int}) == 2 + @test euclideandimension(Tuple{Int,Float64}) == 2 + @test_throws ArgumentError euclideandimension(Vector{Float64}) +end + +function test_prectype() + @test prectype(1.0) == Float64 + @test prectype(big(1.0)) == BigFloat + @test prectype(1) == typeof(float(1)) + @test prectype(SVector(1,2)) == typeof(float(1)) + @test prectype(SVector(1,big(2))) == typeof(float(big(2))) + @test prectype(1.0+2.0im) == Float64 + @test prectype([1.0+2.0im, 3.0]) == Float64 + @test prectype(NTuple{2,Int}) == Float64 + @test prectype(typeof((1.0,))) == Float64 + @test prectype(typeof((1.0, 2.0))) == Float64 + @test prectype(typeof((1.0, 2.0, 3.0))) == Float64 + @test prectype(typeof((1.0, big(2.0), 3.0+im))) == BigFloat + @test prectype(NTuple{4,Int}) == Float64 + @test @inferred(prectype(1, 2.0)) == Float64 + @test @inferred(prectype(typeof((1, 2.0, 3, 40+im)))) == Float64 + + @test convert_prectype(Float64, 2) == 2 + @test convert_prectype(Float64, 2) isa Float64 + @test convert_prectype(BigFloat, 1.0+im) == 1+im + @test convert_prectype(BigFloat, 1.0+im) isa Complex{BigFloat} + @test convert_prectype(Float64, SA[1,2]) == SA[1.0,2.0] + @test convert_prectype(Float64, SA[1,2]) isa SVector{2,Float64} + @test convert_prectype(Float64, MVector(1,2)) == MVector(1.0,2.0) + @test convert_prectype(Float64, MVector(1,2)) isa MVector + @test convert_prectype(Float64, SA[1 2; 3 4]) == SA[1.0 2.0; 3.0 4.0] + @test convert_prectype(Float64, SA[1 2; 3 4]) isa SMatrix{2,2,Float64} + @test convert_prectype(Float64, MMatrix{2,2}(1, 2, 3, 4)) == SA[1.0 3.0; 2.0 4.0] + @test convert_prectype(Float64, MMatrix{2,2}(1, 2, 3, 4)) isa MMatrix{2,2,Float64} + @test convert_prectype(BigFloat, SA[1,2+im]) isa SVector{2,Complex{BigFloat}} + @test_throws ArgumentError convert_prectype(BigFloat, "a") + + @test promote_prectype(2) == 2 + @test promote_prectype(2, 3.0) isa Tuple{Float64,Float64} + @test promote_prectype(2, 3.0+im, big(4)) isa Tuple{BigFloat,Complex{BigFloat},BigFloat} +end + + +function test_numtype() + @test numtype(1.0) == Float64 + @test numtype(big(1.0)) == BigFloat + @test numtype(1) == Int + @test numtype([1,2,3], [4,5,6]) == Int + @test numtype(SVector(1,2)) == Int + @test numtype(SVector(1,big(2))) == BigInt + @test numtype(Array{Float64,2}) == Float64 + @test numtype(1.0+2.0im) == Complex{Float64} + @test numtype([1.0+2.0im, 3.0]) == Complex{Float64} + @test numtype(NTuple{2,Complex{Int}}) == Complex{Int} + @test numtype(typeof((1.0,))) == Float64 + @test numtype(typeof((1.0, 2.0))) == Float64 + @test numtype((1.0, 2.0, 3.0)) == Float64 + @test numtype((1.0, 2.0, 3.0, 4.0)) == Float64 + @test numtype(1.0, big(2.0), 3.0+im) == Complex{BigFloat} + @test numtype(typeof((1.0, big(2.0), 3.0+im))) == Complex{BigFloat} + @test @inferred(numtype(1, 2.0)) == Float64 + @test @inferred(numtype(typeof((1, 2.0, 3, 40+im)))) == Complex{Float64} + + @test convert_numtype(Float64, 2) == 2 + @test convert_numtype(Float64, 2) isa Float64 + @test convert_numtype(Float64, SA[1,2]) == SA[1.0,2.0] + @test convert_numtype(Float64, SA[1,2]) isa SVector{2,Float64} + @test convert_numtype(Float64, MVector(1,2)) == MVector(1.0,2.0) + @test convert_numtype(Float64, MVector(1,2)) isa MVector + @test convert_numtype(Float64, SA[1 2; 3 4]) == SA[1.0 2.0; 3.0 4.0] + @test convert_numtype(Float64, SA[1 2; 3 4]) isa SMatrix{2,2,Float64} + @test convert_numtype(Float64, MMatrix{2,2}(1, 2, 3, 4)) == SA[1.0 3.0; 2.0 4.0] + @test convert_numtype(Float64, MMatrix{2,2}(1, 2, 3, 4)) isa MMatrix{2,2,Float64} + @test_throws ArgumentError convert_numtype(BigFloat, "a") + + @test promote_numtype(2) == 2 + @test promote_numtype(2, 3.0) isa Tuple{Float64,Float64} + @test promote_numtype(2, 3.0+im, big(4)) isa Tuple{Complex{BigFloat},Complex{BigFloat},Complex{BigFloat}} +end + +@testset "common functionality" begin + @testset "dimension" begin + test_dimension() + end + @testset "prectype" begin + test_prectype() + end + @testset "numtype" begin + test_numtype() + end +end diff --git a/test/test_generic.jl b/test/test_generic.jl new file mode 100644 index 0000000..adad1b5 --- /dev/null +++ b/test/test_generic.jl @@ -0,0 +1,145 @@ +function generic_map_tests(T) + for map in maps_to_test(T) + @test prectype(map) == T + test_generic_map(map) + end +end + +# generic map test suite +function test_generic_map(m) + @test convert(Map{domaintype(m)}, m) == m + + x = suitable_point_to_map(m) + @test applymap(m, x) == m(x) + y = m(x) + @test isrealmap(m) == (isreal(x) && isreal(y)) + + S = domaintype(m) + U = codomaintype(m) + @test x isa S + @test y isa U + + if isaffinemap(m) && !isconstantmap(m) && !(prectype(m) == BigFloat) + test_generic_inverse(m) + else + try + # The map may not support an inverse, let's try + test_generic_inverse(m) + catch e + # make sure it failed because the inverse is not defined + @test e isa MethodError || prectype(m) == BigFloat + end + end + + if isaffinemap(m) + M = affinematrix(m) + @test size(M) == mapsize(m) + test_generic_jacobian(m) + else + try # jacobian may not be implemented + test_generic_jacobian(m) + catch e + end + end + + if domaintype(m) == Float64 + @test convert(Map{BigFloat}, m) isa Map{BigFloat} + @test convert(Map{BigFloat}, m) == m + end + if prectype(m) == Float64 + U = widertype(domaintype(m)) + @test convert(Map{U}, m) isa Map{U} + @test convert(Map{U}, m) == m + end + + if isaffinemap(m) + A = affinematrix(m) + b = affinevector(m) + @test m(x) == A*x+b + end + + @test hash(m) == map_hash(m) +end + +function test_generic_inverse(m) + M = mapsize(m,1) + N = mapsize(m,2) + x = suitable_point_to_map(m) + y = m(x) + + if M==N + inverse(m, y) # trigger exception outside test if not implemented + + minv = inverse(m) + @test minv(y) ≈ x + @test inverse(m)(y) ≈ x + @test inverse(m, y) ≈ x + @test m\y ≈ x + @test LazyInverse(m)(y) ≈ inverse(m, y) + end + if M >= N && numtype(m)!=BigFloat + leftinverse(m, y) # trigger exception outside test if not implemented + + mli = leftinverse(m) + @test mli(y) ≈ x + @test leftinverse(m, y) ≈ x + end + if M <= N && numtype(m)!=BigFloat + rightinverse(m, y) # trigger exception outside test if not implemented + + mri = rightinverse(m) + @test m(mri(y)) ≈ y + @test m(rightinverse(m, y)) ≈ y + end +end + +function test_generic_jacobian(m) + x = suitable_point_to_map(m) + # Trigger exception outside of test if jacobian(m, x) is not implemented. + # test_generic_jacobian(m) should be called within a try/catch block + jacobian(m, x) + + δ = sqrt(eps(prectype(m))) + x2 = x .+ δ + # we intentionally test jacobian(m, x) before testing jacobian(m) + if !(m isa ProductMap) + @test norm(m(x2) - m(x) + jacobian(m, x)*(x-x2)) < 100δ + end + jac = jacobian(m) + @test jac(x) ≈ jacobian(m, x) + j = jacobian(m, x) + @test size(j) == mapsize(m) + if j isa AbstractArray{Float64} + y = similar(j) + jacobian!(y, m, x) + @test y ≈ jac(x) + end + if issquarematrix(jac(x)) + @test jacdet(m, x) ≈ det(jacobian(m, x)) + end + if mapsize(m) == () + d = diffvolume(m, x) + s = sqrt(det(transpose(jacobian(m,x))*jacobian(m,x))) + @test d ≈ s || d ≈ -s + else + @test abs(diffvolume(m, x) - sqrt(det(transpose(jacobian(m,x))*jacobian(m,x)))) < 1e-7 + end + if isaffinemap(m) + @test diffvolume(m) isa ConstantMap + end +end + +function test_generic_functionality() + m = LinearMap(2) + @test convert_domaintype(Float64, m) isa ScalarLinearMap{Float64} + @test convert_domaintype(Complex{Float64}, m) isa ScalarLinearMap{Complex{Float64}} + @test convert_domaintype(SVector{2,BigFloat}, m) isa GenericLinearMap{SVector{2, BigFloat}, BigFloat} + @test convert_domaintype(Any, m) === m + @test convert_codomaintype(Float64, m) isa ScalarLinearMap{Float64} + @test convert_codomaintype(Complex{Float64}, m) isa ScalarLinearMap{Complex{Float64}} + @test_throws ArgumentError convert_codomaintype(SVector{2,BigFloat}, m) + + m2 = LinearMap{SVector{2,Float64}}(2) + @test convert_domaintype(SVector{2,BigFloat}, m2) isa GenericLinearMap{SVector{2, BigFloat}, BigFloat} + @test convert_codomaintype(SVector{2,BigFloat}, m2) isa GenericLinearMap{SVector{2, BigFloat}, BigFloat} +end diff --git a/test/test_interface.jl b/test/test_interface.jl new file mode 100644 index 0000000..6920ee2 --- /dev/null +++ b/test/test_interface.jl @@ -0,0 +1,16 @@ + +struct MySimpleMap end +FunctionMaps.MapStyle(::Type{MySimpleMap}) = IsMap() + +@testset "map interface" begin + @test MapStyle(2.0) isa NotMap + @test MapStyle(2.0) == MapStyle(typeof(2.0)) + @test_throws ArgumentError checkmap(2.0) + m = LinearMap(2) + @test MapStyle(m) isa IsMap + @test functionmap(m) === m + @test checkmap(m) === m + m2 = MySimpleMap() + @test MapStyle(m2) isa IsMap + @test checkmap(m2) == m2 +end diff --git a/test/test_maps.jl b/test/test_maps.jl new file mode 100644 index 0000000..264be85 --- /dev/null +++ b/test/test_maps.jl @@ -0,0 +1,191 @@ + +maps_to_test(T) = [ + StaticIdentityMap{T}(), + VectorIdentityMap{T}(10), + ConstantMap{T}(one(T)), + ConstantMap{T}(SVector{2,T}(1,2)), + ZeroMap{T}(), + UnityMap{T}(), + AffineMap(T(1.2), T(2.4)), # scalar map + AffineMap(-T(1.2), T(2.4)), # scalar map with negative A + AffineMap(randvec(T, 2, 2), randvec(T, 2)), # static map + AffineMap(randvec(T, 3, 2), randvec(T, 3)), # static map, rectangular + AffineMap(rand(T, 2, 2), rand(T, 2)), # vector map + AffineMap(rand(T, 3, 2), rand(T, 3)), # vector map, rectangular + AffineMap(LinearAlgebra.I, one(T)/2), # use UniformScaling object as A + AffineMap(LinearAlgebra.I, randvec(T,2)), # use UniformScaling object as A + GenericAffineMap(randvec(T, 2, 2), randvec(T, 2)), + GenericAffineMap(T(1.2), randvec(T, 2)), + GenericAffineMap(randvec(T, 3, 2), randvec(T, 3)), + Translation(randvec(T, 3)), + LinearMap(randvec(T, 2, 2)), + LinearMap(randvec(T, 2)), + LinearMap(randvec(T, 2, 2)) ∘ AffineMap(T(1.2), randvec(T, 2)), + AffineMap(5.0, 2.0) ∘ VectorToComplex{T}() ∘ UnitCircleMap{T}(), + LinearMap(SMatrix{2,2}(1,2,3,T(4))) ∘ CartToPolarMap{T}() ∘ LinearMap(SMatrix{2,2}(1,2,3,T(4))), + # Interval{Any}(0.0, 1.0) +] + +randvec(T,n) = SVector{n,T}(rand(n)) +randvec(T,m,n) = SMatrix{m,n,T}(rand(m,n)) + +suitable_point_to_map(m::Map) = suitable_point_to_map(m, domaintype(m)) + +suitable_point_to_map(m::Map, ::Type{SVector{N,T}}) where {N,T} = SVector{N,T}(rand(N)) +suitable_point_to_map(m::Map, ::Type{T}) where {T<:Number} = rand(T) +suitable_point_to_map(m::Map, ::Type{<:AbstractVector{T}}) where {T} = rand(T, mapsize(m,2)) + +suitable_point_to_map(m::ProductMap) = + map(suitable_point_to_map, components(m)) +suitable_point_to_map(m::VcatMap{T,M,N}) where {T,M,N} = + SVector{N,T}(rand(T,N)) + +suitable_point_to_map(::CartToPolarMap{T}) where {T} = randvec(T,2) +suitable_point_to_map(::PolarToCartMap{T}) where {T} = randvec(T,2) + +widertype(T) = widen(T) +widertype(::Type{SVector{N,T}}) where {N,T} = SVector{N,widen(T)} +widertype(::Type{Vector{T}}) where {T} = Vector{widen(T)} +widertype(::Type{Tuple{A}}) where {A} = Tuple{widen(A)} +widertype(::Type{Tuple{A,B}}) where {A,B} = Tuple{widen(A),widen(B)} +widertype(::Type{Tuple{A,B,C}}) where {A,B,C} = Tuple{widen(A),widen(B),widen(C)} +widertype(::Type{NTuple{N,T}}) where {N,T} = NTuple{N,widen(T)} + +issquarematrix(A) = false +issquarematrix(A::AbstractArray) = size(A,1)==size(A,2) + + +function test_maps() + @testset "generic functionality" begin + test_generic_functionality() + end + + @testset "generic map tests" begin + generic_map_tests(Float64) + generic_map_tests(BigFloat) + end + + # Test special maps + @testset "identity map" begin + test_identity_map(Float64) + test_identity_map(BigFloat) + end + @testset "basic maps" begin + test_basic_maps(Float64) + test_basic_maps(BigFloat) + end + @testset "affine maps" begin + test_affine_maps(Float64) + test_affine_maps(BigFloat) + end + @testset "composite maps" begin + test_composite_map(Float64) + test_composite_map(BigFloat) + end + @testset "product maps" begin + test_product_map(Float64) + test_product_map(BigFloat) + end + @testset "wrapped maps" begin + test_wrapped_maps(Float64) + test_wrapped_maps(BigFloat) + end + @testset "scaling maps" begin + test_scaling_maps(Float64) + test_scaling_maps(BigFloat) + end + @testset "isomorphisms" begin + test_isomorphisms(Float64) + test_isomorphisms(BigFloat) + end + @testset "Mixed maps" begin + test_mixed_maps() + end +end + +function test_composite_map(T) + a = T(0) + b = T(1) + c = T(2) + d = T(3) + ma = StaticIdentityMap{T}() + mb = interval_map(a, b, c, d) + + r = suitable_point_to_map(ma) + m1 = ma∘mb + test_generic_map(m1) + @test m1(r) ≈ ma(mb(r)) + m2 = m1∘mb + test_generic_map(m2) + @test m2(r) ≈ m1(mb(r)) + m3 = mb∘m2 + test_generic_map(m3) + @test m3(r) ≈ mb(m2(r)) + m = m2∘m3 + test_generic_map(m) + @test m(r) ≈ m2(m3(r)) + + 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 mapsize(ComposedMap(LinearMap(2),LinearMap(rand(T,2)),LinearMap(rand(T,2,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 sum_jacobian() == () + @test sum_jacobian(ma) == jacobian(ma) +end + +function test_mixed_maps() + m1 = composedmap(cos, sin) + @test domaintype(m1) == Any + @test m1(0.4) == sin(cos(0.4)) + + m2 = composedmap(AffineMap(2.0, 3.0), cos) + @test domaintype(m2) == Float64 + @test m2(0.4) ≈ cos(2*0.4+3) + @inferred m2(0.4) + @test repr(m2) == "cos ∘ (x -> 2.0 * x + 3.0)" + + m3 = composedmap(sin, AffineMap(2.0, 3.0), cos) + @test m3(0.5) ≈ cos(2*sin(0.5)+3) + @test repr(m3) == "cos ∘ (x -> 2.0 * x + 3.0) ∘ sin" + @test domaintype(m3) == Any + @inferred m3(0.5) + + m4 = productmap(sin, cos) + @test m4 isa TupleProductMap + @test domaintype(m4) == Tuple{Any,Any} + @test m4(0.3,0.5) == (sin(0.3), cos(0.5)) +end + +test_maps() diff --git a/test/test_product.jl b/test/test_product.jl new file mode 100644 index 0000000..db053c0 --- /dev/null +++ b/test/test_product.jl @@ -0,0 +1,48 @@ +function test_product_map(T) + ma = StaticIdentityMap{T}() + mb = interval_map(T(0), T(1), T(2), T(3)) + + r1 = suitable_point_to_map(ma) + r2 = suitable_point_to_map(ma) + r3 = suitable_point_to_map(ma) + r4 = suitable_point_to_map(ma) + r5 = suitable_point_to_map(ma) + + m1 = productmap(ma,mb) + test_generic_map(m1) + @test m1(SVector(r1,r2)) ≈ SVector(ma(r1),mb(r2)) + m2 = productmap(m1,mb) + test_generic_map(m2) + @test m2(SVector(r1,r2,r3)) ≈ SVector(ma(r1),mb(r2),mb(r3)) + m3 = productmap(mb,m2) + test_generic_map(m3) + @test m3(SVector(r1,r2,r3,r4)) ≈ SVector(mb(r1),ma(r2),mb(r3),mb(r4)) + m4 = productmap(m1,m2) + test_generic_map(m4) + @test m4(SVector(r1,r2,r3,r4,r5)) ≈ SVector(m1(SVector(r1,r2))...,m2(SVector(r3,r4,r5))...) + + @test ProductMap(ma, mb) == ProductMap([ma,mb]) + @test hash(ProductMap(ma, mb)) == hash(ProductMap([ma,mb])) + + m5 = productmap(AffineMap(SMatrix{2,2,T}(1.0,2,3,4), SVector{2,T}(1,3)), LinearMap{T}(2.0)) + @test domaintype(component(m5,1)) == SVector{2,T} + @test domaintype(component(m5,2)) == T + test_generic_map(m5) + x = SVector{3,T}(rand(T,3)) + @test m5(x) ≈ SVector(component(m5,1)(SVector(x[1],x[2]))...,component(m5,2)(x[3])) + + m6 = ProductMap([ma,mb]) + @test m6 isa VectorProductMap + @test mapsize(m6) == (2,2) + @test convert(Map{SVector{2,T}}, m6) isa VcatMap + test_generic_map(m6) + + @test components(productmap(LinearMap(1), LinearMap(one(T)))) isa Tuple{<:Map{T},<:Map{T}} + + # test a non-rectangular product map + m_rect = AffineMap(SA[one(T) 2; 3 4; 5 6], SA[one(T),one(T),zero(T)]) + pm = productmap(m_rect, m_rect) + @test pm isa VcatMap{T,6,4,(3,3),(2,2)} + @test jacobian(pm, SA[one(T),0,0,0]) isa SMatrix{6,4,T} + @test diffvolume(pm, SA[one(T),0,0,0]) == 24 +end diff --git a/test/using_fmaps.jl b/test/using_fmaps.jl new file mode 100644 index 0000000..1058146 --- /dev/null +++ b/test/using_fmaps.jl @@ -0,0 +1,53 @@ +using FunctionMaps + +# for test_common.jl +using FunctionMaps: + convert_numtype, + convert_prectype, + promote_numtype, + promote_prectype, + euclideandimension + +# for test_interface.jl +using FunctionMaps: + MapStyle, IsMap, NotMap, + functionmap, checkmap + +# for test_generic.jl +using FunctionMaps: + convert_domaintype, convert_codomaintype, + map_hash, + LazyInverse, jacobian! + +# for test_maps.jl +using FunctionMaps: ScalarAffineMap, + VectorAffineMap, + StaticAffineMap, + GenericAffineMap, + ScalarLinearMap, + VectorLinearMap, + StaticLinearMap, + GenericLinearMap, + ScalarTranslation, + VectorTranslation, + StaticTranslation, + GenericTranslation, + ProductMap, + TupleProductMap, VcatMap, VectorProductMap, + WrappedMap, + interval_map, multiply_map, + SumMap, sum_map, + composedmap, ComposedMap, + composite_jacobian, sum_jacobian, + CartToPolarMap, PolarToCartMap, + UnitCircleMap, AngleMap, UnitDiskMap, + VectorToComplex + +# for test_affine.jl +using FunctionMaps: + to_matrix, to_vector, + matrix_pinv + +# for test_arithmetics.jl +using FunctionMaps: + interval_map, bounded_interval_map