manifolds 0.1.6.3 → 0.2.0.1

raw patch · 12 files changed

+1136/−967 lines, 12 filesdep +trivial-constraintdep ~constrained-categoriesbinary-addedPVP ok

version bump matches the API change (PVP)

Dependencies added: trivial-constraint

Dependency ranges changed: constrained-categories

API changes (from Hackage documentation)

Files

• Data/Function/Affine.hs +330/−37
• Data/Function/Differentiable.hs +126/−92
• Data/Function/Differentiable/Data.hs +26/−13
• Data/LinearMap/Category.hs +111/−1
• Data/LinearMap/HerMetric.hs +120/−71
• Data/Manifold.hs +3/−652
• Data/Manifold/Griddable.hs +4/−4
• Data/Manifold/PseudoAffine.hs +17/−0
• Data/Manifold/TreeCover.hs +393/−94
• Data/Manifold/Types.hs +3/−1
• images/examples/cartesiandisk-2d-ShadeTree.png binary
• manifolds.cabal +3/−2

@@ -19,6 +19,8 @@ {-# LANGUAGE TupleSections            #-} {-# LANGUAGE ConstraintKinds          #-} {-# LANGUAGE PatternGuards            #-}+{-# LANGUAGE PatternSynonyms          #-}+{-# LANGUAGE ViewPatterns             #-} {-# LANGUAGE TypeOperators            #-} {-# LANGUAGE UnicodeSyntax            #-} {-# LANGUAGE MultiWayIf               #-}@@ -28,7 +30,9 @@   module Data.Function.Affine (-              Affine(..)+              Affine+            , linearAffine+            , toOffsetSlope, toOffset'Slope              ) where      @@ -54,7 +58,9 @@ import qualified Prelude import qualified Control.Applicative as Hask +import Data.Constraint.Trivial import Control.Category.Constrained.Prelude hiding ((^))+import Control.Category.Constrained.Reified import Control.Arrow.Constrained import Control.Monad.Constrained import Data.Foldable.Constrained@@ -62,54 +68,343 @@   -data Affine s d c-   = Affine { affineCoOffset :: d-            , affineOffset :: c-            , affineSlope :: Needle d :-* Needle c-            }+data Affine s d c where+   Subtract :: AffineManifold α => Affine s (α,α) (Needle α)+   AddTo :: Affine s (α, Needle α) α+   ScaleWith :: (LinearManifold α, LinearManifold β) => (α:-*β) -> Affine s α β+   ReAffine :: ReWellPointed (Affine s) α β -> Affine s α β -instance (RealDimension s) => EnhancedCat (->) (Affine s) where-  arr (Affine co ao sl) x = ao .+~^ lapply sl (x.-.co)+reAffine :: ReWellPointed (Affine s) α β -> Affine s α β+reAffine (ReWellPointed f) = f+reAffine f = ReAffine f +pattern Specific f = ReWellPointed f+pattern Id = ReAffine WellPointedId+infixr 1 :>>>, :<<<+pattern f :>>> g <- ReAffine (WellPointedCompo (reAffine -> f) (reAffine -> g))+pattern g :<<< f <- ReAffine (WellPointedCompo (reAffine -> f) (reAffine -> g))+pattern Swap = ReAffine WellPointedSwap+pattern AttachUnit = ReAffine WellPointedAttachUnit+pattern DetachUnit = ReAffine WellPointedDetachUnit+pattern Regroup = ReAffine WellPointedRegroup+pattern Regroup' = ReAffine WellPointedRegroup_+pattern Terminal = ReAffine WellPointedTerminal+pattern Fst = ReAffine WellPointedFst+pattern Snd = ReAffine WellPointedSnd+infixr 3 :***, :&&&+pattern f :*** g <- ReAffine (WellPointedPar (reAffine -> f) (reAffine -> g))+pattern f :&&& g <- ReAffine (WellPointedFanout (reAffine -> f) (reAffine -> g))+pattern Const c = ReAffine (WellPointedConst c) -instance (MetricScalar s) => Category (Affine s) where-  type Object (Affine s) o = WithField s LinearManifold o-  id = Affine zeroV zeroV idL-  Affine cof aof slf . Affine cog aog slg-      = Affine cog (aof .+~^ lapply slf (aog.-.cof)) (slf*.*slg) -linearAffine :: ( AdditiveGroup d, AdditiveGroup c-                , HasBasis (Needle d), HasTrie (Basis (Needle d)) )-       => (Needle d -> Needle c) -> Affine s d c-linearAffine = Affine zeroV zeroV . linear+toOffsetSlope :: (MetricScalar s, WithField s LinearManifold d+                                 , WithField s AffineManifold c )+                      => Affine s d c -> (c, Needle d :-* Needle c)+toOffsetSlope f = toOffset'Slope f zeroV +-- | Basically evaluates an affine function as a generic differentiable one,+--   yielding at a given reference point the result and Jacobian. Unlike with+--   'Data.Function.Differentiable.Differentiable', the induced 1st-order Taylor+--   series is equal to the function!+toOffset'Slope :: ( MetricScalar s, WithField s AffineManifold d+                                   , WithField s AffineManifold c )+                      => Affine s d c -> d -> (c, Needle d :-* Needle c)+toOffset'Slope Subtract (a,b) = (a.-.b, linear $ uncurry(^-^))+toOffset'Slope AddTo (p,v) = (p.+^v, linear $ uncurry(^+^))+toOffset'Slope (ScaleWith q) ref = (lapply q ref, q)+toOffset'Slope Id ref = (ref, linear id)+toOffset'Slope (f :>>> g) ref = case toOffset'Slope f ref of+                  (cf,sf) -> case toOffset'Slope g cf of+                     (cg,sg)     -> (cg, sg*.*sf)+toOffset'Slope Swap ref = (swap ref, linear swap)+toOffset'Slope AttachUnit ref = ((ref,Origin), linear (,Origin))+toOffset'Slope DetachUnit ref = (fst ref, linear fst)+toOffset'Slope Regroup ref = (regroup ref, linear regroup)+toOffset'Slope Regroup' ref = (regroup' ref, linear regroup')+toOffset'Slope (f:***g) ref = case ( toOffset'Slope f (fst ref)+                                 , toOffset'Slope g (snd ref) ) of+                  ((cf, sf), (cg, sg)) -> ((cf,cg), linear $ lapply sf *** lapply sg)+toOffset'Slope Terminal ref = (Origin, zeroV)+toOffset'Slope Fst ref = (fst ref, linear fst)+toOffset'Slope Snd ref = (snd ref, linear snd)+toOffset'Slope (f:&&&g) ref = case ( toOffset'Slope (arr f) ref+                                  , toOffset'Slope (arr g) ref ) of+                  ((cf, sf), (cg, sg)) -> ((cf,cg), linear $ lapply sf &&& lapply sg)+toOffset'Slope (Const c) ref = (c, zeroV)+            +coOffsetForm :: ( MetricScalar s, WithField s AffineManifold d+                                , WithField s AffineManifold c )+                      => Affine s d c -> Affine s d c+coOffsetForm (ScaleWith q) = id&&&const zeroV >>> Subtract >>> ScaleWith q+coOffsetForm ((coOffsetForm -> Id:&&&Const cof :>>> Subtract :>>> f) :>>> g)+                    = id&&&const cof >>> Subtract >>> (f >>> g)+coOffsetForm ( (coOffsetForm -> Id:&&&Const cof :>>> Subtract :>>> f)+          :*** (coOffsetForm -> Id:&&&Const cog :>>> Subtract :>>> g) )+     = id&&&const(cof,cog) >>> Subtract >>> (f***g)+coOffsetForm (Id:&&&Const cof :>>> Subtract)+           = (Id&&&Const cof >>> ReAffine (ReWellPointed Subtract`WellPointedCompo`WellPointedId))+coOffsetForm f = f++pattern PreSubtract c f <- (coOffsetForm -> Id:&&&Const c :>>> Subtract :>>> f)++preSubtract :: ( MetricScalar s, WithField s AffineManifold d+                               , WithField s AffineManifold c )+               => c -> Affine s (Diff c) d -> Affine s c d+-- The specialised clauses may not actually be useful here.+preSubtract _ (Const d) = const d+preSubtract _ Terminal = Terminal+preSubtract c (f:>>>g) = preSubtract c f >>>! g+-- preSubtract t (f:***g) | (c,d)<-t = preSubtract c f *** preSubtract d g+preSubtract c (f:&&&g) = preSubtract c f &&& preSubtract c g+preSubtract c f = id&&&const c >>>! Subtract >>>! f+   +pattern PostAdd c f <- f:&&&Const c :>>> AddTo+pattern PostAdd' c f <- Const c:&&&f :>>> AddTo++postAdd :: (MetricScalar s, WithField s AffineManifold d, WithField s AffineManifold c)+               => Diff d -> Affine s c d -> Affine s c d+postAdd c f = f&&&const c >>>! AddTo+postAdd' :: (MetricScalar s, WithField s AffineManifold d, WithField s AffineManifold c)+               => d -> Affine s c (Diff d) -> Affine s c d+postAdd' c f = const c&&&f >>>! AddTo++instance (MetricScalar s) => EnhancedCat (->) (Affine s) where+  arr f = fst . toOffset'Slope f++instance (MetricScalar s) => EnhancedCat (Affine s) (ReWellPointed (Affine s)) where+  arr (Specific c) = c+  arr c = ReAffine c++instance (MetricScalar s, WithField s AffineManifold d, WithField s AffineManifold c)+                  => AffineSpace (Affine s d c) where+  type Diff (Affine s d c) = Affine s d (Diff c)+  +  ScaleWith q .-. ScaleWith r = ScaleWith $ q^-^r+  (PostAdd c (ScaleWith q)) .-. g = let (d, r) = toOffsetSlope g+                                    in postAdd (c.-.d) $ ScaleWith (q^-^r)+  f .-. (PostAdd d (ScaleWith r)) = let (c, q) = toOffsetSlope f+                                    in postAdd (c.-.d) $ ScaleWith (q^-^r)+  (PostAdd' c (ScaleWith q)) .-. g = let (d, r) = toOffsetSlope g+                                     in postAdd (c.-.d) $ ScaleWith (q^-^r)+  f .-. (PostAdd' d (ScaleWith r)) = let (c, q) = toOffsetSlope f+                                     in postAdd (c.-.d) $ ScaleWith (q^-^r)+  +  Id .-. Id = const zeroV+  Fst .-. Fst = const zeroV+  Snd .-. Snd = const zeroV+  Swap .-. Swap = const zeroV+  AttachUnit .-. AttachUnit = const zeroV+  DetachUnit .-. DetachUnit = const zeroV+  Terminal .-. _ = Terminal+  _ .-. Terminal = Terminal+  Subtract .-. Subtract = const zeroV+  AddTo .-. AddTo = const zeroV+  +  Const c .-. Const d = Const $ c.-.d+  +  Fst .-. Snd = Subtract++  (f:***g) .-. (h:***i) = f.-.h *** g.-.i+  (f:***g) .-. Const (c,d) = f.-.const c *** g.-.const d+  ζ .-. (f:***g) | Const (c,d) <- ζ = const c.-.f *** const d.-.g+  (f:&&&g) .-. (h:&&&i) = f.-.h &&& g.-.i+  (f:&&&_) .-. AttachUnit = f.-.id >>>! AttachUnit+  (f:&&&g) .-. Const (c,d) = f.-.const c &&& g.-.const d+  ζ .-. (f:&&&g) | Const (c,d) <- ζ = const c.-.f &&& const d.-.g++  ScaleWith q .-. f = let (c, r) = toOffset'Slope f zeroV+                      in postAdd (negateV c) $ ScaleWith (q^-^r)+  f .-. ScaleWith q = let (c, r) = toOffset'Slope f zeroV+                      in postAdd c $ ScaleWith (r^-^q)+  +  PreSubtract b f .-. g = let (c, q) = toOffsetSlope f+                              (d, r) = toOffset'Slope g b+                          in preSubtract b . postAdd (c.-.d) $ ScaleWith (q^-^r)+      -- f x = q·x + c+      -- g x = r·x + w+      -- d = r·b + w+      -- (q−r)·(x−b) = q·x − q⋅b − r⋅x + r⋅b+      -- s x = f (x−b) − g x+      --     = q⋅(x−b) + c − r⋅x − w+      --     = q⋅x − q⋅b + c − r⋅x − w+      --     = (q−r)·(x−b) + c − r⋅b − w+      --     = (q−r)·(x−b) + c − d+  +  -- According to GHC, this clause overlaps with the above. Hm...+  f .-. PreSubtract b g = let (c, q) = toOffset'Slope f b+                              (d, r) = toOffsetSlope g+                          in preSubtract b $ postAdd (c.-.d) $ ScaleWith (q^-^r)+      -- f x = q·x + v+      -- g x = r·x + d+      -- c = q·b + v+      -- (q−r)·(x−b) = q·x − q⋅b − r⋅x + r⋅b+      -- s x = f x − g (x−b)+      --     = q⋅x + v − r⋅(x−b) − d+      --     = q⋅x + v − r⋅x + r⋅b − d+      --     = (q−r)·(x−b) + q⋅b + v − d+      --     = (q−r)·(x−b) + c − d+  +  f .-. g = f&&&g >>> Subtract+  +  +  ScaleWith q .+^ ScaleWith r = ScaleWith $ q^+^r+  (PostAdd c (ScaleWith q)) .+^ g = let (d, r) = toOffsetSlope g+                                    in postAdd (c.+^d) $ ScaleWith (q^+^r)+  f .+^ (PostAdd d (ScaleWith r)) = let (c, q) = toOffsetSlope f+                                    in postAdd' (c.+^d) $ ScaleWith (q^+^r)+  (PostAdd' c (ScaleWith q)) .+^ g = let (d, r) = toOffsetSlope g+                                     in postAdd' (c.+^d) $ ScaleWith (q^+^r)+  f .+^ (PostAdd' d (ScaleWith r)) = let (c, q) = toOffsetSlope f+                                     in postAdd' (c.+^d) $ ScaleWith (q^+^r)+  (f:***g) .+^ (h:***i) = f.+^h *** g.+^i+  (f:&&&g) .+^ (h:&&&i) = f.+^h &&& g.+^i+  +  Const c .+^ Const c' = const (c.+^c')++  Terminal .+^ _ = Terminal+  Const c .+^ Terminal = Const c+  Const c .+^ f = const c&&&f >>> AddTo+  +  Id .+^ Id = Id >>> ScaleWith (linear (^*2))+  Fst .+^ Fst = Fst >>> ScaleWith (linear (^*2))+  Snd .+^ Snd = Snd >>> ScaleWith (linear (^*2))+  Fst .+^ Snd = AddTo+  Swap .+^ Swap = Swap >>> ScaleWith (linear (^*2))+  +  f .+^ Id = let (c,q) = toOffset'Slope f zeroV+             in const c&&&ScaleWith (q^+^idL) >>>! AddTo+  f .+^ AttachUnit = let (c,q) = toOffset'Slope f zeroV+                     in postAdd' c $ ScaleWith (q^+^linear(,Origin))+  f .+^ DetachUnit = let (c,q) = toOffset'Slope f zeroV+                     in postAdd' c $ ScaleWith (q^+^linear fst)+  f .+^ Swap = let (c,q) = toOffset'Slope f zeroV+               in postAdd' c $ ScaleWith (q^+^linear swap)+  +  PreSubtract b f .+^ g = let (c, q) = toOffsetSlope f+                              (d, r) = toOffset'Slope g b+                          in preSubtract b . postAdd' (c.+^d) $ ScaleWith (q^+^r)+      -- f x = q·x + c+      -- g x = r·x + w+      -- d = r·b + w+      -- (q+r)·(x−b) = q·x − q⋅b + r⋅x − r⋅b+      -- s x = f (x−b) + g x+      --     = q⋅(x−b) + c + r⋅x + w+      --     = q⋅x − q⋅b + c + r⋅x + w+      --     = (q+r)·(x−b) + c + r⋅b + w+      --     = (q−r)·(x−b) + c + d+  +  f .+^ PreSubtract b g = let (c, q) = toOffset'Slope f b+                              (d, r) = toOffsetSlope g+                          in preSubtract b . postAdd' (c.+^d) $ ScaleWith (q^+^r)+      -- f x = q·x + v+      -- g x = r·x + d+      -- c = q·b + v+      -- (q+r)·(x−b) = q·x − q⋅b + r⋅x − r⋅b+      -- s x = f x + g (x−b)+      --     = q⋅x + v + r⋅(x−b) + d+      --     = q⋅x + v + r⋅x − r⋅b + d+      --     = (q+r)·(x−b) + q⋅b + v + d+      --     = (q+r)·(x−b) + c + d+  +  f .+^ g = f&&&g >>> AddTo++++instance (MetricScalar s, WithField s AffineManifold d, WithField s LinearManifold c)+                  => AdditiveGroup (Affine s d c) where+  zeroV = const zeroV+  +  negateV (Const c) = const $ negateV c+  negateV Terminal = Terminal+  negateV (ScaleWith ϕ) = ScaleWith $ negateV ϕ+  negateV (f:***g) = negateV f *** negateV g+  negateV (f:&&&g) = negateV f &&& negateV g+  negateV (f:>>>AddTo) = negateV f >>> AddTo+  negateV (f:>>>Subtract) = (f>>>swap) >>>! Subtract+  negateV (f:>>>ScaleWith ϕ) = negateV f >>>! ScaleWith ϕ+  negateV (f:>>>g) = f >>>! negateV g+  negateV AttachUnit = ScaleWith $ linear (negateV >>> (,Origin))+  negateV Subtract = Swap >>>! Subtract+  negateV f = f >>>! ScaleWith (linear negateV)+  +  (^+^) = (.+^)+  (^-^) = (.-.)+  ++infixr 1 >>>!, <<<!+-- | Affine composition using only the reified skeleton, without trying to be+--   clever in any way.+(>>>!) :: ( MetricScalar s, WithField s AffineManifold α+          , WithField s AffineManifold β, WithField s AffineManifold γ )+      => Affine s α β -> Affine s β γ -> Affine s α γ+ReAffine f >>>! ReAffine g = ReAffine $ f >>> g+f >>>! ReAffine g = ReAffine $ ReWellPointed f >>> g+ReAffine f >>>! g = ReAffine $ f >>> ReWellPointed g+f >>>! g = ReAffine $ ReWellPointed f >>> ReWellPointed g++(<<<!) :: ( MetricScalar s, WithField s AffineManifold α+          , WithField s AffineManifold β, WithField s AffineManifold γ )+      => Affine s β γ -> Affine s α β -> Affine s α γ+(<<<!) = flip (>>>!)++instance (MetricScalar s) => Category (Affine s) where+  type Object (Affine s) o = WithField s AffineManifold o+  +  id = ReAffine id+  +  ScaleWith ϕ . ScaleWith ψ = ScaleWith $ ϕ*.*ψ+  g . ScaleWith ψ = let (d, ϕ) = toOffsetSlope g+                    in postAdd' d $ ScaleWith (ϕ*.*ψ)+  (f:***g) . (h:***i) = f.h *** g.i+  (f:***g) . (h:&&&i) = f.h &&& g.i+  g . (PostAdd' c f) = let (d, ϕ) = toOffset'Slope g c+                      in postAdd' d $ ScaleWith ϕ . f+  +  f . g = f <<<! g+ instance (MetricScalar s) => Cartesian (Affine s) where   type UnitObject (Affine s) = ZeroDim s-  swap = linearAffine swap-  attachUnit = linearAffine (, Origin)-  detachUnit = linearAffine fst-  regroup = linearAffine regroup-  regroup' = linearAffine regroup'+  swap = ReAffine swap+  attachUnit = ReAffine attachUnit+  detachUnit = ReAffine detachUnit+  regroup = ReAffine regroup+  regroup' = ReAffine regroup'  instance (MetricScalar s) => Morphism (Affine s) where-  Affine cof aof slf *** Affine cog aog slg-      = Affine (cof,cog) (aof,aog) (linear $ lapply slf *** lapply slg)+  Const c *** Const c' = const (c,c')+  Terminal *** Terminal = const (mempty, mempty)+  ReAffine f *** ReAffine g = ReAffine $ f *** g+  f *** ReAffine g = ReAffine $ ReWellPointed f *** g+  ReAffine f *** g = ReAffine $ f *** ReWellPointed g+  f *** g = ReAffine $ ReWellPointed f *** ReWellPointed g  instance (MetricScalar s) => PreArrow (Affine s) where-  terminal = linearAffine $ const Origin-  fst = linearAffine fst-  snd = linearAffine snd-  Affine cof aof slf &&& Affine cog aog slg-      = Affine zeroV (aof.-^lapply slf cof, aog.-^lapply slg cog)-                 (linear $ lapply slf &&& lapply slg)+  terminal = ReAffine terminal+  fst = ReAffine fst+  snd = ReAffine snd+  Const c &&& Const c' = const (c,c')+  Terminal &&& Terminal = const (mempty, mempty)+  ReAffine f &&& ReAffine g = ReAffine $ f &&& g+  f &&& ReAffine g = ReAffine $ ReWellPointed f &&& g+  ReAffine f &&& g = ReAffine $ f &&& ReWellPointed g+  f &&& g = ReAffine $ ReWellPointed f &&& ReWellPointed g+        +--   Affine cof aof slf &&& Affine cog aog slg+--       = Affine coh (aof.-^lapply slf rco, aog.+^lapply slg rco)+--                  (linear $ lapply slf &&& lapply slg)+--    where rco = (cog.-.cof)^/2+--          coh = cof .+^ rco  instance (MetricScalar s) => WellPointed (Affine s) where   unit = Tagged Origin-  globalElement x = Affine zeroV x zeroV-  const x = Affine zeroV x zeroV+  const = ReAffine . const  +linearAffine :: (MetricScalar s, WithField s LinearManifold α, WithField s LinearManifold β)+            => (α:-*β) -> Affine s α β+linearAffine = ScaleWith + type AffinFuncValue s = GenericAgent (Affine s)  instance (MetricScalar s) => HasAgent (Affine s) where@@ -127,11 +422,9 @@  instance (WithField s LinearManifold v, WithField s LinearManifold a)     => AdditiveGroup (AffinFuncValue s a v) where-  zeroV = GenericAgent $ Affine zeroV zeroV zeroV-  GenericAgent (Affine cof aof slf) ^+^ GenericAgent (Affine cog aog slg)-       = GenericAgent $ Affine (cof^+^cog) (aof^+^aog) (slf^+^slg)-  negateV (GenericAgent (Affine co ao sl))-      = GenericAgent $ Affine (negateV co) (negateV ao) (negateV sl)+  zeroV = GenericAgent zeroV+  GenericAgent f ^+^ GenericAgent g = GenericAgent $ f ^+^ g+  negateV (GenericAgent f) = GenericAgent $ negateV f   

@@ -55,10 +55,12 @@ import Data.Maybe import Data.Semigroup import Data.Function (on)+import Data.Embedding import Data.Fixed  import Data.VectorSpace import Data.LinearMap+import Data.LinearMap.Category import Data.LinearMap.HerMetric import Data.MemoTrie (HasTrie(..)) import Data.AffineSpace@@ -240,7 +242,8 @@                   = (map (id&&&ivimg) domsL, map (id&&&ivimg) domsR)  where (domsL, domsR) = continuityRanges nLim mx f        ivimg (xl,xr) = go xl 1 i₀ ∪ go xr (-1) i₀-        where (_, Option (Just fdd@(Differentiable fddd))) = fd xc+        where (_, Option (Just fdd@(Differentiable fddd)))+                    = second (fmap genericiseDifferentiable) $ fd xc               xc = (xl+xr)/2               i₀ = minimum&&&maximum $ [fdd$xl, fdd$xc, fdd$xr]               go x dir (a,b)@@ -306,8 +309,9 @@  genericiseDifferentiable :: (LocallyScalable s d, LocallyScalable s c)                     => Differentiable s d c -> Differentiable s d c-genericiseDifferentiable (AffinDiffable (Affine x₀ y₀ f))-     = Differentiable $ \x -> (y₀ .+^ lapply f (x.-.x₀), f, const zeroV)+genericiseDifferentiable (AffinDiffable _ af)+     = Differentiable $ \x -> let (y₀, ϕ) = toOffset'Slope af x+                              in (y₀, ϕ, const zeroV) genericiseDifferentiable f = f  @@ -316,21 +320,23 @@   id = Differentiable $ \x -> (x, idL, const zeroV)   Differentiable f . Differentiable g = Differentiable $      \x -> let (y, g', devg) = g x+               jg = convertLinear $->$ g'                (z, f', devf) = f y-               devfg δz = let δy = transformMetric f' δz+               jf = convertLinear $->$ f'+               devfg δz = let δy = transformMetric jf δz                               εy = devf δz-                          in transformMetric g' εy ^+^ devg δy ^+^ devg εy+                          in transformMetric jg εy ^+^ devg δy ^+^ devg εy            in (z, f'*.*g', devfg)-  AffinDiffable f . AffinDiffable g = AffinDiffable $ f . g+  AffinDiffable ef f . AffinDiffable eg g = AffinDiffable (ef . eg) (f . g)   f . g = genericiseDifferentiable f . genericiseDifferentiable g   -- instance (RealDimension s) => EnhancedCat (Differentiable s) (Affine s) where---   arr (Affine co ao sl) = actuallyAffine (ao .-^ lapply sl co) sl+--   arr (Affine co ao sl) = actuallyAffineEndo (ao .-^ lapply sl co) sl    instance (RealDimension s) => EnhancedCat (->) (Differentiable s) where   arr (Differentiable f) x = let (y,_,_) = f x in y-  arr (AffinDiffable f) x = f $ x+  arr (AffinDiffable _ f) x = f $ x  instance (MetricScalar s) => Cartesian (Differentiable s) where   type UnitObject (Differentiable s) = ZeroDim s@@ -351,14 +357,14 @@    where h (x,y) = ((fx, gy), lPar, devfg)           where (fx, f', devf) = f x                 (gy, g', devg) = g y-                devfg δs = transformMetric lfst δx -                           ^+^ transformMetric lsnd δy-                  where δx = devf $ transformMetric lcofst δs-                        δy = devg $ transformMetric lcosnd δs+                devfg δs = transformMetric fst δx +                           ^+^ transformMetric snd δy+                  where δx = devf $ transformMetric (id&&&zeroV) δs+                        δy = devg $ transformMetric (zeroV&&&id) δs                 lPar = linear $ lapply f'***lapply g'-         lfst = linear fst; lsnd = linear snd-         lcofst = linear (,zeroV); lcosnd = linear (zeroV,)-  AffinDiffable f *** AffinDiffable g = AffinDiffable $ f *** g+  AffinDiffable IsDiffableEndo f *** AffinDiffable IsDiffableEndo g+         = AffinDiffable IsDiffableEndo $ f *** g+  AffinDiffable _ f *** AffinDiffable _ g = AffinDiffable NotDiffableEndo $ f *** g   f *** g = genericiseDifferentiable f *** genericiseDifferentiable g  @@ -372,10 +378,9 @@    where h x = ((fx, gx), lFanout, devfg)           where (fx, f', devf) = f x                 (gx, g', devg) = g x-                devfg δs = (devf $ transformMetric lcofst δs)-                           ^+^ (devg $ transformMetric lcosnd δs)+                devfg δs = (devf $ transformMetric (id&&&zeroV) δs)+                           ^+^ (devg $ transformMetric (zeroV&&&id) δs)                 lFanout = linear $ lapply f'&&&lapply g'-         lcofst = linear (,zeroV); lcosnd = linear (zeroV,)   f &&& g = genericiseDifferentiable f &&& genericiseDifferentiable g  @@ -401,16 +406,22 @@   -actuallyLinear :: ( WithField s LinearManifold x, WithField s LinearManifold y, x~y )+actuallyLinearEndo :: WithField s LinearManifold x+            => (x:-*x) -> Differentiable s x x+actuallyLinearEndo = AffinDiffable IsDiffableEndo . linearAffine++actuallyAffineEndo :: WithField s LinearManifold x+            => x -> (x:-*x) -> Differentiable s x x+actuallyAffineEndo y₀ f = AffinDiffable IsDiffableEndo $ const y₀ .+^ linearAffine f++actuallyLinear :: ( WithField s LinearManifold x, WithField s LinearManifold y )             => (x:-*y) -> Differentiable s x y-actuallyLinear f = actuallyAffine zeroV f+actuallyLinear = AffinDiffable NotDiffableEndo . linearAffine  actuallyAffine :: ( WithField s LinearManifold x-                  , WithField s LinearManifold y -- Really, this should only need `AffineManifold`.-                  , x~y-                  )+                  , WithField s AffineManifold y )             => y -> (x:-*Diff y) -> Differentiable s x y-actuallyAffine y₀ f = AffinDiffable $ Affine zeroV y₀ f+actuallyAffine y₀ f = AffinDiffable NotDiffableEndo $ const y₀ .+^ linearAffine f   -- affinPoint :: (WithField s LinearManifold c, WithField s LinearManifold d)@@ -435,20 +446,21 @@                       (GenericAgent (Differentiable g))      = GenericAgent . Differentiable $         \d -> let (c', f', devf) = f d+                  jf = convertLinear$->$f'                   (c'', g', devg) = g d+                  jg = convertLinear$->$g'                   (c, h', devh) = cmb c' c''-                  h'l = h' *.* lcofst; h'r = h' *.* lcosnd+                  jh = convertLinear$->$h'+                  jhl = jh . (id&&&zeroV); jhr = jh . (zeroV&&&id)               in ( c                  , h' *.* linear (lapply f' &&& lapply g')-                 , \εc -> let εc' = transformMetric h'l εc-                              εc'' = transformMetric h'r εc+                 , \εc -> let εc' = transformMetric jhl εc+                              εc'' = transformMetric jhr εc                               (δc',δc'') = devh εc                            in devf εc' ^+^ devg εc''-                               ^+^ transformMetric f' δc'-                               ^+^ transformMetric g' δc''+                               ^+^ transformMetric jf δc'+                               ^+^ transformMetric jg δc''                  )- where lcofst = linear(,zeroV)-       lcosnd = linear(zeroV,)  dfblFnValsCombine cmb (GenericAgent fa) (GenericAgent ga)           = dfblFnValsCombine cmb (GenericAgent $ genericiseDifferentiable fa)                                  (GenericAgent $ genericiseDifferentiable ga)@@ -460,14 +472,12 @@ instance (WithField s LinearManifold v, LocallyScalable s a, Floating s)     => AdditiveGroup (DfblFuncValue s a v) where   zeroV = point zeroV-  GenericAgent (AffinDiffable f) ^+^ GenericAgent (AffinDiffable g)-       = let (GenericAgent h) = GenericAgent f ^+^ GenericAgent g-         in GenericAgent $ AffinDiffable h+  GenericAgent (AffinDiffable ef f) ^+^ GenericAgent (AffinDiffable eg g)+         = GenericAgent $ AffinDiffable (ef<>eg) (f^+^g)   α^+^β = dfblFnValsCombine (\a b -> (a^+^b, lPlus, const zeroV)) α β       where lPlus = linear $ uncurry (^+^)-  negateV (GenericAgent (AffinDiffable f))-       = let (GenericAgent h) = negateV $ GenericAgent f-         in GenericAgent $ AffinDiffable h+  negateV (GenericAgent (AffinDiffable ef f))+         = GenericAgent $ AffinDiffable ef (negateV f)   negateV α = dfblFnValsFunc (\a -> (negateV a, lNegate, const zeroV)) α       where lNegate = linear negateV   @@ -527,7 +537,7 @@                                ^+^ transformMetric δj d )         where (fx, jf, devf) = f x               (gx, jg, devg) = g x-              δj = jf ^-^ jg+              δj = convertLinear $->$ jf ^-^ jg   postEndo :: ∀ c a b . (HasAgent c, Object c a, Object c b)@@ -620,7 +630,7 @@               xp1² = xp1 ^ 2 negativePreRegion' = PreRegion $ ppr . ngt  where PreRegion ppr = positivePreRegion'-       ngt = actuallyLinear $ linear negate+       ngt = actuallyLinearEndo $ linear negate  preRegionToInfFrom, preRegionFromMinInfTo :: RealDimension s => s -> PreRegion s s preRegionToInfFrom = RealSubray PositiveHalfSphere@@ -629,10 +639,10 @@ preRegionToInfFrom', preRegionFromMinInfTo' :: RealDimension s => s -> PreRegion s s preRegionToInfFrom' xs = PreRegion $ ppr . trl  where PreRegion ppr = positivePreRegion'-       trl = actuallyAffine (-xs) idL+       trl = actuallyAffineEndo (-xs) idL preRegionFromMinInfTo' xe = PreRegion $ ppr . flp  where PreRegion ppr = positivePreRegion'-       flp = actuallyAffine xe (linear negate)+       flp = actuallyAffineEndo xe (linear negate)  intervalPreRegion :: RealDimension s => (s,s) -> PreRegion s s intervalPreRegion (lb,rb) = PreRegion $ Differentiable prr@@ -656,28 +666,29 @@   id = RWDiffable $ \x -> (GlobalRegion, pure id)   RWDiffable f . RWDiffable g = RWDiffable h where    h x₀ = case g x₀ of-           ( rg, Option (Just gr'@(AffinDiffable gr@(Affine cog aog slg))) )-            -> let y₀ = gr $ x₀+           ( rg, Option (Just gr'@(AffinDiffable IsDiffableEndo gr)) )+            -> let (y₀, ϕg) = toOffset'Slope gr x₀                in case f y₀ of-                   (GlobalRegion, Option (Just (AffinDiffable fr)))-                         -> (rg, Option (Just (AffinDiffable (fr.gr))))+                   (GlobalRegion, Option (Just (AffinDiffable fe fr)))+                         -> (rg, Option (Just (AffinDiffable fe (fr.gr))))                    (GlobalRegion, fhr)                          -> (rg, fmap (. gr') fhr)                    (RealSubray diry yl, fhr)                       -> let hhr = fmap (. gr') fhr-                         in case lapply slg 1 of+                         in case lapply ϕg 1 of                               y' | y'>0 -> ( unsafePreRegionIntersect rg-                                                  $ RealSubray diry (cog + (yl-aog)/y')-                                   -- aog + y' * (xl − cog) = yl-                                   -- xl = cog + (yl − aog)/y'+                                                  $ RealSubray diry (x₀ + (yl-y₀)/y')+                                   -- y'⋅(xl−x₀) + y₀ ≝ yl                                            , hhr )                                  | y'<0 -> ( unsafePreRegionIntersect rg                                                   $ RealSubray (otherHalfSphere diry)-                                                               (cog + (yl-aog)/y')+                                                               (x₀ + (yl-y₀)/y')                                            , hhr )                                  | otherwise -> (rg, hhr)                    (PreRegion ry, fhr)                          -> ( PreRegion $ ry . gr', fmap (. gr') fhr )+           ( rg, Option (Just gr'@(AffinDiffable _ gr)) )+            -> error "( rg, Option (Just gr'@(AffinDiffable gr)) )"            (GlobalRegion, Option (Just gr@(Differentiable grd)))             -> let (y₀,_,_) = grd x₀                in case f y₀ of@@ -717,6 +728,7 @@                          -> ( PreRegion $ minDblfuncs (ry . gr) rx                             , notDefinedHere )                    (r, Option (Just fr)) | PreRegion ry <- genericisePreRegion r+                          -> ( PreRegion $ minDblfuncs (ry . gr) rx                             , pure (fr . gr) )            (r, Option Nothing)@@ -816,21 +828,22 @@                       (Option(Just(Differentiable f)), Option(Just(Differentiable g))) ->                         pure . Differentiable $ \d                          -> let (c', f', devf) = f d+                                jf = convertLinear $->$ f'                                 (c'',g', devg) = g d+                                jg = convertLinear $->$ g'                                 (c, h', devh) = cmb c' c''-                                h'l = h' *.* lcofst; h'r = h' *.* lcosnd+                                jh = convertLinear $->$ h'+                                jhl = jh . (id&&&zeroV); jhr = jh . (zeroV&&&id)                             in ( c                                , h' *.* linear (lapply f' &&& lapply g')-                               , \εc -> let εc' = transformMetric h'l εc-                                            εc'' = transformMetric h'r εc+                               , \εc -> let εc' = transformMetric jhl εc+                                            εc'' = transformMetric jhr εc                                             (δc',δc'') = devh εc                                          in devf εc' ^+^ devg εc''-                                             ^+^ transformMetric f' δc'-                                             ^+^ transformMetric g' δc''+                                             ^+^ transformMetric jf δc'+                                             ^+^ transformMetric jg δc''                                )                       _ -> notDefinedHere- where lcofst = linear(,zeroV)-       lcosnd = linear(zeroV,)  grwDfblFnValsCombine cmb fv gv         = grwDfblFnValsCombine cmb (genericiseRWDFV fv) (genericiseRWDFV gv) @@ -849,18 +862,18 @@                  where hd x = (fx^+^gx, jf^+^jg, \ε -> δf(ε^*4) ^+^ δg(ε^*4))                         where (fx, jf, δf) = fd x                               (gx, jg, δg) = gd x-                fgplus (Differentiable fd) (AffinDiffable ga@(Affine cog aog slg))+                fgplus (Differentiable fd) (AffinDiffable _ ga)                                  = Differentiable hd-                 where hd x = (fx^+^gx, jf^+^slg, δf)+                 where hd x = (fx^+^gx, jf^+^ϕg, δf)                         where (fx, jf, δf) = fd x-                              gx = ga $ x-                fgplus (AffinDiffable fa@(Affine cof aof slf)) (Differentiable gd)+                              (gx, ϕg) = toOffset'Slope ga x+                fgplus (AffinDiffable _ fa) (Differentiable gd)                                  = Differentiable hd-                 where hd x = (fx^+^gx, slf^+^jg, δg)+                 where hd x = (fx^+^gx, ϕf^+^jg, δg)                         where (gx, jg, δg) = gd x-                              fx = fa $ x-                fgplus (AffinDiffable fa) (AffinDiffable ga) = AffinDiffable ha-                 where (GenericAgent ha) = GenericAgent fa ^+^ GenericAgent ga+                              (fx, ϕf) = toOffset'Slope fa x+                fgplus (AffinDiffable fe fa) (AffinDiffable ge ga)+                           = AffinDiffable (fe<>ge) (fa^+^ga)  rwDfbl_negateV :: ∀ s a v .         ( WithField s EuclidSpace v, AdditiveGroup v, v ~ Needle (Interior (Needle v))@@ -873,8 +886,7 @@                 fneg (Differentiable fd) = Differentiable hd                  where hd x = (negateV fx, negateV jf, δf)                         where (fx, jf, δf) = fd x-                fneg (AffinDiffable (Affine cof aof slf))-                        = AffinDiffable $ Affine (negateV cof) (negateV aof) (negateV slf)+                fneg (AffinDiffable ef af) = AffinDiffable ef $ negateV af  postCompRW :: ( RealDimension s               , LocallyScalable s a, LocallyScalable s b, LocallyScalable s c )@@ -891,16 +903,17 @@   zeroV = point zeroV   ConstRWDFV c₁ ^+^ ConstRWDFV c₂ = ConstRWDFV (c₁^+^c₂)   ConstRWDFV c₁ ^+^ RWDFV_IdVar = GenericRWDFV $-                               globalDiffable' (actuallyAffine c₁ idL)+                               globalDiffable' (actuallyAffineEndo c₁ idL)   RWDFV_IdVar ^+^ ConstRWDFV c₂ = GenericRWDFV $-                               globalDiffable' (actuallyAffine c₂ idL)+                               globalDiffable' (actuallyAffineEndo c₂ idL)   ConstRWDFV c₁ ^+^ GenericRWDFV g = GenericRWDFV $-                               globalDiffable' (actuallyAffine c₁ idL) . g+                               globalDiffable' (actuallyAffineEndo c₁ idL) . g   GenericRWDFV f ^+^ ConstRWDFV c₂ = GenericRWDFV $-                                  globalDiffable' (actuallyAffine c₂ idL) . f-  GenericRWDFV f ^+^ GenericRWDFV g = GenericRWDFV $ rwDfbl_plus f g+                                  globalDiffable' (actuallyAffineEndo c₂ idL) . f+  fa^+^ga | GenericRWDFV f <- genericiseRWDFV fa+          , GenericRWDFV g <- genericiseRWDFV ga = GenericRWDFV $ rwDfbl_plus f g   negateV (ConstRWDFV c) = ConstRWDFV (negateV c)-  negateV RWDFV_IdVar = GenericRWDFV $ globalDiffable' (actuallyLinear $ linear negateV)+  negateV RWDFV_IdVar = GenericRWDFV $ globalDiffable' (actuallyLinearEndo $ linear negateV)   negateV (GenericRWDFV f) = GenericRWDFV $ rwDfbl_negateV f  instance (RealDimension n, LocallyScalable n a)@@ -909,13 +922,13 @@   (+) = (^+^)   ConstRWDFV c₁ * ConstRWDFV c₂ = ConstRWDFV (c₁*c₂)   ConstRWDFV c₁ * RWDFV_IdVar = GenericRWDFV $-                               globalDiffable' (actuallyLinear $ linear (c₁*))+                               globalDiffable' (actuallyLinearEndo $ linear (c₁*))   RWDFV_IdVar * ConstRWDFV c₂ = GenericRWDFV $-                               globalDiffable' (actuallyLinear $ linear (*c₂))+                               globalDiffable' (actuallyLinearEndo $ linear (*c₂))   ConstRWDFV c₁ * GenericRWDFV g = GenericRWDFV $-                               globalDiffable' (actuallyLinear $ linear (c₁*)) . g+                               globalDiffable' (actuallyLinearEndo $ linear (c₁*)) . g   GenericRWDFV f * ConstRWDFV c₂ = GenericRWDFV $-                                  globalDiffable' (actuallyLinear $ linear (*c₂)) . f+                                  globalDiffable' (actuallyLinearEndo $ linear (*c₂)) . f   f*g = genericiseRWDFV f ⋅ genericiseRWDFV g    where (⋅) :: ∀ n a . (RealDimension n, LocallyScalable n a)            => RWDfblFuncValue n a n -> RWDfblFuncValue n a n -> RWDfblFuncValue n a n @@ -925,19 +938,26 @@                           (rc₂,gmay) = gpcs d₀                       in (unsafePreRegionIntersect rc₁ rc₂, mulDi <$> fmay <*> gmay)           where mulDi :: Differentiable n a n -> Differentiable n a n -> Differentiable n a n-                mulDi (AffinDiffable f@(Affine _ _ slf)) (AffinDiffable g@(Affine _ _ slg))-                   = let f' = lapply slf 1; g' = lapply slg 1+                mulDi f@(AffinDiffable ef af) g@(AffinDiffable eg ag) = case ef<>eg of+                   IsDiffableEndo ->+                  {- let f' = lapply slf 1; g' = lapply slg 1                      in case f'*g' of-                          0 -> AffinDiffable undefined-                          f'g' -> Differentiable $-                           \d -> let c₁ = f $ d; c₂ = g $ d-                                 in ( c₁*c₂-                                    , linear.(*)$ c₁*g' + c₂*f'-                                    , unsafe_dev_ε_δ "*" $ sqrt . (/f'g') )+                          0 -> AffinDiffableEndo $ const (aof*aog)+                          f'g' -> -} Differentiable $+                           \d -> let (fd,ϕf) = toOffset'Slope af d+                                     (gd,ϕg) = toOffset'Slope ag d+                                     f' = lapply ϕf 1; g' = lapply ϕg 1+                                     invf'g' = recip $ f'*g'+                                 in ( fd*gd+                                    , linear.(*)$ fd*g' + gd*f'+                                    , unsafe_dev_ε_δ "*" $ sqrt . (*invf'g') )+                   _ -> mulDi (genericiseDifferentiable f) (genericiseDifferentiable g)                 mulDi (Differentiable f) (Differentiable g)                    = Differentiable $                        \d -> let (c₁, slf, devf) = f d+                                 jf = convertLinear$->$slf                                  (c₂, slg, devg) = g d+                                 jg = convertLinear$->$slg                                  c = c₁*c₂; c₁² = c₁^2; c₂² = c₂^2                                  h' = c₁*^slg ^+^ c₂*^slf                                  in ( c@@ -945,7 +965,7 @@                                     , \εc -> let rε² = metric εc 1                                                  c₁worst² = c₁² + recip(1 + c₂²*rε²)                                                  c₂worst² = c₂² + recip(1 + c₁²*rε²)-                                             in (4*rε²) *^ dualCoCoProduct slf slg+                                             in (4*rε²) *^ dualCoCoProduct jf jg                                                 ^+^ devf (εc^*(4*c₂worst²))                                                 ^+^ devg (εc^*(4*c₁worst²))                     -- TODO: add formal proof for this (or, if necessary, the correct form)@@ -957,8 +977,8 @@    where absPW a₀           | a₀<0       = (negativePreRegion, pure desc)           | otherwise  = (positivePreRegion, pure asc)-         desc = actuallyLinear $ linear negate-         asc = actuallyLinear idL+         desc = actuallyLinearEndo $ linear negate+         asc = actuallyLinearEndo idL   signum = (RWDiffable sgnPW $~)    where sgnPW a₀           | a₀<0       = (negativePreRegion, pure (const $ -1))@@ -1049,7 +1069,7 @@                     -- Safety margins for overlap between quadratic and cubic model                     -- (these aren't naturally compatible to be used both together)                       -  cos = sin . (globalDiffable' (actuallyAffine (pi/2) idL) $~)+  cos = sin . (globalDiffable' (actuallyAffineEndo (pi/2) idL) $~)      sinh x = (exp x - exp (-x))/2     {- = grwDfblFnValsFunc sinhDfb@@ -1179,7 +1199,7 @@  positiveRegionalId :: RealDimension n => RWDiffable n n n positiveRegionalId = RWDiffable $ \x₀ ->-       if x₀ > 0 then (positivePreRegion, pure . AffinDiffable $ id)+       if x₀ > 0 then (positivePreRegion, pure . AffinDiffable IsDiffableEndo $ id)                  else (negativePreRegion, notDefinedHere)  infixl 5 ?> , ?<@@ -1196,10 +1216,12 @@ (?<) :: (RealDimension n, LocallyScalable n a)            => RWDfblFuncValue n a n -> RWDfblFuncValue n a n -> RWDfblFuncValue n a n ConstRWDFV a ?< RWDFV_IdVar = GenericRWDFV . RWDiffable $-       \x₀ -> if a < x₀ then (preRegionToInfFrom a, pure . AffinDiffable $ id)+       \x₀ -> if a < x₀ then ( preRegionToInfFrom a+                             , pure . AffinDiffable IsDiffableEndo $ id)                         else (preRegionFromMinInfTo a, notDefinedHere) RWDFV_IdVar ?< ConstRWDFV a = GenericRWDFV . RWDiffable $-       \x₀ -> if x₀ < a then (preRegionFromMinInfTo a, pure . AffinDiffable $ const a)+       \x₀ -> if x₀ < a then ( preRegionFromMinInfTo a+                             , pure . AffinDiffable IsDiffableEndo $ const a)                         else (preRegionToInfFrom a, notDefinedHere) a ?< b = (positiveRegionalId $~ b-a) ?-> b @@ -1238,6 +1260,18 @@                 (rf, q@(Option (Just _))) -> (rf, q)                 (rf, Option Nothing) | (rg, q) <- g x₀                         -> (unsafePreRegionIntersect rf rg, q)++++++-- | Like 'Data.VectorSpace.lerp', but gives a differentiable function+--   instead of a Hask one.+lerp_diffable :: (WithField s LinearManifold m, RealDimension s)+      => m -> m -> Differentiable s s m+lerp_diffable a b = actuallyAffine a $ linear (*^(b.-.a))++   

@@ -12,8 +12,10 @@ import Data.Manifold.Types.Primitive import Data.Manifold.PseudoAffine +import qualified Control.Category.Constrained as CC  + type LinDevPropag d c = Metric c -> Metric d  @@ -38,18 +40,19 @@ --   and actually be sure you get /all/ solutions correctly, not just /some/ that are --   (hopefully) the closest to some reference point you'd need to laborously define! -- ---   Unfortunately however, this also prevents doing any serious algebra etc. with the---   category, because even something as simple as division necessary introduces singularities---   where the derivatives must diverge.---   Not to speak of many trigonometric e.g. trigonometric functions that---   are undefined on whole regions. The 'PWDiffable' and 'RWDiffable' categories have explicit+--   Unfortunately however, this also prevents doing any serious algebra with the+--   category, because even something as simple as division necessary introduces+--   singularities where the derivatives must diverge.+--   Not to speak of many e.g. trigonometric functions that are undefined+--   on whole regions. The 'PWDiffable' and 'RWDiffable' categories have explicit --   handling for those issues built in; you may simply use these categories even when --   you know the result will be smooth in your relevant domain (or must be, for e.g. --   physics reasons). --   ---   &#xb9;(The implementation does not deal with /&#x3b5;/ and /&#x3b4;/ as difference-bounding---   reals, but rather as metric tensors that define a boundary by prohibiting the---   overlap from exceeding one; this makes the concept actually work on general manifolds.)+--   &#xb9;(The implementation does not deal with /&#x3b5;/ and /&#x3b4;/ as+--   difference-bounding reals, but rather as metric tensors which define a+--   boundary by prohibiting the overlap from exceeding one.+--   This makes the category actually work on general manifolds.) data Differentiable s d c where    Differentiable :: ( d -> ( c   -- function value                             , Needle d :-* Needle c -- Jacobian@@ -59,16 +62,26 @@                                                -- some error margin                               ) )                   -> Differentiable s d c-   AffinDiffable :: LinearManifold d-               => Affine s d d -> Differentiable s d d-                      -- This should ideally map between two general affine spaces,-                      -- but since the special case of affine functions is mostly relevant-                      -- to optimise the propagation of real intervals, we don't do that.+   AffinDiffable :: (AffineManifold d, AffineManifold c)+               => DiffableEndoProof d c -> Affine s d c -> Differentiable s d c    +data DiffableEndoProof d c where+  IsDiffableEndo :: DiffableEndoProof d d+  NotDiffableEndo :: DiffableEndoProof d c +instance Semigroup (DiffableEndoProof d c) where+  IsDiffableEndo <> _ = IsDiffableEndo+  _ <> IsDiffableEndo = IsDiffableEndo+  _ <> _ = NotDiffableEndo+  ++instance CC.Category DiffableEndoProof where+  id = IsDiffableEndo+  IsDiffableEndo . IsDiffableEndo = IsDiffableEndo+  _ . _ = NotDiffableEndo   -- | A pathwise connected subset of a manifold @m@, whose tangent space has scalar @s@.

@@ -35,6 +35,7 @@  import Data.MemoTrie import Data.VectorSpace+import Data.LinearMap import Data.VectorSpace.FiniteDimensional import Data.AffineSpace import Data.Basis@@ -59,12 +60,27 @@       -- | A linear mapping between finite-dimensional spaces, implemeted as a dense matrix.-data Linear s a b = DenseLinear { getDenseMatrix :: HMat.Matrix s }+--  +--   Note that this is equivalent to the tensor product @'DualSpace' a ⊗ b@. One+--   of the types should be deprecated in the future, or either implemented in+--   terms of the other.+newtype Linear s a b = DenseLinear { getDenseMatrix :: HMat.Matrix s }  identMat :: forall v w . FiniteDimensional v => Linear (Scalar v) w v identMat = DenseLinear $ HMat.ident n  where (Tagged n) = dimension :: Tagged v Int +convertLinear :: ∀ v w s . ( FiniteDimensional v, FiniteDimensional w+                           , Scalar v ~ s, Scalar w ~ s )+                   => Isomorphism (->) (v:-*w) (Linear s v w)+convertLinear = Isomorphism (asPackedMatrix >>> DenseLinear)+                            (fromPackedMatrix<<<getDenseMatrix)++denseLinear :: ∀ v w s . (FiniteDimensional v, FiniteDimensional w, Scalar w ~ s)+                   => (v->w) -> Linear s v w+denseLinear f = DenseLinear . HMat.fromColumns $ (asPackedVector . f . basisValue) <$> cbv+ where Tagged cbv = completeBasis :: Tagged v [Basis v]+ instance (SmoothScalar s) => Category (Linear s) where   type Object (Linear s) v = (FiniteDimensional v, Scalar v~s)   id = identMat@@ -118,6 +134,100 @@ instance (SmoothScalar s) => EnhancedCat (->) (Linear s) where   arr (DenseLinear mat) = fromPackedVector . HMat.app mat . asPackedVector +type DenseLinearFuncValue s = GenericAgent (Linear s)++instance (SmoothScalar s) => HasAgent (Linear s) where+  alg = genericAlg+  ($~) = genericAgentMap+instance (SmoothScalar s) => CartesianAgent (Linear s) where+  alg1to2 = genericAlg1to2+  alg2to1 = genericAlg2to1+  alg2to2 = genericAlg2to2+++instance (FiniteDimensional v, Scalar v~s, FiniteDimensional w, Scalar w~s, SmoothScalar s)+                     => AffineSpace (Linear s v w) where+  type Diff (Linear s v w) = Linear s v w+  DenseLinear m.-.DenseLinear n = DenseLinear (m-n)+  DenseLinear m.+^DenseLinear n = DenseLinear (m+n)++instance (FiniteDimensional v, Scalar v~s, FiniteDimensional w, Scalar w~s, SmoothScalar s)+                       => AdditiveGroup (Linear s v w) where+  zeroV = zx+   where zx :: ∀ v w . (FiniteDimensional v, FiniteDimensional w) => Linear s v w+         zx = DenseLinear $ HMat.konst 0 (dw,dv)+          where Tagged dv = dimension :: Tagged v Int+                Tagged dw = dimension :: Tagged w Int+  negateV (DenseLinear m) = DenseLinear $ negate m+  DenseLinear m^+^DenseLinear n = DenseLinear (m+n)+  DenseLinear m^-^DenseLinear n = DenseLinear (m-n)++instance (FiniteDimensional v, Scalar v~s, FiniteDimensional w, Scalar w~s, SmoothScalar s)+             => VectorSpace (Linear s v w) where+  type Scalar (Linear s v w) = s+  μ *^ DenseLinear m = DenseLinear $ HMat.scale μ m++instance (FiniteDimensional v, Scalar v~s, FiniteDimensional w, Scalar w~s, SmoothScalar s)+             => HasBasis (Linear s v w) where+  type Basis (Linear s v w) = (Basis v, Basis w)+  basisValue = bx+   where bx :: ∀ v w . (FiniteDimensional v, FiniteDimensional w)+                          => (Basis v, Basis w)->Linear s v w+         bx = \(bv,bw) -> DenseLinear $ HMat.assoc (dw,dv) 0 [((biw bw, biv bv),1)]+          where Tagged dv = dimension :: Tagged v Int+                Tagged dw = dimension :: Tagged w Int+                Tagged biv = basisIndex :: Tagged v (Basis v->Int)+                Tagged biw = basisIndex :: Tagged w (Basis w->Int)+  decompose = dc+   where dc :: ∀ s v w . ( FiniteDimensional v, Scalar v ~ s+                         , FiniteDimensional w, Scalar w ~ s )+                 => Linear s v w -> [((Basis v, Basis w), s)]+         dc lm = map (id &&& decompose' lm) cb+          where Tagged cb = completeBasis :: Tagged (Linear s v w) [(Basis v, Basis w)]+  decompose' = dc+   where dc :: ∀ s v w . (FiniteDimensional v, FiniteDimensional w, Scalar w ~ s)+               => Linear s v w -> (Basis v, Basis w) -> s+         dc (DenseLinear m) = \(bv,bw) -> m HMat.! biw bw HMat.! biv bv+          where Tagged biv = basisIndex :: Tagged v (Basis v->Int)+                Tagged biw = basisIndex :: Tagged w (Basis w->Int)++instance (FiniteDimensional v, Scalar v ~ s, FiniteDimensional w, Scalar w ~ s)+                => FiniteDimensional (Linear s v w) where+  dimension = d+   where d :: ∀ s v w . (FiniteDimensional v, FiniteDimensional w)+               => Tagged (Linear s v w) Int+         d = Tagged (dv*dw)+          where Tagged dv = dimension::Tagged v Int; Tagged dw = dimension::Tagged w Int+  basisIndex = bi+   where bi :: ∀ s v w . (FiniteDimensional v, FiniteDimensional w)+               => Tagged (Linear s v w) ((Basis v, Basis w) -> Int)+         bi = Tagged $ \(bv,bw) -> dv * biv bv + biw bw where +          Tagged dv=dimension::Tagged v Int; Tagged biv=basisIndex::Tagged v (Basis v->Int)+          Tagged biw = basisIndex :: Tagged w (Basis w -> Int)+  indexBasis = ib+   where ib :: ∀ s v w . (FiniteDimensional v, FiniteDimensional w)+               => Tagged (Linear s v w) (Int -> (Basis v, Basis w))+         ib = Tagged $ (`divMod`dv) >>> \(iv,iw) -> (ibv iv, ibw iw) where+          Tagged dv=dimension::Tagged v Int; Tagged ibv=indexBasis::Tagged v (Int->Basis v)+          Tagged ibw = indexBasis :: Tagged w (Int->Basis w)+  completeBasis = cb+   where cb :: ∀ s v w . (FiniteDimensional v, FiniteDimensional w)+               => Tagged (Linear s v w) [(Basis v, Basis w)]+         cb = Tagged $ liftA2 (,) cbv cbw where+          Tagged cbv = completeBasis :: Tagged v [Basis v]+          Tagged cbw = completeBasis :: Tagged w [Basis w]+  asPackedVector = getDenseMatrix >>> HMat.flatten+  fromPackedVector = fpv+   where fpv :: ∀ s v w . (FiniteDimensional v, Scalar v ~ s, FiniteDimensional w, Scalar w ~ s)+               => HMat.Vector s -> Linear s v w+         fpv = HMat.reshape dv >>> DenseLinear+          where Tagged dv = dimension :: Tagged v Int++instance (FiniteDimensional v, Scalar v ~ s, FiniteDimensional a, Scalar a ~ s)+    => AdditiveGroup (DenseLinearFuncValue s a v) where+  zeroV = GenericAgent zeroV+  GenericAgent f ^+^ GenericAgent g = GenericAgent $ f ^+^ g+  negateV (GenericAgent f) = GenericAgent $ negateV f   

@@ -9,14 +9,16 @@ {-# LANGUAGE ConstraintKinds            #-} {-# LANGUAGE ScopedTypeVariables        #-} {-# LANGUAGE UnicodeSyntax              #-}+{-# LANGUAGE LambdaCase                 #-}     module Data.LinearMap.HerMetric (   -- * Metric operator types-    HerMetric, HerMetric'+    HerMetric(..), HerMetric'(..)   -- * Evaluating metrics+  , toDualWith, fromDualWith   , metricSq, metricSq', metric, metric', metrics, metrics'   -- * Defining metrics   , projector, projector'@@ -39,6 +41,8 @@   , metriScale', metriScale   , adjoint   , extendMetric+  , applyLinMapMetric, applyLinMapMetric'+  , imitateMetricSpanChange   -- * The dual-space class   , HasMetric   , HasMetric'(..)@@ -99,22 +103,22 @@ --   Yet other possible interpretations of this type include /density matrix/ (as in --   quantum mechanics), /standard range of statistical fluctuations/, and /volume element/. newtype HerMetric v = HerMetric {-   -- morally:  @getHerMetric :: v :-* DualSpace v@.-          metricMatrix :: Maybe (HMat.Matrix (Scalar v)) -- @Nothing@ for zero metric.+          metricMatrix :: Maybe (Linear (Scalar v) v (DualSpace v)) -- @Nothing@ for zero metric.                       }  matrixMetric :: HasMetric v => HMat.Matrix (Scalar v) -> HerMetric v-matrixMetric = HerMetric . Just+matrixMetric = HerMetric . Just . DenseLinear +-- | Deprecated (this doesn't preserve positive-definiteness) instance (HasMetric v) => AdditiveGroup (HerMetric v) where   zeroV = HerMetric Nothing-  negateV (HerMetric m) = HerMetric $ negate <$> m+  negateV (HerMetric m) = HerMetric $ negateV <$> m   HerMetric Nothing ^+^ HerMetric n = HerMetric n   HerMetric m ^+^ HerMetric Nothing = HerMetric m-  HerMetric (Just m) ^+^ HerMetric (Just n) = HerMetric . Just $ m + n+  HerMetric (Just m) ^+^ HerMetric (Just n) = HerMetric . Just $ m ^+^ n instance HasMetric v => VectorSpace (HerMetric v) where   type Scalar (HerMetric v) = Scalar v-  s *^ (HerMetric m) = HerMetric $ HMat.scale s <$> m +  s *^ (HerMetric m) = HerMetric $ (s*^) <$> m   -- | A metric on the dual space; equivalent to a linear mapping from the dual space --   to the original vector space.@@ -122,15 +126,15 @@ --   Prime-versions of the functions in this module target those dual-space metrics, so --   we can avoid some explicit handling of double-dual spaces. newtype HerMetric' v = HerMetric' {-          metricMatrix' :: Maybe (HMat.Matrix (Scalar v))+          metricMatrix' :: Maybe (Linear (Scalar v) (DualSpace v) v)                       }  extendMetric :: (HasMetric v, Scalar v~ℝ) => HerMetric v -> v -> HerMetric v extendMetric (HerMetric Nothing) _ = HerMetric Nothing-extendMetric (HerMetric (Just m)) v-      | isInfinite' detm  = HerMetric $ Just m+extendMetric (HerMetric (Just (DenseLinear m))) v+      | isInfinite' detm  = HerMetric . Just $ DenseLinear m       | isInfinite' detmninv  = singularMetric-      | otherwise         = HerMetric $ Just mn+      | otherwise         = HerMetric . Just $ DenseLinear mn  where -- this could probably be done much more efficiently, with only        -- multiplications, no inverses.        (minv, (detm, _)) = HMat.invlndet m@@ -139,17 +143,18 @@                                 matrixMetric' :: HasMetric v => HMat.Matrix (Scalar v) -> HerMetric' v-matrixMetric' = HerMetric' . Just+matrixMetric' = HerMetric' . Just . DenseLinear +-- | Deprecated instance (HasMetric v) => AdditiveGroup (HerMetric' v) where   zeroV = HerMetric' Nothing-  negateV (HerMetric' m) = HerMetric' $ negate <$> m+  negateV (HerMetric' m) = HerMetric' $ negateV <$> m   HerMetric' Nothing ^+^ HerMetric' n = HerMetric' n   HerMetric' m ^+^ HerMetric' Nothing = HerMetric' m-  HerMetric' (Just m) ^+^ HerMetric' (Just n) = matrixMetric' $ m + n+  HerMetric' (Just m) ^+^ HerMetric' (Just n) = HerMetric' . Just $ m ^+^ n instance HasMetric v => VectorSpace (HerMetric' v) where   type Scalar (HerMetric' v) = Scalar v-  s *^ (HerMetric' m) = HerMetric' $ HMat.scale s <$> m +  s *^ (HerMetric' m) = HerMetric' $ (s*^) <$> m        -- | A metric on @v@ that simply yields the squared overlap of a vector with the@@ -183,13 +188,13 @@ --   this will be simply 'magnitudeSq'. metricSq :: HasMetric v => HerMetric v -> v -> Scalar v metricSq (HerMetric Nothing) _ = 0-metricSq (HerMetric (Just m)) v = vDecomp `HMat.dot` HMat.app m vDecomp+metricSq (HerMetric (Just (DenseLinear m))) v = vDecomp `HMat.dot` HMat.app m vDecomp  where vDecomp = asPackedVector v   metricSq' :: HasMetric v => HerMetric' v -> DualSpace v -> Scalar v metricSq' (HerMetric' Nothing) _ = 0-metricSq' (HerMetric' (Just m)) u = uDecomp `HMat.dot` HMat.app m uDecomp+metricSq' (HerMetric' (Just (DenseLinear m))) u = uDecomp `HMat.dot` HMat.app m uDecomp  where uDecomp = asPackedVector u  -- | Evaluate a vector's &#x201c;magnitude&#x201d; through a metric. This assumes an actual@@ -205,8 +210,12 @@  toDualWith :: HasMetric v => HerMetric v -> v -> DualSpace v toDualWith (HerMetric Nothing) = const zeroV-toDualWith (HerMetric (Just m)) = fromPackedVector . HMat.app m . asPackedVector+toDualWith (HerMetric (Just m)) = (m$) +fromDualWith :: HasMetric v => HerMetric' v -> DualSpace v -> v+fromDualWith (HerMetric' Nothing) = const zeroV+fromDualWith (HerMetric' (Just m)) = (m$)+ -- | Divide a vector by its own norm, according to metric, i.e. normalise it --   or &#x201c;project to the metric's boundary&#x201d;. metriNormalise :: (HasMetric v, Floating (Scalar v)) => HerMetric v -> v -> v@@ -237,31 +246,27 @@ metrics' m vs = sqrt . sum $ metricSq' m <$> vs  -transformMetric :: (HasMetric v, HasMetric w, Scalar v ~ Scalar w)-           => (w :-* v) -> HerMetric v -> HerMetric w+transformMetric :: ∀ s v w . (HasMetric v, HasMetric w, Scalar v~s, Scalar w~s)+           => Linear s w v -> HerMetric v -> HerMetric w transformMetric _ (HerMetric Nothing) = HerMetric Nothing-transformMetric t (HerMetric (Just m)) = matrixMetric $ HMat.tr tmat HMat.<> m HMat.<> tmat- where tmat = asPackedMatrix t+transformMetric t (HerMetric (Just m)) = HerMetric . Just $ adjoint t . m . t -transformMetric' :: ( HasMetric v, HasMetric w, Scalar v ~ Scalar w )-           => (v :-* w) -> HerMetric' v -> HerMetric' w+transformMetric' :: ∀ s v w . (HasMetric v, HasMetric w, Scalar v~s, Scalar w~s)+           => Linear s v w -> HerMetric' v -> HerMetric' w transformMetric' _ (HerMetric' Nothing) = HerMetric' Nothing-transformMetric' t (HerMetric' (Just m))-                      = matrixMetric' $ tmat HMat.<> m HMat.<> HMat.tr tmat- where tmat = asPackedMatrix t+transformMetric' t (HerMetric' (Just m)) = HerMetric' . Just $ t . m . adjoint t  -- | This does something vaguely like  @\\s t -> (s⋅t)²@, --   but without actually requiring an inner product on the covectors. --   Used for calculating the superaffine term of multiplications in --   'Differentiable' categories.-dualCoCoProduct :: (HasMetric v, HasMetric w, Scalar v ~ Scalar w)-           => (w :-* v) -> (w :-* v) -> HerMetric w-dualCoCoProduct s t = ( (sArr `HMat.dot` (t²PLUSs² HMat.<\> sArr))+dualCoCoProduct :: (HasMetric v, HasMetric w, Scalar v ~ s, Scalar w ~ s)+           => Linear s w v -> Linear s w v -> HerMetric w+dualCoCoProduct (DenseLinear smat) (DenseLinear tmat)+                  = ( (sArr `HMat.dot` (t²PLUSs² HMat.<\> sArr))                        * (tArr `HMat.dot` (t²PLUSs² HMat.<\> tArr)) )                     *^ matrixMetric t²PLUSs²- where tmat = asPackedMatrix t-       tArr = HMat.flatten tmat-       smat = asPackedMatrix s+ where tArr = HMat.flatten tmat        sArr = HMat.flatten smat        t²PLUSs² = tmat HMat.<> HMat.tr tmat + smat HMat.<> HMat.tr smat @@ -278,16 +283,17 @@ -- | The inverse mapping of a metric tensor. Since a metric maps from --   a space to its dual, the inverse maps from the dual into the --   (double-dual) space &#x2013; i.e., it is a metric on the dual space.+--   Deprecated: the singular case isn't properly handled. recipMetric' :: HasMetric v => HerMetric v -> HerMetric' v recipMetric' (HerMetric Nothing) = singularMetric'-recipMetric' (HerMetric (Just m))+recipMetric' (HerMetric (Just (DenseLinear m)))           | isInfinite' detm  = singularMetric'           | otherwise         = matrixMetric' minv  where (minv, (detm, _)) = HMat.invlndet m  recipMetric :: HasMetric v => HerMetric' v -> HerMetric v recipMetric (HerMetric' Nothing) = singularMetric-recipMetric (HerMetric' (Just m))+recipMetric (HerMetric' (Just (DenseLinear m)))           | isInfinite' detm  = singularMetric           | otherwise         = matrixMetric minv  where (minv, (detm, _)) = HMat.invlndet m@@ -309,24 +315,24 @@ --   &#x201c;scaled length&#x201d; doesn't really makes sense then in the usual way!) eigenSpan :: (HasMetric v, Scalar v ~ ℝ) => HerMetric' v -> [v] eigenSpan (HerMetric' Nothing) = []-eigenSpan (HerMetric' (Just m)) = map fromPackedVector eigSpan+eigenSpan (HerMetric' (Just (DenseLinear m))) = map fromPackedVector eigSpan  where (μs,vsm) = HMat.eigSH' m        eigSpan = zipWith (HMat.scale . sqrt) (HMat.toList μs) (HMat.toColumns vsm)  eigenSpan' :: (HasMetric v, Scalar v ~ ℝ) => HerMetric v -> [DualSpace v] eigenSpan' (HerMetric Nothing) = []-eigenSpan' (HerMetric (Just m)) = map fromPackedVector eigSpan+eigenSpan' (HerMetric (Just (DenseLinear m))) = map fromPackedVector eigSpan  where (μs,vsm) = HMat.eigSH' m        eigSpan = zipWith (HMat.scale . sqrt) (HMat.toList μs) (HMat.toColumns vsm)  eigenCoSpan :: (HasMetric v, Scalar v ~ ℝ) => HerMetric' v -> [DualSpace v] eigenCoSpan (HerMetric' Nothing) = []-eigenCoSpan (HerMetric' (Just m)) = map fromPackedVector eigSpan+eigenCoSpan (HerMetric' (Just (DenseLinear m))) = map fromPackedVector eigSpan  where (μs,vsm) = HMat.eigSH' m        eigSpan = zipWith (HMat.scale . recip . sqrt) (HMat.toList μs) (HMat.toColumns vsm) eigenCoSpan' :: (HasMetric v, Scalar v ~ ℝ) => HerMetric v -> [v] eigenCoSpan' (HerMetric Nothing) = []-eigenCoSpan' (HerMetric (Just m)) = map fromPackedVector eigSpan+eigenCoSpan' (HerMetric (Just (DenseLinear m))) = map fromPackedVector eigSpan  where (μs,vsm) = HMat.eigSH' m        eigSpan = zipWith (HMat.scale . recip . sqrt) (HMat.toList μs) (HMat.toColumns vsm) @@ -344,8 +350,7 @@ --   all about dual spaces. class ( FiniteDimensional v, FiniteDimensional (DualSpace v)       , VectorSpace (DualSpace v), HasBasis (DualSpace v)-      , MetricScalar (Scalar v), Scalar v ~ Scalar (DualSpace v)-      , Basis v ~ Basis (DualSpace v) )+      , MetricScalar (Scalar v), Scalar v ~ Scalar (DualSpace v) )     => HasMetric' v where            -- | @'DualSpace' v@ is isomorphic to the space of linear functionals on @v@, i.e.@@ -374,7 +379,15 @@   doubleDual :: HasMetric' (DualSpace v) => v -> DualSpace (DualSpace v)   doubleDual' :: HasMetric' (DualSpace v) => DualSpace (DualSpace v) -> v   +  basisInDual :: Tagged v (Basis v -> Basis (DualSpace v))+  basisInDual = bid+   where bid :: ∀ v . HasMetric' v => Tagged v (Basis v -> Basis (DualSpace v))+         bid = Tagged $ bi >>> ib'+          where Tagged bi = basisIndex :: Tagged v (Basis v -> Int)+                Tagged ib' = indexBasis :: Tagged (DualSpace v) (Int -> Basis (DualSpace v))+   +    -- | Simple flipped version of '<.>^'. (^<.>) :: HasMetric v => v -> DualSpace v -> Scalar v@@ -382,7 +395,7 @@   euclideanMetric' :: forall v . (HasMetric v, InnerSpace v) => HerMetric v-euclideanMetric' = HerMetric . pure $ HMat.ident n+euclideanMetric' = HerMetric . pure . DenseLinear $ HMat.ident n  where (Tagged n) = dimension :: Tagged v Int  -- -- | Associate a Hilbert space vector canonically with its dual-space counterpart,@@ -401,17 +414,24 @@ instance (MetricScalar k) => HasMetric' (ZeroDim k) where   Origin<.>^Origin = zeroV   functional _ = Origin-  doubleDual = id; doubleDual'= id+  doubleDual = id; doubleDual'= id; basisInDual = pure id instance HasMetric' Double where   (<.>^) = (<.>)   functional f = f 1-  doubleDual = id; doubleDual'= id+  doubleDual = id; doubleDual'= id; basisInDual = pure id instance ( HasMetric v, HasMetric w, Scalar v ~ Scalar w          ) => HasMetric' (v,w) where   type DualSpace (v,w) = (DualSpace v, DualSpace w)   (v,w)<.>^(v',w') = v<.>^v' + w<.>^w'   functional f = (functional $ f . (,zeroV), functional $ f . (zeroV,))   doubleDual = id; doubleDual'= id+  basisInDual = bid+   where bid :: ∀ v w . (HasMetric v, HasMetric w) => Tagged (v,w)+                       (Basis v + Basis w -> Basis (DualSpace v) + Basis (DualSpace w))+         bid = Tagged $ \case Left q -> Left $ bidv q+                              Right q -> Right $ bidw q+          where Tagged bidv = basisInDual :: Tagged v (Basis v -> Basis (DualSpace v))+                Tagged bidw = basisInDual :: Tagged w (Basis w -> Basis (DualSpace w)) instance (SmoothScalar s, Ord s, KnownNat n) => HasMetric' (s^n) where   type DualSpace (s^n) = s^n   (<.>^) = (<.>)@@ -420,7 +440,7 @@          fnal f =     FreeVect . Arr.generate n $             \i -> f . FreeVect . Arr.generate n $ \j -> if i==j then 1 else 0           where Tagged n = theNatN :: Tagged n Int-  doubleDual = id; doubleDual'= id+  doubleDual = id; doubleDual'= id; basisInDual = pure id instance (HasMetric v, s~Scalar v) => HasMetric' (FinVecArrRep t v s) where   type DualSpace (FinVecArrRep t v s) = FinVecArrRep t (DualSpace v) s   FinVecArrRep v <.>^ FinVecArrRep w = HMat.dot v w@@ -432,17 +452,38 @@                      $ (f . FinVecArrRep) <$> HMat.toRows (HMat.ident n)          Tagged n = dimension :: Tagged v Int   doubleDual = id; doubleDual'= id+  basisInDual = bid+   where bid :: ∀ s v t . (HasMetric v, s~Scalar v)+                     => Tagged (FinVecArrRep t v s) (Basis v -> Basis (DualSpace v))+         bid = Tagged bid₀+          where Tagged bid₀ = basisInDual :: Tagged v (Basis v -> Basis (DualSpace v)) +instance (HasMetric v, HasMetric w, s ~ Scalar v, s ~ Scalar w)+               => HasMetric' (Linear s v w) where+  type DualSpace (Linear s v w) = Linear s w v+  DenseLinear bw <.>^ DenseLinear fw+                  = HMat.sumElements (HMat.tr bw * fw) -- trace of product+  functional = completeBasisFunctional+  doubleDual = id; doubleDual' = id +completeBasisFunctional :: ∀ v . HasMetric' v => (v -> Scalar v) -> DualSpace v+completeBasisFunctional f = recompose [ (bid b, f $ basisValue b) | b <- cb ]+          where Tagged cb = completeBasis :: Tagged v [Basis v]+                Tagged bid = basisInDual :: Tagged v (Basis v -> Basis (DualSpace v))   + -- | Transpose a linear operator. Contrary to popular belief, this does not --   just inverse the direction of mapping between the spaces, but also switch to --   their duals.-adjoint :: (HasMetric v, HasMetric w, Scalar w ~ Scalar v)+adjoint :: (HasMetric v, HasMetric w, s~Scalar v, s~Scalar w)+     => (Linear s v w) -> Linear s (DualSpace w) (DualSpace v)+adjoint (DenseLinear m) = DenseLinear $ HMat.tr m++adjoint_fln :: (HasMetric v, HasMetric w, Scalar w ~ Scalar v)      => (v :-* w) -> DualSpace w :-* DualSpace v-adjoint m = linear $ \w -> functional $ \v+adjoint_fln m = linear $ \w -> functional $ \v                      -> w <.>^lapply m v  @@ -458,7 +499,8 @@   negate = negateV               -- | This does /not/ work correctly if the metrics don't share an eigenbasis!-  HerMetric m * HerMetric n = HerMetric $ liftA2 (HMat.<>) m n+  HerMetric m * HerMetric n = HerMetric . fmap DenseLinear+                              $ liftA2 (HMat.<>) (getDenseMatrix<$>m) (getDenseMatrix<$>n)                                  -- | Undefined, though it could actually be done.   abs = error "abs undefined for HerMetric"@@ -468,7 +510,8 @@ metrNumFun :: (HasMetric v, v ~ Scalar v, v ~ DualSpace v, Num v)       => (v -> v) -> HerMetric v -> HerMetric v metrNumFun f (HerMetric Nothing) = matrixMetric . HMat.scalar $ f 0-metrNumFun f (HerMetric (Just m)) = matrixMetric . HMat.scalar . f $ m HMat.! 0 HMat.! 0+metrNumFun f (HerMetric (Just (DenseLinear m)))+              = matrixMetric . HMat.scalar . f $ m HMat.! 0 HMat.! 0  instance (HasMetric v, v ~ Scalar v, v ~ DualSpace v, Fractional v)              => Fractional (HerMetric v) where@@ -530,40 +573,46 @@ productMetric :: ∀ v w . (HasMetric v, HasMetric w, Scalar v ~ ℝ, Scalar w ~ ℝ)                => HerMetric v -> HerMetric w -> HerMetric (v,w) productMetric (HerMetric Nothing) (HerMetric Nothing) = HerMetric Nothing-productMetric (HerMetric (Just mv)) (HerMetric (Just mw))-        = HerMetric . Just $ HMat.diagBlock [mv, mw]-productMetric (HerMetric Nothing) (HerMetric (Just mw))-        = HerMetric . Just $ HMat.diagBlock [HMat.konst 0 (dv,dv), mw]- where (Tagged dv) = dimension :: Tagged v Int-productMetric (HerMetric (Just mv)) (HerMetric Nothing)-        = HerMetric . Just $ HMat.diagBlock [mv, HMat.konst 0 (dw,dw)]- where (Tagged dw) = dimension :: Tagged w Int+productMetric (HerMetric (Just mv)) (HerMetric (Just mw)) = HerMetric . Just $ mv *** mw+productMetric (HerMetric Nothing) (HerMetric (Just mw)) = HerMetric . Just $ zeroV *** mw+productMetric (HerMetric (Just mv)) (HerMetric Nothing) = HerMetric . Just $ mv *** zeroV  productMetric' :: ∀ v w . (HasMetric v, HasMetric w, Scalar v ~ ℝ, Scalar w ~ ℝ)                => HerMetric' v -> HerMetric' w -> HerMetric' (v,w) productMetric' (HerMetric' Nothing) (HerMetric' Nothing) = HerMetric' Nothing-productMetric' (HerMetric' (Just mv)) (HerMetric' (Just mw))-        = HerMetric' . Just $ HMat.diagBlock [mv, mw]-productMetric' (HerMetric' Nothing) (HerMetric' (Just mw))-        = HerMetric' . Just $ HMat.diagBlock [HMat.konst 0 (dv,dv), mw]- where (Tagged dv) = dimension :: Tagged v Int-productMetric' (HerMetric' (Just mv)) (HerMetric' Nothing)-        = HerMetric' . Just $ HMat.diagBlock [mv, HMat.konst 0 (dw,dw)]- where (Tagged dw) = dimension :: Tagged w Int+productMetric' (HerMetric' (Just mv)) (HerMetric' (Just mw)) = HerMetric' . Just $ mv***mw+productMetric' (HerMetric' Nothing) (HerMetric' (Just mw)) = HerMetric' . Just $ zeroV***mw+productMetric' (HerMetric' (Just mv)) (HerMetric' Nothing) = HerMetric' . Just $ mv***zeroV  +applyLinMapMetric :: ∀ v w . (HasMetric v, HasMetric w, Scalar v ~ ℝ, Scalar w ~ ℝ)+               => HerMetric (Linear ℝ v w) -> DualSpace v -> HerMetric w+applyLinMapMetric met v' = transformMetric ap2v met+ where ap2v :: Linear ℝ w (Linear ℝ v w)+       ap2v = denseLinear $ \w -> denseLinear $ \v -> w ^* (v'<.>^v) +applyLinMapMetric' :: ∀ v w . (HasMetric v, HasMetric w, Scalar v ~ ℝ, Scalar w ~ ℝ)+               => HerMetric' (Linear ℝ v w) -> v -> HerMetric' w+applyLinMapMetric' met v = transformMetric' ap2v met+ where ap2v :: Linear ℝ (Linear ℝ v w) w+       ap2v = denseLinear ($v) +++imitateMetricSpanChange :: ∀ v . (HasMetric v, Scalar v ~ ℝ)+                           => HerMetric v -> HerMetric' v -> Linear ℝ v v+imitateMetricSpanChange (HerMetric (Just m)) (HerMetric' (Just n)) = n . m+imitateMetricSpanChange _ _ = zeroV++ covariance :: ∀ v w . (HasMetric v, HasMetric w, Scalar v ~ ℝ, Scalar w ~ ℝ)-          => HerMetric' (v,w) -> Option (v:-*w)+          => HerMetric' (v,w) -> Option (Linear ℝ v w) covariance (HerMetric' Nothing) = pure zeroV covariance (HerMetric' (Just m))     | isInfinite' detvnm  = empty-    | otherwise           = pure . fromPackedMatrix $-                               wmat HMat.<> m HMat.<> vmat HMat.<> vnorml- where wmat = asPackedMatrix (linear snd :: (v,w):-*w)-       vmat = asPackedMatrix (linear (id&&&const zeroV) :: v:-*(v,w))-       (vnorml, (detvnm, _)) = HMat.invlndet (HMat.tr vmat HMat.<> m HMat.<> vmat)+    | otherwise           = return $ snd . m . (id&&&zeroV) . DenseLinear vnorml+ where (vnorml, (detvnm, _))+           = HMat.invlndet . getDenseMatrix $ fst . m . (id&&&zeroV)   metricAsLength :: HerMetric ℝ -> ℝ

@@ -15,661 +15,12 @@  {-# LANGUAGE FlexibleInstances        #-} {-# LANGUAGE UndecidableInstances     #-}--- {-# LANGUAGE OverlappingInstances     #-}-{-# LANGUAGE TypeFamilies             #-}-{-# LANGUAGE FunctionalDependencies   #-}-{-# LANGUAGE FlexibleContexts         #-}-{-# LANGUAGE GADTs                    #-}-{-# LANGUAGE RankNTypes               #-}-{-# LANGUAGE TupleSections            #-}-{-# LANGUAGE ConstraintKinds          #-}-{-# LANGUAGE PatternGuards            #-}-{-# LANGUAGE TypeOperators            #-}-{-# LANGUAGE ScopedTypeVariables      #-}-{-# LANGUAGE RecordWildCards          #-}  -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------ | Continuous mapping.-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)-  -------- | A chart is a homeomorphism from a connected, open subset /Q/ ⊂ /M/ of--- an /n/-manifold /M/ to either the open unit disk /Dⁿ/ ⊂ /V/ ≃ ℝ/ⁿ/, or--- the half-disk /Hⁿ/ = {/x/ ∊ /Dⁿ/: x₀≥0}. In e.g. the former case, 'chartInMap'--- is thus defined ∀ /v/ ∊ /V/ : |/v/| < 1, while 'chartOutMap p' will yield @Just x@--- with /x/ ∊ /Dⁿ/ provided /p/ is in /Q/, and @Nothing@ otherwise.--- Obviously, @fromJust . 'chartOutMap' . 'chartInMap'@ should be equivalent to @id@--- on /Dⁿ/, and @'chartInMap' . fromJust . 'chartOutMap'@ to @id@ on /Q/.-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  -- ^ A /M/ ⇆ /Dⁿ/ chart, for ordinary manifolds-               | RimChart         -- ^ A /M/ ⇆ /Hⁿ/ chart, for manifolds with a rim---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---- rimGuard :: EuclidSpace v => ChartKind -> v -> Maybe v--- rimGuard LandlockedChart v = Just v--- rimGuard RimChart v---  | isInUpperHemi v = Just v---  | otherwise       = Nothing--- --- chartEnv :: Manifold m => Chart m---                -> (TangentSpace m->TangentSpace m)---                -> m -> Maybe m--- chartEnv IdChart f x = Just $ f x--- chartEnv (Chart inMap outMap chKind) f x = do---     vGet <- outMap x---     let v = vGet $ x---     v' <- rimGuard chKind v---     Just $ inMap $ v'--- -  -- --type Atlas m = [Chart m]--class (MetricSpace(TangentSpace m), Metric(TangentSpace m) ~ ℝ) => Manifold m where-  type TangentSpace m :: *-  type TangentSpace m = m   -- For \"flat\", i.e. vector space manifolds.-  -  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   -- "Infinity" in floating-point-                | otherwise            = just $ log (expx + ε) - x--- exp x + ε = exp (x + δ) = exp x * exp δ--- δ = ln ( (exp x + ε)/exp 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 -> (a-b, \ε -> (ε, (ε/2,ε/2), ε))-  -  (*) = cntnFnValsCombine $ \a b -> (a*b, -                             \ε -> ( ε/b-                                   , (ε / (2 * sqrt(2*b^2+ε)), ε / (2 * sqrt(2*a^2+ε)))-                                   , ε/a ))-  --  |δa| < ε / 2·sqrt(2·b² + ε) ∧ |δb| < ε / 2·sqrt(2·a² + ε)-  --  ⇒  | (a+δa) · (b+δb) - a·b | = | a·δb + b·δa + δa·δb | -  --   ≤ | a·δb | + | b·δa | + | δa·δb |-  --   ≤ | a·ε/2·sqrt(2·a² + ε) | + | b·ε/2·sqrt(2·b² + ε) | + | ε² / 4·sqrt(2·b² + ε)·sqrt(2·a² + ε) |-  --   ≤ | a·ε/2·sqrt(2·a²) | + | b·ε/2·sqrt(2·b²) | + | ε² / 4·sqrt(ε)·sqrt(ε) |-  --   ≤ | ε/sqrt(8) | + | ε/sqrt(8) | + | ε / 4 |-  --   ≈ .96·ε < ε--  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¹))-  -- Readily derived from the worst-case of ε = 1 / (|x| – δ) – 1/|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)-------- instance Manifold S2 where---   type TangentSpace S2 = (Double, Double)---   localAtlas (S2 ϑ φ)---    | ϑ<pi-2     = [ Chart (\(x,y)---                              -> S2(2 * sqrt(x^2+y^2)) (atan2 y x) )---                           (\(S2 ϑ' φ')---                              -> let r=ϑ'/2---                                 in guard (r<1) >> Just (r * cos φ', r * sin φ') )---                           LandlockedChart ]---    | ϑ>2        = [ Chart (\(x,y)---                              -> S2(pi - 2*sqrt(x^2+y^2)) (atan2 y x) )---                           (\(S2 ϑ' φ')---                              -> let r=(pi-ϑ')/2---                                 in guard (r<1) >> Just (r * cos φ', r * sin φ') )---                           LandlockedChart ]---    | otherwise  = localAtlas(S2 0 φ) ++ localAtlas(S2 (2*pi) φ)--- --------(.:) :: (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+module Data.Manifold (module Data.Manifold.PseudoAffine, module Data.Manifold.Types) where +import Data.Manifold.PseudoAffine+import Data.Manifold.Types   

@@ -139,12 +139,12 @@                fstGriddingParams :: GriddingParameters m a              , sndGriddingParams :: GriddingParameters n a }   mkGridding (PairGriddingParameters p₁ p₂) n (Shade (c₁,c₂) e₁e₂)-          = gshmap ( uncurry fullShade . (                  (,c₂).(^.shadeCtr)+          = ( gshmap ( uncurry fullShade . (                  (,c₂).(^.shadeCtr)                                          &&& (`productMetric'`e₂).(^.shadeExpanse)) )-              <$> g₁s-         ++ gshmap ( uncurry fullShade . (                  (c₁,).(^.shadeCtr)+              <$> g₁s )+         ++ ( gshmap ( uncurry fullShade . (                  (c₁,).(^.shadeCtr)                                          &&& ( productMetric' e₁).(^.shadeExpanse)) )-              <$> g₂s+              <$> g₂s )    where g₁s = mkGridding p₁ n $ fullShade c₁ e₁          g₂s = mkGridding p₂ n $ fullShade c₂ e₂          (e₁,e₂) = factoriseMetric' e₁e₂ 

@@ -61,6 +61,8 @@             , HilbertSpace             , EuclidSpace             , LocallyScalable+            -- ** Local functions+            , LocalLinear, LocalAffine             -- * Misc             , palerp             ) where@@ -78,6 +80,7 @@ import Data.VectorSpace import Data.LinearMap import Data.LinearMap.HerMetric+import Data.LinearMap.Category import Data.MemoTrie (HasTrie(..)) import Data.AffineSpace import Data.Basis@@ -220,6 +223,9 @@                            , HasMetric (Needle x)                            , s ~ Scalar (Needle x) ) +type LocalLinear x y = Linear (Scalar (Needle x)) (Needle x) (Needle y)+type LocalAffine x y = (Needle y, LocalLinear x y)+ -- | Basically just an “updated” version of the 'VectorSpace' class. --   Every vector space is a manifold, this constraint makes it explicit. --   @@ -400,6 +406,17 @@   p.+~^n = p ^+^ linMapFromTensProd n instance (HasMetric a, FiniteDimensional b, Scalar a~Scalar b) => PseudoAffine (a:-*b) where   a.-~.b = pure . linMapAsTensProd $ a^-^b++instance (HasMetric a, FiniteDimensional b, Scalar a~s, Scalar b~s)+                          => Semimanifold (Linear s a b) where+  type Needle (Linear s a b) = Linear s a b+  fromInterior = id+  toInterior = pure+  translateP = Tagged (.+^)+  (.+~^) = (^+^)+instance (HasMetric a, FiniteDimensional b, Scalar a~s, Scalar b~s)+                          => PseudoAffine (Linear s a b) where+  a.-~.b = pure (a^-^b)  instance Semimanifold S⁰ where   type Needle S⁰ = ℝ⁰

@@ -13,6 +13,7 @@ {-# LANGUAGE DeriveGeneric              #-} {-# LANGUAGE DeriveFunctor              #-} {-# LANGUAGE DeriveFoldable             #-}+{-# LANGUAGE DeriveTraversable          #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE TypeFamilies               #-} {-# LANGUAGE FunctionalDependencies     #-}@@ -34,9 +35,15 @@  module Data.Manifold.TreeCover (        -- * Shades -         Shade(..), Shade'(..)-       -- ** Lenses and constructors-       , shadeCtr, shadeExpanse, shadeNarrowness, fullShade, fullShade', pointsShades+         Shade(..), Shade'(..), IsShade+       -- ** Lenses+       , shadeCtr, shadeExpanse, shadeNarrowness+       -- ** Construction+       , fullShade, fullShade', pointsShades+       -- ** Evaluation+       , occlusion+       -- ** Misc+       , factoriseShade, intersectShade's        -- * Shade trees        , ShadeTree(..), fromLeafPoints        -- * Simple view helpers@@ -45,6 +52,9 @@        , SimpleTree, Trees, NonEmptyTree, GenericTree(..)        -- * Misc        , sShSaw, chainsaw, HasFlatView(..), shadesMerge, smoothInterpolate+       , twigsWithEnvirons, completeTopShading, flexTwigsShading+       , WithAny(..), Shaded, stiAsIntervalMapping, spanShading+       , DifferentialEqn, filterDEqnSolution_static        -- ** Triangulation-builders        , TriangBuild, doTriangBuild, singleFullSimplex, autoglueTriangulation        , AutoTriang, elementaryTriang, breakdownAutoTriang@@ -63,6 +73,7 @@ import Control.DeepSeq  import Data.VectorSpace+import Data.AffineSpace import Data.LinearMap import Data.LinearMap.HerMetric import Data.LinearMap.Category@@ -89,6 +100,7 @@ import Data.Functor.Identity import Control.Monad.Trans.State import Control.Monad.Trans.Writer+import Control.Monad.Trans.Maybe import Control.Monad.Trans.Class import qualified Data.Foldable       as Hask import Data.Foldable (all, elem, toList, sum, foldr1)@@ -98,10 +110,11 @@ import qualified Numeric.LinearAlgebra.HMatrix as HMat  import Control.Category.Constrained.Prelude hiding-     ((^), all, elem, sum, forM, Foldable(..), foldr1, Traversable)+     ((^), all, elem, sum, forM, Foldable(..), foldr1, Traversable, traverse) import Control.Arrow.Constrained import Control.Monad.Constrained hiding (forM) import Data.Foldable.Constrained+import Data.Traversable.Constrained (traverse)  import GHC.Generics (Generic) @@ -122,12 +135,15 @@ --   there is 'Region', whose implementation is vastly more complex. data Shade x = Shade { _shadeCtr :: !(Interior x)                      , _shadeExpanse :: !(Metric' x) }+deriving instance (Show x, Show (Needle x), WithField ℝ Manifold x) => Show (Shade x)  -- | A &#x201c;co-shade&#x201d; can describe ellipsoid regions as well, but unlike --   'Shade' it can be unlimited / infinitely wide in some directions. --   It does OTOH need to have nonzero thickness, which 'Shade' needs not. data Shade' x = Shade' { _shade'Ctr :: !(Interior x)                        , _shade'Narrowness :: !(Metric x) }+deriving instance (Show x, Show (DualSpace (Needle x)), WithField ℝ Manifold x)+             => Show (Shade' x)  class IsShade shade where --  type (*) shade :: *->*@@ -135,15 +151,39 @@   shadeCtr :: Functor f (->) (->) => (Interior x->f (Interior x)) -> shade x -> f (shade x) --  -- | Convert between 'Shade' and 'Shade' (which must be neither singular nor infinite). --  unsafeDualShade :: WithField ℝ Manifold x => shade x -> shade* x+  -- | Check the statistical likelihood-density of a point being within a shade.+  --   This is taken as a normal distribution.+  occlusion :: ( Manifold x, s ~ (Scalar (Needle x)), RealDimension s )+                => shade x -> x -> s+  factoriseShade :: ( Manifold x, RealDimension (Scalar (Needle x))+                    , Manifold y, RealDimension (Scalar (Needle y)) )+                => shade (x,y) -> (shade x, shade y)  instance IsShade Shade where   shadeCtr f (Shade c e) = fmap (`Shade`e) $ f c+  occlusion (Shade p₀ δ) = occ+   where occ p = case p .-~. p₀ of+           Option(Just vd) | mSq <- metricSq δinv vd+                           , mSq == mSq  -- avoid NaN+                           -> exp (negate mSq)+           _               -> zeroV+         δinv = recipMetric δ+  factoriseShade (Shade (x₀,y₀) δxy) = (Shade x₀ δx, Shade y₀ δy)+   where (δx,δy) = factoriseMetric' δxy  shadeExpanse :: Functor f (->) (->) => (Metric' x -> f (Metric' x)) -> Shade x -> f (Shade x) shadeExpanse f (Shade c e) = fmap (Shade c) $ f e  instance IsShade Shade' where   shadeCtr f (Shade' c e) = fmap (`Shade'`e) $ f c+  occlusion (Shade' p₀ δinv) = occ+   where occ p = case p .-~. p₀ of+           Option(Just vd) | mSq <- metricSq δinv vd+                           , mSq == mSq  -- avoid NaN+                           -> exp (negate mSq)+           _               -> zeroV+  factoriseShade (Shade' (x₀,y₀) δxy) = (Shade' x₀ δx, Shade' y₀ δy)+   where (δx,δy) = factoriseMetric δxy  shadeNarrowness :: Functor f (->) (->) => (Metric x -> f (Metric x)) -> Shade' x -> f (Shade' x) shadeNarrowness f (Shade' c e) = fmap (Shade' c) $ f e@@ -186,6 +226,8 @@ pointsShades :: WithField ℝ Manifold x => [x] -> [Shade x] pointsShades = map snd . pointsShades' zeroV +pointsShade's :: WithField ℝ Manifold x => [x] -> [Shade' x]+pointsShade's = map (\(Shade c e) -> Shade' c $ recipMetric e) . pointsShades  pseudoECM :: WithField ℝ Manifold x => NonEmpty x -> (x, ([x],[x])) pseudoECM (p₀ NE.:| psr) = foldl' ( \(acc, (rb,nr)) (i,p)@@ -231,30 +273,34 @@                   = Just $ let cc = c₂ .+~^ v ^/ 2                                Option (Just cv₁) = c₁.-~.cc                                Option (Just cv₂) = c₂.-~.cc-                           in Shade cc . sumV $ [e₁, e₂] ++ projector'<$>[cv₁, cv₂] +                           in Shade cc . sumV $ [e₁, e₂] ++ (projector'<$>[cv₁, cv₂])            | otherwise  = Nothing shadesMerge _ shs = shs -minusLogOcclusion :: ( Manifold x, s ~ (Scalar (Needle x)), RealDimension s )-                => Shade x -> x -> s-minusLogOcclusion (Shade p₀ δ) = occ+-- | Evaluate the shade as a quadratic form; essentially+-- @+-- minusLogOcclusion sh x = x <.>^ (sh^.shadeExpanse $ x - sh^.shadeCtr)+-- @+-- where 'shadeExpanse' gives a metric (matrix) that characterises the+-- width of the shade.+minusLogOcclusion' :: ( Manifold x, s ~ (Scalar (Needle x)), RealDimension s )+              => Shade' x -> x -> s+minusLogOcclusion' (Shade' p₀ δinv) = occ  where occ p = case p .-~. p₀ of          Option(Just vd) | mSq <- metricSq δinv vd                          , mSq == mSq  -- avoid NaN                          -> mSq          _               -> 1/0-       δinv = recipMetric δ-  --- | Check the statistical likelyhood of a point being within a shade.-occlusion :: ( Manifold x, s ~ (Scalar (Needle x)), RealDimension s )-                => Shade x -> x -> s-occlusion (Shade p₀ δ) = occ+minusLogOcclusion :: ( Manifold x, s ~ (Scalar (Needle x)), RealDimension s )+              => Shade x -> x -> s+minusLogOcclusion (Shade p₀ δ) = occ  where occ p = case p .-~. p₀ of          Option(Just vd) | mSq <- metricSq δinv vd                          , mSq == mSq  -- avoid NaN-                         -> exp (negate mSq)-         _               -> zeroV+                         -> mSq+         _               -> 1/0        δinv = recipMetric δ+     @@ -306,8 +352,36 @@   DBranches b1 <> DBranches b2 = DBranches $ NE.zipWith (\(DBranch d1 c1) (DBranch _ c2)                                                               -> DBranch d1 $ c1<>c2 ) b1 b2   +directionChoices :: WithField ℝ Manifold x+               => [DBranch x]+                 -> [ ( (Needle' x, ShadeTree x)+                      ,[(Needle' x, ShadeTree x)] ) ]+directionChoices [] = []+directionChoices (DBranch ѧ (Hourglass t b) : hs)+       =  ( (ѧ,t), (v,b) : map fst uds)+          : ((v,b), (ѧ,t) : map fst uds)+          : map (second $ ((ѧ,t):) . ((v,b):)) uds+ where v = negateV ѧ+       uds = directionChoices hs +traverseDirectionChoices :: (WithField ℝ Manifold x, Hask.Applicative f)+               => (     (Needle' x, ShadeTree x)+                    -> [(Needle' x, ShadeTree x)]+                    -> f (ShadeTree x) )+                 -> [DBranch x]+                 -> f [DBranch x]+traverseDirectionChoices f dbs = td [] (dbs >>=+                            \(DBranch ѧ (Hourglass τ β))+                              -> [(ѧ,τ), (negateV ѧ,β)])+ where td pds ((ѧ,t):(v,b):vds)+         = liftA3 (\t' b' -> (DBranch ѧ (Hourglass t' b') :))+             (f (ѧ,t) $ pds++(v,b):uds)+             (f (v,b) $ pds++(ѧ,t):uds)+             $ td ((ѧ,t):(v,b):pds) vds+        where uds = pds ++ vds+       td _ _ = pure [] + instance (NFData x, NFData (Needle' x)) => NFData (ShadeTree x) where   rnf (PlainLeaves xs) = rnf xs   rnf (DisjointBranches n bs) = n `seq` rnf (NE.toList bs)@@ -345,22 +419,7 @@  -- | Build a quite nicely balanced tree from a cloud of points, on any real manifold. -- ---   Example:--- --- @--- > :m +Graphics.Dynamic.Plot.R2 Data.Manifold.TreeCover Data.VectorSpace Data.AffineSpace --- > import Diagrams.Prelude ((^&), p2, r2, P2, circle, fc, (&), moveTo, opacity)--- --- >   -- Generate sort-of&#x2013;random cloud of lots of points--- > let testPts0 = p2 \<$\> [(0,0), (0,1), (1,1), (1,2), (2,2)] :: [P2 Double]--- > let testPts1 = [p .+^ v^/3 | p\<-testPts0, v \<- r2\<$\>[(0,0), (-1,1), (1,2)]]--- > let testPts2 = [p .+^ v^/4 | p\<-testPts1, v \<- r2\<$\>[(0,0), (-1,1), (1,2)]]--- > let testPts3 = [p .+^ v^/5 | p\<-testPts2, v \<- r2\<$\>[(0,0), (-2,1), (1,2)]]--- > let testPts4 = [p .+^ v^/7 | p\<-testPts3, v \<- r2\<$\>[(0,1), (-1,1), (1,2)]]--- --- > plotWindow [ plot [ shapePlot $ circle 0.06 & moveTo p & opacity 0.3 | p <- testPts4 ]--- >            , plot . onlyNodes $ 'fromLeafPoints' testPts4 ]--- @+--   Example: https://nbviewer.jupyter.org/github/leftaroundabout/manifolds/blob/master/test/generate-ShadeTrees.ipynb#pseudorandomCloudTree --  -- <<images/examples/simple-2d-ShadeTree.png>> fromLeafPoints :: ∀ x. WithField ℝ Manifold x => [x] -> ShadeTree x@@ -375,7 +434,7 @@        fg_sShIdPart (Shade c expa) xs         | b:bs <- [DBranch (v, zeroV) mempty                     | v <- eigenCoSpan-                           (transformMetric' (linear fst) expa :: Metric' x) ]+                           (transformMetric' fst expa :: Metric' x) ]                       = sShIdPartition' c xs $ b:|bs  fromLeafPoints' :: ∀ x. WithField ℝ Manifold x =>@@ -461,78 +520,282 @@ sortByKey = map snd . sortBy (comparing fst)  +trunks :: ∀ x. WithField ℝ Manifold x => ShadeTree x -> [Shade x]+trunks (PlainLeaves lvs) = pointsShades lvs+trunks (DisjointBranches _ brs) = Hask.foldMap trunks brs+trunks (OverlappingBranches _ sh _) = [sh]  +nLeaves :: ShadeTree x -> Int+nLeaves (PlainLeaves lvs) = length lvs+nLeaves (DisjointBranches n _) = n+nLeaves (OverlappingBranches n _ _) = n++overlappingBranches :: Shade x -> NonEmpty (DBranch x) -> ShadeTree x+overlappingBranches shx brs = OverlappingBranches n shx brs+ where n = sum $ fmap (sum . fmap nLeaves) brs++unsafeFmapLeaves :: (x -> x) -> ShadeTree x -> ShadeTree x+unsafeFmapLeaves f (PlainLeaves lvs) = PlainLeaves $ fmap f lvs+unsafeFmapLeaves f (DisjointBranches n brs)+                  = DisjointBranches n $ unsafeFmapLeaves f <$> brs+unsafeFmapLeaves f (OverlappingBranches n sh brs)+                  = OverlappingBranches n sh $ fmap (unsafeFmapLeaves f) <$> brs++unsafeFmapTree :: (NonEmpty x -> NonEmpty y)+               -> (Needle' x -> Needle' y)+               -> (Shade x -> Shade y)+               -> ShadeTree x -> ShadeTree y+unsafeFmapTree _ _ _ (PlainLeaves []) = PlainLeaves []+unsafeFmapTree f _ _ (PlainLeaves lvs) = PlainLeaves . toList . f $ NE.fromList lvs+unsafeFmapTree f fn fs (DisjointBranches n brs)+    = let brs' = unsafeFmapTree f fn fs <$> brs+      in DisjointBranches (sum $ nLeaves<$>brs') brs'+unsafeFmapTree f fn fs (OverlappingBranches n sh brs)+    = let brs' = fmap (\(DBranch dir br)+                        -> DBranch (fn dir) (unsafeFmapTree f fn fs<$>br)+                      ) brs+      in overlappingBranches (fs sh) brs'++ intersectShade's :: ∀ y . WithField ℝ Manifold y => [Shade' y] -> Option (Shade' y) intersectShade's [] = error "Global `Shade'` not implemented, so can't do intersection of zero co-shades." intersectShade's (sh:shs) = Hask.foldrM inter2 sh shs  where inter2 :: Shade' y -> Shade' y -> Option (Shade' y)        inter2 (Shade' c e) (Shade' ζ η)-           | μc > 1 && μζ > 1  = empty-           | otherwise         = return $ Shade' (c.+~^w) (e^+^η)-        where Option (Just c2ζ) = ζ.-~.c-              Option (Just ζ2c) = c.-~.ζ-              ζNearest, cNearest :: y-              ζNearest = c .+~^ metriNormalise e c2ζ-              cNearest = ζ .+~^ metriNormalise η ζ2c-              Option (Just rζ) = ζNearest.-~.ζ-              Option (Just rc) = cNearest.-~.c-              μc = metric e rc-              μζ = metric η rζ-              w = c2ζ ^* (μζ/(μc + μζ))-              -- = (c^*μc + ζ^*μζ)/(μc + μζ) − c-              -- = (c^*μc + ζ^*μζ − c^*(μc+μζ))^/(μc + μζ)-              -- = (ζ^*μζ − c^*μζ)^/(μc + μζ)-              -- = (ζ−c)^*μζ/(μc + μζ)+           | μe < 1 && μη < 1  = return $ Shade' iCtr iExpa+           | otherwise         = empty+        where [c', ζ'] = [ ctr.+~^linearCombo+                                     [ (v, 1 / (1 + metricSq oExpa w))+                                     | v <- (*^) <$> [-1,1] <*> span+                                     , let p = ctr .+~^ v  :: y+                                           Option (Just w) = p.-~.oCtr+                                     ]+                         | ctr                  <- [c,     ζ    ]+                         | span <- eigenCoSpan'<$> [e,     η    ]+                         | (oCtr,oExpa)         <- [(ζ,η), (c,e)]+                         ]+              Option (Just c'2ζ') = ζ'.-~.c'+              Option (Just c2ζ') = ζ'.-~.c+              Option (Just ζ2c') = c'.-~.ζ+              μc = metricSq e c2ζ'+              μζ = metricSq η ζ2c'+              iCtr = c' .+~^ c'2ζ' ^* (μζ/(μc + μζ)) -- weighted mean between c' and ζ'.+              Option (Just rc) = c.-~.iCtr+              Option (Just rζ) = ζ.-~.iCtr+              rcⰰ = toDualWith e rc+              rζⰰ = toDualWith η rζ+              μe = rcⰰ<.>^rc+              μη = rζⰰ<.>^rζ+              iExpa = (e^+^η)^/2 ^+^ projector rcⰰ^/(1-μe) ^+^ projector rζⰰ^/(1-μη)    -type DifferentialEqn x y = RWDiffable ℝ (x,y) (Needle x :-* Needle y)+type DifferentialEqn x y = Shade' (x,y) -> Shade' (LocalLinear x y)   filterDEqnSolution_loc :: ∀ x y . (WithField ℝ Manifold x, WithField ℝ Manifold y)-           => DifferentialEqn x y -> (Shade' (x,y), [Shade' (x,y)]) -> [Shade' (x,y)]-filterDEqnSolution_loc (RWDiffable f) (Shade' (x,y) expa, neighbours) = case f (x,y) of-          (_, Option Nothing) -> []-          (r, Option (Just (Differentiable fl)))-                | (fc, fc', δ) <- fl (x,y)-                   -> let flatMet :: HerMetric (Needle (x,y))-                          flatMet = recipMetric -- this won't work, metric is singular.-                               . transformMetric' (linear $ id &&& lapply fc) -                               $ recipMetric' expax-                          -- fcs = lapply fc' <$> xSpan-                          -- flinRange = δ $ projectors fcs-                          marginδs :: [(Needle x, (Needle y, Metric y))]-                          marginδs = [ (δxm, (δym, expany))-                                     | Shade' (xn, yn) expan <- neighbours-                                     , let (Option (Just δx)) = x.-~.xn-                                           (expanx, expany) = factoriseMetric expan-                                           (Option (Just yc'n))-                                                  = covariance $ recipMetric' expan-                                           xntoMarg = metriNormalise expanx δx-                                           (Option (Just δxm))-                                              = (xn .+~^ xntoMarg :: x) .-~. x-                                           (Option (Just δym))-                                              = (yn .+~^ lapply yc'n xntoMarg :: y-                                                  ) .-~. y-                                     ]-                          ycQuad :: y-                          (Option (Just (Shade' ycQuad _))) = intersectShade's-                                     [ Shade' ycn expany-                                     | (δxm,(δym,expany)) <- marginδs-                                     , let fca :: Needle x:-*Needle y-                                           fca = fc .+~^ lapply fc' ((δxm,δym)^/2)-                                           ycn = y .+~^ (δym ^-^ lapply fca δxm)-                                     ]-                                     :: Option (Shade' y)-                      in [Shade' (x,ycQuad) flatMet]- where (expax, expay) = factoriseMetric expa+           => DifferentialEqn x y -> (Shade' (x,y), [Shade' (x,y)])+                   -> Option (Shade' y, LocalLinear x y)+filterDEqnSolution_loc f (shxy@(Shade' (x,y) expa), neighbours) = (,j₀) <$> yc+ where jShade@(Shade' j₀ jExpa) = f shxy+       marginδs :: [(Needle x, (Needle y, Metric y))]+       marginδs = [ (δxm, (δym, expany))+                  | Shade' (xn, yn) expan <- neighbours+                  , let (Option (Just δx)) = x.-~.xn+                        (expanx, expany) = factoriseMetric expan+                        (Option (Just yc'n))+                               = covariance $ recipMetric' expan+                        xntoMarg = metriNormalise expanx δx+                        (Option (Just δxm))+                           = (xn .+~^ xntoMarg :: x) .-~. x+                        (Option (Just δym))+                           = (yn .+~^ (yc'n $ xntoMarg) :: y+                               ) .-~. y+                  ]+       back2Centre :: (Needle x, (Needle y, Metric y)) -> Shade' y+       back2Centre (δx, (δym, expany))+            = Shade' (y.+~^δyb) . recipMetric+                $ recipMetric' expany+                  ^+^ recipMetric' (applyLinMapMetric jExpa δx')+        where δyb = δym ^-^ (j₀ $ δx)+              δx' = toDualWith expax δx+       yc :: Option (Shade' y)+       yc = intersectShade's $ back2Centre <$> marginδs+       (expax, expay) = factoriseMetric expa        xSpan = eigenCoSpan' expax  +-- Formerly, this was the signature of what has now become 'traverseTwigsWithEnvirons'.+-- The simple list-yielding version (see rev. b4a427d59ec82889bab2fde39225b14a57b694df+-- may well be more efficient than this version via a traversal.+twigsWithEnvirons :: ∀ x. WithField ℝ Manifold x+    => ShadeTree x -> [(ShadeTree x, [ShadeTree x])]+twigsWithEnvirons = execWriter . traverseTwigsWithEnvirons (writer . (fst&&&pure))++data OuterMaybeT f a = OuterNothing | OuterJust (f a) deriving (Hask.Functor)+instance (Hask.Applicative f) => Hask.Applicative (OuterMaybeT f) where+  pure = OuterJust . pure+  OuterJust fs <*> OuterJust xs = OuterJust $ fs <*> xs+  _ <*> _ = OuterNothing++traverseTwigsWithEnvirons :: ∀ x f .+            (WithField ℝ Manifold x, Hask.Applicative f)+    => ((ShadeTree x, [ShadeTree x]) -> f (ShadeTree x))+         -> ShadeTree x -> f (ShadeTree x)+traverseTwigsWithEnvirons f = fst . go []+ where go :: [ShadeTree x] -> ShadeTree x -> (f (ShadeTree x), Bool)+       go _ (DisjointBranches nlvs djbs) = ( fmap (DisjointBranches nlvs)+                                               $ Hask.traverse (fst . go []) djbs+                                           , False )+       go envi ct@(OverlappingBranches nlvs rob@(Shade robc _) brs)+                = ( case descentResult of+                     OuterNothing -> f+                         $ purgeRemotes (ct, Hask.foldMap (twigProximæ robc) envi)+                     OuterJust dR -> fmap (OverlappingBranches nlvs rob . NE.fromList) dR+                  , False )+        where descentResult = traverseDirectionChoices tdc $ NE.toList brs+              tdc (vy, ty) alts = case go envi'' ty of+                                   (_, True) -> OuterNothing+                                   (down, _) -> OuterJust down+               where envi'' = filter (trunks >>> \(Shade ce _:_)+                                         -> let Option (Just δyenv) = ce.-~.robc+                                                qq = vy<.>^δyenv+                                            in qq > -1 && qq < 5+                                       ) envi'+                              ++ map snd alts+              envi' = approach =<< envi+              approach apt@(OverlappingBranches _ (Shade envc _) _)+                  = twigsaveTrim hither apt+               where Option (Just δxenv) = robc .-~. envc+                     hither (DBranch bdir (Hourglass bdc₁ bdc₂))+                       | bdir<.>^δxenv > 0  = [bdc₁]+                       | otherwise          = [bdc₂]+              approach q = [q]+       go envi plvs@(PlainLeaves _) = (f $ purgeRemotes (plvs, envi), True)+       +       twigProximæ :: x -> ShadeTree x -> [ShadeTree x]+       twigProximæ x₀ (DisjointBranches _ djbs) = Hask.foldMap (twigProximæ x₀) djbs+       twigProximæ x₀ ct@(OverlappingBranches _ (Shade xb qb) brs)+                   = twigsaveTrim hither ct+        where Option (Just δxb) = x₀ .-~. xb+              hither (DBranch bdir (Hourglass bdc₁ bdc₂))+                 | bdir<.>^δxb > 0  = twigProximæ x₀ bdc₁+                 | otherwise        = twigProximæ x₀ bdc₂+       twigProximæ _ plainLeaves = [plainLeaves]+       +       twigsaveTrim :: (DBranch x -> [ShadeTree x])+                       -> ShadeTree x -> [ShadeTree x]+       twigsaveTrim f ct@(OverlappingBranches _ _ dbs)+                 = case Hask.mapM (f >>> noLeaf) dbs of+                      Just pqe -> Hask.fold pqe+                      _        -> [ct]+        where noLeaf [PlainLeaves _] = empty+              noLeaf bqs = pure bqs+       +       purgeRemotes :: (ShadeTree x, [ShadeTree x]) -> (ShadeTree x, [ShadeTree x])+       purgeRemotes (ctm@(OverlappingBranches _ sm@(Shade xm _) _), candidates)+                                       = (ctm, filter unobscured closeby)+        where closeby = filter proximate candidates+              proximate (OverlappingBranches _ sh@(Shade xh _) _)+                    = minusLogOcclusion sh xm * minusLogOcclusion sm xh+                       < 1024  -- = (2⋅4²)².  The four-radius occlusion occurs+                               -- if two 𝑟-sized shades have just enough space+                               -- to fit another 𝑟-shade between them; then+                               -- we don't consider the shades neighbours+                               -- anymore. A factor √2 for the discrepancy+                               -- between standard deviation and max distance.+              proximate _ = True+              unobscured ht@(OverlappingBranches _ (Shade xh _) _)+                     = all (don'tObscure (xh, onlyLeaves ht)) closeby+              don'tObscure (xh,lvsh) (OverlappingBranches _ sb@(Shade xb eb) _)+                          = vmc⋅vhc >= 0 || vm⋅vh >= 0+               where Option (Just vm) = pbm .-~. xb+                     Option (Just vh) = pbh .-~. xb+                     Option (Just vmc) = xm .-~. xb+                     Option (Just vhc) = xh .-~. xb+                     [pbm, pbh] = [ maximumBy (comparing $ \l ->+                                               let Option (Just w) = l.-~.xb+                                               in v⋅w ) lvs+                                  | lvs <- [lvsm, lvsh]+                                  | v <- [vhc, vmc] ]+                     (⋅) :: Needle x -> Needle x -> ℝ+                     v⋅w = toDualWith mb v <.>^ w+                     mb = recipMetric eb+              don'tObscure _ _ = True+              lvsm = onlyLeaves ctm+       purgeRemotes xyz = xyz     +    +completeTopShading :: (WithField ℝ Manifold x, WithField ℝ Manifold y)+                   => x`Shaded`y -> [Shade' (x,y)]+completeTopShading (PlainLeaves plvs)+                     = pointsShade's $ (_topological &&& _untopological) <$> plvs+completeTopShading (DisjointBranches _ bqs)+                     = take 1 . completeTopShading =<< NE.toList bqs+completeTopShading t = pointsShade's . map (_topological &&& _untopological) $ onlyLeaves t +flexTopShading :: ∀ x y f . ( WithField ℝ Manifold x, WithField ℝ Manifold y+                            , Applicative f (->) (->) )+                  => (Shade' (x,y) -> f (x, (Shade' y, LocalLinear x y)))+                      -> x`Shaded`y -> f (x`Shaded`y)+flexTopShading f tr = seq (assert_onlyToplevDisjoint tr)+                    $ recst (completeTopShading tr) tr+ where recst qsh@(_:_) (DisjointBranches n bqs)+          = undefined -- DisjointBranches n $ NE.zipWith (recst . (:[])) (NE.fromList qsh) bqs+       recst [sha@(Shade' (_,yc₀) expa₀)] t = fmap fts $ f sha+        where expa'₀ = recipMetric' expa₀+              j₀ :: LocalLinear x y+              Option (Just j₀) = covariance expa'₀+              (_,expay₀) = factoriseMetric expa₀+              fts (xc, (Shade' yc expay, jtg)) = unsafeFmapLeaves applδj t+               where Option (Just δyc) = yc.-~.yc₀+                     tfm = imitateMetricSpanChange expay₀ (recipMetric' expay)+                     applδj (WithAny y x)+                           = WithAny (yc₀ .+~^ ((tfm$δy) ^+^ (jtg$δx) ^+^ δyc)) x+                      where Option (Just δx) = x.-~.xc+                            Option (Just δy) = y.-~.(yc₀.+~^(j₀$δx))+       +       assert_onlyToplevDisjoint, assert_connected :: x`Shaded`y -> ()+       assert_onlyToplevDisjoint (DisjointBranches _ dp) = rnf (assert_connected<$>dp)+       assert_onlyToplevDisjoint t = assert_connected t+       assert_connected (OverlappingBranches _ _ dp)+           = rnf (Hask.foldMap assert_connected<$>dp)+       assert_connected (PlainLeaves _) = ()++flexTwigsShading :: ∀ x y f . ( WithField ℝ Manifold x, WithField ℝ Manifold y+                              , Hask.Applicative f )+                  => (Shade' (x,y) -> f (x, (Shade' y, LocalLinear x y)))+                      -> x`Shaded`y -> f (x`Shaded`y)+flexTwigsShading f = traverseTwigsWithEnvirons locFlex+ where locFlex :: ∀ μ . (x`Shaded`y, μ) -> f (x`Shaded`y)+       locFlex (lsh, _) = flexTopShading f lsh++filterDEqnSolution_static :: ∀ x y . (WithField ℝ Manifold x, WithField ℝ Manifold y)+           => DifferentialEqn x y+               -> x`Shaded`y -> Option (x`Shaded`y)+filterDEqnSolution_static deq tr = traverseTwigsWithEnvirons locSoltn tr+ where locSoltn :: (x`Shaded`y, [x`Shaded`y]) -> Option (x`Shaded`y)+       locSoltn (local, environs) = do+            let enviShades = completeTopShading =<< environs+            flexed <- flexTopShading+                           (\oSh@(Shade' (ox,_) _) -> +                              (ox,) <$> filterDEqnSolution_loc deq (oSh, enviShades)+                           ) local+            top'@(Shade' (top'x,_) top'exp)+                     <- intersectShade's $ completeTopShading =<< [local, flexed]+            let (_, top'ySh) = factoriseShade top'+            j' <- covariance $ recipMetric' top'exp+            flexTopShading (const $ pure (top'x, (top'ySh, j'))) flexed+                +++++++ -- simplexFaces :: forall n x . Simplex (S n) x -> Triangulation n x -- simplexFaces (Simplex p (ZeroSimplex q))    = TriangVertices $ Arr.fromList [p, q] -- simplexFaces splx = carpent splx $ TriangVertices ps@@ -633,7 +896,7 @@ simplexPlane :: forall n x . (KnownNat n, WithField ℝ Manifold x)         => Metric x -> Simplex n x -> Embedding (Linear ℝ) (FreeVect n ℝ) (Needle x) simplexPlane m s = embedding- where bc = barycenter s+ where bc = simplexBarycenter s        spread = init . map ((.-~.bc) >>> \(Option (Just v)) -> v) $ splxVertices s        embedding = case spanHilbertSubspace m spread of                      (Option (Just e)) -> e@@ -641,10 +904,14 @@                                 \ simplex (which cannot span sufficient basis vectors)."  +leavesBarycenter :: WithField ℝ Manifold x => NonEmpty x -> x+leavesBarycenter (x :| xs) = x .+~^ sumV [x'–x | x'<-xs] ^/ (n+1)+ where n = fromIntegral $ length xs+       x' – x = case x'.-~.x of {Option(Just v)->v}  -- simplexShade :: forall x n . (KnownNat n, WithField ℝ Manifold x)-barycenter :: forall x n . (KnownNat n, WithField ℝ Manifold x) => Simplex n x -> x-barycenter = bc +simplexBarycenter :: forall x n . (KnownNat n, WithField ℝ Manifold x) => Simplex n x -> x+simplexBarycenter = bc   where bc (ZS x) = x        bc (x :<| xs') = x .+~^ sumV [x'–x | x'<-splxVertices xs'] ^/ (n+1)        @@ -654,7 +921,7 @@ toISimplex :: forall x n . (KnownNat n, WithField ℝ Manifold x)                  => Metric x -> Simplex n x -> ISimplex n x toISimplex m s = ISimplex $ fromEmbedProject fromBrc toBrc- where bc = barycenter s+ where bc = simplexBarycenter s        (Embedding emb (DenseLinear prj))                          = simplexPlane m s        (r₀:rs) = [ prj HMat.#> asPackedVector v@@ -850,8 +1117,8 @@ partitionsOfFstLength :: Int -> [a] -> [([a],[a])] partitionsOfFstLength 0 l = [([],l)] partitionsOfFstLength n [] = []-partitionsOfFstLength n (x:xs) = first (x:) <$> partitionsOfFstLength (n-1) xs-                              ++ second (x:) <$> partitionsOfFstLength n xs+partitionsOfFstLength n (x:xs) = ( first (x:) <$> partitionsOfFstLength (n-1) xs )+                              ++ ( second (x:) <$> partitionsOfFstLength n xs )  splxVertices :: Simplex n x -> [x] splxVertices (ZS x) = [x]@@ -898,7 +1165,9 @@ type NonEmptyTree = GenericTree NonEmpty []      newtype GenericTree c b x = GenericTree { treeBranches :: c (x,GenericTree b b x) }- deriving (Hask.Functor)+ deriving (Generic, Hask.Functor, Hask.Foldable, Hask.Traversable)+instance (NFData x, Hask.Foldable c, Hask.Foldable b) => NFData (GenericTree c b x) where+  rnf (GenericTree t) = rnf $ toList t instance (Hask.MonadPlus c) => Semigroup (GenericTree c b x) where   GenericTree b1 <> GenericTree b2 = GenericTree $ Hask.mplus b1 b2 instance (Hask.MonadPlus c) => Monoid (GenericTree c b x) where@@ -1028,8 +1297,10 @@ data x`WithAny`y       = WithAny { _untopological :: y                 , _topological :: !x  }- deriving (Hask.Functor)+ deriving (Hask.Functor, Show, Generic) +instance (NFData x, NFData y) => NFData (WithAny x y)+ instance (Semimanifold x) => Semimanifold (x`WithAny`y) where   type Needle (WithAny x y) = Needle x   type Interior (WithAny x y) = Interior x `WithAny` y@@ -1068,6 +1339,9 @@   WithAny y x >>= f = WithAny r $ x^+^q    where WithAny r q = f y +shadeWithAny :: y -> Shade x -> Shade (x`WithAny`y)+shadeWithAny y (Shade x xe) = Shade (WithAny y x) xe+ shadeWithoutAnything :: Shade (x`WithAny`y) -> Shade x shadeWithoutAnything (Shade (WithAny _ b) e) = Shade b e @@ -1104,6 +1378,13 @@                  $ linearCombo [(v, d/dens) | Cℝay d v <- NE.toList contribs]         where dens = sum (hParamCℝay <$> contribs) +stiAsIntervalMapping :: (x ~ ℝ, y ~ ℝ)+            => x`Shaded`y -> [(x, ((y, Diff y), Linear ℝ x y))]+stiAsIntervalMapping = twigsWithEnvirons >=> pure.fst >=> completeTopShading >=> pure.+             \(Shade' (xloc, yloc) shd)+                 -> ( xloc, ( (yloc, recip $ metric shd (0,1))+                            , case covariance (recipMetric' shd) of+                                {Option(Just j)->j} ) )  smoothInterpolate :: (WithField ℝ Manifold x, WithField ℝ LinearManifold y)              => NonEmpty (x,y) -> x -> y@@ -1117,6 +1398,21 @@        ltr = stiWithDensity $ fromLeafPoints l'  +spanShading :: ∀ x y . (WithField ℝ Manifold x, WithField ℝ Manifold y)+          => (Shade x -> Shade y) -> ShadeTree x -> x`Shaded`y+spanShading f = unsafeFmapTree addYs id addYSh+ where addYs :: NonEmpty x -> NonEmpty (x`WithAny`y)+       addYs l = foldr (NE.<|) (fmap ( WithAny ymid) l     )+                               (fmap (`WithAny`xmid) yexamp)+          where [xsh@(Shade xmid _)] = pointsShades $ toList l+                Shade ymid yexpa = f xsh+                yexamp = [ ymid .+~^ σ*^δy+                         | δy <- eigenSpan yexpa, σ <- [-1,1] ]+       addYSh :: Shade x -> Shade (x`WithAny`y)+       addYSh xsh = shadeWithAny (_shadeCtr $ f xsh) xsh+                      ++ coneTip :: (AdditiveGroup v) => Cℝay v coneTip = Cℝay 0 zeroV @@ -1128,6 +1424,9 @@ foci :: [a] -> [(a,[a])] foci [] = [] foci (x:xs) = (x,xs) : fmap (second (x:)) (foci xs)+       +fociNE :: NonEmpty a -> NonEmpty (a,[a])+fociNE (x:|xs) = (x,xs) :| fmap (second (x:)) (foci xs)          (.:) :: (c->d) -> (a->b->c) -> a->b->d 

@@ -50,10 +50,11 @@         , D¹(..), D²(..)         , ℝay         , CD¹(..), Cℝay(..)-        -- * Misc         -- * Cut-planes         , Cutplane(..)         , fathomCutDistance, sideOfCut+        -- * Linear mappings+        , Linear, denseLinear    ) where  @@ -75,6 +76,7 @@ import Data.Manifold.Cone import Data.LinearMap.HerMetric import Data.VectorSpace.FiniteDimensional+import Data.LinearMap.Category (Linear, denseLinear)  import qualified Prelude 

binary file changed (absent → 150317 bytes)

@@ -1,5 +1,5 @@ Name:                manifolds-Version:             0.1.6.3+Version:             0.2.0.1 Category:            Math Synopsis:            Coordinate-free hypersurfaces Description:         Manifolds, a generalisation of the notion of &#x201c;smooth curves&#x201d; or surfaces,@@ -48,7 +48,8 @@                      , void                      , tagged                      , deepseq-                     , constrained-categories >= 0.2 && < 0.3+                     , trivial-constraint >= 0.4+                     , constrained-categories >= 0.2.3 && < 0.3   other-extensions:  FlexibleInstances                      , TypeFamilies                      , FlexibleContexts