-- |
-- Module      : Data.Manifold.Riemannian
-- Copyright   : (c) Justus Sagemüller 2015
-- License     : GPL v3
-- 
-- Maintainer  : (@) jsag $ hvl.no
-- Stability   : experimental
-- Portability : portable
-- 
-- Riemannian manifolds are manifolds equipped with a 'Metric' at each point.
-- That means, these manifolds aren't merely topological objects anymore, but
-- have a geometry as well. This gives, in particular, a notion of distance
-- and shortest paths (geodesics) along which you can interpolate.
-- 
-- Keep in mind that the types in this library are
-- generally defined in an abstract-mathematical spirit, which may not always
-- match the intuition if you think about manifolds as embedded in ℝ³.
-- (For instance, the torus inherits its geometry from the decomposition as
-- @'S¹' × 'S¹'@, not from the “doughnut” embedding; the cone over @S¹@ is
-- simply treated as the unit disk, etc..)

{-# LANGUAGE FlexibleInstances          #-}
{-# LANGUAGE UndecidableInstances       #-}
{-# LANGUAGE StandaloneDeriving         #-}
{-# LANGUAGE DeriveGeneric              #-}
{-# LANGUAGE DeriveFunctor              #-}
{-# LANGUAGE DeriveFoldable             #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE TypeFamilies               #-}
{-# LANGUAGE FunctionalDependencies     #-}
{-# LANGUAGE FlexibleContexts           #-}
{-# LANGUAGE GADTs                      #-}
{-# LANGUAGE RankNTypes                 #-}
{-# LANGUAGE TupleSections              #-}
{-# LANGUAGE ParallelListComp           #-}
{-# LANGUAGE UnicodeSyntax              #-}
{-# LANGUAGE ConstraintKinds            #-}
{-# LANGUAGE PatternGuards              #-}
{-# LANGUAGE TypeOperators              #-}
{-# LANGUAGE ScopedTypeVariables        #-}
{-# LANGUAGE LiberalTypeSynonyms        #-}
{-# LANGUAGE CPP                        #-}
{-# LANGUAGE DataKinds                  #-}
{-# LANGUAGE TypeApplications           #-}
{-# LANGUAGE DefaultSignatures          #-}


module Data.Manifold.Riemannian  where


import Data.Maybe
import Data.List.NonEmpty (NonEmpty(..))
import qualified Data.List.NonEmpty as NE
import qualified Data.Vector as Arr

import Data.VectorSpace
import Data.VectorSpace.Free
import Data.AffineSpace
import Math.LinearMap.Category
import Linear (V0(..), V1(..), V2(..), V3(..), V4(..))

import Data.Manifold.Types
import Data.Manifold.Types.Primitive ( (^), empty, embed, coEmbed )
import Data.Manifold.Types.Stiefel
import Data.Manifold.WithBoundary
import Data.Manifold.WithBoundary.Class
import Data.Manifold.PseudoAffine
import Data.Manifold.Atlas (AffineManifold)
    
import qualified Prelude as Hask hiding(foldl, sum, sequence)
import qualified Control.Applicative as Hask
import qualified Control.Monad       as Hask hiding(forM_, sequence)
import Data.Functor.Identity
import qualified Data.Foldable       as Hask
import qualified Data.Traversable as Hask

import Control.Category.Constrained.Prelude hiding
     ((^), all, elem, sum, forM, Foldable(..), Traversable)
import Control.Arrow.Constrained
import Control.Monad.Constrained hiding (forM)
import Data.Foldable.Constrained


class SemimanifoldWithBoundary x => Geodesic x where
  geodesicBetween ::
          x -- ^ Starting point; the interpolation will yield this at -1.
       -> x -- ^ End point, for +1.
            -- 
            --   If the two points are actually connected by a path...
       -> Maybe (D¹ -> x) -- ^ ...then this is the interpolation function. Attention: 
                          --   the type will change to 'Differentiable' in the future.
  middleBetween :: x -> x -> Maybe x
  middleBetween x
p₀ x
p₁ = ((D¹_ Double -> x) -> D¹_ Double -> x
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ Double -> D¹_ Double
forall r. r -> D¹_ r
D¹ Double
0) ((D¹_ Double -> x) -> x) -> Maybe (D¹_ Double -> x) -> Maybe x
forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> x -> x -> Maybe (D¹_ Double -> x)
forall x. Geodesic x => x -> x -> Maybe (D¹_ Double -> x)
geodesicBetween x
p₀ x
p₁


interpolate :: (Geodesic x, IntervalLike i) => x -> x -> Maybe (i -> x)
interpolate :: x -> x -> Maybe (i -> x)
interpolate x
a x
b = ((D¹_ Double -> x) -> (i -> D¹_ Double) -> i -> x
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. i -> D¹_ Double
forall i. IntervalLike i => i -> D¹_ Double
toClosedInterval) ((D¹_ Double -> x) -> i -> x)
-> Maybe (D¹_ Double -> x) -> Maybe (i -> x)
forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> x -> x -> Maybe (D¹_ Double -> x)
forall x. Geodesic x => x -> x -> Maybe (D¹_ Double -> x)
geodesicBetween x
a x
b




#define deriveAffineGD(x)                                         \
instance Geodesic x where {                                        \
  geodesicBetween a b = return $ alerp a b . (/2) . (+1) . xParamD¹; \
  middleBetween a b = return $ alerp a b (1/2) \
 }

deriveAffineGD (ℝ)

instance (Num' s, OpenManifold s) => Geodesic (ZeroDim s) where
  geodesicBetween :: ZeroDim s -> ZeroDim s -> Maybe (D¹_ Double -> ZeroDim s)
geodesicBetween ZeroDim s
Origin ZeroDim s
Origin = (D¹_ Double -> ZeroDim s) -> Maybe (D¹_ Double -> ZeroDim s)
forall (m :: * -> *) a. Monad m (->) => a -> m a
return ((D¹_ Double -> ZeroDim s) -> Maybe (D¹_ Double -> ZeroDim s))
-> (D¹_ Double -> ZeroDim s) -> Maybe (D¹_ Double -> ZeroDim s)
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ \D¹_ Double
_ -> ZeroDim s
forall s. ZeroDim s
Origin
  middleBetween :: ZeroDim s -> ZeroDim s -> Maybe (ZeroDim s)
middleBetween ZeroDim s
Origin ZeroDim s
Origin = ZeroDim s -> Maybe (ZeroDim s)
forall (m :: * -> *) a. Monad m (->) => a -> m a
return ZeroDim s
forall s. ZeroDim s
Origin

instance ∀ a b . ( Geodesic a, Geodesic b
                 , Scalar (Needle (Interior a)) ~ Scalar (Needle (Interior b))
                 , SemimanifoldWithBoundary (a,b)
                 )
      => Geodesic (a,b) where
  geodesicBetween :: (a, b) -> (a, b) -> Maybe (D¹_ Double -> (a, b))
geodesicBetween (a
a,b
b) (a
α,b
β) = ((D¹_ Double -> a) -> (D¹_ Double -> b) -> D¹_ Double -> (a, b))
-> Maybe (D¹_ Double -> a)
-> Maybe (D¹_ Double -> b)
-> Maybe (D¹_ Double -> (a, b))
forall (f :: * -> *) (r :: * -> * -> *) (t :: * -> * -> *) c b a.
(Applicative f r t, Object r c, ObjectMorphism r b c,
 Object t (f c), ObjectMorphism t (f b) (f c), ObjectPair r a b,
 ObjectPair t (f a) (f b)) =>
r a (r b c) -> t (f a) (t (f b) (f c))
liftA2 (D¹_ Double -> a) -> (D¹_ Double -> b) -> D¹_ Double -> (a, b)
forall (a :: * -> * -> *) b c c'.
(PreArrow a, Object a b, ObjectPair a c c') =>
a b c -> a b c' -> a b (c, c')
(&&&) (a -> a -> Maybe (D¹_ Double -> a)
forall x. Geodesic x => x -> x -> Maybe (D¹_ Double -> x)
geodesicBetween a
a a
α) (b -> b -> Maybe (D¹_ Double -> b)
forall x. Geodesic x => x -> x -> Maybe (D¹_ Double -> x)
geodesicBetween b
b b
β)
  middleBetween :: (a, b) -> (a, b) -> Maybe (a, b)
middleBetween (a
a,b
b) (a
α,b
β) = (Maybe a, Maybe b) -> Maybe (a, b)
forall (f :: * -> *) (r :: * -> * -> *) (t :: * -> * -> *) a b.
(Monoidal f r t, ObjectPair r a b, ObjectPair t (f a) (f b),
 Object t (f (a, b))) =>
t (f a, f b) (f (a, b))
fzip (a -> a -> Maybe a
forall x. Geodesic x => x -> x -> Maybe x
middleBetween a
a a
α, b -> b -> Maybe b
forall x. Geodesic x => x -> x -> Maybe x
middleBetween b
b b
β)

-- instance ∀ a b c . (Geodesic a, Geodesic b, Geodesic c) => Geodesic (a,b,c) where
--   geodesicBetween (a,b,c) (α,β,γ)
--       = liftA3 (\ia ib ic t -> (ia t, ib t, ic t))
--            (geodesicBetween a α) (geodesicBetween b β) (geodesicBetween c γ)

-- instance (KnownNat n) => Geodesic (FreeVect n ℝ) where
--   geodesicBetween (FreeVect v) (FreeVect w)
--       = return $ \(D¹ t) -> let μv = (1-t)/2; μw = (t+1)/2
--                             in FreeVect $ Arr.zipWith (\vi wi -> μv*vi + μw*wi) v w

-- instance ∀ v . ( Geodesic v, FiniteFreeSpace v, FiniteFreeSpace (DualVector v)
--                , LinearSpace v, Scalar v ~ ℝ, Geodesic (DualVector v)
--                , InnerSpace (DualVector v) )
--              => Geodesic (Stiefel1 v) where
--   geodesicBetween = gb dualSpaceWitness
--    where gb :: DualSpaceWitness v -> Stiefel1 v -> Stiefel1 v -> Maybe (D¹ -> Stiefel1 v)
--          gb DualSpaceWitness (Stiefel1 p') (Stiefel1 q')
--            = (\f -> \(D¹ t) -> Stiefel1 . f . D¹ $ g * tan (ϑ*t))
--             <$> geodesicBetween p q
--           where p = normalized p'; q = normalized q'
--                 l = magnitude $ p^-^q
--                 ϑ = asin $ l/2
--                 g = sqrt $ 4/l^2 - 1
--   middleBetween = gb dualSpaceWitness
--    where gb :: DualSpaceWitness v -> Stiefel1 v -> Stiefel1 v -> Maybe (Stiefel1 v)
--          gb DualSpaceWitness  (Stiefel1 p) (Stiefel1 q)
--              = Stiefel1 <$> middleBetween (normalized p) (normalized q)


instance Geodesic S⁰ where
  geodesicBetween :: S⁰ -> S⁰ -> Maybe (D¹_ Double -> S⁰)
geodesicBetween S⁰
PositiveHalfSphere S⁰
PositiveHalfSphere = (D¹_ Double -> S⁰) -> Maybe (D¹_ Double -> S⁰)
forall (m :: * -> *) a. Monad m (->) => a -> m a
return ((D¹_ Double -> S⁰) -> Maybe (D¹_ Double -> S⁰))
-> (D¹_ Double -> S⁰) -> Maybe (D¹_ Double -> S⁰)
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ S⁰ -> D¹_ Double -> S⁰
forall (a :: * -> * -> *) b x.
(WellPointed a, Object a b, ObjectPoint a x) =>
x -> a b x
const S⁰
forall r. S⁰_ r
PositiveHalfSphere
  geodesicBetween S⁰
NegativeHalfSphere S⁰
NegativeHalfSphere = (D¹_ Double -> S⁰) -> Maybe (D¹_ Double -> S⁰)
forall (m :: * -> *) a. Monad m (->) => a -> m a
return ((D¹_ Double -> S⁰) -> Maybe (D¹_ Double -> S⁰))
-> (D¹_ Double -> S⁰) -> Maybe (D¹_ Double -> S⁰)
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ S⁰ -> D¹_ Double -> S⁰
forall (a :: * -> * -> *) b x.
(WellPointed a, Object a b, ObjectPoint a x) =>
x -> a b x
const S⁰
forall r. S⁰_ r
NegativeHalfSphere
  geodesicBetween S⁰
_ S⁰
_ = Maybe (D¹_ Double -> S⁰)
forall (f :: * -> *) a. Alternative f => f a
empty

instance Geodesic S¹ where
  geodesicBetween :: S¹ -> S¹ -> Maybe (D¹_ Double -> S¹)
geodesicBetween (S¹Polar Double
φ) (S¹Polar Double
ϕ)
    | Double -> Double
forall a. Num a => a -> a
abs (Double
φDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
ϕ) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
forall a. Floating a => a
pi  = ((D¹_ Double -> Double) -> (Double -> S¹) -> D¹_ Double -> S¹
forall (k :: * -> * -> *) a b c.
(Category k, Object k a, Object k b, Object k c) =>
k a b -> k b c -> k a c
>>> Double -> S¹
forall r. r -> S¹_ r
S¹Polar) ((D¹_ Double -> Double) -> D¹_ Double -> S¹)
-> Maybe (D¹_ Double -> Double) -> Maybe (D¹_ Double -> S¹)
forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> Double -> Double -> Maybe (D¹_ Double -> Double)
forall x. Geodesic x => x -> x -> Maybe (D¹_ Double -> x)
geodesicBetween Double
φ Double
ϕ
    | Double
φ Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0           = ((D¹_ Double -> Double) -> (Double -> S¹) -> D¹_ Double -> S¹
forall (k :: * -> * -> *) a b c.
(Category k, Object k a, Object k b, Object k c) =>
k a b -> k b c -> k a c
>>> Double -> S¹
forall r. r -> S¹_ r
S¹Polar (Double -> S¹) -> (Double -> Double) -> Double -> S¹
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. \Double
ψ -> Double -> Double
forall a. Num a => a -> a
signum Double
ψDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
ψ)
                        ((D¹_ Double -> Double) -> D¹_ Double -> S¹)
-> Maybe (D¹_ Double -> Double) -> Maybe (D¹_ Double -> S¹)
forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> Double -> Double -> Maybe (D¹_ Double -> Double)
forall x. Geodesic x => x -> x -> Maybe (D¹_ Double -> x)
geodesicBetween (Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
φ) (-Double
ϕDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
forall a. Floating a => a
pi)
    | Bool
otherwise       = ((D¹_ Double -> Double) -> (Double -> S¹) -> D¹_ Double -> S¹
forall (k :: * -> * -> *) a b c.
(Category k, Object k a, Object k b, Object k c) =>
k a b -> k b c -> k a c
>>> Double -> S¹
forall r. r -> S¹_ r
S¹Polar (Double -> S¹) -> (Double -> Double) -> Double -> S¹
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. \Double
ψ -> Double -> Double
forall a. Num a => a -> a
signum Double
ψDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
ψ)
                        ((D¹_ Double -> Double) -> D¹_ Double -> S¹)
-> Maybe (D¹_ Double -> Double) -> Maybe (D¹_ Double -> S¹)
forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> Double -> Double -> Maybe (D¹_ Double -> Double)
forall x. Geodesic x => x -> x -> Maybe (D¹_ Double -> x)
geodesicBetween (-Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
φ) (Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
ϕ)
  middleBetween :: S¹ -> S¹ -> Maybe S¹
middleBetween (S¹Polar Double
φ) (S¹Polar Double
ϕ)
    | Double -> Double
forall a. Num a => a -> a
abs (Double
φDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
ϕ) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
forall a. Floating a => a
pi  = Double -> S¹
forall r. r -> S¹_ r
S¹Polar (Double -> S¹) -> Maybe Double -> Maybe S¹
forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> Double -> Double -> Maybe Double
forall x. Geodesic x => x -> x -> Maybe x
middleBetween Double
φ Double
ϕ
    | Double
φ Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0           = Double -> S¹
forall r. r -> S¹_ r
S¹Polar (Double -> S¹) -> Maybe Double -> Maybe S¹
forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> Double -> Double -> Maybe Double
forall x. Geodesic x => x -> x -> Maybe x
middleBetween (Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
φ) (-Double
ϕDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
forall a. Floating a => a
pi)
    | Bool
otherwise       = Double -> S¹
forall r. r -> S¹_ r
S¹Polar (Double -> S¹) -> Maybe Double -> Maybe S¹
forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> Double -> Double -> Maybe Double
forall x. Geodesic x => x -> x -> Maybe x
middleBetween (-Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
φ) (Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
ϕ)


-- instance Geodesic (Cℝay S⁰) where
--   geodesicBetween p q = (>>> fromℝ) <$> geodesicBetween (toℝ p) (toℝ q)
--    where toℝ (Cℝay h PositiveHalfSphere) = h
--          toℝ (Cℝay h NegativeHalfSphere) = -h
--          fromℝ x | x>0        = Cℝay x PositiveHalfSphere
--                  | otherwise  = Cℝay (-x) NegativeHalfSphere
-- 
-- instance Geodesic (CD¹ S⁰) where
--   geodesicBetween p q = (>>> fromI) <$> geodesicBetween (toI p) (toI q)
--    where toI (CD¹ h PositiveHalfSphere) = h
--          toI (CD¹ h NegativeHalfSphere) = -h
--          fromI x | x>0        = CD¹ x PositiveHalfSphere
--                  | otherwise  = CD¹ (-x) NegativeHalfSphere
-- 
-- instance Geodesic (Cℝay S¹) where
--   geodesicBetween p q = (>>> fromP) <$> geodesicBetween (toP p) (toP q)
--    where fromP = fromInterior
--          toP w = case toInterior w of {Just i -> i}
-- 
-- instance Geodesic (CD¹ S¹) where
--   geodesicBetween p q = (>>> fromI) <$> geodesicBetween (toI p) (toI q)
--    where toI (CD¹ h (S¹ φ)) = (h*cos φ, h*sin φ)
--          fromI (x,y) = CD¹ (sqrt $ x^2+y^2) (S¹ $ atan2 y x)
-- 
-- instance Geodesic (Cℝay S²) where
--   geodesicBetween p q = (>>> fromP) <$> geodesicBetween (toP p) (toP q)
--    where fromP = fromInterior
--          toP w = case toInterior w of {Just i -> i}
-- 
-- instance Geodesic (CD¹ S²) where
--   geodesicBetween p q = (>>> fromI) <$> geodesicBetween (toI p) (toI q :: ℝ³)
--    where toI (CD¹ h sph) = h *^ embed sph
--          fromI v = CD¹ (magnitude v) (coEmbed v)

#define geoVSpCone(c,t)                                               \
instance (c) => Geodesic (Cℝay (t)) where {                            \
  geodesicBetween p q = (>>> fromP) <$> geodesicBetween (toP p) (toP q) \
   where { fromP (x,0) = Cℝay 0 x                                        \
         ; fromP (x,h) = Cℝay h (x^/h)                                    \
         ; toP (Cℝay h w) = ( h*^w, h ) } } ;                              \
instance (c) => Geodesic (CD¹ (t)) where {                                  \
  geodesicBetween p q = (>>> fromP) <$> geodesicBetween (toP p) (toP q)      \
   where { fromP (x,0) = CD¹ 0 x                                              \
         ; fromP (x,h) = CD¹ h (x^/h)                                          \
         ; toP (CD¹ h w) = ( h*^w, h ) } }

-- geoVSpCone ((), ℝ)
-- geoVSpCone ((), ℝ⁰)
-- geoVSpCone ((WithField ℝ HilbertManifold a, WithField ℝ HilbertManifold b
--             , Geodesic (a,b)), (a,b))
-- geoVSpCone (KnownNat n, FreeVect n ℝ)

deriveAffineGD ((V0 ℝ))
deriveAffineGD (ℝ¹)
deriveAffineGD (ℝ²)
deriveAffineGD (ℝ³)
deriveAffineGD (ℝ⁴)

instance (LinearSpace v, Scalar v ~ ℝ, LinearSpace w, Scalar w ~ ℝ)
             => Geodesic (Tensor ℝ v w) where
  geodesicBetween :: Tensor Double v w
-> Tensor Double v w -> Maybe (D¹_ Double -> Tensor Double v w)
geodesicBetween Tensor Double v w
a Tensor Double v w
b = (D¹_ Double -> Tensor Double v w)
-> Maybe (D¹_ Double -> Tensor Double v w)
forall (m :: * -> *) a. Monad m (->) => a -> m a
return ((D¹_ Double -> Tensor Double v w)
 -> Maybe (D¹_ Double -> Tensor Double v w))
-> (D¹_ Double -> Tensor Double v w)
-> Maybe (D¹_ Double -> Tensor Double v w)
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ Tensor Double v w
-> Tensor Double v w
-> Scalar (Diff (Tensor Double v w))
-> Tensor Double v w
forall p.
(AffineSpace p, VectorSpace (Diff p)) =>
p -> p -> Scalar (Diff p) -> p
alerp Tensor Double v w
a Tensor Double v w
b (Double -> Tensor Double v w)
-> (D¹_ Double -> Double) -> D¹_ Double -> Tensor Double v w
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2) (Double -> Double)
-> (D¹_ Double -> Double) -> D¹_ Double -> Double
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1) (Double -> Double)
-> (D¹_ Double -> Double) -> D¹_ Double -> Double
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. D¹_ Double -> Double
forall r. D¹_ r -> r
xParamD¹
instance (LinearSpace v, Scalar v ~ ℝ, LinearSpace w, Scalar w ~ ℝ)
             => Geodesic (LinearMap ℝ v w) where
  geodesicBetween :: LinearMap Double v w
-> LinearMap Double v w
-> Maybe (D¹_ Double -> LinearMap Double v w)
geodesicBetween LinearMap Double v w
a LinearMap Double v w
b = (D¹_ Double -> LinearMap Double v w)
-> Maybe (D¹_ Double -> LinearMap Double v w)
forall (m :: * -> *) a. Monad m (->) => a -> m a
return ((D¹_ Double -> LinearMap Double v w)
 -> Maybe (D¹_ Double -> LinearMap Double v w))
-> (D¹_ Double -> LinearMap Double v w)
-> Maybe (D¹_ Double -> LinearMap Double v w)
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ LinearMap Double v w
-> LinearMap Double v w
-> Scalar (Diff (LinearMap Double v w))
-> LinearMap Double v w
forall p.
(AffineSpace p, VectorSpace (Diff p)) =>
p -> p -> Scalar (Diff p) -> p
alerp LinearMap Double v w
a LinearMap Double v w
b (Double -> LinearMap Double v w)
-> (D¹_ Double -> Double) -> D¹_ Double -> LinearMap Double v w
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2) (Double -> Double)
-> (D¹_ Double -> Double) -> D¹_ Double -> Double
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1) (Double -> Double)
-> (D¹_ Double -> Double) -> D¹_ Double -> Double
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. D¹_ Double -> Double
forall r. D¹_ r -> r
xParamD¹
instance (LinearSpace v, Scalar v ~ ℝ, LinearSpace w, Scalar w ~ ℝ)
             => Geodesic (LinearFunction ℝ v w) where
  geodesicBetween :: LinearFunction Double v w
-> LinearFunction Double v w
-> Maybe (D¹_ Double -> LinearFunction Double v w)
geodesicBetween LinearFunction Double v w
a LinearFunction Double v w
b = (D¹_ Double -> LinearFunction Double v w)
-> Maybe (D¹_ Double -> LinearFunction Double v w)
forall (m :: * -> *) a. Monad m (->) => a -> m a
return ((D¹_ Double -> LinearFunction Double v w)
 -> Maybe (D¹_ Double -> LinearFunction Double v w))
-> (D¹_ Double -> LinearFunction Double v w)
-> Maybe (D¹_ Double -> LinearFunction Double v w)
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ LinearFunction Double v w
-> LinearFunction Double v w
-> Scalar (Diff (LinearFunction Double v w))
-> LinearFunction Double v w
forall p.
(AffineSpace p, VectorSpace (Diff p)) =>
p -> p -> Scalar (Diff p) -> p
alerp LinearFunction Double v w
a LinearFunction Double v w
b (Double -> LinearFunction Double v w)
-> (D¹_ Double -> Double)
-> D¹_ Double
-> LinearFunction Double v w
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2) (Double -> Double)
-> (D¹_ Double -> Double) -> D¹_ Double -> Double
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1) (Double -> Double)
-> (D¹_ Double -> Double) -> D¹_ Double -> Double
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. D¹_ Double -> Double
forall r. D¹_ r -> r
xParamD¹


-- | One-dimensional manifolds, whose closure is homeomorpic to the unit interval.
class WithField ℝ PseudoAffine (Interior i) => IntervalLike i where
  toClosedInterval :: i -> D¹ -- Differentiable ℝ i D¹

instance IntervalLike D¹ where
  toClosedInterval :: D¹_ Double -> D¹_ Double
toClosedInterval = D¹_ Double -> D¹_ Double
forall κ (k :: κ -> κ -> *) (a :: κ).
(Category k, Object k a) =>
k a a
id
-- instance IntervalLike (CD¹ S⁰) where
--   toClosedInterval (CD¹ h PositiveHalfSphere) = D¹ h
--   toClosedInterval (CD¹ h NegativeHalfSphere) = D¹ (-h)
-- instance IntervalLike (Cℝay S⁰) where
--   toClosedInterval (Cℝay h PositiveHalfSphere) = D¹ $ tanh h
--   toClosedInterval (Cℝay h NegativeHalfSphere) = D¹ $ -tanh h
-- instance IntervalLike (CD¹ ℝ⁰) where
--   toClosedInterval (CD¹ h Origin) = D¹ $ h*2 - 1
-- instance IntervalLike (Cℝay ℝ⁰) where
--   toClosedInterval (Cℝay h Origin) = D¹ $ 1 - 2/(h+1)
instance IntervalLike ℝ where
  toClosedInterval :: Double -> D¹_ Double
toClosedInterval Double
x = Double -> D¹_ Double
forall r. r -> D¹_ r
D¹ (Double -> D¹_ Double) -> Double -> D¹_ Double
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ Double -> Double
forall a. Floating a => a -> a
tanh Double
x





class Geodesic m => Riemannian m where
  rieMetric :: RieMetric m

instance Riemannian ℝ where
  rieMetric :: RieMetric Double
rieMetric = Norm Double -> Double -> Norm Double
forall (a :: * -> * -> *) b x.
(WellPointed a, Object a b, ObjectPoint a x) =>
x -> a b x
const Norm Double
forall v. HilbertSpace v => Norm v
euclideanNorm





pointsBarycenter :: Geodesic m => NonEmpty m -> Maybe m
pointsBarycenter :: NonEmpty m -> Maybe m
pointsBarycenter (m
p:|[]) = m -> Maybe m
forall a. a -> Maybe a
Just m
p
pointsBarycenter NonEmpty m
ps = case ( NonEmpty m -> Maybe m
forall m. Geodesic m => NonEmpty m -> Maybe m
pointsBarycenter ([m] -> NonEmpty m
forall a. [a] -> NonEmpty a
NE.fromList [m]
group₀)
                           , NonEmpty m -> Maybe m
forall m. Geodesic m => NonEmpty m -> Maybe m
pointsBarycenter ([m] -> NonEmpty m
forall a. [a] -> NonEmpty a
NE.fromList [m]
group₁) ) of
            (Just m
bc₀, Just m
bc₁)
                | Int
δn Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0      -> m -> m -> Maybe m
forall x. Geodesic x => x -> x -> Maybe x
middleBetween m
bc₀ m
bc₁
                | Bool
otherwise    -> ((D¹_ Double -> m) -> D¹_ Double -> m
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ Double -> D¹_ Double
forall r. r -> D¹_ r
D¹ (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
δnDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
ntot))
                                    ((D¹_ Double -> m) -> m) -> Maybe (D¹_ Double -> m) -> Maybe m
forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> m -> m -> Maybe (D¹_ Double -> m)
forall x. Geodesic x => x -> x -> Maybe (D¹_ Double -> x)
geodesicBetween m
bc₀ m
bc₁
            (Maybe m, Maybe m)
_                  -> Maybe m
forall a. Maybe a
Nothing
 where psl :: [m]
psl = NonEmpty m -> [m]
forall (t :: * -> *) a. Foldable t => t a -> [a]
Hask.toList NonEmpty m
ps
       ([m]
group₀, [m]
group₁) = Int -> [m] -> ([m], [m])
forall a. Int -> [a] -> ([a], [a])
splitAt Int
nl [m]
psl
       ntot :: Int
ntot = [m] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [m]
psl
       nl :: Int
nl = Int
ntotInt -> Int -> Int
forall a. Integral a => a -> a -> a
`quot`Int
2
       δn :: Int
δn = Int
ntot  Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
nlInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
2


type FlatSpace x = (AffineManifold x, Geodesic x, SimpleSpace x)