module Data.Manifold (module Data.Manifold, module Data.Manifold.Types.Primitive) where
import Data.List
import Data.Maybe
import Data.Semigroup
import Data.Function (on)
import Data.VectorSpace
import Data.AffineSpace
import Data.Basis
import Data.Complex hiding (magnitude)
import Data.Void
import Data.Manifold.Types.Primitive
import qualified Prelude
import Control.Category.Constrained.Prelude hiding ((^), Foldable(..))
import Control.Arrow.Constrained
import Control.Monad.Constrained
import Data.Foldable.Constrained
data domain :--> codomain where
Continuous :: ( Manifold d, Manifold c
, v ~ TangentSpace d, u ~ TangentSpace c
, δ ~ Metric v, ε ~ Metric u ) =>
{ runContinuous :: Chart d -> v -> (Chart c, u, ε->Option δ) }
-> d :--> c
continuous_id' :: Manifold m => m :--> m
continuous_id' = Continuous id'
where id' chart v = (chart, v, return)
const__ :: (Manifold c, Manifold d)
=> c -> d:-->c
const__ x = Continuous f
where f _ _ = (tgtChart, w, const mzero)
tgtChart = head $ localAtlas x
w = case tgtChart of
IdChart -> x
Chart _ tchOut _ -> fromJust (tchOut x) $ x
flatContinuous :: ( FlatManifold v, FlatManifold w, δ~Metric v, ε~Metric w )
=> (v -> (w, ε -> Option δ)) -> (v:-->w)
flatContinuous f = Continuous cnt
where cnt IdChart v = let (w, postEps) = f v
in (IdChart, w, postEps)
cnt (Chart inMap _ _) v = let (v', preEps) = runFlatContinuous inMap v
(w, postEps) = f v'
in (IdChart, w, preEps>=>postEps)
runFlatContinuous :: ( FlatManifold v, FlatManifold w, δ~Metric v, ε~Metric w )
=> (v:-->w) -> v -> (w, ε -> Option δ)
runFlatContinuous (Continuous cnf) v = (w, preEps>=>postEps)
where (cc', v', preEps) = cnf IdChart v
(w, postEps) = case cc' of
IdChart -> (v', return)
Chart inMap _ _ -> runFlatContinuous inMap v'
instance Category (:-->) where
type Object (:-->) t = Manifold t
id = Continuous $ \c v -> (c, v, just)
Continuous f . Continuous g = Continuous h
where h srcChart u = (tgtChart, w, q>=>p)
where (interChart, v, p) = g srcChart u
(tgtChart, w, q) = f interChart v
instance EnhancedCat (->) (:-->) where
Continuous f `arr` x = y
where (tch, v, _) = f sch u
y = case tch of Chart tchIn _ _ -> tchIn $ v
IdChart -> v
u = case sch of Chart _ schOut _ -> fromJust (schOut x) $ x
IdChart -> x
sch = head $ localAtlas x
instance Cartesian (:-->) where
type PairObjects (:-->) a b = ( FlatManifold a, FlatManifold b, Manifold(a,b) )
swap = Continuous $ \c t -> case c of
IdChart -> let (v,w) = t in (IdChart, (w,v), return)
Chart inMap _ _ -> let ((v,w), epsP) = runFlatContinuous inMap t
in (IdChart, (w,v), epsP)
attachUnit = Continuous $ \c v -> case c of
IdChart -> (IdChart, (v,()), return)
Chart inMap _ _ -> let (v', epsP) = runFlatContinuous inMap v
in (IdChart, (v',()), epsP)
detachUnit = Continuous $ \c t -> case c of
IdChart -> let (v,()) = t in (IdChart, v, return)
Chart inMap _ _ -> let ((v,()), epsP) = runFlatContinuous inMap t
in (IdChart, v, epsP)
regroup = Continuous $ \c t -> case c of
IdChart -> let (u,(v,w)) = t in (IdChart, ((u,v),w), return)
Chart inMap _ _ -> let ((u,(v,w)), epsP) = runFlatContinuous inMap t
in (IdChart, ((u,v),w), epsP)
regroup' = Continuous $ \c t -> case c of
IdChart -> let ((u,v),w) = t in (IdChart, (u,(v,w)), return)
Chart inMap _ _ -> let (((u,v),w), epsP) = runFlatContinuous inMap t
in (IdChart, (u,(v,w)), epsP)
instance Morphism (:-->) where
first (Continuous f) = Continuous $ \c t -> case c of
IdChart -> let (v,w) = t
(IdChart, v', epsP) = f IdChart v
in (IdChart, (v',w), (/ sqrt 2) >>>
\ε -> fmap getMin $ (fmap Min $ epsP ε)
<>(just $ Min ε) )
second (Continuous g) = Continuous $ \c t -> case c of
IdChart -> let (v,w) = t
(IdChart, w', epsP) = g IdChart w
in (IdChart, (v,w'), (/ sqrt 2) >>>
\ε -> fmap getMin $ (just $ Min ε)
<>(fmap Min $ epsP ε) )
Continuous f *** Continuous g = Continuous $ \c t -> case c of
IdChart -> let (v,w) = t
(IdChart, v', epsPv) = f IdChart v
(IdChart, w', epsPw) = g IdChart w
in (IdChart, (v',w'), (/ sqrt 2) >>>
\ε -> fmap getMin $ (fmap Min $ epsPv ε)
<>(fmap Min $ epsPw ε) )
instance PreArrow (:-->) where
terminal = const__ ()
Continuous f &&& Continuous g = Continuous $ \c v -> case c of
IdChart -> let (IdChart, v', epsPv) = f IdChart v
(IdChart, w', epsPw) = g IdChart v
in (IdChart, (v',w'), (/ sqrt 2) >>>
\ε -> fmap getMin $ (fmap Min $ epsPv ε)
<>(fmap Min $ epsPw ε) )
fst = Continuous $ \c t -> case c of
IdChart -> let (v,_) = t
in (IdChart, v, return)
snd = Continuous $ \c t -> case c of
IdChart -> let (_,v) = t
in (IdChart, v, return)
data Chart :: * -> * where
IdChart :: (FlatManifold v) => Chart v
Chart :: (Manifold m, v ~ TangentSpace m, FlatManifold v) =>
{ chartInMap :: v :--> m
, chartOutMap :: m -> Maybe (m:-->v)
, chartKind :: ChartKind } -> Chart m
data ChartKind = LandlockedChart
| RimChart
type FlatManifold v = (MetricSpace v, Manifold v, v~TangentSpace v)
type EuclidSpace v = (HasBasis v, EqFloating(Scalar v), Eq v)
isInUpperHemi :: EuclidSpace v => v -> Bool
isInUpperHemi v = (snd . head) (decompose v) >= 0
type Atlas m = [Chart m]
class (MetricSpace(TangentSpace m), Metric(TangentSpace m) ~ ℝ) => Manifold m where
type TangentSpace m :: *
type TangentSpace m = m
localAtlas :: m -> Atlas m
vectorSpaceAtlas :: FlatManifold v => v -> Atlas v
vectorSpaceAtlas _ = [IdChart]
instance Manifold () where
type TangentSpace () = ()
localAtlas = vectorSpaceAtlas
instance Manifold Double where
localAtlas = vectorSpaceAtlas
instance ( FlatManifold v₁, FlatManifold v₂, Scalar v₁~Scalar v₂
, MetricSpace (Scalar v₁), Metric (Scalar v₁)~ℝ
, VectorSpace (v₁,v₂), Scalar (v₁,v₂) ~ Scalar v₁
) => Manifold (v₁,v₂) where
localAtlas = vectorSpaceAtlas
type Representsℝ r = (EqFloating r, FlatManifold r, r~Scalar r, r~Metric r)
continuousFlatFunction :: ( FlatManifold d, FlatManifold c, ε~Metric c, δ~Metric d )
=> (d -> (c, ε->Option δ)) -> d:-->c
continuousFlatFunction f = Continuous f'
where f' IdChart x = (IdChart, y, eps2Delta)
where (y, eps2Delta) = f x
f' (Chart inMap _ _) v = (IdChart, y, postEps>=>preEps)
where (v', preEps) = runFlatContinuous inMap v
(y, postEps) = f v'
type CntnRealFunction = Representsℝ r => r :--> r
sin__, cos__, atan__ , exp__ , sinh__, cosh__, tanh__, asinh__ :: CntnRealFunction
sin__ = continuousFlatFunction sin'
where sin' x = (sinx, eps2Delta)
where eps2Delta ε
| ε > 1 + abs sinx = nothing
| otherwise = just $ ε / (dsinx + sqrt ε)
dsinx = abs $ cos x
sinx = sin x
cos__ = continuousFlatFunction cos'
where cos' x = (cosx, eps2Delta)
where eps2Delta ε
| ε > 1 + abs cosx = nothing
| otherwise = just $ ε / (dcosx + sqrt ε)
dcosx = abs $ sin x
cosx = cos x
atan__ = continuousFlatFunction atan'
where atan' x = (atanx, eps2Delta)
where eps2Delta ε
| ε >= pi/2 + abs atanx = nothing
| otherwise = just $ abs x tan (abs atanx ε)
atanx = atan x
exp__ = continuousFlatFunction exp'
where exp' x = (expx, eps2Delta)
where expx = exp x
eps2Delta ε
| x>0, expx*2 == expx = just 0
| otherwise = just $ log (expx + ε) x
sinh__ = continuousFlatFunction sinh'
where sinh' x = (sinhx, eps2Delta)
where eps2Delta ε = just $ asinh (abs sinhx + ε) abs x
sinhx = sinh x
cosh__ = continuousFlatFunction cosh'
where cosh' x = (coshx, eps2Delta)
where eps2Delta ε = just $ acosh (coshx + ε) abs x
coshx = cosh x
tanh__ = continuousFlatFunction tanh'
where tanh' x = (tanhx, eps2Delta)
where eps2Delta ε
| ε >= 1 + abs tanhx = nothing
| otherwise = just $ abs x atanh (abs tanhx ε)
tanhx = tanh x
asinh__ = continuousFlatFunction asinh'
where asinh' x = (asinhx, eps2Delta)
where eps2Delta ε = just $ abs x sinh (abs asinhx ε)
asinhx = asinh x
cntnFuncsCombine :: forall d v c c' c'' ε ε' ε''.
( FlatManifold c, FlatManifold c', FlatManifold c''
, ε ~ Metric c , ε' ~ Metric c' , ε'' ~ Metric c'', ε~ε', ε~ε'' )
=> (c'->c''->(c, ε->(ε',ε''))) -> (d:-->c') -> (d:-->c'') -> d:-->c
cntnFuncsCombine cmb (Continuous f) (Continuous g) = Continuous h
where h ζd u = case (ζc', ζc'') of
(IdChart, IdChart)
-> let (y, epsSplit) = cmb fu gu
fullEps ε = fmap getMin $ (fmap Min $ fEps ε')
<>(fmap Min $ gEps ε'')
where (ε', ε'') = epsSplit ε
in (IdChart, y, fullEps)
(IdChart, Chart c''In _ _)
-> let (y'', c''Eps) = runFlatContinuous c''In gu
(y, epsSplit) = cmb fu y''
fullEps ε = fmap getMin $ (fmap Min $ fEps ε')
<>(fmap Min $ gEps =<< c''Eps ε'')
where (ε', ε'') = epsSplit ε
in (IdChart, y, fullEps)
(Chart c'In _ _, IdChart)
-> let (y', c'Eps) = runFlatContinuous c'In fu
(y, epsSplit) = cmb y' gu
fullEps ε = fmap getMin $ (fmap Min $ fEps =<< c'Eps ε')
<>(fmap Min $ gEps ε'')
where (ε', ε'') = epsSplit ε
in (IdChart, y, fullEps)
(Chart c'In _ _, Chart c''In _ _)
-> let (y', c'Eps) = runFlatContinuous c'In fu
(y'', c''Eps) = runFlatContinuous c''In gu
(y, epsSplit) = cmb y' y''
fullEps ε = fmap getMin $ (fmap Min $ fEps =<< c'Eps ε')
<>(fmap Min $ gEps =<< c''Eps ε'')
where (ε', ε'') = epsSplit ε
in (IdChart, y, fullEps)
where (ζc', fu, fEps) = f ζd u
(ζc'',gu, gEps) = g ζd u
data CntnFuncValue d c = CntnFuncValue { runCntnFuncValue :: d :--> c }
| CntnFuncConst c
instance HasAgent (:-->) where
type AgentVal (:-->) d c = CntnFuncValue d c
alg f = case f $ CntnFuncValue id of
CntnFuncValue q -> q
CntnFuncConst c -> const__ c
f $~ CntnFuncValue g = CntnFuncValue $ f . g
f $~ CntnFuncConst c = CntnFuncConst $ f $ c
instance PointAgent CntnFuncValue (:-->) d c where
point = CntnFuncConst
instance CartesianAgent (:-->) where
alg1to2 f = case f $ CntnFuncValue id of
(CntnFuncConst c₁, CntnFuncConst c₂) -> const__ (c₁, c₂)
(CntnFuncConst c₁, CntnFuncValue f₂)
-> Continuous $ \IdChart x -> let (fx, epsP) = runFlatContinuous f₂ x
in (IdChart, (c₁, fx), epsP)
(CntnFuncValue f₁, CntnFuncConst c₂)
-> Continuous $ \IdChart x -> let (fx, epsP) = runFlatContinuous f₁ x
in (IdChart, (fx, c₂), epsP)
(CntnFuncValue f₁, CntnFuncValue f₂) -> f₁ &&& f₂
alg2to1 f = case f (CntnFuncValue fst) (CntnFuncValue snd) of
CntnFuncConst c -> const__ c
CntnFuncValue f -> f
alg2to2 f = case f (CntnFuncValue fst) (CntnFuncValue snd) of
(CntnFuncConst c₁, CntnFuncConst c₂) -> const__ (c₁, c₂)
(CntnFuncConst c₁, CntnFuncValue f₂)
-> Continuous $ \IdChart x -> let (fx, epsP) = runFlatContinuous f₂ x
in (IdChart, (c₁, fx), epsP)
(CntnFuncValue f₁, CntnFuncConst c₂)
-> Continuous $ \IdChart x -> let (fx, epsP) = runFlatContinuous f₁ x
in (IdChart, (fx, c₂), epsP)
(CntnFuncValue f₁, CntnFuncValue f₂) -> f₁ &&& f₂
cntnFnValsFunc :: ( FlatManifold c, FlatManifold c', Manifold d
, ε~Metric c, ε~Metric c' )
=> (c' -> (c, ε->Option ε)) -> CntnFuncValue d c' -> CntnFuncValue d c
cntnFnValsFunc = ($~) . continuousFlatFunction
cntnFnValsCombine :: forall d c c' c'' ε ε' ε''.
( FlatManifold c, FlatManifold c', FlatManifold c'', Manifold d
, ε ~ Metric c , ε' ~ Metric c' , ε'' ~ Metric c'', ε~ε', ε~ε'' )
=> ( c' -> c'' -> (c, ε -> (ε',(ε',ε''),ε'')) )
-> CntnFuncValue d c' -> CntnFuncValue d c'' -> CntnFuncValue d c
cntnFnValsCombine cmb (CntnFuncValue f) (CntnFuncValue g)
= CntnFuncValue $ cntnFuncsCombine (second (>>> \(_,splε,_)->splε) .: cmb) f g
cntnFnValsCombine cmb (CntnFuncConst p) (CntnFuncConst q)
= CntnFuncConst . fst $ cmb p q
cntnFnValsCombine cmb f (CntnFuncConst q)
= cntnFnValsFunc (\c' -> second (>>> \(ε',_,_)->return ε') $ cmb c' q) f
cntnFnValsCombine cmb (CntnFuncConst p) g
= cntnFnValsFunc (second (>>> \(_,_,ε'')->return ε'') . cmb p) g
instance (Representsℝ r, Manifold d) => Num (CntnFuncValue d r) where
fromInteger = point . fromInteger
(+) = cntnFnValsCombine $ \a b -> (a+b, \ε -> (ε, (ε/2,ε/2), ε))
() = cntnFnValsCombine $ \a b -> (ab, \ε -> (ε, (ε/2,ε/2), ε))
(*) = cntnFnValsCombine $ \a b -> (a*b,
\ε -> ( ε/b
, (ε / (2 * sqrt(2*b^2+ε)), ε / (2 * sqrt(2*a^2+ε)))
, ε/a ))
negate = cntnFnValsFunc $ \x -> (negate x, return)
abs = cntnFnValsFunc $ \x -> (abs x, return)
signum = cntnFnValsFunc $ \x -> (signum x, \ε -> if ε>2 then nothing else just $ abs x)
instance (Representsℝ r, Manifold d) => Fractional (CntnFuncValue d r) where
fromRational = point . fromRational
recip = cntnFnValsFunc $ \x -> let x¹ = recip x
in (x¹, \ε -> just $ abs x recip(ε + abs x¹))
instance (Representsℝ r, Manifold d) => Floating (CntnFuncValue d r) where
pi = point pi
exp x = exp__$~ x
sin x = sin__$~ x
cos x = cos__$~ x
atan x = atan__$~ x
sinh x = sinh__$~ x
cosh x = cosh__$~ x
tanh x = tanh__$~ x
asinh x = asinh__$~ x
log x = continuousFlatFunction ln' $~ x
where ln' x = (lnx, eps2Delta)
where lnx = log x
eps2Delta ε = just $ x exp (lnx ε)
asin x = continuousFlatFunction asin' $~ x
where asin' x = (asinx, eps2Delta)
where asinx = asin x
eps2Delta ε = just $
if ε > pi/2 abs asinx
then 1 abs x
else sin (abs asinx + ε) abs x
acos x = continuousFlatFunction acos' $~ x
where acos' x = (acosx, eps2Delta)
where acosx = acos x
eps2Delta ε = just $
if ε > pi/2 abs (acosx pi/2)
then 1 abs x
else cos (abs acosx + ε) abs x
acosh x = continuousFlatFunction acosh' $~ x
where acosh' x = (acoshx, eps2Delta)
where acoshx = acosh x
eps2Delta ε = just $
if ε > acoshx
then x 1
else x cosh (acoshx ε)
atanh x = continuousFlatFunction atanh' $~ x
where atanh' x = (atanhx, eps2Delta)
where atanhx = atanh x
eps2Delta ε = just $ tanh (abs atanhx + ε) abs x
instance (FlatManifold v, Manifold d) => AdditiveGroup (CntnFuncValue d v) where
zeroV = point zeroV
(^+^) = cntnFnValsCombine $ \a b -> (a^+^b, \ε -> (ε, (ε/2,ε/2), ε))
negateV = cntnFnValsFunc $ \x -> (negateV x, return)
instance ( FlatManifold v, MetricSpace v, Metric v~ℝ, FlatManifold (Scalar v)
, MetricSpace (Scalar v), Metric (Scalar v) ~ ℝ, Manifold d )
=> VectorSpace (CntnFuncValue d v) where
type Scalar (CntnFuncValue d v) = CntnFuncValue d (Scalar v)
(*^) = cntnFnValsCombine
$ \λ v -> ( λ*^v
, \ε -> let l = metric v
λ' = metric λ
in ( ε/l
, ( ε / (2 * sqrt(2 * l^2 + ε))
, ε / (2 * sqrt(2 * λ'^2 + ε)))
, ε / λ' ))
finiteGraphContinℝtoℝ :: GraphWindowSpec -> (Double:-->Double) -> [(Double, Double)]
finiteGraphContinℝtoℝ (GraphWindowSpec{..}) fc
= connect [(x, f x, δyG) | x<-[lBound, rBound] ] [(rBound, fst (f rBound))]
where connect [(x₁, (y₁, eps₁), ε₁), (x₂, (y₂, eps₂), ε₂)]
= case (getOption $ eps₁ ε₁, getOption $ eps₂ ε₂) of
(Nothing, Nothing) -> done
(Just δ₁, Nothing) | δ₁>δxS -> done
| otherwise -> refine
(Nothing, Just δ₂) | δ₂>δxS -> done
| otherwise -> refine
(Just δ₁, Just δ₂) | δ₁>δxS, δ₂>δxS -> done
| otherwise -> refine
where δxS = x₂x₁
m = x₁ + δxS/2
fm@(ym, _) = f m
done = ((x₁, y₁) :)
refine = connect [(x₁, (y₁, eps₁), ε₁), (m, fm, ε')]
. connect [(m, fm, ε'), (x₂, (y₂, eps₂), ε₂)]
ε' = (if δxS < δxG then max (min (abs $ ym y₁) (abs $ ym y₂)) else id)
$ max ε₁ ε₂
f = runFlatContinuous fc
δxG = (rBound lBound) / fromIntegral xResolution
δyG = (tBound bBound) / fromIntegral yResolution
finiteGraphContinℝtoℝ² :: GraphWindowSpec -> (Double:-->(Double, Double)) -> [[(Double, Double)]]
finiteGraphContinℝtoℝ² (GraphWindowSpec{..}) fc
= map (\(tl, tu) -> reCoarsen $ connect (tl, f tl) (tu, f tu) [fst (f tu)]) segments
where connect n₁@(t₁, (p₁, eps₁)) n₂@(t₂, (p₂, eps₂))
| and . catMaybes $ map (getOption . fmap( > t₂ t₁ ) . ($reso)) [eps₁, eps₂]
= (p₁ : )
| m <- (id &&& f) $ midBetween [t₁, t₂] = connect n₁ m . connect m n₂
segments = do
(start, dir) <- [ (Just 0 , 1)
, (go (\_ -> not . inRange) reasonable 1 0, 1 ) ]
foldMap (`explore`dir) start
where explore t₀ dir
| Just ti <- go (\_ -> inRange) reasonable dir t₀
, Just tb <- exitWindow (dir) ti
, Just te <- exitWindow dir ti
= (if dir > 0 then (tb, te) else (te, tb)) : explore te dir
| otherwise = []
where exitWindow = go (\t p -> not $ reasonable t && inRange p) (const True)
go isDone hasHope dir t
| not $ hasHope t = Nothing
| isDone t p = Just t
| Just s <- getOption(epsP $ mobility p)
= go isDone hasHope dir $ t + dir * s
| otherwise = Nothing
where (p, epsP) = f t
f = runFlatContinuous fc
inRange (x, y) = x > lBound && x < rBound && y > bBound && y < tBound
reasonable = (< 1e+250) . abs
mobility = \p -> sqrt $ max (distanceSq p cp₁) (distanceSq p cp₂)
where cp₁ = ( midBetween[lBound, rBound, rBound], midBetween[bBound, tBound, tBound] )
cp₂ = ( midBetween[lBound, lBound, rBound], midBetween[bBound, bBound, tBound] )
resoSq = reso ^ 2
reso = min ( (rBound lBound) / fromIntegral xResolution )
( (tBound bBound) / fromIntegral yResolution ) * 2
firstJust = head . catMaybes
reCoarsen (p₁ : p₂ : ps)
| distanceSq p₁ p₂ > resoSq = p₁ : reCoarsen (p₂ : ps)
| otherwise = reCoarsen (p₁ : ps)
reCoarsen ps = ps
midBetween :: (VectorSpace v, Fractional(Scalar v)) => [v] -> v
midBetween vs = sumV vs ^/ (fromIntegral $ Prelude.length vs)
(.:) :: (c->d) -> (a->b->c) -> a->b->d
(.:) = (.) . (.)
just = Option . Just
nothing = Option Nothing
class (RealFloat (Metric v), InnerSpace v) => MetricSpace v where
type Metric v :: *
type Metric v = ℝ
metric :: v -> Metric v
metric = sqrt . metricSq
metricSq :: v -> Metric v
metricSq = (^2) . metric
(|*^) :: Metric v -> v -> v
μ |*^ v = metricToScalar v μ *^ v
metricToScalar :: v -> Metric v -> Scalar v
instance MetricSpace () where
metric = const 0
metricToScalar = const id
instance MetricSpace ℝ where
metric = id
metricToScalar = const id
instance ( RealFloat r, MetricSpace r, Scalar (Complex r)~Metric r )
=> MetricSpace (Complex r) where
type Metric (Complex r) = Metric r
metricSq (a :+ b) = metricSq a + metricSq b
metricToScalar = const id
instance ( MetricSpace v, MetricSpace (Scalar v)
, MetricSpace w, Scalar v~Scalar w
, Metric v~Metric (Scalar v), Metric w~Metric v
, Metric(Scalar w)~Metric v, RealFloat (Metric v)
) => MetricSpace (v,w) where
type Metric (v,w) = Metric v
metricSq (v,w) = metric (magnitudeSq v) + metric (magnitudeSq w)
metricToScalar (v,_) = metricToScalar v