diff --git a/Data/Function/Affine.hs b/Data/Function/Affine.hs --- a/Data/Function/Affine.hs +++ b/Data/Function/Affine.hs @@ -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 >>>!, <<>>!) :: ( 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 + +(<< Affine s β γ -> Affine s α β -> Affine s α γ +(<<>>!) + +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 << 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 diff --git a/Data/Function/Differentiable.hs b/Data/Function/Differentiable.hs --- a/Data/Function/Differentiable.hs +++ b/Data/Function/Differentiable.hs @@ -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)) + + diff --git a/Data/Function/Differentiable/Data.hs b/Data/Function/Differentiable/Data.hs --- a/Data/Function/Differentiable/Data.hs +++ b/Data/Function/Differentiable/Data.hs @@ -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). -- --- ¹(The implementation does not deal with /ε/ and /δ/ 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.) +-- ¹(The implementation does not deal with /ε/ and /δ/ 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@. diff --git a/Data/LinearMap/Category.hs b/Data/LinearMap/Category.hs --- a/Data/LinearMap/Category.hs +++ b/Data/LinearMap/Category.hs @@ -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<< (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 diff --git a/Data/LinearMap/HerMetric.hs b/Data/LinearMap/HerMetric.hs --- a/Data/LinearMap/HerMetric.hs +++ b/Data/LinearMap/HerMetric.hs @@ -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 “magnitude” 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 “project to the metric's boundary”. 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 – 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 @@ -- “scaled length” 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 ℝ -> ℝ diff --git a/Data/Manifold.hs b/Data/Manifold.hs --- a/Data/Manifold.hs +++ b/Data/Manifold.hs @@ -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 ϑ φ) --- | ϑ 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 diff --git a/Data/Manifold/Griddable.hs b/Data/Manifold/Griddable.hs --- a/Data/Manifold/Griddable.hs +++ b/Data/Manifold/Griddable.hs @@ -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₂ diff --git a/Data/Manifold/PseudoAffine.hs b/Data/Manifold/PseudoAffine.hs --- a/Data/Manifold/PseudoAffine.hs +++ b/Data/Manifold/PseudoAffine.hs @@ -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⁰ = ℝ⁰ diff --git a/Data/Manifold/TreeCover.hs b/Data/Manifold/TreeCover.hs --- a/Data/Manifold/TreeCover.hs +++ b/Data/Manifold/TreeCover.hs @@ -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 “co-shade” 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–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 -- -- <> 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 diff --git a/Data/Manifold/Types.hs b/Data/Manifold/Types.hs --- a/Data/Manifold/Types.hs +++ b/Data/Manifold/Types.hs @@ -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 diff --git a/images/examples/cartesiandisk-2d-ShadeTree.png b/images/examples/cartesiandisk-2d-ShadeTree.png new file mode 100644 Binary files /dev/null and b/images/examples/cartesiandisk-2d-ShadeTree.png differ diff --git a/manifolds.cabal b/manifolds.cabal --- a/manifolds.cabal +++ b/manifolds.cabal @@ -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 “smooth curves” 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