From 4376b8b654a2ede3185522fdef05bb37289ac237 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Fri, 24 Sep 2021 07:28:11 +0200 Subject: [PATCH 01/20] Squashed commit of the following: MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit commit 26ad8a51be8497e2aa18b210d03bd114b3c9bd03 Author: Ronny Bergmann Date: Thu Sep 23 20:17:01 2021 +0200 runs formatter. commit 99bcb23c981308d58a2721bd90fe1773d53728bd Merge: b027f4c82 2a03fb147 Author: Ronny Bergmann Date: Thu Sep 23 18:33:25 2021 +0200 Merge branch 'master' into kellertuer/polish-diff-interfaces commit 2a03fb14795c216768b6212b2ed170ab3e4e9a40 Author: Mateusz Baran Date: Thu Sep 23 18:08:23 2021 +0200 Support for Zygote and ReverseDiff gradients (#427) * Support for Zygote and ReverseDiff gradients * Add Zygote test dependency * bump ambiguity limit because of Zygote * fix tests and Zygote backend * bump Julia to 1.5 * fixing some issues on Julia 1.7-rc1 * more fixing for Julia 1.7 * formatting * reduce tangent vector length in a test since the approximation only works very locally (and they changed the default random number generator) * Update Project.toml Co-authored-by: Ronny Bergmann * bump atol on Rotations * reduce tangent vector size even further. * adapt one more tolerance Co-authored-by: Ronny Bergmann commit b027f4c82e190ce657631ab6e39efae4c571f7de Author: Ronny Bergmann Date: Thu Sep 23 17:11:19 2021 +0200 fix a few typos. commit dcc471e273732185acd0670459667b29996924fa Author: Ronny Bergmann Date: Thu Sep 23 13:53:07 2021 +0200 unify naming. commit 3de9a309fa9b3f29d3e368275a192f1c85e07f01 Author: Ronny Bergmann Date: Thu Sep 23 13:33:59 2021 +0200 Reduce dependencies, wort a little bit further on polishing. commit 87029029a80f1ccb8562c5feaf0115a91fa48185 Author: Ronny Bergmann Date: Thu Sep 23 10:45:36 2021 +0200 runs formatter. commit 378bbbafa121cf00739dc04c9abe1d44907278a6 Author: Ronny Bergmann Date: Thu Sep 23 10:45:15 2021 +0200 fix the other jacobian. commit a7d1edbd2b18b9650ec6c7583ab89f077bedc05e Author: Ronny Bergmann Date: Thu Sep 23 10:29:37 2021 +0200 renaming and dependencies. commit 02f97d69904a6f9ae3c652a03c2e4c49d1096084 Merge: aca6d5c2c 9fc772a00 Author: Ronny Bergmann Date: Wed Sep 22 19:48:36 2021 +0200 Merge branch 'master' into kellertuer/polish-diff-interfaces # Conflicts: # Project.toml # src/Manifolds.jl commit aca6d5c2cf66119873b1727c6d1241847602decc Author: Ronny Bergmann Date: Wed Sep 22 19:46:59 2021 +0200 Add a proper working default. commit 9fc772a00e256db7f0ad51dfe3f063849128942c Author: Ronny Bergmann Date: Wed Sep 22 19:46:01 2021 +0200 Sketch change metric. (#423) * Adds a change_metric and a change_representer * Apply suggestions from code review * bump version. Co-authored-by: Mateusz Baran commit 14620189df19c526dfdad89469fcb49e1bb79b44 Merge: 923a718f0 4501a271f Author: Ronny Bergmann Date: Wed Sep 22 19:20:52 2021 +0200 Merge branch 'kellertuer/gradient-conversion' into kellertuer/polish-diff-interfaces # Conflicts: # Project.toml # src/differentiation/riemannian_diff.jl commit 923a718f01c3f1a95965990477067cab99716b34 Merge: d39148376 07c905eab Author: Ronny Bergmann Date: Wed Sep 22 18:20:39 2021 +0200 Merge branch 'mbaran/reverse-diff-zygote' into kellertuer/polish-diff-interfaces commit 07c905eab7de9aa0ca827ce6977ce7c5fc4ae980 Author: Mateusz Baran Date: Wed Sep 22 17:27:10 2021 +0200 Update Project.toml Co-authored-by: Ronny Bergmann commit 4501a271ffe1bf94d28d5a1389b27642a414c21e Author: Ronny Bergmann Date: Wed Sep 22 17:19:42 2021 +0200 Apply suggestions from code review Co-authored-by: Mateusz Baran commit 3c47350fd3ce87550757b81020aeb00b6dc5b7d7 Author: Ronny Bergmann Date: Wed Sep 22 08:59:28 2021 +0200 reduce tangent vector length in a test since the approximation only works very locally (and they changed the default random number generator) commit 3b7396d3329073a31328220f87833960f700f0bc Author: Mateusz Baran Date: Tue Sep 21 19:19:09 2021 +0200 formatting commit 5d75e6ebc6b05a63c7c54937bf478e41bdecc651 Author: Mateusz Baran Date: Tue Sep 21 19:17:28 2021 +0200 more fixing for Julia 1.7 commit c723345258687723479167c18680a683b29e4c5e Author: Mateusz Baran Date: Tue Sep 21 12:31:30 2021 +0200 fixing some issues on Julia 1.7-rc1 commit cf00e220a9d1c1c73b9b4f255e5c42be27b5ba7b Author: Mateusz Baran Date: Tue Sep 21 08:49:40 2021 +0200 bump Julia to 1.5 commit 1ed46db14f8152c3235cfd40e789703ca57ff032 Author: Mateusz Baran Date: Mon Sep 20 23:02:03 2021 +0200 fix tests and Zygote backend commit 7849fc0a3057c590df189f184bc2ac3f6fa67cde Author: Mateusz Baran Date: Mon Sep 20 19:45:41 2021 +0200 bump ambiguity limit because of Zygote commit 91634807265d1da5fe38b3486122bcc1d1b7b043 Author: Mateusz Baran Date: Mon Sep 20 14:40:50 2021 +0200 Add Zygote test dependency commit 7ed57dac28cecd8f662da1bb85180751d90118fa Author: Mateusz Baran Date: Mon Sep 20 12:48:09 2021 +0200 Support for Zygote and ReverseDiff gradients commit 2f498dace948ca30f13a072de5b03bc9dab31230 Author: Ronny Bergmann Date: Sat Sep 18 23:17:32 2021 +0200 add a test that double metric wrapping is avoided; remove unnecessary dispatch. commit e180bc4d0cf1b6376e5ff99ffcf0075c5bd60d17 Author: Ronny Bergmann Date: Sat Sep 18 21:58:49 2021 +0200 Simplify default implementation. Extend coverage. Order functions alphabetically. commit 420047d837ee6f30f21ff422845af53be9f2d09a Author: Ronny Bergmann Date: Sat Sep 18 20:07:21 2021 +0200 Fix the error message for hyperbolic. commit ed2871c2a81681e4768d6142366487d5712ea9a1 Author: Ronny Bergmann Date: Sat Sep 18 19:54:19 2021 +0200 maybe fix the new test by slightly overtyping the basis. commit 49dc500ca5cdd44b350128e00886942dd43f7a5e Author: Ronny Bergmann Date: Fri Sep 17 18:44:59 2021 +0200 format this interims non-working test. commit a78ace43b0f9099b65963e6f5ff27f9aa43fd2e5 Author: Ronny Bergmann Date: Fri Sep 17 18:40:57 2021 +0200 Maybe. just maybe fixed local metric, but now get vector and get coordinates do not work for the test. commit 97417aea2259cbdb2ed3e3d787bcf078d7db9b86 Author: Ronny Bergmann Date: Thu Sep 16 14:31:49 2021 +0200 Improve docs. commit d34ef6c63d5247c1cf09839d0eeab4b2c324bf15 Merge: 95a3c8bcf d1dbb7cd7 Author: Ronny Bergmann Date: Thu Sep 16 11:09:24 2021 +0200 Merge branch 'kellertuer/gradient-conversion' of github.com:JuliaManifolds/Manifolds.jl into kellertuer/gradient-conversion commit 95a3c8bcfe15e1dfe4c02f8ed008cf5fa1f33ce4 Author: Ronny Bergmann Date: Thu Sep 16 11:09:16 2021 +0200 adds a test for positive vectors and starts testing the default Metric fallbacks using the local metric. commit d1dbb7cd73e5f374b92528a349a10100ffd79b27 Author: Ronny Bergmann Date: Thu Sep 16 10:47:57 2021 +0200 Apply suggestions from code review Co-authored-by: Mateusz Baran commit 14fdb06d2020fb34298fc48d7f837865de1ca21e Author: Ronny Bergmann Date: Wed Sep 15 22:39:38 2021 +0200 bump version. commit cbcef945dd8e2c23e04f4a8f730487d607801335 Merge: da98b5f86 bd2b42599 Author: Ronny Bergmann Date: Wed Sep 15 22:38:45 2021 +0200 Merge branch 'master' into kellertuer/gradient-conversion commit da98b5f860da76b24d1193766d8c1e2c82a368cd Author: Ronny Bergmann Date: Wed Sep 15 22:37:18 2021 +0200 append this step to riemannian differentiation. commit c1545164b7cf92c3ef4d302496a56ff766aa5ab4 Author: Ronny Bergmann Date: Wed Sep 15 22:37:01 2021 +0200 finish tests for product and power – extend features to also work on the sphere. commit 844f8c4e79aa66cb4877625151212791a223ac4f Author: Ronny Bergmann Date: Wed Sep 15 21:42:59 2021 +0200 Fixes a few typos in the docs. commit 5005965700e7e0d13f59a3823ad529d7ae19f6dc Author: Ronny Bergmann Date: Wed Sep 15 20:40:11 2021 +0200 Document and test SPDs commit d2cec4d838689a2ed9d4619987911dea53386c03 Author: Ronny Bergmann Date: Wed Sep 15 20:30:00 2021 +0200 Finish positive Numbers. commit 7737c924d669ce64ba4339db7c1a1d31e1820ca2 Author: Ronny Bergmann Date: Wed Sep 15 20:03:49 2021 +0200 Extend ProbabilitySimples to also cover change_metric. commit 6300f8025c7ddfb4a915b340bea402367a2aa924 Author: Ronny Bergmann Date: Wed Sep 15 20:03:33 2021 +0200 delete change_representer since it should be inherited from ProbabilitySimplex combined with AbstractPowerManifold and EmbeddedManifold. commit 95b8bc6b3064917d08a6bdbd1413b2fc5f013782 Author: Ronny Bergmann Date: Wed Sep 15 19:48:20 2021 +0200 finishes testing for hyperbolic. commit 6f1fb647d13a634bc86d929acad171887b4e2e07 Author: Ronny Bergmann Date: Wed Sep 15 19:47:37 2021 +0200 fix a small bug in Hyperbolic copyto! which returned the wrong value. commit a62e134bb961193acda299b4eed26eecbe5ec3e2 Author: Ronny Bergmann Date: Wed Sep 15 10:31:24 2021 +0200 Documentation and Test for Poincaré Ball conversion. commit d84a3cc6b835d51042ece7cdcfb75b284b28f16d Author: Ronny Bergmann Date: Tue Sep 14 19:52:52 2021 +0200 adds two further small tests. commit ea76128bead6cf3c58855608ed2174cbbbca1fdf Author: Ronny Bergmann Date: Tue Sep 14 19:21:08 2021 +0200 Document and test hyperbolic. commit 257b2233b03953ec7a9ba55559f01135e70fd2c8 Author: Ronny Bergmann Date: Tue Sep 14 18:20:55 2021 +0200 Finish documentation and test for Generalised Grassmann. commit 152e9a164fc64836e5caed4a01c7e4556c99abbc Author: Ronny Bergmann Date: Tue Sep 14 18:20:34 2021 +0200 fix transparency of nonmutating to parent on a Manifold level. commit 1b79065ec525a81156ffbc19de248ad6fece7f50 Author: Ronny Bergmann Date: Mon Sep 13 17:49:20 2021 +0200 Fix two typos. commit 2eac3ee3e9d3dd304bb2c2458db930145924df8c Author: Ronny Bergmann Date: Mon Sep 13 16:30:35 2021 +0200 Change `change_gradient` to `change_represender` and introduce mutating variant commit 5539f4b470b2eb31234f97cf0c6a63fb4324b2f8 Author: Ronny Bergmann Date: Mon Sep 13 12:11:52 2021 +0200 introduce mutating variants and decorator transparency. commit 714862e41f051f6cfa6a6469f131ead24783ff6d Merge: 4912e4feb 05968dc42 Author: Ronny Bergmann Date: Sun Sep 12 20:09:17 2021 +0200 Merge branch 'kellertuer/gradient-conversion' of github.com:JuliaManifolds/Manifolds.jl into kellertuer/gradient-conversion commit 4912e4feb353df8d8073c659a8ce243f652d404a Author: Ronny Bergmann Date: Sun Sep 12 20:09:11 2021 +0200 Adress points from code review. commit 05968dc4288ab5498f67aeb800309e9a8df88dfa Author: Ronny Bergmann Date: Sun Sep 12 20:03:18 2021 +0200 Apply suggestions from code review Co-authored-by: Mateusz Baran commit 9769a1a2c6424a43ae445a88114e7cbc16383941 Author: Ronny Bergmann Date: Sun Sep 12 18:17:45 2021 +0200 add the generic case for doubly stochastic (includes symmetric, too). commit 42875a6605fa332e91502f442671c42aa0e05b57 Author: Ronny Bergmann Date: Sun Sep 12 18:13:42 2021 +0200 implement the remaining conversions. commit fb304c27e1bc17c9e201ccc5c203216590cca49a Merge: cb2f93278 2ee207293 Author: Ronny Bergmann Date: Sat Sep 11 16:51:29 2021 +0200 Merge branch 'kellertuer/gradient-conversion' of github.com:JuliaManifolds/Manifolds.jl into kellertuer/gradient-conversion commit cb2f93278c51e0bb5c30932e14ea9f854ef0518d Author: Ronny Bergmann Date: Sat Sep 11 16:51:17 2021 +0200 adds a few conversion, where the embedding is not isometric. commit 2ee207293e19c83e4224426081f3db26832ba83d Author: Ronny Bergmann Date: Sat Sep 11 13:04:01 2021 +0200 Update src/manifolds/MetricManifold.jl Co-authored-by: Mateusz Baran commit c742ea529dfb798ed48550a7e333f4b5776e9547 Author: Ronny Bergmann Date: Sat Sep 11 12:55:35 2021 +0200 Fix two typos. commit a17b4ca64ac12946e0d1796aba6763fd5ae9658f Merge: 540f3747b 056e24ca1 Author: Ronny Bergmann Date: Sat Sep 11 12:47:27 2021 +0200 Merge branch 'kellertuer/gradient-conversion' of github.com:JuliaManifolds/Manifolds.jl into kellertuer/gradient-conversion # Conflicts: # src/manifolds/MetricManifold.jl commit 540f3747bcff7ca639190019437ed50e58f9ef39 Author: Ronny Bergmann Date: Sat Sep 11 12:46:50 2021 +0200 Specify Reisz representer, change adjoint to ^H. commit 056e24ca1c8be5bcd4b9671bcb153a17e7fa8c1f Author: Ronny Bergmann Date: Sat Sep 11 12:38:40 2021 +0200 Update src/manifolds/MetricManifold.jl Co-authored-by: Mateusz Baran commit 158bbf221f7a593a83dbccc2867ba3a69f10ebe4 Author: Ronny Bergmann Date: Sat Sep 11 12:33:10 2021 +0200 Apply suggestions from code review Co-authored-by: Mateusz Baran commit 616db5d62ee6fdeb39e488c20644a879293c8c60 Author: Ronny Bergmann Date: Sat Sep 11 10:39:47 2021 +0200 Fix a few typos in docs. commit c9737b3d05589352899a84d9a30cb1bcc9a7595c Author: Ronny Bergmann Date: Sat Sep 11 10:21:02 2021 +0200 adds a change_gradient. commit 0857b2149f990ccd698bdbe2562700ac908335bf Author: Ronny Bergmann Date: Sat Sep 11 09:25:54 2021 +0200 fix a few typos. commit a3027332e9ce52bd47d55a77397d82d1ce3cc259 Author: Ronny Bergmann Date: Sat Sep 11 09:04:07 2021 +0200 rfF. commit 334793066ecb5b647c4d6f79df0a2ef660611d2a Author: Ronny Bergmann Date: Sat Sep 11 09:01:30 2021 +0200 Document insights from yesterday evening and – getting stuck at the next point. commit ae0cdd1031f005c9a36e7d0a05c35d35c19ee489 Author: Ronny Bergmann Date: Fri Sep 10 20:38:56 2021 +0200 Sketch change metric. --- .github/workflows/ci.yml | 2 +- Project.toml | 10 +- src/Manifolds.jl | 37 ++- src/differentiation.jl | 144 ---------- src/differentiation/differentiation.jl | 117 +++++++++ src/differentiation/finite_diff.jl | 38 +++ src/differentiation/finite_differences.jl | 28 ++ src/{ => differentiation}/forward_diff.jl | 4 +- src/differentiation/reverse_diff.jl | 13 + src/differentiation/riemannian_diff.jl | 212 +++++++++++++++ src/differentiation/zygote.jl | 13 + src/finite_diff.jl | 25 -- src/groups/group.jl | 2 +- src/manifolds/CenteredMatrices.jl | 4 +- src/manifolds/ConnectionManifold.jl | 38 ++- src/manifolds/GeneralizedGrassmann.jl | 33 +++ src/manifolds/Hyperbolic.jl | 5 +- src/manifolds/HyperbolicHyperboloid.jl | 28 +- src/manifolds/HyperbolicPoincareBall.jl | 45 ++++ src/manifolds/MetricManifold.jl | 224 ++++++++++++++-- src/manifolds/Multinomial.jl | 1 + src/manifolds/PositiveNumbers.jl | 49 ++++ src/manifolds/PowerManifold.jl | 42 +++ src/manifolds/ProbabilitySimplex.jl | 38 +++ src/manifolds/ProductManifold.jl | 36 +++ src/manifolds/Spectrahedron.jl | 2 +- src/manifolds/Sphere.jl | 6 +- .../SymmetricPositiveDefiniteLinearAffine.jl | 47 ++++ src/ode.jl | 2 +- src/riemannian_diff.jl | 248 ------------------ test/ambiguities.jl | 18 +- test/approx_inverse_retraction.jl | 4 +- test/differentiation.jl | 157 +++++------ test/groups/special_euclidean.jl | 2 + test/manifolds/generalized_grassmann.jl | 75 +++--- test/manifolds/hyperbolic.jl | 29 +- test/manifolds/positive_numbers.jl | 6 + test/manifolds/power_manifold.jl | 21 +- test/manifolds/probability_simplex.jl | 6 + test/manifolds/product_manifold.jl | 22 ++ test/manifolds/rotations.jl | 2 +- test/manifolds/sphere.jl | 9 + test/manifolds/stiefel.jl | 2 +- test/manifolds/symmetric_positive_definite.jl | 33 ++- test/metric.jl | 69 ++++- 45 files changed, 1338 insertions(+), 610 deletions(-) delete mode 100644 src/differentiation.jl create mode 100644 src/differentiation/differentiation.jl create mode 100644 src/differentiation/finite_diff.jl create mode 100644 src/differentiation/finite_differences.jl rename src/{ => differentiation}/forward_diff.jl (81%) create mode 100644 src/differentiation/reverse_diff.jl create mode 100644 src/differentiation/riemannian_diff.jl create mode 100644 src/differentiation/zygote.jl delete mode 100644 src/finite_diff.jl delete mode 100644 src/riemannian_diff.jl diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index efd106f843..2e3a0a6108 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -11,7 +11,7 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - julia-version: ["1.4", "1.6"] + julia-version: ["1.5", "1.6", "~1.7.0-0"] os: [ubuntu-latest, macOS-latest] steps: - uses: actions/checkout@v2 diff --git a/Project.toml b/Project.toml index b1cca42bea..d731092df7 100644 --- a/Project.toml +++ b/Project.toml @@ -1,13 +1,12 @@ name = "Manifolds" uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" authors = ["Seth Axen ", "Mateusz Baran ", "Ronny Bergmann ", "Antoine Levitt "] -version = "0.6.7" +version = "0.6.9" [deps] Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Einsum = "b7d42ee7-0b51-5a75-98ca-779d3107e4c0" -FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" HybridArrays = "1baab800-613f-4b0a-84e4-9cd3431bfbb9" Kronecker = "2c470bb0-bcc8-11e8-3dad-c9649493f05e" LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d" @@ -29,7 +28,6 @@ Colors = "0.12" Distributions = "0.22.6, 0.23, 0.24, 0.25" Einsum = "0.4" FiniteDiff = "2" -FiniteDifferences = "0.12" HybridArrays = "0.4" Kronecker = "0.4, 0.5" LightGraphs = "1" @@ -42,12 +40,13 @@ SimpleWeightedGraphs = "1" SpecialFunctions = "0.8, 0.9, 0.10, 1.0" StaticArrays = "1.0" StatsBase = "0.32, 0.33" -julia = "1.4" +julia = "1.5" [extras] Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" +FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Gtk = "4c0ca9eb-093a-5379-98c5-f87ac0bbbf44" ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19" @@ -62,6 +61,7 @@ RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" VisualRegressionTests = "34922c18-7c2a-561c-bac1-01e79b2c4c92" +Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [targets] -test = ["Test", "Colors", "DoubleFloats", "FiniteDiff", "ForwardDiff", "Gtk", "ImageIO", "ImageMagick", "OrdinaryDiffEq", "NLsolve", "Plots", "PyPlot", "Quaternions", "QuartzImageIO", "RecipesBase", "ReverseDiff"] +test = ["Test", "Colors", "DoubleFloats", "FiniteDiff", "FiniteDifferences", "ForwardDiff", "Gtk", "ImageIO", "ImageMagick", "OrdinaryDiffEq", "NLsolve", "Plots", "PyPlot", "Quaternions", "QuartzImageIO", "RecipesBase", "ReverseDiff", "Zygote"] diff --git a/src/Manifolds.jl b/src/Manifolds.jl index ddf474c2ef..98022f2ad3 100644 --- a/src/Manifolds.jl +++ b/src/Manifolds.jl @@ -86,7 +86,6 @@ import Base: using Base.Iterators: repeated using Distributions using Einsum: @einsum -using FiniteDifferences using HybridArrays using Kronecker using LightGraphs @@ -108,7 +107,6 @@ using ManifoldsBase: AbstractLinearVectorTransportMethod, ApproximateInverseRetraction, ApproximateRetraction, - DifferentiatedRetractionVectorTransport, ComponentManifoldError, CompositeManifoldError, CotangentSpaceType, @@ -158,8 +156,8 @@ using RecursiveArrayTools: ArrayPartition include("utils.jl") include("product_representations.jl") -include("differentiation.jl") -include("riemannian_diff.jl") +include("differentiation/differentiation.jl") +include("differentiation/riemannian_diff.jl") # Main Meta Manifolds include("manifolds/ConnectionManifold.jl") @@ -285,12 +283,17 @@ end function __init__() @require FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" begin using .FiniteDiff - include("finite_diff.jl") + include("differentiation/finite_diff.jl") + end + + @require FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" begin + using .FiniteDifferences + include("differentiation/finite_differences.jl") end @require ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" begin using .ForwardDiff - include("forward_diff.jl") + include("differentiation/forward_diff.jl") end @require OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" begin @@ -303,6 +306,11 @@ function __init__() include("nlsolve.jl") end + @require ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" begin + using .ReverseDiff: ReverseDiff + include("differentiation/reverse_diff.jl") + end + @require Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" begin using .Test: Test include("tests/tests_general.jl") @@ -333,6 +341,12 @@ function __init__() include("recipes.jl") end end + + @require Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" begin + using .Zygote: Zygote + include("differentiation/zygote.jl") + end + return nothing end @@ -460,6 +474,10 @@ export ×, allocate_result, base_manifold, bundle_projection, + change_metric, + change_metric!, + change_representer, + change_representer!, check_point, check_vector, christoffel_symbols_first, @@ -648,8 +666,11 @@ export get_basis, get_coordinates, get_coordinates!, get_vector, get_vector!, get_vectors, number_system # differentiation export AbstractDiffBackend, - AbstractRiemannianDiffBackend, FiniteDifferencesBackend, RiemannianONBDiffBackend -export diff_backend, diff_backend!, diff_backends + AbstractRiemannianDiffBackend, + FiniteDifferencesBackend, + TangentDiffBackend, + RiemannianProjectionBackend +export default_differential_backend, set_default_differential_backend! # atlases and charts export get_point, get_point!, get_parameters, get_parameters! diff --git a/src/differentiation.jl b/src/differentiation.jl deleted file mode 100644 index 17ed083186..0000000000 --- a/src/differentiation.jl +++ /dev/null @@ -1,144 +0,0 @@ - -""" - AbstractDiffBackend - -An abstract type for diff backends. See [`FiniteDifferencesBackend`](@ref) for -an example. -""" -abstract type AbstractDiffBackend end - -struct NoneDiffBackend <: AbstractDiffBackend end - -""" - _derivative(f, t[, backend::AbstractDiffBackend]) - -Compute the derivative of a callable `f` at time `t` computed using the given `backend`, -an object of type [`Manifolds.AbstractDiffBackend`](@ref). If the backend is not explicitly -specified, it is obtained using the function [`Manifolds.diff_backend`](@ref). - -This function calculates plain Euclidean derivatives, for Riemannian differentiation see -for example [`differential`](@ref Manifolds.differential(::AbstractManifold, ::Any, ::Real, ::AbstractRiemannianDiffBackend)). - -!!! note - - Not specifying the backend explicitly will usually result in a type instability - and decreased performance. -""" -function _derivative end - -function _derivative!(f, X, t, backend::AbstractDiffBackend) - return copyto!(X, _derivative(f, t, backend)) -end - -""" - _gradient(f, p[, backend::AbstractDiffBackend]) - -Compute the gradient of a callable `f` at point `p` computed using the given `backend`, -an object of type [`AbstractDiffBackend`](@ref). If the backend is not explicitly -specified, it is obtained using the function [`diff_backend`](@ref). - -This function calculates plain Euclidean gradients, for Riemannian gradient calculation see -for example [`gradient`](@ref Manifolds.gradient(::AbstractManifold, ::Any, ::Any, ::AbstractRiemannianDiffBackend)). - -!!! note - - Not specifying the backend explicitly will usually result in a type instability - and decreased performance. -""" -function _gradient end - -function _gradient!(f, X, p, backend::AbstractDiffBackend) - return copyto!(X, _gradient(f, p, backend)) -end - -""" - CurrentDiffBackend(backend::AbstractDiffBackend) - -A mutable struct for storing the current differentiation backend in a global -constant [`_current_diff_backend`](@ref). - -# See also - -[`AbstractDiffBackend`](@ref), [`diff_backend`](@ref), [`diff_backend!`](@ref) -""" -mutable struct CurrentDiffBackend - backend::AbstractDiffBackend -end - -""" - _current_diff_backend - -The instance of [`Manifolds.CurrentDiffBackend`](@ref) that stores the globally default -differentiation backend. -""" -const _current_diff_backend = CurrentDiffBackend(NoneDiffBackend()) - -""" - _diff_backends - -A vector of valid [`Manifolds.AbstractDiffBackend`](@ref). -""" -const _diff_backends = AbstractDiffBackend[] - -""" - diff_backend() -> AbstractDiffBackend - -Get the current differentiation backend. -""" -diff_backend() = _current_diff_backend.backend - -""" - diff_backend!(backend::AbstractDiffBackend) - -Set current backend for differentiation to `backend`. -""" -function diff_backend!(backend::AbstractDiffBackend) - _current_diff_backend.backend = backend - return backend -end - -""" - diff_backends() -> Vector{AbstractDiffBackend} - -Get vector of currently valid differentiation backends. -""" -diff_backends() = _diff_backends - -_derivative(f, t) = _derivative(f, t, diff_backend()) - -_derivative!(f, X, t) = _derivative!(f, X, t, diff_backend()) - -_gradient(f, p) = _gradient(f, p, diff_backend()) - -_gradient!(f, X, p) = _gradient!(f, X, p, diff_backend()) - -# Finite differences - -""" - FiniteDifferencesBackend(method::FiniteDifferenceMethod = central_fdm(5, 1)) - -Differentiation backend based on the FiniteDifferences package. -""" -struct FiniteDifferencesBackend{TM<:FiniteDifferenceMethod} <: AbstractDiffBackend - method::TM -end - -function FiniteDifferencesBackend() - return FiniteDifferencesBackend(central_fdm(5, 1)) -end - -push!(_diff_backends, FiniteDifferencesBackend()) - -diff_backend!(_diff_backends[end]) - -function _derivative(f, t, backend::FiniteDifferencesBackend) - return backend.method(f, t) -end - -function _gradient(f, p, backend::FiniteDifferencesBackend) - return FiniteDifferences.grad(backend.method, f, p)[1] -end - -function _jacobian(f, p, backend::FiniteDifferencesBackend) - return FiniteDifferences.jacobian(backend.method, f, p)[1] -end diff --git a/src/differentiation/differentiation.jl b/src/differentiation/differentiation.jl new file mode 100644 index 0000000000..67e64ca1db --- /dev/null +++ b/src/differentiation/differentiation.jl @@ -0,0 +1,117 @@ + +""" + AbstractDiffBackend + +An abstract type for diff backends. See [`FiniteDifferencesBackend`](@ref) for +an example. +""" +abstract type AbstractDiffBackend end + +struct NoneDiffBackend <: AbstractDiffBackend end + +""" + _derivative(f, t[, backend::AbstractDiffBackend]) + +Compute the derivative of a callable `f` at time `t` computed using the given `backend`, +an object of type [`Manifolds.AbstractDiffBackend`](@ref). If the backend is not explicitly +specified, it is obtained using the function [`default_differential_backend`](@ref). + +This function calculates plain Euclidean derivatives, for Riemannian differentiation see +for example [`differential`](@ref Manifolds.differential(::AbstractManifold, ::Any, ::Real, ::AbstractRiemannianDiffBackend)). + +!!! note + + Not specifying the backend explicitly will usually result in a type instability + and decreased performance. +""" +function _derivative end + +_derivative(f, t) = _derivative(f, t, default_differential_backend()) + +function _derivative!(f, X, t, backend::AbstractDiffBackend=default_differential_backend()) + return copyto!(X, _derivative(f, t, backend)) +end + +""" + _gradient(f, p[, backend::AbstractDiffBackend]) + +Compute the gradient of a callable `f` at point `p` computed using the given `backend`, +an object of type [`AbstractDiffBackend`](@ref). If the backend is not explicitly +specified, it is obtained using the function [`default_differential_backend`](@ref). + +This function calculates plain Euclidean gradients, for Riemannian gradient calculation see +for example [`gradient`](@ref Manifolds.gradient(::AbstractManifold, ::Any, ::Any, ::AbstractRiemannianDiffBackend)). + +!!! note + + Not specifying the backend explicitly will usually result in a type instability + and decreased performance. +""" +function _gradient end + +_gradient(f, p) = _gradient(f, p, default_differential_backend()) + +function _gradient!(f, X, p, backend::AbstractDiffBackend=default_differential_backend()) + return copyto!(X, _gradient(f, p, backend)) +end + +""" + _jacobian(f, p[, backend::AbstractDiffBackend]) + +Compute the jacobian of a callable `f` at point `p` computed using the given `backend`, +an object of type [`AbstractDiffBackend`](@ref). If the backend is not explicitly +specified, it is obtained using the function [`default_differential_backend`](@ref). + +This function calculates plain Euclidean gradients, for Riemannian gradient calculation see +for example [`gradient`](@ref Manifolds.gradient(::AbstractManifold, ::Any, ::Any, ::AbstractRiemannianDiffBackend)). + +!!! note + + Not specifying the backend explicitly will usually result in a type instability + and decreased performance. +""" +function _jacobian end + +_jacobian(f, p) = _jacobian(f, p, default_differential_backend()) + +function _jacobian!(f, X, p, backend::AbstractDiffBackend=default_differential_backend()) + return copyto!(X, _jacobian(f, p, backend)) +end + +""" + CurrentDiffBackend(backend::AbstractDiffBackend) + +A mutable struct for storing the current differentiation backend in a global +constant [`_current_default_differential_backend`](@ref). + +# See also + +[`AbstractDiffBackend`](@ref), [`default_differential_backend`](@ref), [`set_default_differential_backend`](@ref) +""" +mutable struct CurrentDiffBackend + backend::AbstractDiffBackend +end + +""" + _current_default_differential_backend + +The instance of [`Manifolds.CurrentDiffBackend`](@ref) that stores the globally default +differentiation backend. +""" +const _current_default_differential_backend = CurrentDiffBackend(NoneDiffBackend()) +""" + default_differential_backend() -> AbstractDiffBackend + +Get the default differentiation backend. +""" +default_differential_backend() = _current_default_differential_backend.backend + +""" + set_default_differential_backend!(backend::AbstractDiffBackend) + +Set current backend for differentiation to `backend`. +""" +function set_default_differential_backend!(backend::AbstractDiffBackend) + _current_default_differential_backend.backend = backend + return backend +end diff --git a/src/differentiation/finite_diff.jl b/src/differentiation/finite_diff.jl new file mode 100644 index 0000000000..bd5403d60a --- /dev/null +++ b/src/differentiation/finite_diff.jl @@ -0,0 +1,38 @@ + +""" + FiniteDiffBackend <: AbstractDiffBackend + +A type to specify / use differentiation backend based on FiniteDiff package. + +# Constructor + FiniteDiffBackend(method::Val{Symbol} = Val{:central}) +""" +struct FiniteDiffBackend{TM<:Val} <: AbstractDiffBackend + method::TM +end + +FiniteDiffBackend() = FiniteDiffBackend(Val(:central)) + +function _derivative(f, p, ::FiniteDiffBackend{Method}) where {Method} + return FiniteDiff.finite_difference_derivative(f, p, Method) +end + +function _gradient(f, p, ::FiniteDiffBackend{Method}) where {Method} + return FiniteDiff.finite_difference_gradient(f, p, Method) +end + +function _gradient!(f, X, p, ::FiniteDiffBackend{Method}) where {Method} + return FiniteDiff.finite_difference_gradient!(X, f, p, Method) +end + +function _jacobian(f, p, ::FiniteDiffBackend{Method}) where {Method} + return FiniteDiff.finite_difference_jacobian(f, p, Method) +end + +function _jacobian!(f, X, p, ::FiniteDiffBackend{Method}) where {Method} + return FiniteDiff.finite_difference_jacobian!(X, f, p, Method) +end + +if default_differential_backend() === NoneDiffBackend() + set_default_differential_backend!(FiniteDiffBackend()) +end diff --git a/src/differentiation/finite_differences.jl b/src/differentiation/finite_differences.jl new file mode 100644 index 0000000000..6129218c28 --- /dev/null +++ b/src/differentiation/finite_differences.jl @@ -0,0 +1,28 @@ +""" + FiniteDifferencesBackend(method::FiniteDifferenceMethod = central_fdm(5, 1)) + +Differentiation backend based on the FiniteDifferences package. +""" +struct FiniteDifferencesBackend{TM<:FiniteDifferenceMethod} <: AbstractDiffBackend + method::TM +end + +function FiniteDifferencesBackend() + return FiniteDifferencesBackend(central_fdm(5, 1)) +end + +function _derivative(f, t, backend::FiniteDifferencesBackend) + return backend.method(f, t) +end + +function _gradient(f, p, backend::FiniteDifferencesBackend) + return FiniteDifferences.grad(backend.method, f, p)[1] +end + +function _jacobian(f, p, backend::FiniteDifferencesBackend) + return FiniteDifferences.jacobian(backend.method, f, p)[1] +end + +if default_differential_backend() === NoneDiffBackend() + set_default_differential_backend!(FiniteDifferencesBackend()) +end diff --git a/src/forward_diff.jl b/src/differentiation/forward_diff.jl similarity index 81% rename from src/forward_diff.jl rename to src/differentiation/forward_diff.jl index d6aabe3583..4a06b927ef 100644 --- a/src/forward_diff.jl +++ b/src/differentiation/forward_diff.jl @@ -21,4 +21,6 @@ function _jacobian(f, p, ::ForwardDiffBackend) return ForwardDiff.jacobian(f, p) end -push!(_diff_backends, ForwardDiffBackend()) +if default_differential_backend() === NoneDiffBackend() + set_default_differential_backend!(ForwardDiffBackend()) +end diff --git a/src/differentiation/reverse_diff.jl b/src/differentiation/reverse_diff.jl new file mode 100644 index 0000000000..4ec0081c01 --- /dev/null +++ b/src/differentiation/reverse_diff.jl @@ -0,0 +1,13 @@ +struct ReverseDiffBackend <: AbstractDiffBackend end + +function Manifolds._gradient(f, p, ::ReverseDiffBackend) + return ReverseDiff.gradient(f, p) +end + +function Manifolds._gradient!(f, X, p, ::ReverseDiffBackend) + return ReverseDiff.gradient!(X, f, p) +end + +if default_differential_backend() === NoneDiffBackend() + set_default_differential_backend!(ReverseDiffBackend()) +end diff --git a/src/differentiation/riemannian_diff.jl b/src/differentiation/riemannian_diff.jl new file mode 100644 index 0000000000..63395992ae --- /dev/null +++ b/src/differentiation/riemannian_diff.jl @@ -0,0 +1,212 @@ + +""" + AbstractRiemannianDiffBackend + +An abstract type for backends for differentiation. +""" +abstract type AbstractRiemannianDiffBackend end + +@doc raw""" + differential(M::AbstractManifold, f, t::Real, backend::AbstractDiffBackend) + differential!(M::AbstractManifold, f, X, t::Real, backend::AbstractDiffBackend) + +Compute the Riemannian differential of a curve $f: ℝ\to M$ on a manifold `M` +represented by function `f` at time `t` using the given backend. +It is calculated as the tangent vector equal to $\mathrm{d}f_t(t)[1]$. + +The mutating variant computes the differential in place of `X`. +""" +differential(::AbstractManifold, ::Any, ::Real, ::AbstractRiemannianDiffBackend) + +@doc raw""" + gradient(M::AbstractManifold, f, p, backend::AbstractRiemannianDiffBackend) + gradient!(M::AbstractManifold, f, X, p, backend::AbstractRiemannianDiffBackend) + +Compute the Riemannian gradient ``∇f(p)`` of a real-valued function ``f:\mathcal M \to ℝ`` +at at point `p` on the manifold `M` using the specified [`AbstractRiemannianDiffBackend`](@ref). + +The mutating variant computes the gradient in place of `X`. +""" +gradient(::AbstractManifold, ::Any, ::Any, ::AbstractRiemannianDiffBackend) + +function differential!( + M::AbstractManifold, + f::Any, + X, + t, + backend::AbstractRiemannianDiffBackend, +) + return copyto!(X, differential(M, f, t, backend)) +end + +function gradient!(M::AbstractManifold, f, X, p, backend::AbstractRiemannianDiffBackend) + return copyto!(X, gradient(M, f, p, backend)) +end + +@doc raw""" + TangentDiffBackend <: AbstractRiemannianDiffBackend + +A backend that uses a tangent space and a basis therein to derive an +(intrinsic, approximate) differentiation scheme. + +Since it works in a tangent space, methods might require a retraction and an +inverse retraction as well as a basis. + +In the tangent space itself, this backend then employs an (Euclidean) +[`AbstractDiffBackend`](@ref) + +# Constructor + + TangentDiffBackend(diff_backend) + +where `diff_backend` is an [`AbstractDiffBackend`](@ref) to be used on the tangent space. + +With the keyword arguments + +* `retraction` an [`AbstractRetractionMethod`](@ref) ([`ExponentialRetraction`]('ref) by default) +* `inverse_retraction` an [`AbstractInverseRetractionMethod`](@ref) ([`LogarithmicInverseRetraction`]('ref) by default) +* `basis` an [`AbstractBasis`](@ref) ([`DefaultOrthogonalBasis`]('ref) by default) +""" +struct TangentDiffBackend{ + TAD<:AbstractDiffBackend, + TR<:AbstractRetractionMethod, + TIR<:AbstractInverseRetractionMethod, + TB<:AbstractBasis, +} <: AbstractRiemannianDiffBackend + diff_backend::TAD + retraction::TR + inverse_retraction::TIR + basis::TB +end +function TangentDiffBackend( + diff_backend::TAD; + retraction::TR=ExponentialRetraction(), + inverse_retraction::TIR=LogarithmicInverseRetraction(), + basis::TB=DefaultOrthonormalBasis(), +) where { + TAD<:AbstractDiffBackend, + TR<:AbstractRetractionMethod, + TIR<:AbstractInverseRetractionMethod, + TB<:AbstractBasis, +} + return TangentDiffBackend{TAD,TR,TIR,TB}( + diff_backend, + retraction, + inverse_retraction, + basis, + ) +end + +function differential(M::AbstractManifold, f, t::Real, backend::TangentDiffBackend) + p = f(t) + onb_coords = _derivative(zero(number_eltype(p)), backend.diff_backend) do h + return get_coordinates( + M, + p, + inverse_retract(M, p, f(t + h), backend.inverse_retraction), + backend.basis, + ) + end + return get_vector(M, p, onb_coords, backend.basis) +end + +function differential!(M::AbstractManifold, f, X, t::Real, backend::TangentDiffBackend) + p = f(t) + onb_coords = _derivative(zero(number_eltype(p)), backend.diff_backend) do h + return get_coordinates( + M, + p, + inverse_retract(M, p, f(t + h), backend.inverse_retraction), + backend.basis, + ) + end + return get_vector!(M, X, p, onb_coords, backend.basis) +end + +@doc raw""" + gradient(M, f, p, backend::TangentDiffBackend) + +This method uses the internal `backend.diff_backend` (Euclidean) on the function + +```math + f(\retr_p(\cdot)) +``` + +which is given on the tangent space. In detail, the gradient can be written in +terms of the `backend.basis`. We illustrate it here for an [`AbstractOrthonormalBasis`](@ref), +since that simplifies notations: + +```math +\operatorname{grad}f(p) = \operatorname{grad}f(p) = \sum_{i=1}^{d} g_p(\operatorname{grad}f(p),X_i)X_i + = \sum_{i=1}^{d} Df(p)[X_i]X_i +``` + +where the last equality is due to the definition of the gradient as the Riesz representer of the differential. + +If the backend is a forward (or backward) finite difference, these coefficients in the sum can be approximates as + +```math +DF(p)[Y] ≈ \frac{1}{h}\bigl( f(\exp_p(hY)) - f(p) \bigr) +``` +writing ``p=\exp_p(0)`` we see that this is a finite difference of ``f\circ\exp_p``, i.e. for +a function on the tangent space, so we can also use other (Euclidean) backends +""" +function gradient(M::AbstractManifold, f, p, backend::TangentDiffBackend) + X = get_coordinates(M, p, zero_vector(M, p), backend.basis) + onb_coords = _gradient(X, backend.diff_backend) do Y + return f(retract(M, p, get_vector(M, p, Y, backend.basis), backend.retraction)) + end + return get_vector(M, p, onb_coords, backend.basis) +end + +function gradient!(M::AbstractManifold, f, X, p, backend::TangentDiffBackend) + X2 = get_coordinates(M, p, zero_vector(M, p), backend.basis) + onb_coords = _gradient(X2, backend.diff_backend) do Y + return f(retract(M, p, get_vector(M, p, Y, backend.basis), backend.retraction)) + end + return get_vector!(M, X, p, onb_coords, backend.basis) +end + +@doc raw""" + RiemannianProjectionBackend <: AbstractRiemannianDiffBackend + +This backend computes the differentiation in the embedding, which is currently limited +to the gradient. Let ``mathcal M`` denote a manifold embedded in some ``R^m``, where ``m`` +is usually (much) larger than the manifold dimension. +Then we require three tools + +* A function ``f̃: ℝ^m → ℝ`` such that its restriction to the manifold yields the cost + function ``f`` of interest. +* A [`project`](@ref) function to project tangent vectors from the embedding (at `T_pℝ^m``) + back onto the tangent space ``T_p\mathcal M``. This also includes possible changes + of the representation of the tangent vector (e.g. in the Lie algebra or in a different data format). +* A [`change_representer`](@ref) for non-isometrically embedded manifolds, + i.e. where the tangent space ``T_p\mathcal M`` of the manifold does not inherit + the inner product from restriction of the inner product from the tangent space ``T_pℝ^m`` + of the embedding + +For more details see [^AbsilMahonySepulchre2008], Section 3.6.1 for a derivation on submanifolds. + +[^AbsilMahonySepulchre2008]: + > Absil, P.-A., Mahony, R. and Sepulchre R., + > _Optimization Algorithms on Matrix Manifolds_ + > Princeton University Press, 2008, + > doi: [10.1515/9781400830244](https://doi.org/10.1515/9781400830244) + > [open access](http://press.princeton.edu/chapters/absil/) +""" +struct RiemannianProjectionBackend{TADBackend<:AbstractDiffBackend} <: + AbstractRiemannianDiffBackend + diff_backend::TADBackend +end + +function gradient(M::AbstractManifold, f, p, backend::RiemannianProjectionBackend) + amb_grad = _gradient(f, p, backend.diff_backend) + return change_representer(M, EuclideanMetric(), p, project(M, p, amb_grad)) +end + +function gradient!(M::AbstractManifold, f, X, p, backend::RiemannianProjectionBackend) + amb_grad = embed(M, p, X) + _gradient!(f, amb_grad, p, backend.diff_backend) + project!(M, X, p, amb_grad) + return change_representer!(M, X, EuclideanMetric(), p, X) +end diff --git a/src/differentiation/zygote.jl b/src/differentiation/zygote.jl new file mode 100644 index 0000000000..883e027e30 --- /dev/null +++ b/src/differentiation/zygote.jl @@ -0,0 +1,13 @@ +struct ZygoteDiffBackend <: AbstractDiffBackend end + +function Manifolds._gradient(f, p, ::ZygoteDiffBackend) + return Zygote.gradient(f, p)[1] +end + +function Manifolds._gradient!(f, X, p, ::ZygoteDiffBackend) + return copyto!(X, Zygote.gradient(f, p)[1]) +end + +if default_differential_backend() === NoneDiffBackend() + set_default_differential_backend!(ZygoteDiffBackend()) +end diff --git a/src/finite_diff.jl b/src/finite_diff.jl deleted file mode 100644 index 3599522a4d..0000000000 --- a/src/finite_diff.jl +++ /dev/null @@ -1,25 +0,0 @@ - -""" - FiniteDiffBackend(method::Val{Symbol} = Val{:central}) - -Differentiation backend based on FiniteDiff package. -""" -struct FiniteDiffBackend{TM<:Val} <: AbstractDiffBackend - method::TM -end - -FiniteDiffBackend() = FiniteDiffBackend(Val(:central)) - -push!(_diff_backends, FiniteDiffBackend()) - -function _derivative(f, p, backend::FiniteDiffBackend{Method}) where {Method} - return FiniteDiff.finite_difference_derivative(f, p, Method) -end - -function _gradient(f, p, backend::FiniteDiffBackend{Method}) where {Method} - return FiniteDiff.finite_difference_gradient(f, p, Method) -end - -function _gradient!(f, X, p, backend::FiniteDiffBackend{Method}) where {Method} - return FiniteDiff.finite_difference_gradient!(X, f, p, Method) -end diff --git a/src/groups/group.jl b/src/groups/group.jl index 7b74c0edc4..75b6dc6770 100644 --- a/src/groups/group.jl +++ b/src/groups/group.jl @@ -404,7 +404,7 @@ end Compose elements ``p,q ∈ \mathcal{G}`` using the group operation ``p \circ q``. -For implementing composition on a new group manifold, please overload [`_compose`](@ref) +For implementing composition on a new group manifold, please overload `_compose` instead so that methods with [`Identity`](@ref) arguments are not ambiguous. """ compose(::AbstractGroupManifold, ::Any...) diff --git a/src/manifolds/CenteredMatrices.jl b/src/manifolds/CenteredMatrices.jl index 1f6757c2e3..33d922b727 100644 --- a/src/manifolds/CenteredMatrices.jl +++ b/src/manifolds/CenteredMatrices.jl @@ -101,7 +101,7 @@ where $c_i = \frac{1}{m}\sum_{j=1}^m p_{j,i}$ for $i = 1, \dots, n$. """ project(::CenteredMatrices, ::Any) -project!(M::CenteredMatrices, q, p) = copyto!(q, p .- mean(p, dims=1)) +project!(::CenteredMatrices, q, p) = copyto!(q, p .- mean(p, dims=1)) @doc raw""" project(M::CenteredMatrices, p, X) @@ -119,7 +119,7 @@ where $c_i = \frac{1}{m}\sum_{j=1}^m x_{j,i}$ for $i = 1, \dots, n$. """ project(::CenteredMatrices, ::Any, ::Any) -project!(M::CenteredMatrices, Y, p, X) = (Y .= X .- mean(X, dims=1)) +project!(::CenteredMatrices, Y, p, X) = (Y .= X .- mean(X, dims=1)) @generated representation_size(::CenteredMatrices{m,n,𝔽}) where {m,n,𝔽} = (m, n) diff --git a/src/manifolds/ConnectionManifold.jl b/src/manifolds/ConnectionManifold.jl index 38dd212a28..7ee5d4f876 100644 --- a/src/manifolds/ConnectionManifold.jl +++ b/src/manifolds/ConnectionManifold.jl @@ -59,7 +59,7 @@ end M::AbstractManifold, p, B::AbstractBasis; - backend::AbstractDiffBackend = diff_backend(), + backend::AbstractDiffBackend = default_differential_backend(), ) Compute the Christoffel symbols of the second kind in local coordinates of basis `B`. @@ -89,7 +89,7 @@ christoffel_symbols_second(::AbstractManifold, ::Any, ::AbstractBasis) M::AbstractManifold, p, B::AbstractBasis; - backend::AbstractDiffBackend = diff_backend(), + backend::AbstractDiffBackend = default_differential_backend(), ) Get partial derivatives of the Christoffel symbols of the second kind @@ -106,15 +106,25 @@ function christoffel_symbols_second_jacobian( M::AbstractManifold, p, B::AbstractBasis; - backend::AbstractDiffBackend=diff_backend(), + backend::AbstractDiffBackend=default_differential_backend(), + retraction::AbstractRetractionMethod=ManifoldsBase.default_retraction_method(M), ) - n = size(p, 1) + d = manifold_dimension(M) ∂Γ = reshape( - _jacobian(q -> christoffel_symbols_second(M, q, B; backend=backend), p, backend), - n, - n, - n, - n, + _jacobian( + c -> christoffel_symbols_second( + M, + retract(M, p, get_vector(M, p, c, B), retraction), + B; + backend=backend, + ), + p, + backend, + ), + d, + d, + d, + d, ) return ∂Γ end @@ -162,7 +172,7 @@ exp(::AbstractConnectionManifold, ::Any...) end """ - gaussian_curvature(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend()) + gaussian_curvature(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = default_differential_backend()) Compute the Gaussian curvature of the manifold `M` at the point `p` using basis `B`. This is equal to half of the scalar Ricci curvature, see [`ricci_curvature`](@ref). @@ -195,7 +205,7 @@ function injectivity_radius(M::AbstractConnectionManifold, p, m::ExponentialRetr end """ - ricci_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend()) + ricci_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = default_differential_backend()) Compute the Ricci tensor, also known as the Ricci curvature tensor, of the manifold `M` at the point `p` using basis `B`, @@ -217,7 +227,7 @@ end ) @doc raw""" - riemann_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend=diff_backend()) + riemann_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend=default_differential_backend()) Compute the Riemann tensor ``R^l_{ijk}``, also known as the Riemann curvature tensor, at the point `p` in local coordinates defined by `B`. The dimensions of the @@ -236,7 +246,7 @@ function riemann_tensor( M::AbstractManifold, p, B::AbstractBasis; - backend::AbstractDiffBackend=diff_backend(), + backend::AbstractDiffBackend=default_differential_backend(), ) n = size(p, 1) Γ = christoffel_symbols_second(M, p, B; backend=backend) @@ -260,7 +270,7 @@ end X, tspan, B::AbstractBasis; - backend::AbstractDiffBackend = diff_backend(), + backend::AbstractDiffBackend = default_differential_backend(), solver = AutoVern9(Rodas5()), kwargs..., ) diff --git a/src/manifolds/GeneralizedGrassmann.jl b/src/manifolds/GeneralizedGrassmann.jl index cd46d6b43b..41112e5dda 100644 --- a/src/manifolds/GeneralizedGrassmann.jl +++ b/src/manifolds/GeneralizedGrassmann.jl @@ -57,6 +57,39 @@ function GeneralizedGrassmann( return GeneralizedGrassmann{n,k,field,typeof(B)}(B) end +@doc raw""" + change_representer(M::GeneralizedGrassmann, ::EuclideanMetric, p, X) + +Change `X` to the corresponding representer of a cotangent vector at `p` with respect to the scaled metric +of the [`GeneralizedGrassmann`](@ref) `M`, i.e, since + +```math +g_p(X,Y) = \operatorname{tr}(Y^{\mathrm{H}}BZ) = \operatorname{tr}(X^{\mathrm{H}}Z) = ⟨X,Z⟩ +``` + +has to hold for all ``Z``, where the repreenter `X` is given, the resulting representer with +respect to the metric on the [`GeneralizedGrassmann`](@ref) is given by ``Y = B^{-1}X``. +""" +change_representer(::GeneralizedGrassmann, ::EuclideanMetric, ::Any, ::Any) + +function change_representer!(M::GeneralizedGrassmann, Y, ::EuclideanMetric, p, X) + return copyto!(M, Y, p, M.B \ X) +end + +@doc raw""" + change_metric(M::GeneralizedGrassmann, ::EuclideanMetric, p X) + +Change `X` to the corresponding vector with respect to the metric of the [`GeneralizedGrassmann`](@ref) `M`, +i.e. let ``B=LL'`` be the Cholesky decomposition of the matrix `M.B`, then the corresponding vector is ``L\X``. +""" +change_metric(M::GeneralizedGrassmann, ::EuclideanMetric, ::Any, ::Any) + +function change_metric!(M::GeneralizedGrassmann, Y, ::EuclideanMetric, p, X) + C2 = cholesky(M.B).L + Y .= C2 \ X + return Y +end + @doc raw""" check_point(M::GeneralizedGrassmann{n,k,𝔽}, p) diff --git a/src/manifolds/Hyperbolic.jl b/src/manifolds/Hyperbolic.jl index 8e44b43ad8..a9e5af25f3 100644 --- a/src/manifolds/Hyperbolic.jl +++ b/src/manifolds/Hyperbolic.jl @@ -133,7 +133,10 @@ for T in _HyperbolicTypes allocate(p::$T, ::Type{P}, dims::Tuple) where {P} = $T(allocate(p.value, P, dims)) @inline Base.copy(p::$T) = $T(copy(p.value)) - Base.copyto!(q::$T, p::$T) = copyto!(q.value, p.value) + function Base.copyto!(q::$T, p::$T) + copyto!(q.value, p.value) + return q + end Base.similar(p::$T) = $T(similar(p.value)) diff --git a/src/manifolds/HyperbolicHyperboloid.jl b/src/manifolds/HyperbolicHyperboloid.jl index e1411c99a6..798545c43f 100644 --- a/src/manifolds/HyperbolicHyperboloid.jl +++ b/src/manifolds/HyperbolicHyperboloid.jl @@ -1,3 +1,29 @@ +@doc raw""" + change_representer(M::Hyperbolic{n}, ::EuclideanMetric, p, X) + +Change the Eucliden representer `X` of a cotangent vector at point `p`. +We only have to correct for the metric, which means that the sign of the last entry changes, since +for the result ``Y`` we are looking for a tangent vector such that + +```math + g_p(Y,Z) = -y_{n+1}z_{n+1} + \sum_{i=1}^n y_iz_i = \sum_{i=1}^{n+1} z_ix_i +``` + +holds, which directly yields ``y_i=x_i`` for ``i=1,\ldots,n`` and ``y_{n+1}=-x_{n+1}``. +""" +change_representer(::Hyperbolic, ::EuclideanMetric, ::Any, ::Any) + +function change_representer!(M::Hyperbolic, Y, ::EuclideanMetric, p, X) + copyto!(M, Y, p, X) + Y[end] *= -1 + return Y +end + +function change_metric!(::Hyperbolic, ::Any, ::EuclideanMetric, ::Any, ::Any) + return error( + "Changing metric from Euclidean to Hyperbolic is not possible (see Sylvester's law of inertia).", + ) +end function check_point(M::Hyperbolic, p; kwargs...) mpv = invoke(check_point, Tuple{supertype(typeof(M)),typeof(p)}, M, p; kwargs...) @@ -325,7 +351,7 @@ component such that for the resulting `p` we have that its [`minkowski_metric`](@ref) is $⟨p,X⟩_{\mathrm{M}} = 0$, i.e. $X_{n+1} = \frac{⟨\tilde p, Y⟩}{p_{n+1}}$, where $\tilde p = (p_1,\ldots,p_n)$. """ -_hyperbolize(M::Hyperbolic, p, Y) = vcat(Y, dot(p[1:(end - 1)], Y) / p[end]) +_hyperbolize(::Hyperbolic, p, Y) = vcat(Y, dot(p[1:(end - 1)], Y) / p[end]) @doc raw""" inner(M::Hyperbolic{n}, p, X, Y) diff --git a/src/manifolds/HyperbolicPoincareBall.jl b/src/manifolds/HyperbolicPoincareBall.jl index 651e208648..0499b02f2a 100644 --- a/src/manifolds/HyperbolicPoincareBall.jl +++ b/src/manifolds/HyperbolicPoincareBall.jl @@ -1,3 +1,48 @@ +@doc raw""" + change_representer(M::Hyperbolic{n}, ::EuclideanMetric, p::PoincareBallPoint, X::PoincareBallTVector) + +Since in the metric we have the term `` α = \frac{2}{1-\sum_{i=1}^n p_i^2}`` per element, +the correction for the gradient reads `` Y = \frac{1}{α^2}X``. +""" +change_representer( + ::Hyperbolic, + ::EuclideanMetric, + ::PoincareBallPoint, + ::PoincareBallTVector, +) + +function change_representer!( + ::Hyperbolic, + Y::PoincareBallTVector, + ::EuclideanMetric, + p::PoincareBallPoint, + X::PoincareBallTVector, +) + α = 2 / (1 - norm(p.value)^2) + Y.value .= X.value ./ α^2 + return Y +end + +@doc raw""" + change_metric(M::Hyperbolic{n}, ::EuclideanMetric, p::PoincareBallPoint, X::PoincareBallTVector) + +Since in the metric we always have the term `` α = \frac{2}{1-\sum_{i=1}^n p_i^2}`` per element, +the correction for the metric reads ``Z = \frac{1}{α}X``. +""" +change_metric(::Hyperbolic, ::EuclideanMetric, ::PoincareBallPoint, ::PoincareBallTVector) + +function change_metric!( + ::Hyperbolic, + Y::PoincareBallTVector, + ::EuclideanMetric, + p::PoincareBallPoint, + X::PoincareBallTVector, +) + α = 2 / (1 - norm(p.value)^2) + Y.value .= X.value ./ α + return Y +end + function check_point(M::Hyperbolic{N}, p::PoincareBallPoint; kwargs...) where {N} mpv = check_point(Euclidean(N), p.value; kwargs...) mpv === nothing || return mpv diff --git a/src/manifolds/MetricManifold.jl b/src/manifolds/MetricManifold.jl index 066b00760c..f3c91073ba 100644 --- a/src/manifolds/MetricManifold.jl +++ b/src/manifolds/MetricManifold.jl @@ -6,11 +6,13 @@ varying inner products on the tangent space. See [`inner`](@ref). # Functor - (metric::Metric)(M::Manifold) + (metric::Metric)(M::AbstractManifold) + (metric::Metric)(M::MetricManifold) Generate the `MetricManifold` that wraps the manifold `M` with given `metric`. This works for both a variable containing the metric as well as a subtype `T<:AbstractMetric`, where a zero parameter constructor `T()` is availabe. +If `M` is already a metric manifold, the inner manifold with the new `metric` is returned. """ abstract type AbstractMetric end @@ -41,6 +43,9 @@ struct MetricManifold{𝔽,M<:AbstractManifold{𝔽},G<:AbstractMetric} <: manifold::M metric::G end +# remetricise instead of double-decorating +(metric::AbstractMetric)(M::MetricManifold) = MetricManifold(M.manifold, metric) +(::Type{T})(M::MetricManifold) where {T<:AbstractMetric} = MetricManifold(M.manifold, T()) @doc raw""" RiemannianMetric <: AbstractMetric @@ -51,20 +56,159 @@ inner product ``g(X, X) > 0`` whenever ``X`` is not the zero vector. """ abstract type RiemannianMetric <: AbstractMetric end +@doc raw""" + change_metric(M::AbstractcManifold, G2::AbstractMetric, p, X) + +On the [`AbstractManifold`](@ref) `M` with implicitly given metric ``g_1`` +and a second [`AbstractMetric`](@ref) ``g_2`` this function performs a change of metric in the +sense that it returns the tangent vector ``Z=BX`` such that the linear map ``B`` fulfills + +````math +g_2(Y_1,Y_2) = g_1(BY_1,BY_2) \quad \text{for all } Y_1, Y_2 ∈ T_p\mathcal M. +```` + +If both metrics are given in their [`local_metric`](@ref) (symmetric positive defintie) matrix +representations ``G_1 = C_1C_1^{\mathrm{H}}`` and ``G_2 = C_2C_2^{\mathrm{H}}``, where ``C_1,C_2`` denote their respective +Cholesky factors, then solving ``C_2C_2^{\mathrm{H}} = G_2 = B^{\mathrm{H}}G_1B = B^{\mathrm{H}}C_1C_1^{\mathrm{H}}B`` yields ``B = (C_1 \backslash C_2)^{\mathrm{H}}``, +where ``\cdot^{\mathrm{H}}`` denotes the conjugate transpose. + +This function returns `Z = BX`. + +# Examples + + change_metric(Sphere(2), EuclideanMetric(), p, X) + +Since the metric in ``T_p\mathbb S^2`` is the Euclidean metric from the embedding restricted to ``T_p\mathbb S^2``, this just returns `X` + + change_metric(SymmetricPOsitiveDefinite(3), EuclideanMetric, p, X) + +Here, the default metric in ``\mathcal P(3)`` is the [`LinearAffineMetric`](@ref) and the transformation can be computed as ``B=p`` +""" +change_metric(::AbstractManifold, ::AbstractMetric, ::Any, ::Any) + +function change_metric(M::AbstractManifold, G::AbstractMetric, p, X) + Y = allocate_result(M, change_metric, X, p) # this way we allocate a tangent + return change_metric!(M, Y, G, p, X) +end +function change_metric!(M::AbstractManifold, Y, G::AbstractMetric, p, X) + is_default_metric(M, G) && return copyto!(M, Y, p, X) + M.metric === G && return copyto!(M, Y, p, X) # no metric change + # TODO: For local metric, inverse_local metric, det_local_metric: Introduce a default basis? + B = DefaultOrthogonalBasis() + G1 = local_metric(M, p, B) + G2 = local_metric(G(M), p, B) + x = get_coordinates(M, p, X, B) + C1 = cholesky(G1).L + C2 = cholesky(G2).L + z = (C1 \ C2)'x + return get_vector!(M, Y, p, z, B) +end + +@decorator_transparent_signature change_metric( + M::AbstractDecoratorManifold, + G::AbstractMetric, + X, + p, +) +@decorator_transparent_signature change_metric!( + M::AbstractDecoratorManifold, + Y, + G::AbstractMetric, + X, + p, +) + +@doc raw""" + change_representer(M::AbstractManifold, G2::AbstractMetric, p, X) + +Convert the representer `X` of a linear function (in other words a cotangent vector at `p`) +in the tangent space at `p` on the [`AbstractManifold`](@ref) `M` given with respect to the +[`AbstractMetric`](@ref) `G2` into the representer with respect to the (implicit) metric of `M`. + +In order to convert `X` into the representer with respect to the (implicitly given) metric ``g_1`` of `M`, +we have to find the conversion function ``c: T_p\mathcal M \to T_p\mathcal M`` such that + +```math + g_2(X,Y) = g_1(c(X),Y) +``` + +If both metrics are given in their [`local_metric`](@ref) (symmetric positive defintie) matrix +representations ``G_1`` and ``G_2`` and ``x,y`` are the local coordinates with respect to +the same basis of the tangent space, the equation reads + +```math + x^{\mathrm{H}}G_2y = c(x)^{\mathrm{H}}G_1 y \quad \text{for all } y \in ℝ^d, +``` +where ``\cdot^{\mathrm{H}}`` denotes the conjugate transpose. +We obtain ``c(X) = (G_1\backslash G_2)^{\mathrm{H}}X`` + +For example `X` could be the gradient ``\operatorname{grad}f`` of a real-valued function +``f: \mathcal M \to ℝ``, i.e. + +```math + g_2(X,Y) = Df(p)[Y] \quad \text{for all } Y ∈ T_p\mathcal M. +``` + +and we would change the Riesz representer `X` to the representer with respect to the metric ``g_1``. + +# Examples + + change_representer(Sphere(2), EuclideanMetric(), p, X) + +Since the metric in ``T_p\mathbb S^2`` is the Euclidean metric from the embedding restricted to ``T_p\mathbb S^2``, this just returns `X` + + change_representer(SymmetricPositiveDefinite(3), EuclideanMetric(), p, X) + +Here, the default metric in ``\mathcal P(3)`` is the [`LinearAffineMetric`](@ref) and the transformation can be computed as ``pXp`` +""" +change_representer(::AbstractManifold, ::AbstractMetric, ::Any, ::Any) + +function change_representer(M::AbstractManifold, G::AbstractMetric, p, X) + Y = allocate_result(M, change_representer, X, p) # this way we allocate a tangent + return change_representer!(M, Y, G, p, X) +end + +@decorator_transparent_signature change_representer( + M::AbstractDecoratorManifold, + G::AbstractMetric, + X, + p, +) +@decorator_transparent_signature change_representer!( + M::AbstractDecoratorManifold, + Y, + G::AbstractMetric, + X, + p, +) + +# Default fallback I: compute in local metric representations +function change_representer!(M::AbstractManifold, Y, G::AbstractMetric, p, X) + is_default_metric(M, G) && return copyto!(M, Y, p, X) + M.metric === G && return copyto!(M, Y, p, X) # no metric change + # TODO: For local metric, inverse_local metric, det_local_metric: Introduce a default basis? + B = DefaultOrthogonalBasis() + G1 = local_metric(M, p, B) + G2 = local_metric(G(M), p, B) + x = get_coordinates(M, p, X, B) + z = (G1 \ G2)'x + return get_vector!(M, Y, p, z, B) +end + @doc raw""" christoffel_symbols_first( M::MetricManifold, p, B::AbstractBasis; - backend::AbstractDiffBackend = diff_backend(), + backend::AbstractDiffBackend = default_differential_backend(), ) Compute the Christoffel symbols of the first kind in local coordinates of basis `B`. The Christoffel symbols are (in Einstein summation convention) -```math +````math Γ_{ijk} = \frac{1}{2} \Bigl[g_{kj,i} + g_{ik,j} - g_{ij,k}\Bigr], -``` +```` where ``g_{ij,k}=\frac{∂}{∂ p^k} g_{ij}`` is the coordinate derivative of the local representation of the metric tensor. The dimensions of @@ -75,7 +219,7 @@ function christoffel_symbols_first( M::AbstractManifold, p, B::AbstractBasis; - backend::AbstractDiffBackend=diff_backend(), + backend::AbstractDiffBackend=default_differential_backend(), ) ∂g = local_metric_jacobian(M, p, B; backend=backend) n = size(∂g, 1) @@ -94,7 +238,7 @@ function christoffel_symbols_second( M::AbstractManifold, p, B::AbstractBasis; - backend::AbstractDiffBackend=diff_backend(), + backend::AbstractDiffBackend=default_differential_backend(), ) Ginv = inverse_local_metric(M, p, B) Γ₁ = christoffel_symbols_first(M, p, B; backend=backend) @@ -128,7 +272,7 @@ end B::AbstractBasis, ) """ - einstein_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend()) + einstein_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_badefault_differential_backendckend()) Compute the Einstein tensor of the manifold `M` at the point `p`, see [https://en.wikipedia.org/wiki/Einstein_tensor](https://en.wikipedia.org/wiki/Einstein_tensor) """ @@ -137,7 +281,7 @@ function einstein_tensor( M::AbstractManifold, p, B::AbstractBasis; - backend::AbstractDiffBackend=diff_backend(), + backend::AbstractDiffBackend=default_differential_backend(), ) Ric = ricci_tensor(M, p, B; backend=backend) g = local_metric(M, p, B) @@ -179,12 +323,16 @@ flat(::MetricManifold, ::Any...) end @doc raw""" - inverse_local_metric(M::AbstractcManifold, p, B::AbstractBasis) + inverse_local_metric(M::AbstractcManifold{𝔽}, p, B::AbstractBasis) -Return the local matrix representation of the inverse metric (cometric) tensor, usually -written ``g^{ij}``. +Return the local matrix representation of the inverse metric (cometric) tensor +of the tangent space at `p` on the [`AbstractManifold`](@ref) `M` with respect +to the [`AbstractBasis`](@ref) basis `B`. -See also [`local_metric`](@ref) +The metric tensor (see [`local_metric`](@ref)) is usually denoted by ``G = (g_{ij}) ∈ 𝔽^{d×d}``, +where ``d`` is the dimension of the manifold. + +Then the inverse local metric is denoted by ``G^{-1} = g^{ij}``. """ inverse_local_metric(::AbstractManifold, ::Any, ::AbstractBasis) function inverse_local_metric(M::AbstractManifold, p, B::AbstractBasis) @@ -280,13 +428,16 @@ inner(::MetricManifold, ::Any, ::Any, ::Any) end @doc raw""" - local_metric(M::AbstractManifold, p, B::AbstractBasis) + local_metric(M::AbstractManifold{𝔽}, p, B::AbstractBasis) Return the local matrix representation at the point `p` of the metric tensor ``g`` with -respect to the [`AbstractBasis`](@ref) `B` on the [`AbstractManifold`](@ref) `M`, usually written ``g_{ij}``. -The matrix has the property that ``g(X, Y)=X^\mathrm{T} [g_{ij}] Y = g_{ij} X^i Y^j``, -where the latter expression uses Einstein summation convention. -The metric tensor is such that the formula works for the given [`AbstractBasis`](@ref) `B`. +respect to the [`AbstractBasis`](@ref) `B` on the [`AbstractManifold`](@ref) `M`. +Let ``d``denote the dimension of the manifold and $b_1,\ldots,b_d$ the basis vectors. +Then the local matrix representation is a matrix ``G\in 𝔽^{n\times n}`` whose entries are +given by ``g_{ij} = g_p(b_i,b_j), i,j\in\{1,…,d\}``. + +This yields the property for two tangent vectors (using Einstein summation convention) +``X = X^ib_i, Y=Y^ib_i \in T_p\mathcal M`` we get ``g_p(X, Y) = g_{ij} X^i Y^j``. """ local_metric(::AbstractManifold, ::Any, ::AbstractBasis) @decorator_transparent_signature local_metric( @@ -301,22 +452,32 @@ local_metric(::AbstractManifold, ::Any, ::AbstractBasis) M::AbstractManifold, p, B::AbstractBasis; - backend::AbstractDiffBackend = diff_backend(), + backend::AbstractDiffBackend, ) Get partial derivatives of the local metric of `M` at `p` in basis `B` with respect to the coordinates of `p`, ``\frac{∂}{∂ p^k} g_{ij} = g_{ij,k}``. The dimensions of the resulting multi-dimensional array are ordered ``(i,j,k)``. """ -local_metric_jacobian(::AbstractManifold, ::Any, B::AbstractBasis) +local_metric_jacobian(::AbstractManifold, ::Any, B::AbstractBasis, ::AbstractDiffBackend) function local_metric_jacobian( M::AbstractManifold, p, B::AbstractBasis; - backend::AbstractDiffBackend=diff_backend(), + backend::AbstractDiffBackend=default_differential_backend(), + retraction::AbstractRetractionMethod=ManifoldsBase.default_retraction_method(M), ) - n = size(p, 1) - ∂g = reshape(_jacobian(q -> local_metric(M, q, B), p, backend), n, n, n) + d = manifold_dimension(M) + ∂g = reshape( + _jacobian( + c -> local_metric(M, retract(M, p, get_vector(M, p, c, B), retraction), B), + p, + backend, + ), + d, + d, + d, + ) return ∂g end @decorator_transparent_signature local_metric_jacobian( @@ -365,7 +526,7 @@ function metric(M::MetricManifold) end @doc raw""" - ricci_curvature(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend()) + ricci_curvature(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = default_differential_backend()) Compute the Ricci scalar curvature of the manifold `M` at the point `p` using basis `B`. The curvature is computed as the trace of the Ricci curvature tensor with respect to @@ -380,7 +541,7 @@ function ricci_curvature( M::AbstractManifold, p, B::AbstractBasis; - backend::AbstractDiffBackend=diff_backend(), + backend::AbstractDiffBackend=default_differential_backend(), ) Ginv = inverse_local_metric(M, p, B) Ric = ricci_tensor(M, p, B; backend=backend) @@ -436,7 +597,7 @@ for f in [ quote function decorator_transparent_dispatch( ::typeof($f), - ::AbstractConnectionManifold, + M::AbstractConnectionManifold, args..., ) return Val(:parent) @@ -445,6 +606,19 @@ for f in [ ) end +for f in [change_metric, change_representer, change_metric!, change_representer!] + eval( + quote + function decorator_transparent_dispatch( + ::typeof($f), + ::AbstractManifold, + args..., + ) + return Val(:parent) + end + end, + ) +end function decorator_transparent_dispatch( ::typeof(christoffel_symbols_second), ::MetricManifold, diff --git a/src/manifolds/Multinomial.jl b/src/manifolds/Multinomial.jl index 36f86eef64..8db4e6e63d 100644 --- a/src/manifolds/Multinomial.jl +++ b/src/manifolds/Multinomial.jl @@ -53,6 +53,7 @@ function check_point(M::MultinomialMatrices{n,m}, p; kwargs...) where {n,m} end return check_point(PowerManifold(M.manifold, m), p; kwargs...) end + @doc raw""" check_vector(M::MultinomialMatrices p, X; kwargs...) diff --git a/src/manifolds/PositiveNumbers.jl b/src/manifolds/PositiveNumbers.jl index 05f16c3c05..c3b22362fa 100644 --- a/src/manifolds/PositiveNumbers.jl +++ b/src/manifolds/PositiveNumbers.jl @@ -38,6 +38,55 @@ This manifold is modeled as a [`PowerManifold`](@ref) of [`PositiveNumbers`](@re """ PositiveArrays(n::Vararg{Int,I}) where {I} = PositiveNumbers()^(n) +@doc raw""" + change_representer(M::PositiveNumbers, E::EuclideanMetric, p, X) + +Given a tangent vector ``X ∈ T_p\mathcal M`` representing a linear function with respect +to the [`EuclideanMetric`](@ref) `g_E`, this function changes the representer into the one +with respect to the positivity metric representation of +[`PositiveNumbers`](@ref) `M`. + +For all tangent vectors ``Y`` the result ``Z`` has to fulfill + +```math + ⟨X,Y⟩ = XY = \frac{ZY}{p^2} = g_p(YZ) +``` + +and hence ``Z = p^2X`` + +""" +change_representer(::PositiveNumbers, ::EuclideanMetric, ::Any, ::Any) +change_representer(::PositiveNumbers, ::EuclideanMetric, p::Real, X::Real) = p * X * p + +function change_representer!(::PositiveNumbers, Y, ::EuclideanMetric, p, X) + Y .= p .* X .* p + return Y +end + +@doc raw""" + change_metric(M::PositiveNumbers, E::EuclideanMetric, p, X) + +Given a tangent vector ``X ∈ T_p\mathcal M`` representing a linear function with respect to +the [`EuclideanMetric`](@ref) `g_E`, +this function changes the representer into the one with respect to the positivity metric +of [`PositiveNumbers`](@ref) `M`. + +For all ``Z,Y`` we are looking for the function ``c`` on the tangent space at ``p`` such that + +```math + ⟨Z,Y⟩ = XY = \frac{c(Z)c(Y)}{p^2} = g_p(c(Y),c(Z)) +``` + +and hence ``C(X) = pX``. +""" +change_metric(::PositiveNumbers, ::EuclideanMetric, ::Any, ::Any) +change_metric(::PositiveNumbers, ::EuclideanMetric, p::Real, X::Real) = p * X + +function change_metric!(::PositiveNumbers, Y, ::EuclideanMetric, p, X) + Y .= p .* X + return Y +end + @doc raw""" check_point(M::PositiveNumbers, p) diff --git a/src/manifolds/PowerManifold.jl b/src/manifolds/PowerManifold.jl index 833a952d04..e5cc8c6ee8 100644 --- a/src/manifolds/PowerManifold.jl +++ b/src/manifolds/PowerManifold.jl @@ -78,6 +78,48 @@ end default_metric_dispatch(::AbstractPowerManifold, ::PowerMetric) = Val(true) +""" + change_representer(M::AbstractPowerManifold, ::AbstractMetric, p, X) + +Since the metric on a power manifold decouples, the change of a representer can be done elementwise +""" +change_representer(::AbstractPowerManifold, ::AbstractMetric, ::Any, ::Any) + +function change_representer!(M::AbstractPowerManifold, Y, G::AbstractMetric, p, X) + rep_size = representation_size(M.manifold) + for i in get_iterator(M) + change_representer!( + M.manifold, + _write(M, rep_size, Y, i), + G, + _read(M, rep_size, p, i), + _read(M, rep_size, X, i), + ) + end + return Y +end + +""" + change_metric(M::AbstractPowerManifold, ::AbstractMetric, p, X) + +Since the metric on a power manifold decouples, the change of metric can be done elementwise. +""" +change_metric(M::AbstractPowerManifold, ::AbstractMetric, ::Any, ::Any) + +function change_metric!(M::AbstractPowerManifold, Y, G::AbstractMetric, p, X) + rep_size = representation_size(M.manifold) + for i in get_iterator(M) + change_metric!( + M.manifold, + _write(M, rep_size, Y, i), + G, + _read(M, rep_size, p, i), + _read(M, rep_size, X, i), + ) + end + return Y +end + @doc raw""" flat(M::AbstractPowerManifold, p, X) diff --git a/src/manifolds/ProbabilitySimplex.jl b/src/manifolds/ProbabilitySimplex.jl index 557c344c25..c66aebff99 100644 --- a/src/manifolds/ProbabilitySimplex.jl +++ b/src/manifolds/ProbabilitySimplex.jl @@ -58,6 +58,44 @@ See for example the [`ProbabilitySimplex`](@ref). """ struct FisherRaoMetric <: AbstractMetric end +@doc raw""" + change_representer(M::ProbabilitySimplex, ::EuclideanMetric, p, X) + +Given a tangent vector with respect to the metric from the embedding, the [`EuclideanMetric`](@ref), +the representer of a linear functional on the tangent space is adapted as ``Z = p .* X``, since +this “compensates” for the divsion by ``p`` in the Riemannian metric on the [`ProbabilitySimplex`](@ref). + +To be precise for any ``Y ∈ T_pΔ^n`` we are looking for ``Z ∈ T_pΔ^n`` such that + +```math + ⟨X,Y⟩ = X^\mathrm{T}Y = \sum_{i=1}^{n+1}\frac{Z_iY_i}{p_i} = g_p(Z,Y) +``` + +and hence ``Z_i = X_ip_i, i=1,…,n+1``. +""" +change_representer(::ProbabilitySimplex, ::EuclideanMetric, ::Any, ::Any) + +function change_representer!(::ProbabilitySimplex, Y, ::EuclideanMetric, p, X) + return Y .= p .* X +end + +@doc raw""" + change_metric(M::ProbabilitySimplex, ::EuclideanMetric, p, X) + +To change the metric, we are looking for a function ``c\colon T_pΔ^n \to T_pΔ^n`` such that for all ``X,Y ∈ T_pΔ^n`` + +```math + ⟨X,Y⟩ = X^\mathrm{T}Y = \sum_{i=1}^{n+1}\frac{c(X)_ic(Y)_i}{p_i} = g_p(X,Y) +``` + +and hence ``C(X)_i = X_i\sqrt{p_i}, i=1,…,n+1``. +""" +change_metric(::ProbabilitySimplex, ::EuclideanMetric, ::Any, ::Any) + +function change_metric!(::ProbabilitySimplex, Y, ::EuclideanMetric, p, X) + return Y .= sqrt.(p) .* X +end + """ check_point(M::ProbabilitySimplex, p; kwargs...) diff --git a/src/manifolds/ProductManifold.jl b/src/manifolds/ProductManifold.jl index 4aa0e60833..aa26540412 100644 --- a/src/manifolds/ProductManifold.jl +++ b/src/manifolds/ProductManifold.jl @@ -135,6 +135,42 @@ function ProductVectorTransport(methods::AbstractVectorTransportMethod...) return ProductVectorTransport{typeof(methods)}(methods) end +""" + change_representer(M::ProductManifold, ::AbstractMetric, p, X) + +Since the metric on a product manifold decouples, the change of a representer can be done elementwise +""" +change_representer(::ProductManifold, ::AbstractMetric, ::Any, ::Any) + +function change_representer!(M::ProductManifold, Y, G::AbstractMetric, p, X) + map( + (m, y, P, x) -> change_representer!(m, y, G, P, x), + M.manifolds, + submanifold_components(M, Y), + submanifold_components(M, p), + submanifold_components(M, X), + ) + return Y +end + +""" + change_metric(M::ProductManifold, ::AbstractMetric, p, X) + +Since the metric on a product manifold decouples, the change of metric can be done elementwise. +""" +change_metric(::ProductManifold, ::AbstractMetric, ::Any, ::Any) + +function change_metric!(M::ProductManifold, Y, G::AbstractMetric, p, X) + map( + (m, y, P, x) -> change_metric!(m, y, G, P, x), + M.manifolds, + submanifold_components(M, Y), + submanifold_components(M, p), + submanifold_components(M, X), + ) + return Y +end + """ check_point(M::ProductManifold, p; kwargs...) diff --git a/src/manifolds/Spectrahedron.jl b/src/manifolds/Spectrahedron.jl index 53289fccc8..efd14f7f12 100644 --- a/src/manifolds/Spectrahedron.jl +++ b/src/manifolds/Spectrahedron.jl @@ -138,7 +138,7 @@ project!(::Spectrahedron, r, q) = copyto!(r, q ./ norm(q)) """ project(M::Spectrahedron, q, Y) -Project `Y` onto the tangent space at `q`, i.e. row-wise onto the oblique manifold. +Project `Y` onto the tangent space at `q`, i.e. row-wise onto the Spectrahedron manifold. """ project(::Spectrahedron, ::Any...) diff --git a/src/manifolds/Sphere.jl b/src/manifolds/Sphere.jl index 594e05bc6c..92f07d1a25 100644 --- a/src/manifolds/Sphere.jl +++ b/src/manifolds/Sphere.jl @@ -144,7 +144,9 @@ end function decorated_manifold(M::AbstractSphere{𝔽}) where {𝔽} return Euclidean(representation_size(M)...; field=𝔽) end -get_embedding(M::AbstractSphere{𝔽}) where {𝔽} = decorated_manifold(M) + +# Since on every tangent space the Euclidean matric (restricted to this space) is used, this should be fine +default_metric_dispatch(::AbstractSphere, ::EuclideanMetric) = Val(true) @doc raw""" distance(M::AbstractSphere, p, q) @@ -227,6 +229,8 @@ function get_coordinates!( return Y end +get_embedding(M::AbstractSphere{𝔽}) where {𝔽} = decorated_manifold(M) + @doc raw""" get_vector(M::AbstractSphere{ℝ}, p, X, B::DefaultOrthonormalBasis) diff --git a/src/manifolds/SymmetricPositiveDefiniteLinearAffine.jl b/src/manifolds/SymmetricPositiveDefiniteLinearAffine.jl index 237dce32ce..e072a86676 100644 --- a/src/manifolds/SymmetricPositiveDefiniteLinearAffine.jl +++ b/src/manifolds/SymmetricPositiveDefiniteLinearAffine.jl @@ -6,6 +6,53 @@ matrix logarithms and exponentials, which yields a linear and affine metric. """ struct LinearAffineMetric <: RiemannianMetric end +@doc raw""" + change_representer(M::SymmetricPositiveDefinite, E::EuclideanMetric, p, X) + +Given a tangent vector ``X ∈ T_p\mathcal M`` representing a linear function on the tangent +space at `p` with respect to the [`EuclideanMetric`](@ref) `g_E`, +this is turned into the representer with respect to the (default) metric, +the [`LinearAffineMetric`](@ref) on the [`SymmetricPositiveDefinite`](@ref) `M`. + +To be precise we are looking for ``Z∈T_p\mathcal P(n)`` such that for all ``Y∈T_p\mathcal P(n)``` +it holds + +```math +⟨X,Y⟩ = \operatorname{tr}(XY) = \operatorname{tr}(p^{-1}Zp^{-1}Y) = g_p(Z,Y) +``` + +and hence ``Z = pXp``. +""" +change_representer(::SymmetricPositiveDefinite, ::EuclideanMetric, ::Any, ::Any) + +function change_representer!(::SymmetricPositiveDefinite, Y, ::EuclideanMetric, p, X) + Y .= p * X * p + return Y +end + +@doc raw""" + change_metric(M::SymmetricPositiveDefinite{n}, E::EuclideanMetric, p, X) + +Given a tangent vector ``X ∈ T_p\mathcal P(n)`` with respect to the [`EuclideanMetric`](@ref) `g_E`, +this function changes into the [`LinearAffineMetric`](@ref) (default) metric on the +[`SymmetricPositiveDefinite`](@ref) `M`. + +To be precise we are looking for ``c\colon T_p\mathcal P(n) \to T_p\mathcal P(n) `` +such that for all ``Y,Z ∈ T_p\mathcal P(n)``` it holds + +```math +⟨Y,Z⟩ = \operatorname{tr}(YZ) = \operatorname{tr}(p^{-1}c(Y)p^{-1}c(Z)) = g_p(c(Z),c(Y)) +``` + +and hence ``c(X) = pX`` is computed. +""" +change_metric(::SymmetricPositiveDefinite, ::EuclideanMetric, ::Any, ::Any) + +function change_metric!(::SymmetricPositiveDefinite, Y, ::EuclideanMetric, p, X) + Y .= p * X + return Y +end + default_metric_dispatch(::SymmetricPositiveDefinite, ::LinearAffineMetric) = Val(true) @doc raw""" diff --git a/src/ode.jl b/src/ode.jl index 48433b05fc..5e7c783b04 100644 --- a/src/ode.jl +++ b/src/ode.jl @@ -5,7 +5,7 @@ function solve_exp_ode( tspan, B::AbstractBasis; solver=AutoVern9(Rodas5()), - backend=diff_backend(), + backend=default_differential_backend(), kwargs..., ) n = length(x) diff --git a/src/riemannian_diff.jl b/src/riemannian_diff.jl deleted file mode 100644 index de8654bf98..0000000000 --- a/src/riemannian_diff.jl +++ /dev/null @@ -1,248 +0,0 @@ - -""" - AbstractRiemannianDiffBackend - -An abstract type for diff backends. See [`RiemannianONBDiffBackend`](@ref) for -an example. -""" -abstract type AbstractRiemannianDiffBackend end - -@doc raw""" - differential(M::AbstractManifold, f, t::Real, backend::AbstractDiffBackend = rdifferential_backend()) - -Compute the Riemannian differential of a curve $f: ℝ\to M$ on a manifold `M` -represented by function `f` at time `t` using the given backend. -It is calculated as the tangent vector equal to $\mathrm{d}f_t(t)[1]$. -""" -differential(::AbstractManifold, ::Any, ::Real, ::AbstractRiemannianDiffBackend) - -@doc raw""" - gradient(M::AbstractManifold, f, p, backend::AbstractRiemannianDiffBackend = rgradient_backend()) - -Compute the Riemannian gradient $∇f(p)$ of a real field on manifold `M` represented by -function `f` at point `p` using the given backend. -""" -gradient(::AbstractManifold, ::Any, ::Any, ::AbstractRiemannianDiffBackend) - -function differential!( - M::AbstractManifold, - f::Any, - X, - t, - backend::AbstractRiemannianDiffBackend, -) - return copyto!(X, differential(M, f, t, backend)) -end - -function gradient!(M::AbstractManifold, f, X, p, backend::AbstractRiemannianDiffBackend) - return copyto!(X, gradient(M, f, p, backend)) -end - -differential(M::AbstractManifold, f, p) = differential(M, f, p, rdifferential_backend()) - -function differential!(M::AbstractManifold, f, X, p) - return differential!(M, f, X, p, rdifferential_backend()) -end - -gradient(M::AbstractManifold, f, p) = gradient(M, f, p, rgradient_backend()) - -gradient!(M::AbstractManifold, f, X, p) = gradient!(M, f, X, p, rgradient_backend()) - -""" - RiemannianONBDiffBackend( - diff_backend::AbstractDiffBackend - retraction::AbstractRetractionMethod - inverse_retraction::AbstractInverseRetractionMethod - basis::Union{AbstractOrthonormalBasis,CachedBasis{<:AbstractOrthonormalBasis}}, - ) <: AbstractRiemannianDiffBackend - -Riemannian differentiation based on differentiation in an [`AbstractOrthonormalBasis`](@ref) -`basis` using specified `retraction`, `inverse_retraction` and using backend `diff_backend`. -""" -struct RiemannianONBDiffBackend{ - TADBackend<:AbstractDiffBackend, - TRetr<:AbstractRetractionMethod, - TInvRetr<:AbstractInverseRetractionMethod, - TBasis<:Union{ - AbstractOrthonormalBasis, - CachedBasis{𝔽,<:AbstractOrthonormalBasis{𝔽}} where {𝔽}, - }, -} <: AbstractRiemannianDiffBackend - diff_backend::TADBackend - retraction::TRetr - inverse_retraction::TInvRetr - basis::TBasis -end - -function differential(M::AbstractManifold, f, t::Real, backend::RiemannianONBDiffBackend) - p = f(t) - onb_coords = _derivative(zero(number_eltype(p)), backend.diff_backend) do h - return get_coordinates( - M, - p, - inverse_retract(M, p, f(t + h), backend.inverse_retraction), - backend.basis, - ) - end - return get_vector(M, p, onb_coords, backend.basis) -end - -function differential!( - M::AbstractManifold, - f, - X, - t::Real, - backend::RiemannianONBDiffBackend, -) - p = f(t) - onb_coords = _derivative(zero(number_eltype(p)), backend.diff_backend) do h - return get_coordinates( - M, - p, - inverse_retract(M, p, f(t + h), backend.inverse_retraction), - backend.basis, - ) - end - return get_vector!(M, X, p, onb_coords, backend.basis) -end - -function gradient(M::AbstractManifold, f, p, backend::RiemannianONBDiffBackend) - X = get_coordinates(M, p, zero_vector(M, p), backend.basis) - onb_coords = _gradient(X, backend.diff_backend) do Y - return f(retract(M, p, get_vector(M, p, Y, backend.basis), backend.retraction)) - end - return get_vector(M, p, onb_coords, backend.basis) -end - -function gradient!(M::AbstractManifold, f, X, p, backend::RiemannianONBDiffBackend) - X2 = get_coordinates(M, p, zero_vector(M, p), backend.basis) - onb_coords = _gradient(X2, backend.diff_backend) do Y - return f(retract(M, p, get_vector(M, p, Y, backend.basis), backend.retraction)) - end - return get_vector!(M, X, p, onb_coords, backend.basis) -end - -""" - CurrentRiemannianDiffBackend(backend::AbstractRiemannianDiffBackend) - -A mutable struct for storing the current Riemannian differentiation backend in global -constants [`Manifolds._current_rgradient_backend`](@ref) and -[`Manifolds._current_rdifferential_backend`](@ref). - -# See also - -[`AbstractRiemannianDiffBackend`](@ref), [`rdifferential_backend`](@ref), -[`rdifferential_backend!`](@ref) -""" -mutable struct CurrentRiemannianDiffBackend - backend::AbstractRiemannianDiffBackend -end - -""" - _current_rgradient_backend - -The instance of [`Manifolds.CurrentRiemannianDiffBackend`](@ref) that stores the -globally default differentiation backend for calculating gradients. - -# See also - -[`Manifolds.gradient(::AbstractManifold, ::Any, ::Any, ::AbstractRiemannianDiffBackend)`](@ref) -""" -const _current_rgradient_backend = CurrentRiemannianDiffBackend( - RiemannianONBDiffBackend( - diff_backend(), - ExponentialRetraction(), - LogarithmicInverseRetraction(), - DefaultOrthonormalBasis(), - ), -) - -""" - _current_rdifferential_backend - -The instance of [`Manifolds.CurrentRiemannianDiffBackend`](@ref) that stores the -globally default differentiation backend for calculating differentials. - -# See also - -[`Manifolds.differential`](@ref) -""" -const _current_rdifferential_backend = CurrentRiemannianDiffBackend( - RiemannianONBDiffBackend( - diff_backend(), - ExponentialRetraction(), - LogarithmicInverseRetraction(), - DefaultOrthonormalBasis(), - ), -) - -""" - rgradient_backend() -> AbstractRiemannianDiffBackend - -Get the current differentiation backend for Riemannian gradients. -""" -rgradient_backend() = _current_rgradient_backend.backend - -""" - rgradient_backend!(backend::AbstractRiemannianDiffBackend) - -Set current Riemannian gradient backend for differentiation to `backend`. -""" -function rgradient_backend!(backend::AbstractRiemannianDiffBackend) - _current_rgradient_backend.backend = backend - return backend -end - -""" - rdifferential_backend() -> AbstractRiemannianDiffBackend - -Get the current differentiation backend for Riemannian differentials. -""" -rdifferential_backend() = _current_rdifferential_backend.backend - -""" - rdifferential_backend!(backend::AbstractRiemannianDiffBackend) - -Set current Riemannian differential backend for differentiation to `backend`. -""" -function rdifferential_backend!(backend::AbstractRiemannianDiffBackend) - _current_rdifferential_backend.backend = backend - return backend -end - -""" - RiemannianProjectionGradientBackend( - diff_backend::AbstractDiffBackend - ) <: AbstractRiemannianDiffBackend - -Riemannian differentiation based on differentiation in the ambient space and projection to -the given manifold. Differentiation in the ambient space is performed using -the backend `diff_backend`. - -Only valid for manifolds that are embedded in a special way in the Euclidean space. -See [^Absil2008], Section 3.6.1 for details. - -[^Absil2008]: - > Absil, P. A., et al. Optimization Algorithms on Matrix Manifolds. 2008. -""" -struct RiemannianProjectionGradientBackend{TADBackend<:AbstractDiffBackend} <: - AbstractRiemannianDiffBackend - diff_backend::TADBackend -end - -function gradient(M::AbstractManifold, f, p, backend::RiemannianProjectionGradientBackend) - amb_grad = _gradient(f, p, backend.diff_backend) - return project(M, p, amb_grad) -end - -function gradient!( - M::AbstractManifold, - f, - X, - p, - backend::RiemannianProjectionGradientBackend, -) - amb_grad = embed(M, p, X) - _gradient!(f, amb_grad, p, backend.diff_backend) - return project!(M, X, p, amb_grad) -end diff --git a/test/ambiguities.jl b/test/ambiguities.jl index 418d801665..57988178d2 100644 --- a/test/ambiguities.jl +++ b/test/ambiguities.jl @@ -4,12 +4,26 @@ # Interims solution until we follow what was proposed in # https://discourse.julialang.org/t/avoid-ambiguities-with-individual-number-element-identity/62465/2 fmbs = filter(x -> !any(has_type_in_signature.(x, Identity)), mbs) - @test length(fmbs) <= 20 + FMBS_LIMIT = 20 + @test length(fmbs) <= FMBS_LIMIT + if length(fmbs) > FMBS_LIMIT + for amb in fmbs + println(amb) + println() + end + end ms = Test.detect_ambiguities(Manifolds) # Interims solution until we follow what was proposed in # https://discourse.julialang.org/t/avoid-ambiguities-with-individual-number-element-identity/62465/2 fms = filter(x -> !any(has_type_in_signature.(x, Identity)), ms) - @test length(fms) <= 17 + FMS_LIMIT = 21 + if length(fms) > FMS_LIMIT + for amb in fms + println(amb) + println() + end + end + @test length(fms) <= FMS_LIMIT # this test takes way too long to perform regularly # @test length(our_base_ambiguities()) <= 4 else diff --git a/test/approx_inverse_retraction.jl b/test/approx_inverse_retraction.jl index 00c41360e5..83fa282874 100644 --- a/test/approx_inverse_retraction.jl +++ b/test/approx_inverse_retraction.jl @@ -3,6 +3,8 @@ using LinearAlgebra include("utils.jl") +Random.seed!(10) + @testset "approximate inverse retractions" begin @testset "NLsolveInverseRetraction" begin @testset "constructor" begin @@ -62,7 +64,7 @@ include("utils.jl") NLsolveInverseRetraction(ProjectionRetraction(), X0; project_point=true) X = inverse_retract(M, p, q, inv_retr_method) @test is_vector(M, p, X; atol=1e-9) - @test X ≈ X_exp + @test X ≈ X_exp atol = 1e-8 @test_throws OutOfInjectivityRadiusError inverse_retract( M, p, diff --git a/test/differentiation.jl b/test/differentiation.jl index f6e2bd2992..5bc49d8c67 100644 --- a/test/differentiation.jl +++ b/test/differentiation.jl @@ -10,63 +10,89 @@ using Manifolds: _gradient, _gradient! -using FiniteDifferences +using FiniteDifferences, FiniteDiff using LinearAlgebra: Diagonal, dot @testset "Differentiation backend" begin fd51 = Manifolds.FiniteDifferencesBackend() - @testset "diff_backend" begin - @test diff_backend() isa Manifolds.FiniteDifferencesBackend - @test length(diff_backends()) == 2 - @test diff_backends()[1] isa Manifolds.FiniteDifferencesBackend + @testset "default_differential_backend" begin + #ForwardDiff is loaded first in utils. + @test default_differential_backend() === Manifolds.ForwardDiffBackend() @test length(fd51.method.grid) == 5 # check method order @test typeof(fd51.method).parameters[2] == 1 fd71 = Manifolds.FiniteDifferencesBackend(central_fdm(7, 1)) - @test diff_backend!(fd71) == fd71 - @test diff_backend() == fd71 + @test set_default_differential_backend!(fd71) == fd71 + @test default_differential_backend() == fd71 end using ForwardDiff fwd_diff = Manifolds.ForwardDiffBackend() @testset "ForwardDiff" begin - @test diff_backend() isa Manifolds.FiniteDifferencesBackend - @test length(diff_backends()) == 2 - @test diff_backends()[1] isa Manifolds.FiniteDifferencesBackend - @test diff_backends()[2] == fwd_diff - - @test diff_backend!(fwd_diff) == fwd_diff - @test diff_backend() == fwd_diff - @test diff_backend!(fd51) isa Manifolds.FiniteDifferencesBackend - @test diff_backend() isa Manifolds.FiniteDifferencesBackend - - diff_backend!(fwd_diff) - @test diff_backend() == fwd_diff - diff_backend!(fd51) + @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend + + @test set_default_differential_backend!(fwd_diff) == fwd_diff + @test default_differential_backend() == fwd_diff + @test set_default_differential_backend!(fd51) isa Manifolds.FiniteDifferencesBackend + @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend + + set_default_differential_backend!(fwd_diff) + @test default_differential_backend() == fwd_diff + set_default_differential_backend!(fd51) end using FiniteDiff finite_diff = Manifolds.FiniteDiffBackend() @testset "FiniteDiff" begin - @test diff_backend() isa Manifolds.FiniteDifferencesBackend - @test length(diff_backends()) == 3 - @test diff_backends()[3] == finite_diff - - @test diff_backend!(finite_diff) == finite_diff - @test diff_backend() == finite_diff - @test diff_backend!(fd51) isa Manifolds.FiniteDifferencesBackend - @test diff_backend() isa Manifolds.FiniteDifferencesBackend - - diff_backend!(finite_diff) - @test diff_backend() == finite_diff - diff_backend!(fd51) + @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend + + @test set_default_differential_backend!(finite_diff) == finite_diff + @test default_differential_backend() == finite_diff + @test set_default_differential_backend!(fd51) isa Manifolds.FiniteDifferencesBackend + @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend + + set_default_differential_backend!(finite_diff) + @test default_differential_backend() == finite_diff + set_default_differential_backend!(fd51) + end + + using ReverseDiff + + reverse_diff = Manifolds.ReverseDiffBackend() + @testset "ReverseDiff" begin + @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend + + @test set_default_differential_backend!(reverse_diff) == reverse_diff + @test default_differential_backend() == reverse_diff + @test set_default_differential_backend!(fd51) isa Manifolds.FiniteDifferencesBackend + @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend + + set_default_differential_backend!(reverse_diff) + @test default_differential_backend() == reverse_diff + set_default_differential_backend!(fd51) + end + + using Zygote: Zygote + + zygote_diff = Manifolds.ZygoteDiffBackend() + @testset "Zygote" begin + @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend + + @test set_default_differential_backend!(zygote_diff) == zygote_diff + @test default_differential_backend() == zygote_diff + @test set_default_differential_backend!(fd51) isa Manifolds.FiniteDifferencesBackend + @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend + + set_default_differential_backend!(zygote_diff) + @test default_differential_backend() == zygote_diff + set_default_differential_backend!(fd51) end @testset "gradient" begin - diff_backend!(fd51) + set_default_differential_backend!(fd51) r2 = Euclidean(2) c1(t) = [sin(t), cos(t)] @@ -74,66 +100,70 @@ using LinearAlgebra: Diagonal, dot f2(x) = 3 * x[1] * x[2] + x[2]^3 @testset "Inference" begin - v = [-1.0, -1.0] + X = [-1.0, -1.0] @test (@inferred _derivative(c1, 0.0, Manifolds.ForwardDiffBackend())) ≈ [1.0, 0.0] - @test (@inferred _derivative!(c1, v, 0.0, Manifolds.ForwardDiffBackend())) === v - @test v ≈ [1.0, 0.0] + @test (@inferred _derivative!(c1, X, 0.0, Manifolds.ForwardDiffBackend())) === X + @test X ≈ [1.0, 0.0] @test (@inferred _derivative(c1, 0.0, finite_diff)) ≈ [1.0, 0.0] @test (@inferred _gradient(f1, [1.0, -1.0], finite_diff)) ≈ [1.0, -2.0] end @testset for backend in [fd51, fwd_diff, finite_diff] - diff_backend!(backend) + set_default_differential_backend!(backend) @test _derivative(c1, 0.0) ≈ [1.0, 0.0] - v = [-1.0, -1.0] - @test _derivative!(c1, v, 0.0) === v - @test isapprox(v, [1.0, 0.0]) + X = [-1.0, -1.0] + @test _derivative!(c1, X, 0.0) === X + @test isapprox(X, [1.0, 0.0]) + end + @testset for backend in [fd51, fwd_diff, finite_diff, reverse_diff, zygote_diff] + set_default_differential_backend!(backend) + X = [-1.0, -1.0] @test _gradient(f1, [1.0, -1.0]) ≈ [1.0, -2.0] - @test _gradient!(f1, v, [1.0, -1.0]) === v - @test v ≈ [1.0, -2.0] + @test _gradient!(f1, X, [1.0, -1.0]) === X + @test X ≈ [1.0, -2.0] end - diff_backend!(Manifolds.NoneDiffBackend()) + set_default_differential_backend!(Manifolds.NoneDiffBackend()) @testset for backend in [fd51, Manifolds.ForwardDiffBackend()] @test _derivative(c1, 0.0, backend) ≈ [1.0, 0.0] @test _gradient(f1, [1.0, -1.0], backend) ≈ [1.0, -2.0] end - diff_backend!(fd51) + set_default_differential_backend!(fd51) end end -rb_onb_default = RiemannianONBDiffBackend( - diff_backend(), +rb_onb_default = TangentDiffBackend( + default_differential_backend(), Manifolds.ExponentialRetraction(), Manifolds.LogarithmicInverseRetraction(), DefaultOrthonormalBasis(), ) -rb_onb_fd51 = RiemannianONBDiffBackend( +rb_onb_fd51 = TangentDiffBackend( Manifolds.FiniteDifferencesBackend(), Manifolds.ExponentialRetraction(), Manifolds.LogarithmicInverseRetraction(), DefaultOrthonormalBasis(), ) -rb_onb_fwd_diff = RiemannianONBDiffBackend( +rb_onb_fwd_diff = TangentDiffBackend( Manifolds.ForwardDiffBackend(), Manifolds.ExponentialRetraction(), Manifolds.LogarithmicInverseRetraction(), DefaultOrthonormalBasis(), ) -rb_onb_finite_diff = RiemannianONBDiffBackend( +rb_onb_finite_diff = TangentDiffBackend( Manifolds.FiniteDiffBackend(), Manifolds.ExponentialRetraction(), Manifolds.LogarithmicInverseRetraction(), DefaultOrthonormalBasis(), ) -rb_onb_default2 = RiemannianONBDiffBackend( - diff_backend(), +rb_onb_default2 = TangentDiffBackend( + default_differential_backend(), Manifolds.ExponentialRetraction(), Manifolds.LogarithmicInverseRetraction(), CachedBasis( @@ -142,23 +172,7 @@ rb_onb_default2 = RiemannianONBDiffBackend( ), ) -rb_proj = Manifolds.RiemannianProjectionGradientBackend(diff_backend()) - -@testset "rdiff_ functions" begin - @test Manifolds.rdifferential_backend() === - Manifolds._current_rdifferential_backend.backend - @test Manifolds.rgradient_backend() === Manifolds._current_rgradient_backend.backend - - tmp_diff = Manifolds.rdifferential_backend() - Manifolds.rdifferential_backend!(rb_onb_finite_diff) - @test Manifolds.rdifferential_backend() === rb_onb_finite_diff - Manifolds.rdifferential_backend!(tmp_diff) - - tmp_grad = Manifolds.rgradient_backend() - Manifolds.rgradient_backend!(rb_onb_finite_diff) - @test Manifolds.rgradient_backend() === rb_onb_finite_diff - Manifolds.rgradient_backend!(tmp_grad) -end +rb_proj = Manifolds.RiemannianProjectionBackend(default_differential_backend()) @testset "Riemannian differentials" begin s2 = Sphere(2) @@ -167,9 +181,9 @@ end c1(t) = geodesic(s2, q, p, t) Xval = [-sqrt(2) / 2, 0.0, sqrt(2) / 2] - @test isapprox(s2, c1(π / 4), differential(s2, c1, π / 4), Xval) + @test isapprox(s2, c1(π / 4), differential(s2, c1, π / 4, rb_onb_default), Xval) X = similar(p) - differential!(s2, c1, X, π / 4) + differential!(s2, c1, X, π / 4, rb_onb_default) @test isapprox(s2, c1(π / 4), X, Xval) @testset for backend in [rb_onb_fd51, rb_onb_fwd_diff, rb_onb_finite_diff] @@ -185,13 +199,10 @@ end f1(p) = p[1] q = [sqrt(2) / 2, 0, sqrt(2) / 2] - @test isapprox(s2, q, gradient(s2, f1, q), [0.5, 0.0, -0.5]) for backend in [rb_onb_default, rb_proj] @test isapprox(s2, q, gradient(s2, f1, q, backend), [0.5, 0.0, -0.5]) end X = similar(q) - gradient!(s2, f1, X, q) - @test isapprox(s2, q, X, [0.5, 0.0, -0.5]) for backend in [rb_onb_default, rb_proj] gradient!(s2, f1, X, q, backend) @test isapprox(s2, q, X, [0.5, 0.0, -0.5]) diff --git a/test/groups/special_euclidean.jl b/test/groups/special_euclidean.jl index eedbab4575..09793dd56a 100644 --- a/test/groups/special_euclidean.jl +++ b/test/groups/special_euclidean.jl @@ -78,6 +78,7 @@ Random.seed!(10) X_pts; test_diff=true, diff_convs=[(), (LeftAction(),), (RightAction(),)], + atol=1e-9, ) end end @@ -128,6 +129,7 @@ Random.seed!(10) test_diff=true, test_lie_bracket=true, diff_convs=[(), (LeftAction(),), (RightAction(),)], + atol=1e-9, ) # specific affine tests p = copy(G, pts[1]) diff --git a/test/manifolds/generalized_grassmann.jl b/test/manifolds/generalized_grassmann.jl index 1e81a7ca41..b8c6905321 100644 --- a/test/manifolds/generalized_grassmann.jl +++ b/test/manifolds/generalized_grassmann.jl @@ -4,7 +4,7 @@ include("../utils.jl") @testset "Real" begin B = [1.0 0.0 0.0; 0.0 4.0 0.0; 0.0 0.0 1.0] M = GeneralizedGrassmann(3, 2, B) - x = [1.0 0.0; 0.0 0.5; 0.0 0.0] + p = [1.0 0.0; 0.0 0.5; 0.0 0.0] @testset "Basics" begin @test repr(M) == "GeneralizedGrassmann(3, 2, [1.0 0.0 0.0; 0.0 4.0 0.0; 0.0 0.0 1.0], ℝ)" @@ -13,62 +13,69 @@ include("../utils.jl") @test base_manifold(M) === M @test_throws DomainError is_point(M, [1.0, 0.0, 0.0, 0.0], true) @test_throws DomainError is_point(M, 1im * [1.0 0.0; 0.0 1.0; 0.0 0.0], true) - @test !is_vector(M, x, [0.0, 0.0, 1.0, 0.0]) - @test_throws DomainError is_vector(M, x, [0.0, 0.0, 1.0, 0.0], true) - @test_throws DomainError is_vector(M, x, 1 * im * zero_vector(M, x), true) + @test !is_vector(M, p, [0.0, 0.0, 1.0, 0.0]) + @test_throws DomainError is_vector(M, p, [0.0, 0.0, 1.0, 0.0], true) + @test_throws DomainError is_vector(M, p, 1 * im * zero_vector(M, p), true) @test injectivity_radius(M) == π / 2 @test injectivity_radius(M, ExponentialRetraction()) == π / 2 - @test injectivity_radius(M, x) == π / 2 - @test injectivity_radius(M, x, ExponentialRetraction()) == π / 2 - @test mean(M, [x, x, x]) == x + @test injectivity_radius(M, p) == π / 2 + @test injectivity_radius(M, p, ExponentialRetraction()) == π / 2 + @test mean(M, [p, p, p]) == p end @testset "Embedding and Projection" begin - y = similar(x) - z = embed(M, x) - @test z == x - embed!(M, y, x) - @test y == z + q = similar(p) + p2 = embed(M, p) + @test p2 == p + embed!(M, q, p) + @test q == p2 a = [1.0 0.0; 0.0 2.0; 0.0 0.0] @test !is_point(M, a) b = similar(a) c = project(M, a) - @test c == x + @test c == p project!(M, b, a) - @test b == x + @test b == p X = [0.0 0.0; 0.0 0.0; -1.0 1.0] Y = similar(X) - Z = embed(M, x, X) - embed!(M, Y, x, X) + Z = embed(M, p, X) + embed!(M, Y, p, X) @test Y == X @test Z == X end - + @testset "gradient and metric conversion" begin + L = cholesky(B).L + X = [0.0 0.0; 0.0 0.0; 1.0 -1.0] + Y = change_metric(M, EuclideanMetric(), p, X) + @test Y == L \ X + Z = change_representer(M, EuclideanMetric(), p, X) + @test Z == B \ X + end types = [Matrix{Float64}] TEST_STATIC_SIZED && push!(types, MMatrix{3,2,Float64,6}) X = [0.0 0.0; 1.0 0.0; 0.0 2.0] Y = [0.0 0.0; -1.0 0.0; 0.0 2.0] - @test inner(M, x, X, Y) == 0 - y = retract(M, x, X) - z = retract(M, x, Y) - @test is_point(M, y) - @test is_point(M, z) - @test retract(M, x, X) == exp(M, x, X) + @test inner(M, p, X, Y) == 0 + q = retract(M, p, X) + r = retract(M, p, Y) + @test is_point(M, q) + @test is_point(M, r) + @test retract(M, p, X) == exp(M, p, X) - a = project(M, x + X) - c = retract(M, x, X, ProjectionRetraction()) - d = retract(M, x, X, PolarRetraction()) + a = project(M, p + X) + c = retract(M, p, X, ProjectionRetraction()) + d = retract(M, p, X, PolarRetraction()) @test a == c @test c == d e = similar(a) - retract!(M, e, x, X) - @test e == exp(M, x, X) - @test vector_transport_to(M, x, X, y, ProjectionTransport()) == project(M, y, X) + retract!(M, e, p, X) + @test e == exp(M, p, X) + @test vector_transport_to(M, p, X, q, ProjectionTransport()) == project(M, q, X) @testset "Type $T" for T in types - pts = convert.(T, [x, y, z]) - @test !is_point(M, 2 * x) - @test_throws DomainError !is_point(M, 2 * x, true) - @test !is_vector(M, x, y) - @test_throws DomainError is_vector(M, x, y, true) + pts = convert.(T, [p, q, r]) + @test !is_point(M, 2 * p) + @test_throws DomainError !is_point(M, 2 * r, true) + @test !is_vector(M, p, q) + @test_throws DomainError is_vector(M, p, q, true) test_manifold( M, pts, diff --git a/test/manifolds/hyperbolic.jl b/test/manifolds/hyperbolic.jl index af1272510e..c0caa4782c 100644 --- a/test/manifolds/hyperbolic.jl +++ b/test/manifolds/hyperbolic.jl @@ -44,8 +44,9 @@ include("../utils.jl") copyto!(M, pC, p) @test pC.value == p.value XC = allocate(X) - copyto!(M, XC, p, X) - @test XC.value == X.value + @test copyto!(M, XC, p, X) == X # does copyto return the right value? + @test XC == X # does copyto store the right value? + @test XC.value == X.value # another check end end @testset "Hyperbolic Representation Conversion I" begin @@ -261,4 +262,28 @@ include("../utils.jl") @test isa(Z2, PoincareHalfSpaceTVector) @test Z2 == X2 end + @testset "Metric conversion on Hyperboloid" begin + M = Hyperbolic(2) + p = [1.0, 1.0, sqrt(3)] + X = [1.0, 2.0, sqrt(3)] + Y = change_representer(M, EuclideanMetric(), p, X) + @test inner(M, p, X, Y) == inner(Euclidean(3), p, X, X) + # change metric not possible from Euclidean, since the embedding is Lorenzian + @test_throws ErrorException change_metric(M, EuclideanMetric(), p, X) + # but if we come from the same metric, we have the identity + @test change_metric(M, MinkowskiMetric(), p, X) == X + end + @testset "Metric conversion on Poincare Ball" begin + M = Hyperbolic(2) + p = convert(PoincareBallPoint, [1.0, 1.0, sqrt(3)]) + X = convert(PoincareBallTVector, [1.0, 1.0, sqrt(3)], [1.0, 2.0, sqrt(3)]) + Y = change_representer(M, EuclideanMetric(), p, X) + @test inner(M, p, X, Y) == inner(Euclidean(3), p, X.value, X.value) + α = 2 / (1 - norm(p.value)^2) + @test Y.value == X.value ./ α^2 + Z = change_metric(M, EuclideanMetric(), p, X) + @test Z.value == X.value ./ α + A = change_metric(M, MinkowskiMetric(), p, X) + @test A == X + end end diff --git a/test/manifolds/positive_numbers.jl b/test/manifolds/positive_numbers.jl index 6a37883a7e..3e4b00ee1a 100644 --- a/test/manifolds/positive_numbers.jl +++ b/test/manifolds/positive_numbers.jl @@ -23,6 +23,12 @@ include("../utils.jl") X = similar([1.0]) zero_vector!(M, X, 1.0) @test X == [0.0] + + @test change_metric(M, EuclideanMetric(), 2, 3) == 3 * 2 + @test change_representer(M, EuclideanMetric(), 2, 3) == 3 * 2^2 + N = PositiveVectors(2) + @test change_metric(M, EuclideanMetric(), [1, 2], [3, 4]) == [3, 4 * 2] + @test change_representer(M, EuclideanMetric(), [1, 2], [3, 4]) == [3, 4 * 2^2] end types = [Float64] TEST_FLOAT32 && push!(types, Float32) diff --git a/test/manifolds/power_manifold.jl b/test/manifolds/power_manifold.jl index bc7d03b95b..a66e1f448f 100644 --- a/test/manifolds/power_manifold.jl +++ b/test/manifolds/power_manifold.jl @@ -240,7 +240,7 @@ end rand_tvector_atol_multiplier=6.0, retraction_atol_multiplier=12, is_tangent_atol_multiplier=12.0, - exp_log_atol_multiplier=2 * prod(power_dimensions(Ms2)), + exp_log_atol_multiplier=3 * prod(power_dimensions(Ms2)), test_inplace=true, ) end @@ -413,4 +413,23 @@ end @test isapprox(a, a2) @test_throws ErrorException get_point(M, A2, p, a2) end + + @testset "metric conversion" begin + M = SymmetricPositiveDefinite(3) + N = PowerManifold(M, NestedPowerRepresentation(), 2) + e = EuclideanMetric() + p = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1] + q = [2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1] + P = [p, q] + X = [log(M, p, q), log(M, q, p)] + Y = change_metric(N, e, P, X) + Yc = [change_metric(M, e, p, log(M, p, q)), change_metric(M, e, q, log(M, q, p))] + @test norm(N, P, Y .- Yc) ≈ 0 + Z = change_representer(N, e, P, X) + Zc = [ + change_representer(M, e, p, log(M, p, q)), + change_representer(M, e, q, log(M, q, p)), + ] + @test norm(N, P, Z .- Zc) ≈ 0 + end end diff --git a/test/manifolds/probability_simplex.jl b/test/manifolds/probability_simplex.jl index 071903636b..5d2e09dc70 100644 --- a/test/manifolds/probability_simplex.jl +++ b/test/manifolds/probability_simplex.jl @@ -72,6 +72,12 @@ include("../utils.jl") Y = project(M, p, X2) @test isapprox(M, p, X, Y) + # Check adaption of metric and representer + Y1 = change_metric(M, EuclideanMetric(), p, X) + @test Y1 == X .* sqrt.(p) + Y2 = change_representer(M, EuclideanMetric(), p, X) + @test Y2 == X .* p + X = log(M, q, p) X2 = X + [1, 2, 3] Y = project(M, q, X2) diff --git a/test/manifolds/product_manifold.jl b/test/manifolds/product_manifold.jl index 5d9a742174..07dc94e807 100644 --- a/test/manifolds/product_manifold.jl +++ b/test/manifolds/product_manifold.jl @@ -693,4 +693,26 @@ end p2 = get_point(M, A, p, a) @test all(p2.parts .== p.parts) end + + @testset "metric conversion" begin + M = SymmetricPositiveDefinite(3) + N = ProductManifold(M, M) + e = EuclideanMetric() + p = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1] + q = [2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1] + P = ProductRepr(p, q) + X = ProductRepr(log(M, p, q), log(M, q, p)) + Y = change_metric(N, e, P, X) + Yc = ProductRepr( + change_metric(M, e, p, log(M, p, q)), + change_metric(M, e, q, log(M, q, p)), + ) + @test norm(N, P, Y - Yc) ≈ 0 + Z = change_representer(N, e, P, X) + Zc = ProductRepr( + change_representer(M, e, p, log(M, p, q)), + change_representer(M, e, q, log(M, q, p)), + ) + @test norm(N, P, Z - Zc) ≈ 0 + end end diff --git a/test/manifolds/rotations.jl b/test/manifolds/rotations.jl index 66701f17b0..2eb11efef6 100644 --- a/test/manifolds/rotations.jl +++ b/test/manifolds/rotations.jl @@ -116,7 +116,7 @@ include("../utils.jl") point_distributions=[ptd], tvector_distributions=[tvd], basis_types_to_from=basis_types, - exp_log_atol_multiplier=20, + exp_log_atol_multiplier=250, retraction_atol_multiplier=12, test_inplace=true, ) diff --git a/test/manifolds/sphere.jl b/test/manifolds/sphere.jl index ff533f8bf1..f5889e9c26 100644 --- a/test/manifolds/sphere.jl +++ b/test/manifolds/sphere.jl @@ -214,4 +214,13 @@ using ManifoldsBase: TFVector end end end + + @testset "Metric conversion is the identity" begin + p = [1.0, 0.0, 0.0] + X = [0.0, 1.0, 1.0] + Y = change_representer(M, EuclideanMetric(), p, X) + @test Y == X + Z = change_metric(M, EuclideanMetric(), p, X) + @test Z == X + end end diff --git a/test/manifolds/stiefel.jl b/test/manifolds/stiefel.jl index e9d89ee7c6..b5c172efd2 100644 --- a/test/manifolds/stiefel.jl +++ b/test/manifolds/stiefel.jl @@ -311,7 +311,7 @@ using Manifolds: default_metric_dispatch M4 = MetricManifold(Stiefel(10, 2), CanonicalMetric()) p = Matrix{Float64}(I, 10, 2) Random.seed!(42) - Z = project(base_manifold(M4), p, randn(size(p))) + Z = project(base_manifold(M4), p, 0.2 .* randn(size(p))) s = exp(M4, p, Z) Z2 = log(M4, p, s) @test isapprox(M4, p, Z, Z2) diff --git a/test/manifolds/symmetric_positive_definite.jl b/test/manifolds/symmetric_positive_definite.jl index eeeb447cb1..0b134d378e 100644 --- a/test/manifolds/symmetric_positive_definite.jl +++ b/test/manifolds/symmetric_positive_definite.jl @@ -87,35 +87,35 @@ using Manifolds: default_metric_dispatch end end end - x = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1] - y = [2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1] + p = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1] + q = [2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1] @testset "Convert SPD to Cholesky" begin - v = log(M1, x, y) - (l, w) = Manifolds.spd_to_cholesky(x, v) + v = log(M1, p, q) + (l, w) = Manifolds.spd_to_cholesky(p, v) (xs, vs) = Manifolds.cholesky_to_spd(l, w) - @test isapprox(xs, x) + @test isapprox(xs, p) @test isapprox(vs, v) end @testset "Preliminary tests for LogEuclidean" begin @test representation_size(M4) == (3, 3) - @test isapprox(distance(M4, x, y), sqrt(2) * log(2)) + @test isapprox(distance(M4, p, q), sqrt(2) * log(2)) @test manifold_dimension(M4) == manifold_dimension(M1) end @testset "Test for tangent ONB on LinearAffineMetric" begin - v = log(M2, x, y) - donb = get_basis(base_manifold(M2), x, DiagonalizingOrthonormalBasis(v)) - X = get_vectors(base_manifold(M2), x, donb) + v = log(M2, p, q) + donb = get_basis(base_manifold(M2), p, DiagonalizingOrthonormalBasis(v)) + X = get_vectors(base_manifold(M2), p, donb) k = donb.data.eigenvalues @test isapprox(0.0, first(k)) for i in 1:length(X) - @test isapprox(1.0, norm(M2, x, X[i])) + @test isapprox(1.0, norm(M2, p, X[i])) for j in (i + 1):length(X) - @test isapprox(0.0, inner(M2, x, X[i], X[j])) + @test isapprox(0.0, inner(M2, p, X[i], X[j])) end end - d2onb = get_basis(M2, x, DiagonalizingOrthonormalBasis(v)) + d2onb = get_basis(M2, p, DiagonalizingOrthonormalBasis(v)) @test donb.data.eigenvalues == d2onb.data.eigenvalues - @test get_vectors(base_manifold(M2), x, donb) == get_vectors(M2, x, d2onb) + @test get_vectors(base_manifold(M2), p, donb) == get_vectors(M2, p, d2onb) end @testset "Vector transport and transport along with Schild and Pole ladder" begin A(α) = [1.0 0.0 0.0; 0.0 cos(α) sin(α); 0.0 -sin(α) cos(α)] @@ -141,4 +141,11 @@ using Manifolds: default_metric_dispatch @test inner(M, p1, X1, X2) ≈ inner(M, p2, Y1, Y4) # parallel transport isometric @test inner(M, p1, X1, X2) ≈ inner(M, p2, Y2, Y5) # pole ladder transport isometric end + @testset "Metric change for Linear Affine Metric" begin + X = log(M1, p, q) + Y = change_metric(M1, EuclideanMetric(), p, X) + @test Y == p * X + Z = change_representer(M1, EuclideanMetric(), p, X) + @test Z == p * X * p + end end diff --git a/test/metric.jl b/test/metric.jl index ae7d7bb6d6..a1562bfe2c 100644 --- a/test/metric.jl +++ b/test/metric.jl @@ -7,6 +7,7 @@ include("utils.jl") struct TestEuclidean{N} <: AbstractManifold{ℝ} end struct TestEuclideanMetric <: AbstractMetric end +struct TestScaledEuclideanMetric <: AbstractMetric end Manifolds.manifold_dimension(::TestEuclidean{N}) where {N} = N function Manifolds.local_metric( @@ -19,11 +20,57 @@ end function Manifolds.local_metric( M::MetricManifold{ℝ,<:TestEuclidean,<:TestEuclideanMetric}, ::Any, - ::InducedBasis, -) + ::T, +) where {T<:ManifoldsBase.AbstractOrthogonalBasis} return Diagonal(1.0:manifold_dimension(M)) end - +function Manifolds.local_metric( + M::MetricManifold{ℝ,<:TestEuclidean,<:TestScaledEuclideanMetric}, + ::Any, + ::T, +) where {T<:ManifoldsBase.AbstractOrthogonalBasis} + return 2 .* Diagonal(1.0:manifold_dimension(M)) +end +function Manifolds.get_coordinates!( + M::MetricManifold{ℝ,<:TestEuclidean,<:TestEuclideanMetric}, + c, + ::Any, + X, + ::DefaultOrthogonalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, +) + c .= 1 ./ [1.0:manifold_dimension(M)...] .* X + return c +end +function Manifolds.get_vector!( + M::MetricManifold{ℝ,<:TestEuclidean,<:TestEuclideanMetric}, + X, + ::Any, + c, + ::DefaultOrthogonalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, +) + X .= [1.0:manifold_dimension(M)...] .* c + return X +end +function Manifolds.get_coordinates!( + M::MetricManifold{ℝ,<:TestEuclidean,<:TestScaledEuclideanMetric}, + c, + ::Any, + X, + ::DefaultOrthogonalBasis, +) + c .= 1 ./ (2 .* [1.0:manifold_dimension(M)...]) .* X + return c +end +function Manifolds.get_vector!( + M::MetricManifold{ℝ,<:TestEuclidean,<:TestScaledEuclideanMetric}, + X, + ::Any, + c, + ::DefaultOrthogonalBasis, +) + X .= 2 .* [1.0:manifold_dimension(M)...] .* c + return X +end struct TestSphere{N,T} <: AbstractManifold{ℝ} r::T end @@ -603,4 +650,20 @@ end ExponentialRetraction(), ) === Val{:parent}() end + + @testset "change metric and representer" begin + M = MetricManifold(TestEuclidean{2}(), TestEuclideanMetric()) + G = TestScaledEuclideanMetric() + M2 = TestScaledEuclideanMetric(M) + @test M2.manifold === M.manifold + @test M2.metric == G + p = ones(2) + X = 2 * ones(2) + @test change_metric(M, TestEuclideanMetric(), p, X) == X + Y = change_metric(M, G, p, X) + @test Y ≈ sqrt(2) .* X #scaled metric has a factor 2, removing introduces this factor + @test change_representer(M, TestEuclideanMetric(), p, X) == X + Y2 = change_representer(M, G, p, X) + @test Y2 ≈ 2 .* X #scaled metric has a factor 2, removing introduces this factor + end end From 59c0467508ac8c49cdb914a9f68c5cd8c70cac01 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sat, 25 Sep 2021 12:18:48 +0200 Subject: [PATCH 02/20] first sketch of the ODEExponenitalRetraction --- src/manifolds/ConnectionManifold.jl | 77 +++++++++++++++++++++++++---- src/manifolds/MetricManifold.jl | 9 ++-- src/ode.jl | 6 ++- test/differentiation.jl | 52 ------------------- 4 files changed, 79 insertions(+), 65 deletions(-) diff --git a/src/manifolds/ConnectionManifold.jl b/src/manifolds/ConnectionManifold.jl index 7ee5d4f876..81b03da12d 100644 --- a/src/manifolds/ConnectionManifold.jl +++ b/src/manifolds/ConnectionManifold.jl @@ -54,6 +54,41 @@ struct ConnectionManifold{𝔽,M<:AbstractManifold{𝔽},C<:AbstractAffineConnec connection::C end +@doc raw""" + ODEExponentialRetraction{T<:AbstractRetractionMethod, B<: AbstractBasis} <: AbstractRetractionMethod + +This retraction approximates the exponential map by solving the correspondig ODE. +Let ``p\in\mathal M`` and ``X\in T_p\mathcal M`` denote the input for the exponential map +and ``d`` denote the [`manifold_dimension`](@ref) of `M`. + +This the ODE is formulated in a chart constructed using an [`AbstractBasis`](@ref) `B` and an +[`AbstractRetractionMethod`](@ref) `R` as follows. +Given some coordinates ``c\in ℝ^d`` - these can be used to form a tangent vector +with restect to th basis `B, i.e. ``c \mapsto Y=``[`get_vector`](@ref)`(M, p, c, B)`. +Further, using the retraction we can map ``Y`` to a point on the manifold +``Y \mapsto q =``[`retract`](@ref)`(M, p, X, R)`. + +Hence the ODE can be formulated in a curve ``c(t)`` in parameter space. +This is – for sure – only possible locally as fas as the retraction is well-defined. +""" +struct ODEExponentialRetraction{T<:AbstractRetractionMethod,B<:AbstractBasis} <: + AbstractRetractionMethod + retraction::T + basis::B +end +function ODEExponentialRetraction(r::T) where {T<:AbstractRetractionMethod} + return ODEExponentialRetraction{T,ODefaultOrthonormalBasis}( + r, + ODefaultOrthonormalBasis(), + ) +end +function ODEExponentialRetraction(r::ExponentialRetraction, ::AbstractBasis) + return DomainError( + r, + "You can not use the exponential map as an inner method to solve the ode for the exponential map.", + ) +end + @doc raw""" christoffel_symbols_second( M::AbstractManifold, @@ -162,13 +197,13 @@ in an embedded space. exp(::AbstractConnectionManifold, ::Any...) @decorator_transparent_fallback function exp!(M::AbstractConnectionManifold, q, p, X) - tspan = (0.0, 1.0) - A = get_default_atlas(M) - i = get_chart_index(M, A, p) - B = induced_basis(M, A, i, TangentSpace) - sol = solve_exp_ode(M, p, X, tspan, B; dense=false, saveat=[1.0]) - n = length(p) - return copyto!(q, sol.u[1][(n + 1):end]) + return retract!( + M, + q, + p, + X, + ODEXxponentialRetration(ManifoldsBase.default_retraction_method(M)), + ) end """ @@ -247,10 +282,18 @@ function riemann_tensor( p, B::AbstractBasis; backend::AbstractDiffBackend=default_differential_backend(), + retraction::AbstractRetractionMethod=ManifoldsBase.default_retraction_method(M), ) n = size(p, 1) - Γ = christoffel_symbols_second(M, p, B; backend=backend) - ∂Γ = christoffel_symbols_second_jacobian(M, p, B; backend=backend) ./ n + Γ = christoffel_symbols_second(M, p, B; backend=backend, retraction=retraction) + ∂Γ = + christoffel_symbols_second_jacobian( + M, + p, + B; + backend=backend, + retraction=retraction, + ) ./ n R = allocate(∂Γ, Size(n, n, n, n)) @einsum R[l, i, j, k] = ∂Γ[l, i, k, j] - ∂Γ[l, i, j, k] + Γ[s, i, k] * Γ[l, s, j] - Γ[s, i, j] * Γ[l, s, k] @@ -263,6 +306,22 @@ end kwargs..., ) +function retract(::AbstractConnectionManifold, q, p, X, r::ODEExponentialRetraction) + tspan = (0.0, 1.0) + sol = solve_exp_ode( + M, + p, + X, + tspan, + r.basis; + retraction=r.retraction, + dense=false, + saveat=[1.0], + ) + n = length(p) + return copyto!(q, sol.u[1][(n + 1):end]) +end + @doc raw""" solve_exp_ode( M::AbstractConnectionManifold, diff --git a/src/manifolds/MetricManifold.jl b/src/manifolds/MetricManifold.jl index f3c91073ba..56f6cdec86 100644 --- a/src/manifolds/MetricManifold.jl +++ b/src/manifolds/MetricManifold.jl @@ -220,8 +220,9 @@ function christoffel_symbols_first( p, B::AbstractBasis; backend::AbstractDiffBackend=default_differential_backend(), + retraction::AbstractRetractionMethod=ManifoldsBase.default_retraction_method(M), ) - ∂g = local_metric_jacobian(M, p, B; backend=backend) + ∂g = local_metric_jacobian(M, p, B; backend=backend, retraction=retraction) n = size(∂g, 1) Γ = allocate(∂g, Size(n, n, n)) @einsum Γ[i, j, k] = 1 / 2 * (∂g[k, j, i] + ∂g[i, k, j] - ∂g[i, j, k]) @@ -239,9 +240,10 @@ function christoffel_symbols_second( p, B::AbstractBasis; backend::AbstractDiffBackend=default_differential_backend(), + retraction::AbstractRetractionMethod=ManifoldsBase.default_retraction_method(M), ) Ginv = inverse_local_metric(M, p, B) - Γ₁ = christoffel_symbols_first(M, p, B; backend=backend) + Γ₁ = christoffel_symbols_first(M, p, B; backend=backend, retraction=retraction) Γ₂ = allocate(Γ₁) @einsum Γ₂[l, i, j] = Ginv[k, l] * Γ₁[i, j, k] return Γ₂ @@ -542,9 +544,10 @@ function ricci_curvature( p, B::AbstractBasis; backend::AbstractDiffBackend=default_differential_backend(), + retraction::AbstractRetractionMethod=ManifoldsBase.default_retraction_method(M), ) Ginv = inverse_local_metric(M, p, B) - Ric = ricci_tensor(M, p, B; backend=backend) + Ric = ricci_tensor(M, p, B; backend=backend, retraction=retraction) S = sum(Ginv .* Ric) return S end diff --git a/src/ode.jl b/src/ode.jl index 5e7c783b04..4cbe2dc579 100644 --- a/src/ode.jl +++ b/src/ode.jl @@ -1,3 +1,6 @@ +@doc """ + + """ function solve_exp_ode( M::MetricManifold, x, @@ -6,6 +9,7 @@ function solve_exp_ode( B::AbstractBasis; solver=AutoVern9(Rodas5()), backend=default_differential_backend(), + retraction::AbstractRetractionMethod=ManifoldsBase.default_retraction_method(M), kwargs..., ) n = length(x) @@ -21,7 +25,7 @@ function solve_exp_ode( x = u[ix] ddx = allocate(u, Size(n)) du = allocate(u) - Γ = christoffel_symbols_second(M, x, B; backend=backend) + Γ = christoffel_symbols_second(M, x, B; backend=backend, retraction=retraction) @einsum ddx[k] = -Γ[k, i, j] * dx[i] * dx[j] du[iv] .= ddx du[ix] .= dx diff --git a/test/differentiation.jl b/test/differentiation.jl index e3e1999325..2372d90d2a 100644 --- a/test/differentiation.jl +++ b/test/differentiation.jl @@ -75,58 +75,6 @@ using LinearAlgebra: Diagonal, dot set_default_differential_backend!(fd51) end - using Zygote: Zygote - - zygote_diff = Manifolds.ZygoteDiffBackend() - @testset "Zygote" begin - @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend - - @test set_default_differential_backend!(zygote_diff) == zygote_diff - @test default_differential_backend() == zygote_diff - @test set_default_differential_backend!(fd51) isa Manifolds.FiniteDifferencesBackend - @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend - - set_default_differential_backend!(zygote_diff) - @test default_differential_backend() == zygote_diff - set_default_differential_backend!(fd51) - end - - using ReverseDiff - - reverse_diff = Manifolds.ReverseDiffBackend() - @testset "ReverseDiff" begin - @test diff_backend() isa Manifolds.FiniteDifferencesBackend - @test length(diff_backends()) == 4 - @test diff_backends()[3] == reverse_diff - - @test diff_backend!(reverse_diff) == reverse_diff - @test diff_backend() == reverse_diff - @test diff_backend!(fd51) isa Manifolds.FiniteDifferencesBackend - @test diff_backend() isa Manifolds.FiniteDifferencesBackend - - diff_backend!(reverse_diff) - @test diff_backend() == reverse_diff - diff_backend!(fd51) - end - - using Zygote: Zygote - - zygote_diff = Manifolds.ZygoteDiffBackend() - @testset "Zygote" begin - @test diff_backend() isa Manifolds.FiniteDifferencesBackend - @test length(diff_backends()) == 5 - @test diff_backends()[5] == zygote_diff - - @test diff_backend!(zygote_diff) == zygote_diff - @test diff_backend() == zygote_diff - @test diff_backend!(fd51) isa Manifolds.FiniteDifferencesBackend - @test diff_backend() isa Manifolds.FiniteDifferencesBackend - - diff_backend!(zygote_diff) - @test diff_backend() == zygote_diff - diff_backend!(fd51) - end - @testset "gradient" begin set_default_differential_backend!(fd51) r2 = Euclidean(2) From f6bede48cdb478f32462de5f5e6b536f1ac945cc Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sat, 25 Sep 2021 13:02:29 +0200 Subject: [PATCH 03/20] Fixes two typos. --- src/manifolds/ConnectionManifold.jl | 2 +- test/differentiation.jl | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/manifolds/ConnectionManifold.jl b/src/manifolds/ConnectionManifold.jl index 81b03da12d..04b035cb66 100644 --- a/src/manifolds/ConnectionManifold.jl +++ b/src/manifolds/ConnectionManifold.jl @@ -202,7 +202,7 @@ exp(::AbstractConnectionManifold, ::Any...) q, p, X, - ODEXxponentialRetration(ManifoldsBase.default_retraction_method(M)), + ODEExponentialRetration(ManifoldsBase.default_retraction_method(M)), ) end diff --git a/test/differentiation.jl b/test/differentiation.jl index 2372d90d2a..588ec3de9b 100644 --- a/test/differentiation.jl +++ b/test/differentiation.jl @@ -75,6 +75,8 @@ using LinearAlgebra: Diagonal, dot set_default_differential_backend!(fd51) end + zygote_diff = Manifold.ZygoteDiffBackend() + @testset "gradient" begin set_default_differential_backend!(fd51) r2 = Euclidean(2) From 475140ef375aadd936742f028821d3b814777de0 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sat, 25 Sep 2021 14:40:18 +0200 Subject: [PATCH 04/20] Fix more typos. --- src/Manifolds.jl | 1 + src/manifolds/ConnectionManifold.jl | 2 +- test/ambiguities.jl | 2 +- test/differentiation.jl | 3 ++- 4 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/Manifolds.jl b/src/Manifolds.jl index 98022f2ad3..681c523084 100644 --- a/src/Manifolds.jl +++ b/src/Manifolds.jl @@ -444,6 +444,7 @@ export AbstractRetractionMethod, PolarRetraction, ProjectionRetraction, SoftmaxRetraction, + ODEExponentialRetraction, PadeRetraction, ProductRetraction, PowerRetraction diff --git a/src/manifolds/ConnectionManifold.jl b/src/manifolds/ConnectionManifold.jl index 04b035cb66..0fc6ac83b2 100644 --- a/src/manifolds/ConnectionManifold.jl +++ b/src/manifolds/ConnectionManifold.jl @@ -202,7 +202,7 @@ exp(::AbstractConnectionManifold, ::Any...) q, p, X, - ODEExponentialRetration(ManifoldsBase.default_retraction_method(M)), + ODEExponentialRetraction(ManifoldsBase.default_retraction_method(M)), ) end diff --git a/test/ambiguities.jl b/test/ambiguities.jl index 57988178d2..0f2ba5716d 100644 --- a/test/ambiguities.jl +++ b/test/ambiguities.jl @@ -4,7 +4,7 @@ # Interims solution until we follow what was proposed in # https://discourse.julialang.org/t/avoid-ambiguities-with-individual-number-element-identity/62465/2 fmbs = filter(x -> !any(has_type_in_signature.(x, Identity)), mbs) - FMBS_LIMIT = 20 + FMBS_LIMIT = 21 @test length(fmbs) <= FMBS_LIMIT if length(fmbs) > FMBS_LIMIT for amb in fmbs diff --git a/test/differentiation.jl b/test/differentiation.jl index 588ec3de9b..9366dd2f00 100644 --- a/test/differentiation.jl +++ b/test/differentiation.jl @@ -75,7 +75,8 @@ using LinearAlgebra: Diagonal, dot set_default_differential_backend!(fd51) end - zygote_diff = Manifold.ZygoteDiffBackend() + using Zygote + zygote_diff = Manifolds.ZygoteDiffBackend() @testset "gradient" begin set_default_differential_backend!(fd51) From 61b4dcbffe91642a13a01f16e4ab5b68917858e8 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sat, 25 Sep 2021 19:05:34 +0200 Subject: [PATCH 05/20] Move solve_exp_ode with OridnaryDiffEq to differentiation and change the return value to be just `p`. --- src/Manifolds.jl | 2 +- src/differentiation/ode.jl | 35 ++++++++++++++++++++++++++ src/manifolds/ConnectionManifold.jl | 39 ++++++++++++++--------------- src/manifolds/Euclidean.jl | 23 ++++++++++++++--- src/ode.jl | 39 ----------------------------- test/ambiguities.jl | 2 +- 6 files changed, 76 insertions(+), 64 deletions(-) create mode 100644 src/differentiation/ode.jl delete mode 100644 src/ode.jl diff --git a/src/Manifolds.jl b/src/Manifolds.jl index 681c523084..f89f3246cc 100644 --- a/src/Manifolds.jl +++ b/src/Manifolds.jl @@ -298,7 +298,7 @@ function __init__() @require OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" begin using .OrdinaryDiffEq: ODEProblem, AutoVern9, Rodas5, solve - include("ode.jl") + include("differentiation/ode.jl") end @require NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" begin diff --git a/src/differentiation/ode.jl b/src/differentiation/ode.jl new file mode 100644 index 0000000000..d2fdbfdf01 --- /dev/null +++ b/src/differentiation/ode.jl @@ -0,0 +1,35 @@ +function solve_exp_ode( + M::AbstractConnectionManifold, + p, + X; + basis::AbstractBasis=DefaultOrthonormalBasis(), + solver=AutoVern9(Rodas5()), + backend=default_differential_backend(), + retraction::AbstractRetractionMethod=ManifoldsBase.default_retraction_method(M), + kwargs..., +) + d = manifold_dimension(M) + iv = SVector{n}(1:d) + ix = SVector{n}((n + 1):(2 * d)) + u0 = allocate(p, 2 * d) + u0[iv] .= X + u0[ix] .= p + + function exp_problem(u, params, t) + M = params[1] + dx = u[iv] + x = u[ix] + ddx = allocate(u, Size(d)) + du = allocate(u) + Γ = christoffel_symbols_second(M, x, basis; backend=backend, retraction=retraction) + @einsum ddx[k] = -Γ[k, i, j] * dx[i] * dx[j] + du[iv] .= ddx + du[ix] .= dx + return Base.convert(typeof(u), du) + end + params = (M,) + prob = ODEProblem(exp_problem, u0, (0.0, 1.0), params) + sol = solve(prob, solver; kwargs...) + q = sol.u[1][(d + 1):(2 * d)] + return q +end diff --git a/src/manifolds/ConnectionManifold.jl b/src/manifolds/ConnectionManifold.jl index 0fc6ac83b2..dd22df5a79 100644 --- a/src/manifolds/ConnectionManifold.jl +++ b/src/manifolds/ConnectionManifold.jl @@ -77,9 +77,12 @@ struct ODEExponentialRetraction{T<:AbstractRetractionMethod,B<:AbstractBasis} <: basis::B end function ODEExponentialRetraction(r::T) where {T<:AbstractRetractionMethod} - return ODEExponentialRetraction{T,ODefaultOrthonormalBasis}( + return ODEExponentialRetraction{T,DefaultOrthonormalBasis}(r, DefaultOrthonormalBasis()) +end +function ODEExponentialRetraction(r::T, ::CachedBasis) where {T<:AbstractRetractionMethod} + return DomainError( r, - ODefaultOrthonormalBasis(), + "Cached Bases are currently not supported, since the basis has to be implemented in a surrounding of the start point as well.", ) end function ODEExponentialRetraction(r::ExponentialRetraction, ::AbstractBasis) @@ -88,6 +91,12 @@ function ODEExponentialRetraction(r::ExponentialRetraction, ::AbstractBasis) "You can not use the exponential map as an inner method to solve the ode for the exponential map.", ) end +function ODEExponentialRetraction(r::ExponentialRetraction, ::CachedBasis) + return DomainError( + r, + "Neither the exponential map nor a Cached Basis can be used with this retraction type.", + ) +end @doc raw""" christoffel_symbols_second( @@ -109,6 +118,7 @@ where ``Γ_{ijk}`` are the Christoffel symbols of the first kind (see [`christoffel_symbols_first`](@ref)), and ``g^{kl}`` is the inverse of the local representation of the metric tensor. The dimensions of the resulting multi-dimensional array are ordered ``(l,i,j)``. + """ christoffel_symbols_second(::AbstractManifold, ::Any, ::AbstractBasis) @@ -307,19 +317,9 @@ end ) function retract(::AbstractConnectionManifold, q, p, X, r::ODEExponentialRetraction) - tspan = (0.0, 1.0) - sol = solve_exp_ode( - M, - p, - X, - tspan, - r.basis; - retraction=r.retraction, - dense=false, - saveat=[1.0], - ) - n = length(p) - return copyto!(q, sol.u[1][(n + 1):end]) + sol = + solve_exp_ode(M, p, X, r.basis; retraction=r.retraction, dense=false, saveat=[1.0]) + return copyto!(q, sol) end @doc raw""" @@ -327,14 +327,14 @@ end M::AbstractConnectionManifold, p, X, - tspan, B::AbstractBasis; backend::AbstractDiffBackend = default_differential_backend(), solver = AutoVern9(Rodas5()), + retraction = ManifoldsBase.default_retraction_method(M), kwargs..., ) -Approximate the exponential map on the manifold over the provided timespan +Approximate the exponential map on the manifold by evaluating the ODE descripting the geodesic at 1, assuming the default connection of the given manifold by solving the ordinary differential equation @@ -347,9 +347,8 @@ the Einstein summation convention is assumed. The arguments `tspan` and `solver` follow the `OrdinaryDiffEq` conventions. `kwargs...` specify keyword arguments that will be passed to `OrdinaryDiffEq.solve`. -Currently, the numerical integration is only accurate when using a single -coordinate chart that covers the entire manifold. This excludes coordinates -in an embedded space. +Using a retraction and a basis, the solver works on any intrinsic manifold that +provides a basis for tangent spaces around the start point `p`. !!! note This function only works when diff --git a/src/manifolds/Euclidean.jl b/src/manifolds/Euclidean.jl index 71888e6133..f542ee79fc 100644 --- a/src/manifolds/Euclidean.jl +++ b/src/manifolds/Euclidean.jl @@ -199,7 +199,16 @@ function get_basis(M::Euclidean, p, B::DiagonalizingOrthonormalBasis) return CachedBasis(B, DiagonalizingBasisData(B.frame_direction, eigenvalues, vecs)) end -function get_coordinates!(M::Euclidean, Y, p, X, ::DefaultOrDiagonalizingBasis{ℝ}) +function get_coordinates!( + M::Euclidean, + Y, + p, + X, + ::Union{ + DefaultOrDiagonalizingBasis{ℝ}, + InducedBasis{ℝ,TangentSpaceType,<:RetractionAtlas}, + }, +) S = representation_size(M) PS = prod(S) copyto!(Y, reshape(X, PS)) @@ -218,7 +227,16 @@ function get_coordinates!( return Y end -function get_vector!(M::Euclidean, Y, ::Any, X, ::DefaultOrDiagonalizingBasis{ℝ}) +function get_vector!( + M::Euclidean, + Y, + ::Any, + X, + ::Union{ + DefaultOrDiagonalizingBasis{ℝ}, + InducedBasis{ℝ,TangentSpaceType,<:RetractionAtlas}, + }, +) S = representation_size(M) copyto!(Y, reshape(X, S)) return Y @@ -239,7 +257,6 @@ function get_vector!(M::Euclidean{<:Tuple,ℂ}, Y, ::Any, X, ::DefaultOrDiagonal copyto!(Y, reshape(X[1:N] + im * X[(N + 1):end], S)) return Y end - @doc raw""" injectivity_radius(M::Euclidean) diff --git a/src/ode.jl b/src/ode.jl deleted file mode 100644 index 4cbe2dc579..0000000000 --- a/src/ode.jl +++ /dev/null @@ -1,39 +0,0 @@ -@doc """ - - """ -function solve_exp_ode( - M::MetricManifold, - x, - v, - tspan, - B::AbstractBasis; - solver=AutoVern9(Rodas5()), - backend=default_differential_backend(), - retraction::AbstractRetractionMethod=ManifoldsBase.default_retraction_method(M), - kwargs..., -) - n = length(x) - iv = SVector{n}(1:n) - ix = SVector{n}((n + 1):(2n)) - u0 = allocate(x, 2n) - u0[iv] .= v - u0[ix] .= x - - function exp_problem(u, p, t) - M = p[1] - dx = u[iv] - x = u[ix] - ddx = allocate(u, Size(n)) - du = allocate(u) - Γ = christoffel_symbols_second(M, x, B; backend=backend, retraction=retraction) - @einsum ddx[k] = -Γ[k, i, j] * dx[i] * dx[j] - du[iv] .= ddx - du[ix] .= dx - return Base.convert(typeof(u), du) - end - - p = (M,) - prob = ODEProblem(exp_problem, u0, tspan, p) - sol = solve(prob, solver; kwargs...) - return sol -end diff --git a/test/ambiguities.jl b/test/ambiguities.jl index 0f2ba5716d..6377f04ccf 100644 --- a/test/ambiguities.jl +++ b/test/ambiguities.jl @@ -16,7 +16,7 @@ # Interims solution until we follow what was proposed in # https://discourse.julialang.org/t/avoid-ambiguities-with-individual-number-element-identity/62465/2 fms = filter(x -> !any(has_type_in_signature.(x, Identity)), ms) - FMS_LIMIT = 21 + FMS_LIMIT = 22 if length(fms) > FMS_LIMIT for amb in fms println(amb) From 49ecd3cbcd8466a1755393fa443b50aa2773d50c Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sat, 25 Sep 2021 19:42:02 +0200 Subject: [PATCH 06/20] fixes a first test. --- src/differentiation/ode.jl | 8 ++++---- src/manifolds/ConnectionManifold.jl | 22 +++++++++++----------- test/metric.jl | 20 +++++++++++++++----- 3 files changed, 30 insertions(+), 20 deletions(-) diff --git a/src/differentiation/ode.jl b/src/differentiation/ode.jl index d2fdbfdf01..847be6116b 100644 --- a/src/differentiation/ode.jl +++ b/src/differentiation/ode.jl @@ -9,8 +9,8 @@ function solve_exp_ode( kwargs..., ) d = manifold_dimension(M) - iv = SVector{n}(1:d) - ix = SVector{n}((n + 1):(2 * d)) + iv = SVector{d}(1:d) + ix = SVector{d}((d + 1):(2 * d)) u0 = allocate(p, 2 * d) u0[iv] .= X u0[ix] .= p @@ -18,10 +18,10 @@ function solve_exp_ode( function exp_problem(u, params, t) M = params[1] dx = u[iv] - x = u[ix] + q = u[ix] ddx = allocate(u, Size(d)) du = allocate(u) - Γ = christoffel_symbols_second(M, x, basis; backend=backend, retraction=retraction) + Γ = christoffel_symbols_second(M, q, basis; backend=backend, retraction=retraction) @einsum ddx[k] = -Γ[k, i, j] * dx[i] * dx[j] du[iv] .= ddx du[ix] .= dx diff --git a/src/manifolds/ConnectionManifold.jl b/src/manifolds/ConnectionManifold.jl index dd22df5a79..42d06b6b7b 100644 --- a/src/manifolds/ConnectionManifold.jl +++ b/src/manifolds/ConnectionManifold.jl @@ -77,25 +77,25 @@ struct ODEExponentialRetraction{T<:AbstractRetractionMethod,B<:AbstractBasis} <: basis::B end function ODEExponentialRetraction(r::T) where {T<:AbstractRetractionMethod} - return ODEExponentialRetraction{T,DefaultOrthonormalBasis}(r, DefaultOrthonormalBasis()) + return ODEExponentialRetraction(r, DefaultOrthonormalBasis()) end function ODEExponentialRetraction(r::T, ::CachedBasis) where {T<:AbstractRetractionMethod} - return DomainError( + throw(DomainError( r, "Cached Bases are currently not supported, since the basis has to be implemented in a surrounding of the start point as well.", - ) + )) end function ODEExponentialRetraction(r::ExponentialRetraction, ::AbstractBasis) - return DomainError( + throw(DomainError( r, "You can not use the exponential map as an inner method to solve the ode for the exponential map.", - ) + )) end function ODEExponentialRetraction(r::ExponentialRetraction, ::CachedBasis) - return DomainError( + throw(DomainError( r, "Neither the exponential map nor a Cached Basis can be used with this retraction type.", - ) + )) end @doc raw""" @@ -316,9 +316,9 @@ end kwargs..., ) -function retract(::AbstractConnectionManifold, q, p, X, r::ODEExponentialRetraction) +function retract!(M::AbstractConnectionManifold, q, p, X, r::ODEExponentialRetraction) sol = - solve_exp_ode(M, p, X, r.basis; retraction=r.retraction, dense=false, saveat=[1.0]) + solve_exp_ode(M, p, X; basis=r.basis, retraction=r.retraction, dense=false, saveat=[1.0]) return copyto!(q, sol) end @@ -357,9 +357,9 @@ provides a basis for tangent spaces around the start point `p`. using OrdinaryDiffEq ``` """ -function solve_exp_ode(M, p, X, tspan, B::AbstractBasis; kwargs...) +function solve_exp_ode(M, p, X; kwargs...) return error( - "solve_exp_ode not implemented on $(typeof(M)) for point $(typeof(p)), vector $(typeof(X)), and timespan $(typeof(tspan)). For a suitable default, enter `using OrdinaryDiffEq` on Julia 1.1 or greater.", + "solve_exp_ode not implemented on $(typeof(M)) for point $(typeof(p)), vector $(typeof(X)). For a suitable default, enter `using OrdinaryDiffEq` on Julia 1.1 or greater.", ) end diff --git a/test/metric.jl b/test/metric.jl index a1562bfe2c..fcaf0b4822 100644 --- a/test/metric.jl +++ b/test/metric.jl @@ -1,13 +1,21 @@ using FiniteDifferences, ForwardDiff using LinearAlgebra: I using StatsBase: AbstractWeights, pweights -import Manifolds: mean!, median!, InducedBasis, induced_basis, get_chart_index, connection +import Manifolds: mean!, median!, InducedBasis, induced_basis, get_chart_index, connection, retract! include("utils.jl") struct TestEuclidean{N} <: AbstractManifold{ℝ} end struct TestEuclideanMetric <: AbstractMetric end struct TestScaledEuclideanMetric <: AbstractMetric end +struct TestRetraction <: AbstractRetractionMethod end + +ManifoldsBase.default_retraction_method(::TestEuclidean) = TestRetraction() +function ManifoldsBase.default_retraction_method( + ::MetricManifold{ℝ,<:TestEuclidean,<:TestEuclideanMetric}, +) + return TestRetraction() +end Manifolds.manifold_dimension(::TestEuclidean{N}) where {N} = N function Manifolds.local_metric( @@ -36,7 +44,7 @@ function Manifolds.get_coordinates!( c, ::Any, X, - ::DefaultOrthogonalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, + ::DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, ) c .= 1 ./ [1.0:manifold_dimension(M)...] .* X return c @@ -46,7 +54,7 @@ function Manifolds.get_vector!( X, ::Any, c, - ::DefaultOrthogonalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, + ::DefaultOrthonormalgonalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, ) X .= [1.0:manifold_dimension(M)...] .* c return X @@ -56,7 +64,7 @@ function Manifolds.get_coordinates!( c, ::Any, X, - ::DefaultOrthogonalBasis, + ::DefaultOrthonormalBasis, ) c .= 1 ./ (2 .* [1.0:manifold_dimension(M)...]) .* X return c @@ -66,11 +74,13 @@ function Manifolds.get_vector!( X, ::Any, c, - ::DefaultOrthogonalBasis, + ::DefaultOrthonormalBasis, ) X .= 2 .* [1.0:manifold_dimension(M)...] .* c return X end +retract!(::TestEuclidean, q, p, X, ::TestRetraction) = p + X + struct TestSphere{N,T} <: AbstractManifold{ℝ} r::T end From 188597350b1f99ed3a76edfbe82d4fe2d2444402 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sat, 25 Sep 2021 19:44:50 +0200 Subject: [PATCH 07/20] Ran Juliaformatter 23.564 times today, just not before the last commit. --- test/metric.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/metric.jl b/test/metric.jl index fcaf0b4822..d3887a1fdc 100644 --- a/test/metric.jl +++ b/test/metric.jl @@ -1,7 +1,8 @@ using FiniteDifferences, ForwardDiff using LinearAlgebra: I using StatsBase: AbstractWeights, pweights -import Manifolds: mean!, median!, InducedBasis, induced_basis, get_chart_index, connection, retract! +import Manifolds: + mean!, median!, InducedBasis, induced_basis, get_chart_index, connection, retract! include("utils.jl") From 2181cff2888b85d2da474dfcc44b25da3308fc8d Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sat, 25 Sep 2021 21:05:34 +0200 Subject: [PATCH 08/20] trying to improve tests? --- src/manifolds/ConnectionManifold.jl | 41 +++++++++++++++++++---------- test/metric.jl | 12 +++++++-- 2 files changed, 37 insertions(+), 16 deletions(-) diff --git a/src/manifolds/ConnectionManifold.jl b/src/manifolds/ConnectionManifold.jl index 42d06b6b7b..bb2db85ab6 100644 --- a/src/manifolds/ConnectionManifold.jl +++ b/src/manifolds/ConnectionManifold.jl @@ -80,22 +80,28 @@ function ODEExponentialRetraction(r::T) where {T<:AbstractRetractionMethod} return ODEExponentialRetraction(r, DefaultOrthonormalBasis()) end function ODEExponentialRetraction(r::T, ::CachedBasis) where {T<:AbstractRetractionMethod} - throw(DomainError( - r, - "Cached Bases are currently not supported, since the basis has to be implemented in a surrounding of the start point as well.", - )) + return throw( + DomainError( + r, + "Cached Bases are currently not supported, since the basis has to be implemented in a surrounding of the start point as well.", + ), + ) end function ODEExponentialRetraction(r::ExponentialRetraction, ::AbstractBasis) - throw(DomainError( - r, - "You can not use the exponential map as an inner method to solve the ode for the exponential map.", - )) + return throw( + DomainError( + r, + "You can not use the exponential map as an inner method to solve the ode for the exponential map.", + ), + ) end function ODEExponentialRetraction(r::ExponentialRetraction, ::CachedBasis) - throw(DomainError( - r, - "Neither the exponential map nor a Cached Basis can be used with this retraction type.", - )) + return throw( + DomainError( + r, + "Neither the exponential map nor a Cached Basis can be used with this retraction type.", + ), + ) end @doc raw""" @@ -317,8 +323,15 @@ end ) function retract!(M::AbstractConnectionManifold, q, p, X, r::ODEExponentialRetraction) - sol = - solve_exp_ode(M, p, X; basis=r.basis, retraction=r.retraction, dense=false, saveat=[1.0]) + sol = solve_exp_ode( + M, + p, + X; + basis=r.basis, + retraction=r.retraction, + dense=false, + saveat=[1.0], + ) return copyto!(q, sol) end diff --git a/test/metric.jl b/test/metric.jl index d3887a1fdc..d145c42c97 100644 --- a/test/metric.jl +++ b/test/metric.jl @@ -55,7 +55,7 @@ function Manifolds.get_vector!( X, ::Any, c, - ::DefaultOrthonormalgonalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, + ::DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, ) X .= [1.0:manifold_dimension(M)...] .* c return X @@ -80,7 +80,15 @@ function Manifolds.get_vector!( X .= 2 .* [1.0:manifold_dimension(M)...] .* c return X end -retract!(::TestEuclidean, q, p, X, ::TestRetraction) = p + X +function retract!( + ::MetricManifold{ℝ,<:TestEuclidean,<:TestEuclideanMetric}, + q, + p, + X, + ::TestRetraction, +) + return p + X +end struct TestSphere{N,T} <: AbstractManifold{ℝ} r::T From cde3c8ebbe5e74b373b84d915dc44011d63c8723 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sun, 26 Sep 2021 16:59:36 +0200 Subject: [PATCH 09/20] Fixes most errors, just some default retractions need to be fixed (still returning ExponentialRetraction). --- test/metric.jl | 31 +++++++++++++++++++++++++++---- 1 file changed, 27 insertions(+), 4 deletions(-) diff --git a/test/metric.jl b/test/metric.jl index d145c42c97..42d5e8c6ca 100644 --- a/test/metric.jl +++ b/test/metric.jl @@ -55,7 +55,10 @@ function Manifolds.get_vector!( X, ::Any, c, - ::DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, + ::Union{ + DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, + InducedBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, + }, ) X .= [1.0:manifold_dimension(M)...] .* c return X @@ -87,7 +90,8 @@ function retract!( X, ::TestRetraction, ) - return p + X + copyto!(q, p + X) + return q end struct TestSphere{N,T} <: AbstractManifold{ℝ} @@ -96,6 +100,20 @@ end struct TestSphericalMetric <: AbstractMetric end +function retract!( + ::MetricManifold{ℝ,<:TestSphere,<:TestSphericalMetric}, + q, + p, + X, + ::ProjectionRetraction, +) + copyto(q, (p + X) ./ norm(p + X)) + return q +end +function default_retraction_method(::MetricManifold{ℝ,<:TestSphere,<:TestSphericalMetric}) + return ProjectionRetraction() +end + Manifolds.manifold_dimension(::TestSphere{N}) where {N} = N function Manifolds.local_metric( M::MetricManifold{ℝ,<:TestSphere,<:TestSphericalMetric}, @@ -275,8 +293,13 @@ end i_zeros = get_chart_index(M, A, zeros(3)) B_i_zeros = induced_basis(M, A, i_zeros, TangentSpace) - @test_throws MethodError local_metric_jacobian(E, zeros(3), B_i_zeros) - @test_throws MethodError christoffel_symbols_second_jacobian(E, zeros(3), B_i_zeros) + # get_vector! not implemented -> ErrorException + @test_throws ErrorException local_metric_jacobian(E, zeros(3), B_i_zeros) + @test_throws ErrorException christoffel_symbols_second_jacobian( + E, + zeros(3), + B_i_zeros, + ) for vtype in (Vector, MVector{n}) p, X, Y = vtype(randn(n)), vtype(randn(n)), vtype(randn(n)) From e1bf764bc4ec6abf65e4e5580c214737f705264e Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sun, 26 Sep 2021 17:25:45 +0200 Subject: [PATCH 10/20] reduce errors further. --- test/metric.jl | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/test/metric.jl b/test/metric.jl index 42d5e8c6ca..28ab54522f 100644 --- a/test/metric.jl +++ b/test/metric.jl @@ -3,6 +3,7 @@ using LinearAlgebra: I using StatsBase: AbstractWeights, pweights import Manifolds: mean!, median!, InducedBasis, induced_basis, get_chart_index, connection, retract! +import ManifoldsBase: default_retraction_method include("utils.jl") @@ -45,7 +46,10 @@ function Manifolds.get_coordinates!( c, ::Any, X, - ::DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, + ::Union{ + DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, + DefaultOrthogonalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, + }, ) c .= 1 ./ [1.0:manifold_dimension(M)...] .* X return c @@ -57,6 +61,7 @@ function Manifolds.get_vector!( c, ::Union{ DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, + DefaultOrthogonalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, InducedBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, }, ) @@ -110,7 +115,9 @@ function retract!( copyto(q, (p + X) ./ norm(p + X)) return q end -function default_retraction_method(::MetricManifold{ℝ,<:TestSphere,<:TestSphericalMetric}) +function ManifoldsBase.default_retraction_method( + ::MetricManifold{ℝ,<:TestSphere,<:TestSphericalMetric}, +) return ProjectionRetraction() end From b1124999ff6d4e06f8dd582f73e33c186456b766 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Mon, 27 Sep 2021 11:25:08 +0200 Subject: [PATCH 11/20] reduices errors as far as I am able to do. --- test/metric.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/metric.jl b/test/metric.jl index 28ab54522f..db8e68fe84 100644 --- a/test/metric.jl +++ b/test/metric.jl @@ -549,8 +549,9 @@ end chart_p = get_chart_index(MM2, A, p) B_p = induced_basis(MM2, A, chart_p, TangentSpace) @test_throws MethodError local_metric(MM2, p, B_p) - @test_throws MethodError local_metric_jacobian(MM2, p, B_p) - @test_throws MethodError christoffel_symbols_second_jacobian(MM2, p, B_p) + # the following two default to solving the ODE, but fail since get_vector is not implemented + @test_throws ErrorException local_metric_jacobian(MM2, p, B_p) + @test_throws ErrorException christoffel_symbols_second_jacobian(MM2, p, B_p) # MM falls back to nondefault error @test_throws MethodError projected_distribution(MM, 1, p) @test_throws MethodError projected_distribution(MM, 1) From f2625fb8534fdf18dc1b5aad9b531f0533ae0092 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Mon, 27 Sep 2021 22:47:47 +0200 Subject: [PATCH 12/20] =?UTF-8?q?Implements=20my=20understanding=20of=20wh?= =?UTF-8?q?at=20the=20methods=20should=20be=20=20=E2=80=93=20error=20persi?= =?UTF-8?q?sts.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- test/metric.jl | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/test/metric.jl b/test/metric.jl index db8e68fe84..d691fdef07 100644 --- a/test/metric.jl +++ b/test/metric.jl @@ -112,7 +112,7 @@ function retract!( X, ::ProjectionRetraction, ) - copyto(q, (p + X) ./ norm(p + X)) + copyto!(q, (p + X) ./ norm(p + X)) return q end function ManifoldsBase.default_retraction_method( @@ -133,6 +133,26 @@ function Manifolds.local_metric( d[2] = d[1] * sin(p[1])^2 return Diagonal(d) end +function Manifolds.inverse_local_metric( + ::MetricManifold{ℝ,<:TestSphere{N},<:TestSphericalMetric}, + p, + ::DefaultOrthonormalBasis, +) where {N} + return Diagonal(ones(N)) +end +function Manifolds.get_vector!( + ::MetricManifold{ℝ,<:TestSphere{N},<:TestSphericalMetric}, + Y, + p, + c, + ::Union{ + DefaultOrthonormalBasis{ℝ, <:ManifoldsBase.TangentSpaceType}, + InducedBasis{ℝ, <:ManifoldsBase.TangentSpaceType}, + } +) where {N} + Y .= 1 + return Y # this is just a dummy to check that dispatch works +end sph_to_cart(θ, ϕ) = [cos(ϕ) * sin(θ), sin(ϕ) * sin(θ), cos(θ)] struct BaseManifold{N} <: AbstractManifold{ℝ} end From a6702f62f46a043d2ce489e36f1d6b2c3b6cdb5b Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Mon, 27 Sep 2021 23:11:28 +0200 Subject: [PATCH 13/20] fix some tests --- test/metric.jl | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/test/metric.jl b/test/metric.jl index d691fdef07..c5d81f04e8 100644 --- a/test/metric.jl +++ b/test/metric.jl @@ -125,7 +125,7 @@ Manifolds.manifold_dimension(::TestSphere{N}) where {N} = N function Manifolds.local_metric( M::MetricManifold{ℝ,<:TestSphere,<:TestSphericalMetric}, p, - ::InducedBasis, + ::Union{InducedBasis,DefaultOrthonormalBasis}, ) r = base_manifold(M).r d = allocate(p) @@ -133,21 +133,14 @@ function Manifolds.local_metric( d[2] = d[1] * sin(p[1])^2 return Diagonal(d) end -function Manifolds.inverse_local_metric( - ::MetricManifold{ℝ,<:TestSphere{N},<:TestSphericalMetric}, - p, - ::DefaultOrthonormalBasis, -) where {N} - return Diagonal(ones(N)) -end function Manifolds.get_vector!( ::MetricManifold{ℝ,<:TestSphere{N},<:TestSphericalMetric}, Y, p, c, ::Union{ - DefaultOrthonormalBasis{ℝ, <:ManifoldsBase.TangentSpaceType}, - InducedBasis{ℝ, <:ManifoldsBase.TangentSpaceType}, + DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, + InducedBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, } ) where {N} Y .= 1 From a4a2716c55ca587dd11690adc6404bca3e07f182 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Tue, 28 Sep 2021 14:00:49 +0200 Subject: [PATCH 14/20] =?UTF-8?q?starts=20seperate=20ODE=20tests=20?= =?UTF-8?q?=E2=80=93=C2=A0but=20facing=20a=20new=20problem.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- test/metric.jl | 2 +- test/ode.jl | 26 ++++++++++++++++++++++++++ test/runtests.jl | 1 + 3 files changed, 28 insertions(+), 1 deletion(-) create mode 100644 test/ode.jl diff --git a/test/metric.jl b/test/metric.jl index c5d81f04e8..212a71f9c0 100644 --- a/test/metric.jl +++ b/test/metric.jl @@ -141,7 +141,7 @@ function Manifolds.get_vector!( ::Union{ DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, InducedBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, - } + }, ) where {N} Y .= 1 return Y # this is just a dummy to check that dispatch works diff --git a/test/ode.jl b/test/ode.jl new file mode 100644 index 0000000000..270366fd81 --- /dev/null +++ b/test/ode.jl @@ -0,0 +1,26 @@ +include("utils.jl") + +struct TestSphere{N} <: AbstractManifold{ℝ} end +struct TestSphericalMetric <: AbstractMetric end + +using FiniteDifferences, ForwardDiff +using LinearAlgebra: I +import Manifolds: retract! +import ManifoldsBase: manifold_dimension, default_retraction_method + +function default_retraction_method(::MetricManifold{ℝ,<:TestSphere,<:TestSphericalMetric}) + return ProjectionRetraction() +end +manifold_dimension(::TestSphere{n}) where {n} = n +retract!(::TestSphere{n}, q, p, X) where {n} = retract!(Sphere(n), q, p, X) + +@testset "Test ODE setup for computing geodesics" begin + M = TestSphere{2}() + p = [0.0, 0.0, 1.0] + X = π / (2 * sqrt(2)) .* [0.0, 1.0, 1.0] + M2 = MetricManifold(M, TestSphericalMetric()) + @test_throws ErrorException exp(M, p, X) + @test_throws ErrorException exp(M2, p, X) + using OrdinaryDiffEq + exp(M2, p, X) +end diff --git a/test/runtests.jl b/test/runtests.jl index 25e841c377..27115de4e4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -147,6 +147,7 @@ include("utils.jl") include_test("manifolds/graph.jl") include_test("metric.jl") + include_test("ode.jl") include_test("statistics.jl") include_test("approx_inverse_retraction.jl") From 5ed5feb6612e706964bfb68b60c1ab861a299cc2 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Tue, 28 Sep 2021 14:17:18 +0200 Subject: [PATCH 15/20] getting again stuck on metrics. --- src/differentiation/ode.jl | 3 ++- test/metric.jl | 11 ----------- test/ode.jl | 8 +++++++- 3 files changed, 9 insertions(+), 13 deletions(-) diff --git a/src/differentiation/ode.jl b/src/differentiation/ode.jl index 847be6116b..2cf58df370 100644 --- a/src/differentiation/ode.jl +++ b/src/differentiation/ode.jl @@ -8,7 +8,8 @@ function solve_exp_ode( retraction::AbstractRetractionMethod=ManifoldsBase.default_retraction_method(M), kwargs..., ) - d = manifold_dimension(M) + # d = manifold_dimension(M) + d = length(p) iv = SVector{d}(1:d) ix = SVector{d}((d + 1):(2 * d)) u0 = allocate(p, 2 * d) diff --git a/test/metric.jl b/test/metric.jl index 212a71f9c0..6c84a2d527 100644 --- a/test/metric.jl +++ b/test/metric.jl @@ -273,17 +273,6 @@ end @test length(methods(is_default_metric)) == 2 end - @testset "solve_exp_ode error message" begin - E = TestEuclidean{3}() - g = TestEuclideanMetric() - M = MetricManifold(E, g) - - p = [1.0, 2.0, 3.0] - X = [2.0, 3.0, 4.0] - @test_throws ErrorException exp(M, p, X) - using OrdinaryDiffEq - exp(M, p, X) - end @testset "Local Metric Error message" begin M = MetricManifold(BaseManifold{2}(), NotImplementedMetric()) A = Manifolds.get_default_atlas(M) diff --git a/test/ode.jl b/test/ode.jl index 270366fd81..3000db8b10 100644 --- a/test/ode.jl +++ b/test/ode.jl @@ -13,7 +13,13 @@ function default_retraction_method(::MetricManifold{ℝ,<:TestSphere,<:TestSpher end manifold_dimension(::TestSphere{n}) where {n} = n retract!(::TestSphere{n}, q, p, X) where {n} = retract!(Sphere(n), q, p, X) - +function local_metric( + ::TestSphere{n}, + p, + B::DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, +) where {n} + return local_metric(Sphere(n), p, B) +end @testset "Test ODE setup for computing geodesics" begin M = TestSphere{2}() p = [0.0, 0.0, 1.0] From d544ce9fccfc2575602c89edf4b6f56c586cecae Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Tue, 28 Sep 2021 14:37:01 +0200 Subject: [PATCH 16/20] fixing ode tests a bit --- test/ode.jl | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/test/ode.jl b/test/ode.jl index 3000db8b10..33c5e87693 100644 --- a/test/ode.jl +++ b/test/ode.jl @@ -12,13 +12,31 @@ function default_retraction_method(::MetricManifold{ℝ,<:TestSphere,<:TestSpher return ProjectionRetraction() end manifold_dimension(::TestSphere{n}) where {n} = n -retract!(::TestSphere{n}, q, p, X) where {n} = retract!(Sphere(n), q, p, X) -function local_metric( - ::TestSphere{n}, +function Manifolds.retract!( + ::MetricManifold{ℝ,TestSphere{n},<:TestSphericalMetric}, + q, + p, + X, + ::ProjectionRetraction, +) where {n} + return retract!(Sphere(n), q, p, X) +end +function Manifolds.local_metric( + ::MetricManifold{ℝ,TestSphere{n},<:TestSphericalMetric}, p, B::DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, ) where {n} - return local_metric(Sphere(n), p, B) + return Diagonal(ones(n)) # TODO: fix? +end +function Manifolds.get_vector!( + ::MetricManifold{ℝ,<:TestSphere{N},<:TestSphericalMetric}, + Y, + p, + c, + ::DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, +) where {N} + Y .= 1 + return Y # this is just a dummy to check that dispatch works end @testset "Test ODE setup for computing geodesics" begin M = TestSphere{2}() From 26ab2ec3eaaa331e6e405e9ac1fc24e78750af48 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Tue, 28 Sep 2021 17:42:39 +0200 Subject: [PATCH 17/20] Fighting dimensions within the ODE solver. --- src/manifolds/Sphere.jl | 11 +++++++++++ test/ode.jl | 27 +++++++++++++-------------- 2 files changed, 24 insertions(+), 14 deletions(-) diff --git a/src/manifolds/Sphere.jl b/src/manifolds/Sphere.jl index 92f07d1a25..669e5886a0 100644 --- a/src/manifolds/Sphere.jl +++ b/src/manifolds/Sphere.jl @@ -308,6 +308,17 @@ function inverse_retract!(::AbstractSphere, X, p, q, ::ProjectionInverseRetracti return (X .= q ./ real(dot(p, q)) .- p) end +@doc raw""" + local_metric(M::Sphere{n}, p, ::DefaultOrthonormalBasis) + +return the local representation of the metric in a [`DefaultOrthonormalBasis`](@ref), namely +the diagonal matrix of size ``n×n`` with ones on the diagonal, since the metric is obtained +from the embedding by restriction to the tangent space ``T_p\mathcal M`` at ``p``. +""" +function local_metric(::Sphere{n,ℝ}, p, B::DefaultOrthonormalBasis) where {n} + return Diagonal(ones(SVector{n,eltype(p)})) +end + @doc raw""" log(M::AbstractSphere, p, q) diff --git a/test/ode.jl b/test/ode.jl index 33c5e87693..e9d600abfc 100644 --- a/test/ode.jl +++ b/test/ode.jl @@ -26,25 +26,24 @@ function Manifolds.local_metric( p, B::DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, ) where {n} - return Diagonal(ones(n)) # TODO: fix? + return Manifolds.local_metric(MetricManifold(Sphere(n), EuclideanMetric()), p, B) end function Manifolds.get_vector!( ::MetricManifold{ℝ,<:TestSphere{N},<:TestSphericalMetric}, Y, p, c, - ::DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, + B::DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, ) where {N} - Y .= 1 - return Y # this is just a dummy to check that dispatch works -end -@testset "Test ODE setup for computing geodesics" begin - M = TestSphere{2}() - p = [0.0, 0.0, 1.0] - X = π / (2 * sqrt(2)) .* [0.0, 1.0, 1.0] - M2 = MetricManifold(M, TestSphericalMetric()) - @test_throws ErrorException exp(M, p, X) - @test_throws ErrorException exp(M2, p, X) - using OrdinaryDiffEq - exp(M2, p, X) + return get_vector!(Sphere(N), Y, p, c, B) end +#@testset "Test ODE setup for computing geodesics" begin +M = TestSphere{2}() +p = [0.0, 0.0, 1.0] +X = π / (2 * sqrt(2)) .* [0.0, 1.0, 1.0] +M2 = MetricManifold(M, TestSphericalMetric()) +# @test_throws ErrorException exp(M, p, X) +# @test_throws ErrorException exp(M2, p, X) +using OrdinaryDiffEq +exp(M2, p, X) +#end From bd39060c2f3f2f60a661e0d65501f4f4c6156756 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Tue, 28 Sep 2021 19:07:05 +0200 Subject: [PATCH 18/20] Dimensions fixed, now just the result is wrong :D --- src/differentiation/ode.jl | 18 +++++++++--------- src/manifolds/MetricManifold.jl | 2 +- test/ode.jl | 9 +++++++++ 3 files changed, 19 insertions(+), 10 deletions(-) diff --git a/src/differentiation/ode.jl b/src/differentiation/ode.jl index 2cf58df370..89ec7d14db 100644 --- a/src/differentiation/ode.jl +++ b/src/differentiation/ode.jl @@ -8,29 +8,29 @@ function solve_exp_ode( retraction::AbstractRetractionMethod=ManifoldsBase.default_retraction_method(M), kwargs..., ) - # d = manifold_dimension(M) - d = length(p) - iv = SVector{d}(1:d) - ix = SVector{d}((d + 1):(2 * d)) - u0 = allocate(p, 2 * d) + d = manifold_dimension(M) + n = length(p) + iv = SVector{n}(1:n) + ix = SVector{n}((n + 1):(2 * n)) + u0 = allocate(p, 2 * n) u0[iv] .= X u0[ix] .= p function exp_problem(u, params, t) M = params[1] - dx = u[iv] q = u[ix] + dx = get_coordinates(M, q, u[iv], basis) ddx = allocate(u, Size(d)) du = allocate(u) Γ = christoffel_symbols_second(M, q, basis; backend=backend, retraction=retraction) @einsum ddx[k] = -Γ[k, i, j] * dx[i] * dx[j] - du[iv] .= ddx - du[ix] .= dx + du[iv] .= get_vector(M, q, ddx, basis) + du[ix] .= u[iv] return Base.convert(typeof(u), du) end params = (M,) prob = ODEProblem(exp_problem, u0, (0.0, 1.0), params) sol = solve(prob, solver; kwargs...) - q = sol.u[1][(d + 1):(2 * d)] + q = sol.u[1][(n + 1):(2 * n)] return q end diff --git a/src/manifolds/MetricManifold.jl b/src/manifolds/MetricManifold.jl index 56f6cdec86..95bacabafd 100644 --- a/src/manifolds/MetricManifold.jl +++ b/src/manifolds/MetricManifold.jl @@ -473,7 +473,7 @@ function local_metric_jacobian( ∂g = reshape( _jacobian( c -> local_metric(M, retract(M, p, get_vector(M, p, c, B), retraction), B), - p, + zeros(d), backend, ), d, diff --git a/test/ode.jl b/test/ode.jl index e9d600abfc..4fc61b9c4e 100644 --- a/test/ode.jl +++ b/test/ode.jl @@ -28,6 +28,15 @@ function Manifolds.local_metric( ) where {n} return Manifolds.local_metric(MetricManifold(Sphere(n), EuclideanMetric()), p, B) end +function Manifolds.get_coordinates!( + ::MetricManifold{ℝ,<:TestSphere{N},<:TestSphericalMetric}, + c, + p, + X, + B::DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, +) where {N} + return get_coordinates!(Sphere(N), c, p, X, B) +end function Manifolds.get_vector!( ::MetricManifold{ℝ,<:TestSphere{N},<:TestSphericalMetric}, Y, From 44ceb3e5ff689b2e7c86f86ff091296aed2895c4 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Wed, 29 Sep 2021 18:51:05 +0200 Subject: [PATCH 19/20] reorganise to get the ODE stuff (a) defined on arbitrary manifolds and (b) into its own test file (away from metric). --- src/differentiation/differentiation.jl | 69 ++++++++++++++++++++++++++ src/differentiation/ode.jl | 12 ++--- src/manifolds/ConnectionManifold.jl | 65 +----------------------- test/metric.jl | 17 ------- test/ode.jl | 41 +++++++++------ 5 files changed, 102 insertions(+), 102 deletions(-) diff --git a/src/differentiation/differentiation.jl b/src/differentiation/differentiation.jl index 67e64ca1db..ff198337bf 100644 --- a/src/differentiation/differentiation.jl +++ b/src/differentiation/differentiation.jl @@ -115,3 +115,72 @@ function set_default_differential_backend!(backend::AbstractDiffBackend) _current_default_differential_backend.backend = backend return backend end + +@doc raw""" + ODEExponentialRetraction{T<:AbstractRetractionMethod, B<: AbstractBasis} <: AbstractRetractionMethod + +This retraction approximates the exponential map by solving the correspondig ODE. +Let ``p\in\mathal M`` and ``X\in T_p\mathcal M`` denote the input for the exponential map +and ``d`` denote the [`manifold_dimension`](@ref) of `M`. + +This the ODE is formulated in a chart constructed using an [`AbstractBasis`](@ref) `B` and an +[`AbstractRetractionMethod`](@ref) `R` as follows. +Given some coordinates ``c\in ℝ^d`` - these can be used to form a tangent vector +with restect to th basis `B, i.e. ``c \mapsto Y=``[`get_vector`](@ref)`(M, p, c, B)`. +Further, using the retraction we can map ``Y`` to a point on the manifold +``Y \mapsto q =``[`retract`](@ref)`(M, p, X, R)`. + +Hence the ODE can be formulated in a curve ``c(t)`` in parameter space. +This is – for sure – only possible locally as fas as the retraction is well-defined. +""" +struct ODEExponentialRetraction{T<:AbstractRetractionMethod,B<:AbstractBasis} <: + AbstractRetractionMethod + retraction::T + basis::B +end +function ODEExponentialRetraction(r::T) where {T<:AbstractRetractionMethod} + return ODEExponentialRetraction(r, DefaultOrthonormalBasis()) +end +function ODEExponentialRetraction(r::T, ::CachedBasis) where {T<:AbstractRetractionMethod} + return throw( + DomainError( + r, + "Cached Bases are currently not supported, since the basis has to be implemented in a surrounding of the start point as well.", + ), + ) +end +function ODEExponentialRetraction(r::ExponentialRetraction, ::AbstractBasis) + return throw( + DomainError( + r, + "You can not use the exponential map as an inner method to solve the ode for the exponential map.", + ), + ) +end +function ODEExponentialRetraction(r::ExponentialRetraction, ::CachedBasis) + return throw( + DomainError( + r, + "Neither the exponential map nor a Cached Basis can be used with this retraction type.", + ), + ) +end + +function retract!( + M::AbstractManifold, + q, + p, + X, + r::ODEExponentialRetraction, +) + sol = solve_exp_ode( + M, + p, + X; + basis=r.basis, + retraction=r.retraction, + dense=false, + saveat=[1.0], + ) + return copyto!(q, sol) +end diff --git a/src/differentiation/ode.jl b/src/differentiation/ode.jl index 89ec7d14db..820cfa41f9 100644 --- a/src/differentiation/ode.jl +++ b/src/differentiation/ode.jl @@ -1,5 +1,5 @@ function solve_exp_ode( - M::AbstractConnectionManifold, + M::AbstractManifold, p, X; basis::AbstractBasis=DefaultOrthonormalBasis(), @@ -9,11 +9,10 @@ function solve_exp_ode( kwargs..., ) d = manifold_dimension(M) - n = length(p) - iv = SVector{n}(1:n) - ix = SVector{n}((n + 1):(2 * n)) - u0 = allocate(p, 2 * n) - u0[iv] .= X + iv = SVector{n}(1:d) + ix = SVector{n}((d + 1):(2 * d)) + u0 = allocate(p, 2 * d) + u0[iv] .= get_get_coordinates!(M, u0[iv]) u0[ix] .= p function exp_problem(u, params, t) @@ -26,6 +25,7 @@ function solve_exp_ode( @einsum ddx[k] = -Γ[k, i, j] * dx[i] * dx[j] du[iv] .= get_vector(M, q, ddx, basis) du[ix] .= u[iv] + 1 / 0 return Base.convert(typeof(u), du) end params = (M,) diff --git a/src/manifolds/ConnectionManifold.jl b/src/manifolds/ConnectionManifold.jl index bb2db85ab6..57d9c39bcc 100644 --- a/src/manifolds/ConnectionManifold.jl +++ b/src/manifolds/ConnectionManifold.jl @@ -54,56 +54,6 @@ struct ConnectionManifold{𝔽,M<:AbstractManifold{𝔽},C<:AbstractAffineConnec connection::C end -@doc raw""" - ODEExponentialRetraction{T<:AbstractRetractionMethod, B<: AbstractBasis} <: AbstractRetractionMethod - -This retraction approximates the exponential map by solving the correspondig ODE. -Let ``p\in\mathal M`` and ``X\in T_p\mathcal M`` denote the input for the exponential map -and ``d`` denote the [`manifold_dimension`](@ref) of `M`. - -This the ODE is formulated in a chart constructed using an [`AbstractBasis`](@ref) `B` and an -[`AbstractRetractionMethod`](@ref) `R` as follows. -Given some coordinates ``c\in ℝ^d`` - these can be used to form a tangent vector -with restect to th basis `B, i.e. ``c \mapsto Y=``[`get_vector`](@ref)`(M, p, c, B)`. -Further, using the retraction we can map ``Y`` to a point on the manifold -``Y \mapsto q =``[`retract`](@ref)`(M, p, X, R)`. - -Hence the ODE can be formulated in a curve ``c(t)`` in parameter space. -This is – for sure – only possible locally as fas as the retraction is well-defined. -""" -struct ODEExponentialRetraction{T<:AbstractRetractionMethod,B<:AbstractBasis} <: - AbstractRetractionMethod - retraction::T - basis::B -end -function ODEExponentialRetraction(r::T) where {T<:AbstractRetractionMethod} - return ODEExponentialRetraction(r, DefaultOrthonormalBasis()) -end -function ODEExponentialRetraction(r::T, ::CachedBasis) where {T<:AbstractRetractionMethod} - return throw( - DomainError( - r, - "Cached Bases are currently not supported, since the basis has to be implemented in a surrounding of the start point as well.", - ), - ) -end -function ODEExponentialRetraction(r::ExponentialRetraction, ::AbstractBasis) - return throw( - DomainError( - r, - "You can not use the exponential map as an inner method to solve the ode for the exponential map.", - ), - ) -end -function ODEExponentialRetraction(r::ExponentialRetraction, ::CachedBasis) - return throw( - DomainError( - r, - "Neither the exponential map nor a Cached Basis can be used with this retraction type.", - ), - ) -end - @doc raw""" christoffel_symbols_second( M::AbstractManifold, @@ -322,22 +272,9 @@ end kwargs..., ) -function retract!(M::AbstractConnectionManifold, q, p, X, r::ODEExponentialRetraction) - sol = solve_exp_ode( - M, - p, - X; - basis=r.basis, - retraction=r.retraction, - dense=false, - saveat=[1.0], - ) - return copyto!(q, sol) -end - @doc raw""" solve_exp_ode( - M::AbstractConnectionManifold, + M::AbstractManifold, p, X, B::AbstractBasis; diff --git a/test/metric.jl b/test/metric.jl index 6c84a2d527..a08101da42 100644 --- a/test/metric.jl +++ b/test/metric.jl @@ -329,11 +329,6 @@ end @test inner(M, p, fX, fY) ≈ dot(X, G * Y) atol = 1e-6 @test norm(M, p, fX) ≈ sqrt(dot(X, G * X)) atol = 1e-6 - if VERSION ≥ v"1.1" - T = 0:0.5:10 - @test geodesic(M, p, X, T) ≈ [p + t * X for t in T] atol = 1e-6 - end - @test christoffel_symbols_first(M, p, B_chart_p) ≈ zeros(n, n, n) atol = 1e-6 @test christoffel_symbols_second(M, p, B_chart_p) ≈ zeros(n, n, n) atol = 1e-6 @test riemann_tensor(M, p, B_chart_p) ≈ zeros(n, n, n, n) atol = 1e-6 @@ -407,18 +402,6 @@ end -sin(θ) 0 ] * X - if !Sys.iswindows() || Sys.ARCH == :x86_64 - @testset "numerically integrated geodesics for $vtype" begin - T = 0:0.1:1 - @test isapprox( - [sph_to_cart(yi...) for yi in geodesic(M, p, X, T)], - geodesic(S, pcart, Xcart, T); - atol=1e-3, - rtol=1e-3, - ) - end - end - Γ₁ = christoffel_symbols_first(M, p, B_p) for i in 1:n, j in 1:n, k in 1:n if (i, j, k) == (1, 2, 2) || (i, j, k) == (2, 1, 2) diff --git a/test/ode.jl b/test/ode.jl index 4fc61b9c4e..49897336a8 100644 --- a/test/ode.jl +++ b/test/ode.jl @@ -1,35 +1,46 @@ include("utils.jl") -struct TestSphere{N} <: AbstractManifold{ℝ} end -struct TestSphericalMetric <: AbstractMetric end - using FiniteDifferences, ForwardDiff using LinearAlgebra: I import Manifolds: retract! import ManifoldsBase: manifold_dimension, default_retraction_method -function default_retraction_method(::MetricManifold{ℝ,<:TestSphere,<:TestSphericalMetric}) +# +# Part I: Euclidean +# +struct TestODEEuclidean{N} <: AbstractManifold{ℝ} end +struct TestODEEuclideanMetric <: AbstractMetric end +# +# Part II Spherical +# +struct TestODESphere{N} <: AbstractManifold{ℝ} end +struct TestODESphericalMetric <: AbstractMetric end + +function default_retraction_method( + ::MetricManifold{ℝ,<:TestODESphere,<:TestODESphericalMetric}, +) return ProjectionRetraction() end -manifold_dimension(::TestSphere{n}) where {n} = n + +manifold_dimension(::TestODESphere{N}) where {N} = N function Manifolds.retract!( - ::MetricManifold{ℝ,TestSphere{n},<:TestSphericalMetric}, + ::MetricManifold{ℝ,<:TestODESphere{N},<:TestODESphericalMetric}, q, p, X, ::ProjectionRetraction, -) where {n} - return retract!(Sphere(n), q, p, X) +) where {N} + return retract!(Sphere(N), q, p, X) end function Manifolds.local_metric( - ::MetricManifold{ℝ,TestSphere{n},<:TestSphericalMetric}, + ::MetricManifold{ℝ,<:TestODESphere{N},<:TestODESphericalMetric}, p, B::DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, -) where {n} - return Manifolds.local_metric(MetricManifold(Sphere(n), EuclideanMetric()), p, B) +) where {N} + return Manifolds.local_metric(MetricManifold(Sphere(N), EuclideanMetric()), p, B) end function Manifolds.get_coordinates!( - ::MetricManifold{ℝ,<:TestSphere{N},<:TestSphericalMetric}, + ::MetricManifold{ℝ,<:TestODESphere{N},<:TestODESphericalMetric}, c, p, X, @@ -38,7 +49,7 @@ function Manifolds.get_coordinates!( return get_coordinates!(Sphere(N), c, p, X, B) end function Manifolds.get_vector!( - ::MetricManifold{ℝ,<:TestSphere{N},<:TestSphericalMetric}, + ::MetricManifold{ℝ,<:TestODESphere{N},<:TestODESphericalMetric}, Y, p, c, @@ -47,10 +58,10 @@ function Manifolds.get_vector!( return get_vector!(Sphere(N), Y, p, c, B) end #@testset "Test ODE setup for computing geodesics" begin -M = TestSphere{2}() +M = TestODESphere{2}() p = [0.0, 0.0, 1.0] X = π / (2 * sqrt(2)) .* [0.0, 1.0, 1.0] -M2 = MetricManifold(M, TestSphericalMetric()) +M2 = MetricManifold(M, TestODESphericalMetric()) # @test_throws ErrorException exp(M, p, X) # @test_throws ErrorException exp(M2, p, X) using OrdinaryDiffEq From 6ac5d7555618b42e6bb0a9dfc9e817dccc8d7a52 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Fri, 1 Oct 2021 20:10:57 +0200 Subject: [PATCH 20/20] Revert to an embedded approach. Maybe. --- src/differentiation/differentiation.jl | 8 +------- src/differentiation/ode.jl | 17 ++++++++--------- 2 files changed, 9 insertions(+), 16 deletions(-) diff --git a/src/differentiation/differentiation.jl b/src/differentiation/differentiation.jl index ff198337bf..45d12c177a 100644 --- a/src/differentiation/differentiation.jl +++ b/src/differentiation/differentiation.jl @@ -166,13 +166,7 @@ function ODEExponentialRetraction(r::ExponentialRetraction, ::CachedBasis) ) end -function retract!( - M::AbstractManifold, - q, - p, - X, - r::ODEExponentialRetraction, -) +function retract!(M::AbstractManifold, q, p, X, r::ODEExponentialRetraction) sol = solve_exp_ode( M, p, diff --git a/src/differentiation/ode.jl b/src/differentiation/ode.jl index 820cfa41f9..4b6dcdb412 100644 --- a/src/differentiation/ode.jl +++ b/src/differentiation/ode.jl @@ -8,24 +8,23 @@ function solve_exp_ode( retraction::AbstractRetractionMethod=ManifoldsBase.default_retraction_method(M), kwargs..., ) - d = manifold_dimension(M) - iv = SVector{n}(1:d) - ix = SVector{n}((d + 1):(2 * d)) - u0 = allocate(p, 2 * d) + n = length(p) + iv = SVector{n}(1:n) + ix = SVector{n}((n + 1):(2 * n)) + u0 = allocate(p, 2 * n) u0[iv] .= get_get_coordinates!(M, u0[iv]) u0[ix] .= p function exp_problem(u, params, t) M = params[1] q = u[ix] - dx = get_coordinates(M, q, u[iv], basis) - ddx = allocate(u, Size(d)) + dx = u[iv] + ddx = allocate(u, Size(n)) du = allocate(u) Γ = christoffel_symbols_second(M, q, basis; backend=backend, retraction=retraction) @einsum ddx[k] = -Γ[k, i, j] * dx[i] * dx[j] - du[iv] .= get_vector(M, q, ddx, basis) - du[ix] .= u[iv] - 1 / 0 + du[iv] .= ddx + du[ix] .= dx return Base.convert(typeof(u), du) end params = (M,)