-- |
-- 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₁ = (forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall r. r -> D¹_ r
D¹ ℝ
0) forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> forall x. Geodesic x => x -> x -> Maybe (D¹ -> x)
geodesicBetween x
p₀ x
p₁


interpolate :: (Geodesic x, IntervalLike i) => x -> x -> Maybe (i -> x)
interpolate :: forall x i.
(Geodesic x, IntervalLike i) =>
x -> x -> Maybe (i -> x)
interpolate x
a x
b = (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
. forall i. IntervalLike i => i -> D¹
toClosedInterval) forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> forall x. Geodesic x => x -> x -> Maybe (D¹ -> 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¹ -> ZeroDim s)
geodesicBetween ZeroDim s
Origin ZeroDim s
Origin = forall (m :: * -> *) a. Monad m (->) => a -> m a
return forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ \D¹
_ -> forall s. ZeroDim s
Origin
  middleBetween :: ZeroDim s -> ZeroDim s -> Maybe (ZeroDim s)
middleBetween ZeroDim s
Origin ZeroDim s
Origin = forall (m :: * -> *) a. Monad m (->) => a -> m a
return 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¹ -> (a, b))
geodesicBetween (a
a,b
b) (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 forall (a :: * -> * -> *) b c c'.
(PreArrow a, Object a b, ObjectPair a c c') =>
a b c -> a b c' -> a b (c, c')
(&&&) (forall x. Geodesic x => x -> x -> Maybe (D¹ -> x)
geodesicBetween a
a a
α) (forall x. Geodesic x => x -> x -> Maybe (D¹ -> x)
geodesicBetween b
b b
β)
  middleBetween :: (a, b) -> (a, b) -> Maybe (a, b)
middleBetween (a
a,b
b) (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 (forall x. Geodesic x => x -> x -> Maybe x
middleBetween a
a a
α, 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¹ -> S⁰)
geodesicBetween S⁰
PositiveHalfSphere S⁰
PositiveHalfSphere = forall (m :: * -> *) a. Monad m (->) => a -> m a
return forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall (a :: * -> * -> *) b x.
(WellPointed a, Object a b, ObjectPoint a x) =>
x -> a b x
const forall r. S⁰_ r
PositiveHalfSphere
  geodesicBetween S⁰
NegativeHalfSphere S⁰
NegativeHalfSphere = forall (m :: * -> *) a. Monad m (->) => a -> m a
return forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall (a :: * -> * -> *) b x.
(WellPointed a, Object a b, ObjectPoint a x) =>
x -> a b x
const forall r. S⁰_ r
NegativeHalfSphere
  geodesicBetween S⁰
_ S⁰
_ = forall (f :: * -> *) a. Alternative f => f a
empty

instance Geodesic S¹ where
  geodesicBetween :: S¹ -> S¹ -> Maybe (D¹ -> S¹)
geodesicBetween (S¹Polar ℝ
φ) (S¹Polar ℝ
ϕ)
    | forall a. Num a => a -> a
abs (ℝ
φforall a. Num a => a -> a -> a
-ℝ
ϕ) forall a. Ord a => a -> a -> Bool
< forall a. Floating a => a
pi  = (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
>>> forall r. r -> S¹_ r
S¹Polar) forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> forall x. Geodesic x => x -> x -> Maybe (D¹ -> x)
geodesicBetween ℝ
φ ℝ
ϕ
    | ℝ
φ forall a. Ord a => a -> a -> Bool
> ℝ
0           = (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
>>> forall r. r -> S¹_ r
S¹Polar 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
. \ℝ
ψ -> forall a. Num a => a -> a
signum ℝ
ψforall a. Num a => a -> a -> a
*forall a. Floating a => a
pi forall a. Num a => a -> a -> a
- ℝ
ψ)
                        forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> forall x. Geodesic x => x -> x -> Maybe (D¹ -> x)
geodesicBetween (forall a. Floating a => a
piforall a. Num a => a -> a -> a
-ℝ
φ) (-ℝ
ϕforall a. Num a => a -> a -> a
-forall a. Floating a => a
pi)
    | Bool
otherwise       = (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
>>> forall r. r -> S¹_ r
S¹Polar 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
. \ℝ
ψ -> forall a. Num a => a -> a
signum ℝ
ψforall a. Num a => a -> a -> a
*forall a. Floating a => a
pi forall a. Num a => a -> a -> a
- ℝ
ψ)
                        forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> forall x. Geodesic x => x -> x -> Maybe (D¹ -> x)
geodesicBetween (-forall a. Floating a => a
piforall a. Num a => a -> a -> a
-ℝ
φ) (forall a. Floating a => a
piforall a. Num a => a -> a -> a
-ℝ
ϕ)
  middleBetween :: S¹ -> S¹ -> Maybe S¹
middleBetween (S¹Polar ℝ
φ) (S¹Polar ℝ
ϕ)
    | forall a. Num a => a -> a
abs (ℝ
φforall a. Num a => a -> a -> a
-ℝ
ϕ) forall a. Ord a => a -> a -> Bool
< forall a. Floating a => a
pi  = forall r. r -> S¹_ r
S¹Polar forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> forall x. Geodesic x => x -> x -> Maybe x
middleBetween ℝ
φ ℝ
ϕ
    | ℝ
φ forall a. Ord a => a -> a -> Bool
> ℝ
0           = forall r. r -> S¹_ r
S¹Polar forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> forall x. Geodesic x => x -> x -> Maybe x
middleBetween (forall a. Floating a => a
piforall a. Num a => a -> a -> a
-ℝ
φ) (-ℝ
ϕforall a. Num a => a -> a -> a
-forall a. Floating a => a
pi)
    | Bool
otherwise       = forall r. r -> S¹_ r
S¹Polar forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> forall x. Geodesic x => x -> x -> Maybe x
middleBetween (-forall a. Floating a => a
piforall a. Num a => a -> a -> a
-ℝ
φ) (forall a. Floating a => a
piforall a. Num a => a -> a -> a
-ℝ
ϕ)


-- 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 ℝ v w -> Tensor ℝ v w -> Maybe (D¹ -> Tensor ℝ v w)
geodesicBetween Tensor ℝ v w
a Tensor ℝ v w
b = forall (m :: * -> *) a. Monad m (->) => a -> m a
return forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall p.
(AffineSpace p, VectorSpace (Diff p)) =>
p -> p -> Scalar (Diff p) -> p
alerp Tensor ℝ v w
a Tensor ℝ v w
b 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
. (forall a. Fractional a => a -> a -> a
/ℝ
2) 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
. (forall a. Num a => a -> a -> a
+ℝ
1) 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
. forall r. D¹_ r -> r
xParamD¹
instance (LinearSpace v, Scalar v ~ ℝ, LinearSpace w, Scalar w ~ ℝ)
             => Geodesic (LinearMap ℝ v w) where
  geodesicBetween :: LinearMap ℝ v w -> LinearMap ℝ v w -> Maybe (D¹ -> LinearMap ℝ v w)
geodesicBetween LinearMap ℝ v w
a LinearMap ℝ v w
b = forall (m :: * -> *) a. Monad m (->) => a -> m a
return forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall p.
(AffineSpace p, VectorSpace (Diff p)) =>
p -> p -> Scalar (Diff p) -> p
alerp LinearMap ℝ v w
a LinearMap ℝ v w
b 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
. (forall a. Fractional a => a -> a -> a
/ℝ
2) 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
. (forall a. Num a => a -> a -> a
+ℝ
1) 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
. forall r. D¹_ r -> r
xParamD¹
instance (LinearSpace v, Scalar v ~ ℝ, LinearSpace w, Scalar w ~ ℝ)
             => Geodesic (LinearFunction ℝ v w) where
  geodesicBetween :: LinearFunction ℝ v w
-> LinearFunction ℝ v w -> Maybe (D¹ -> LinearFunction ℝ v w)
geodesicBetween LinearFunction ℝ v w
a LinearFunction ℝ v w
b = forall (m :: * -> *) a. Monad m (->) => a -> m a
return forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall p.
(AffineSpace p, VectorSpace (Diff p)) =>
p -> p -> Scalar (Diff p) -> p
alerp LinearFunction ℝ v w
a LinearFunction ℝ v w
b 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
. (forall a. Fractional a => a -> a -> a
/ℝ
2) 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
. (forall a. Num a => a -> a -> a
+ℝ
1) 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
. 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¹ -> D¹
toClosedInterval = 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 :: ℝ -> D¹
toClosedInterval ℝ
x = forall r. r -> D¹_ r
D¹ forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall a. Floating a => a -> a
tanh ℝ
x





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

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





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


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