-- |
-- Module      : Data.Manifold.Atlas
-- Copyright   : (c) Justus Sagemüller 2015
-- License     : GPL v3
-- 
-- Maintainer  : (@) jsag $ hvl.no
-- Stability   : experimental
-- Portability : portable
-- 

{-# LANGUAGE TypeFamilies              #-}
{-# LANGUAGE ConstraintKinds           #-}
{-# LANGUAGE FlexibleContexts          #-}
{-# LANGUAGE FlexibleInstances         #-}
{-# LANGUAGE UndecidableInstances      #-}
{-# LANGUAGE EmptyDataDecls, EmptyCase #-}
{-# LANGUAGE CPP                       #-}
{-# LANGUAGE ScopedTypeVariables       #-}

module Data.Manifold.Atlas where

import Prelude as Hask

import Data.VectorSpace
import Data.Manifold.PseudoAffine
import Data.Manifold.Types.Primitive
import Data.Manifold.WithBoundary
import Data.Manifold.WithBoundary.Class

import Data.Void

import Data.VectorSpace.Free
import Math.LinearMap.Category

import Control.Arrow

import Data.MemoTrie (HasTrie)

import qualified Linear.Affine as LinAff

class SemimanifoldWithBoundary m => Atlas m where
  type ChartIndex m :: *
  chartReferencePoint :: ChartIndex m -> m
  lookupAtlas :: m -> ChartIndex m

#define VectorSpaceAtlas(c,v)              \
instance (c) => Atlas (v) where {           \
  type ChartIndex (v) = ();                  \
  chartReferencePoint () = zeroV;              \
  lookupAtlas _ = () }

type NumPrime s = (Num' s, Eq s, OpenManifold s, ProjectableBoundary s)

VectorSpaceAtlas(NumPrime s, ZeroDim s)
VectorSpaceAtlas((), )
VectorSpaceAtlas(NumPrime s, V0 s)
VectorSpaceAtlas(NumPrime s, V1 s)
VectorSpaceAtlas(NumPrime s, V2 s)
VectorSpaceAtlas(NumPrime s, V3 s)
VectorSpaceAtlas(NumPrime s, V4 s)
VectorSpaceAtlas((NumPrime s, LinearSpace v, Scalar v ~ s, LinearSpace w, Scalar w ~ s), LinearMap s v w)
VectorSpaceAtlas((NumPrime s, LinearSpace v, Scalar v ~ s, LinearSpace w, Scalar w ~ s), Tensor s v w)

instance (Atlas x, Atlas y, SemimanifoldWithBoundary (x,y)) => Atlas (x,y) where
  type ChartIndex (x,y) = (ChartIndex x, ChartIndex y)
  chartReferencePoint :: ChartIndex (x, y) -> (x, y)
chartReferencePoint = ChartIndex x -> x
forall m. Atlas m => ChartIndex m -> m
chartReferencePoint (ChartIndex x -> x)
-> (ChartIndex y -> y) -> (ChartIndex x, ChartIndex y) -> (x, y)
forall (a :: * -> * -> *) b c b' c'.
Arrow a =>
a b c -> a b' c' -> a (b, b') (c, c')
*** ChartIndex y -> y
forall m. Atlas m => ChartIndex m -> m
chartReferencePoint
  lookupAtlas :: (x, y) -> ChartIndex (x, y)
lookupAtlas = x -> ChartIndex x
forall m. Atlas m => m -> ChartIndex m
lookupAtlas (x -> ChartIndex x)
-> (y -> ChartIndex y) -> (x, y) -> (ChartIndex x, ChartIndex y)
forall (a :: * -> * -> *) b c b' c'.
Arrow a =>
a b c -> a b' c' -> a (b, b') (c, c')
*** y -> ChartIndex y
forall m. Atlas m => m -> ChartIndex m
lookupAtlas

instance Atlas S⁰ where
  type ChartIndex S⁰ = S⁰
  chartReferencePoint :: ChartIndex S⁰ -> S⁰
chartReferencePoint = ChartIndex S⁰ -> S⁰
forall a. a -> a
id
  lookupAtlas :: S⁰ -> ChartIndex S⁰
lookupAtlas = S⁰ -> ChartIndex S⁰
forall a. a -> a
id
instance Atlas  where
  type ChartIndex  = S⁰
  chartReferencePoint :: ChartIndex S¹ -> S¹
chartReferencePoint ChartIndex S¹
NegativeHalfSphere = ℝ -> S¹
forall r. r -> S¹_ r
S¹Polar (ℝ -> S¹) -> ℝ -> S¹
forall a b. (a -> b) -> a -> b
$ -ℝ
forall a. Floating a => a
piℝ -> ℝ -> ℝ
forall a. Fractional a => a -> a -> a
/2
  chartReferencePoint ChartIndex S¹
PositiveHalfSphere = ℝ -> S¹
forall r. r -> S¹_ r
S¹Polar (ℝ -> S¹) -> ℝ -> S¹
forall a b. (a -> b) -> a -> b
$ ℝ
forall a. Floating a => a
piℝ -> ℝ -> ℝ
forall a. Fractional a => a -> a -> a
/2
  lookupAtlas :: S¹ -> ChartIndex S¹
lookupAtlas (S¹Polar φ) | φℝ -> ℝ -> Bool
forall a. Ord a => a -> a -> Bool
<0        = ChartIndex S¹
forall r. S⁰_ r
NegativeHalfSphere
                     | Bool
otherwise  = ChartIndex S¹
forall r. S⁰_ r
PositiveHalfSphere
instance Atlas  where
  type ChartIndex  = S⁰
  chartReferencePoint :: ChartIndex S² -> S²
chartReferencePoint ChartIndex S²
PositiveHalfSphere = ℝ -> ℝ -> S²
forall r. r -> r -> S²_ r
S²Polar 0 0
  chartReferencePoint ChartIndex S²
NegativeHalfSphere = ℝ -> ℝ -> S²
forall r. r -> r -> S²_ r
S²Polar ℝ
forall a. Floating a => a
pi 0
  lookupAtlas :: S² -> ChartIndex S²
lookupAtlas (S²Polar ϑ _) | ϑℝ -> ℝ -> Bool
forall a. Ord a => a -> a -> Bool
<ℝ
forall a. Floating a => a
piℝ -> ℝ -> ℝ
forall a. Fractional a => a -> a -> a
/2     = ChartIndex S²
forall r. S⁰_ r
PositiveHalfSphere
                            | Bool
otherwise  = ChartIndex S²
forall r. S⁰_ r
NegativeHalfSphere

instance (Num'' n, LinearManifold (a n), Scalar (a n) ~ n, Needle (a n) ~ a n)
              => Atlas (LinAff.Point a n) where
  type ChartIndex (LinAff.Point a n) = ()
  chartReferencePoint :: ChartIndex (Point a n) -> Point a n
chartReferencePoint () = a n -> Point a n
forall (f :: * -> *) a. f a -> Point f a
LinAff.P a n
forall v. AdditiveGroup v => v
zeroV
  lookupAtlas :: Point a n -> ChartIndex (Point a n)
lookupAtlas Point a n
_ = ()

type Atlas' m = (Atlas m, HasTrie (ChartIndex m))


-- | The 'AffineSpace' class plus manifold constraints.
type AffineManifold m = ( Atlas' m, Manifold m, AffineSpace m
                        , Needle m ~ Diff m )

-- | An euclidean space is a real affine space whose tangent space is a Hilbert space.
type EuclidSpace x = ( AffineManifold x, InnerSpace (Diff x)
                     , DualVector (Diff x) ~ Diff x, Floating (Scalar (Diff x)) )

euclideanMetric :: EuclidSpace x => proxy x -> Metric x
euclideanMetric :: proxy x -> Metric x
euclideanMetric proxy x
_ = Metric x
forall v. HilbertSpace v => Norm v
euclideanNorm