diff --git a/Data/Function/Affine.hs b/Data/Function/Affine.hs --- a/Data/Function/Affine.hs +++ b/Data/Function/Affine.hs @@ -40,8 +40,6 @@ import Data.Semigroup import Data.VectorSpace -import Data.LinearMap -import Data.LinearMap.HerMetric import Data.AffineSpace import Data.Tagged import Data.Manifold.Types.Primitive @@ -56,13 +54,14 @@ import Control.Monad.Constrained import Data.Foldable.Constrained +import Math.LinearMap.Category data Affine s d c where Subtract :: AffineManifold α => Affine s (α,α) (Needle α) AddTo :: Affine s (α, Needle α) α - ScaleWith :: (LinearManifold α, LinearManifold β) => (α:-*β) -> Affine s α β + ScaleWith :: (LinearManifold α, LinearManifold β) => (α+>β) -> Affine s α β ReAffine :: ReWellPointed (Affine s) α β -> Affine s α β reAffine :: ReWellPointed (Affine s) α β -> Affine s α β @@ -90,23 +89,29 @@ toOffsetSlope :: (MetricScalar s, WithField s LinearManifold d , WithField s AffineManifold c ) - => Affine s d c -> (c, Needle d :-* Needle c) + => Affine s d c -> (c, Needle d +> Needle c) toOffsetSlope f = toOffset'Slope f zeroV +type MetricScalar s = (Num''' s, LSpace (ZeroDim s)) + +linear :: (LSpace a, LSpace b, Scalar a ~ Scalar b) + => (a -> b) -> (a+>b) +linear = arr . LinearFunction + -- | 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) + => 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 (ScaleWith q) ref = (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) + (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) @@ -114,13 +119,13 @@ 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) + ((cf, sf), (cg, sg)) -> ((cf,cg), sf *** 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) + ((cf, sf), (cg, sg)) -> ((cf,cg), sf &&& sg) toOffset'Slope (Const c) ref = (c, zeroV) coOffsetForm :: ( MetricScalar s, WithField s AffineManifold d @@ -263,7 +268,7 @@ Swap .+^ Swap = Swap >>> ScaleWith (linear (^*2)) f .+^ Id = let (c,q) = toOffset'Slope f zeroV - in const c&&&ScaleWith (q^+^idL) >>>! AddTo + in const c&&&ScaleWith (q^+^id) >>>! 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 @@ -343,9 +348,9 @@ id = ReAffine id - ScaleWith ϕ . ScaleWith ψ = ScaleWith $ ϕ*.*ψ + ScaleWith ϕ . ScaleWith ψ = ScaleWith $ ϕ . ψ g . ScaleWith ψ = let (d, ϕ) = toOffsetSlope g - in postAdd' d $ ScaleWith (ϕ*.*ψ) + 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 @@ -392,7 +397,7 @@ linearAffine :: (MetricScalar s, WithField s LinearManifold α, WithField s LinearManifold β) - => (α:-*β) -> Affine s α β + => (α+>β) -> Affine s α β linearAffine = ScaleWith @@ -411,7 +416,7 @@ -instance (WithField s LinearManifold v, WithField s LinearManifold a) +instance (MetricScalar s, WithField s LinearManifold v, WithField s LinearManifold a) => AdditiveGroup (AffinFuncValue s a v) where zeroV = GenericAgent zeroV GenericAgent f ^+^ GenericAgent g = GenericAgent $ f ^+^ g diff --git a/Data/Function/Differentiable.hs b/Data/Function/Differentiable.hs --- a/Data/Function/Differentiable.hs +++ b/Data/Function/Differentiable.hs @@ -55,9 +55,7 @@ import Data.Embedding import Data.VectorSpace -import Data.LinearMap -import Data.LinearMap.Category -import Data.LinearMap.HerMetric +import Math.LinearMap.Category import Data.AffineSpace import Data.Function.Differentiable.Data import Data.Function.Affine @@ -78,7 +76,7 @@ -discretisePathIn :: WithField ℝ Manifold y +discretisePathIn :: (WithField ℝ Manifold y, SimpleSpace (Needle y)) => Int -- ^ Limit the number of steps taken in either direction. Note this will not cap the resolution but /length/ of the discretised path. -> ℝInterval -- ^ Parameter interval of interest. -> (RieMetric ℝ, RieMetric y) -- ^ Inaccuracy allowance /ε/. @@ -91,8 +89,8 @@ | signum (x₀-xlim) == signum dir = [(xlim, fxlim)] | otherwise = (x₀, fx₀) : traceFwd xlim (x₀+xstep) dir where (fx₀, jf, δx²) = f x₀ - εx = my fx₀ `extendMetric` lapply jf (metricAsLength $ mx x₀) - χ = metric (δx² εx) 1 + εx = my fx₀ `relaxNorm` [jf $ normalLength $ mx x₀] + χ = δx² εx |$| 1 xstep = dir * min (abs x₀+1) (recip χ) (fxlim, _, _) = f xlim xm = (xr + xl) / 2 @@ -109,10 +107,10 @@ = ([], [(-huge,huge)]) | otherwise = glueMid (go xc (-1)) (go xc 1) where go x₀ dir - | yq₀ <= abs (lapply jq₀ 1 * step₀) + | yq₀ <= abs ((jq₀$1) * step₀) = go (x₀ + step₀/2) dir | RealSubray PositiveHalfSphere xl' <- rangeHere - = let stepl' = dir/metric (δbf xl') 2 + = let stepl' = dir/(δbf xl'|$| 2) in if dir>0 then if definedHere then [(max (xl'+stepl') x₀, huge)] else [] @@ -120,7 +118,7 @@ then (xl'+stepl',x₀) : go (xl'-stepl') dir else go (xl'-stepl') dir | RealSubray NegativeHalfSphere xr' <- rangeHere - = let stepr' = dir/metric (δbf xr') 2 + = let stepr' = dir/(δbf xr'|$| 2) in if dir<0 then if definedHere then [(-huge, min (xr'-stepr') x₀)] else [] @@ -131,7 +129,7 @@ where (rangeHere, fq₀) = f x₀ (PreRegion (Differentiable r₀)) = genericisePreRegion rangeHere (yq₀, jq₀, δyq₀) = r₀ x₀ - step₀ = dir/metric (δbf x₀) 1 + step₀ = dir/(δbf x₀|$| 1) exit 0 _ xq | not definedHere = [] | xq < xc = [(xq,x₀)] @@ -150,10 +148,10 @@ xq₂ = xq₁ + stepp yq₁ = yq + f'x*stepp yq₂ = yq₁ + f'x*stepp - f'x = lapply jq 1 + f'x = jq $ 1 stepp | f'x*dir < 0 = -0.9 * abs dir' * yq/f'x | otherwise = dir' * as_devεδ δyq yq -- TODO: memoise in `exit` recursion - resoHere = metricSq $ δbf xq + resoHere = normSq $ δbf xq resoStep = dir/sqrt(resoHere 1) definedHere = case fq₀ of Option (Just _) -> True @@ -163,7 +161,7 @@ huge = exp $ fromIntegral nLim xc = 0 -discretisePathSegs :: WithField ℝ Manifold y +discretisePathSegs :: (WithField ℝ Manifold y, SimpleSpace (Needle y)) => Int -- ^ Maximum number of path segments and/or points per segment. -> ( RieMetric ℝ , RieMetric y ) -- ^ Inaccuracy allowance /δ/ for arguments @@ -195,11 +193,11 @@ | inRegion r x₀ -> return $ let (fx, j, δf) = fd x₀ epsprop ε - | ε>0 = case metric (δf $ metricFromLength ε) 1 of + | ε>0 = case (δf $ spanNorm [recip ε])|$| 1 of 0 -> empty δ' -> return $ recip δ' | otherwise = pure 0 - in ((fx, lapply j 1), epsprop) + in ((fx, j $ 1), epsprop) _ -> empty where -- This check shouldn't really be necessary, -- because the initial value lies by definition @@ -242,10 +240,10 @@ | y > b-resoHere = go (x + dir/χ) dir (a,y) | otherwise = go (x + safeStep stepOut₀) dir (a,b) where (y, j, δε) = fddd x - y' = lapply j 1 + y' = j $ 1 εx = my y - resoHere = metricAsLength εx - χ = metric (δε εx) 1 + resoHere = normalLength εx + χ = δε εx|$| 1 safeStep s₀ | as_devεδ δε (safetyMarg s₀) > abs s₀ = s₀ | otherwise = safeStep (s₀*0.5) @@ -267,29 +265,29 @@ unsafe_dev_ε_δ :: RealDimension a => String -> (a -> a) -> LinDevPropag a a unsafe_dev_ε_δ errHint f d - = let ε'² = metricSq d 1 + = let ε'² = normSq d 1 in if ε'²>0 then let δ = f . sqrt $ recip ε'² in if δ > 0 - then projector $ recip δ + then spanNorm [recip δ] else error $ "ε-δ propagator function for " ++errHint++", with ε=" ++show(sqrt $ recip ε'²) ++ " gives non-positive δ="++show δ++"." - else zeroV + else mempty dev_ε_δ :: RealDimension a => (a -> a) -> Metric a -> Option (Metric a) -dev_ε_δ f d = let ε'² = metricSq d 1 +dev_ε_δ f d = let ε'² = normSq d 1 in if ε'²>0 then let δ = f . sqrt $ recip ε'² in if δ > 0 - then pure . projector $ recip δ + then pure (spanNorm [recip δ]) else empty - else pure zeroV + else pure mempty as_devεδ :: RealDimension a => LinDevPropag a a -> a -> a as_devεδ ldp ε | ε>0 - , δ'² <- metricSq (ldp . projector $ recip ε) 1 + , δ'² <- normSq (ldp $ spanNorm [recip ε]) 1 , δ'² > 0 = sqrt $ recip δ'² | otherwise = 0 @@ -299,22 +297,20 @@ => Differentiable s d c -> Differentiable s d c genericiseDifferentiable (AffinDiffable _ af) = Differentiable $ \x -> let (y₀, ϕ) = toOffset'Slope af x - in (y₀, ϕ, const zeroV) + in (y₀, ϕ, const mempty) genericiseDifferentiable f = f -instance (MetricScalar s) => Category (Differentiable s) where +instance RealFrac' s => Category (Differentiable s) where type Object (Differentiable s) o = LocallyScalable s o - id = Differentiable $ \x -> (x, idL, const zeroV) + id = Differentiable $ \x -> (x, id, const mempty) Differentiable f . Differentiable g = Differentiable $ \x -> let (y, g', devg) = g x - jg = convertLinear $->$ g' (z, f', devf) = f y - jf = convertLinear $->$ f' - devfg δz = let δy = transformMetric jf δz + devfg δz = let δy = transformNorm f' δz εy = devf δz - in transformMetric jg εy ^+^ devg δy ^+^ devg εy - in (z, f'*.*g', devfg) + in transformNorm g' εy <> devg δy <> devg εy + in (z, f' . g', devfg) AffinDiffable ef f . AffinDiffable eg g = AffinDiffable (ef . eg) (f . g) f . g = genericiseDifferentiable f . genericiseDifferentiable g @@ -326,89 +322,80 @@ arr (Differentiable f) x = let (y,_,_) = f x in y arr (AffinDiffable _ f) x = f $ x -instance (MetricScalar s) => Cartesian (Differentiable s) where +instance (RealFrac' s) => Cartesian (Differentiable s) where type UnitObject (Differentiable s) = ZeroDim s - swap = Differentiable $ \(x,y) -> ((y,x), lSwap, const zeroV) - where lSwap = linear swap - attachUnit = Differentiable $ \x -> ((x, Origin), lAttachUnit, const zeroV) - where lAttachUnit = linear $ \x -> (x, Origin) - detachUnit = Differentiable $ \(x, Origin) -> (x, lDetachUnit, const zeroV) - where lDetachUnit = linear $ \(x, Origin) -> x - regroup = Differentiable $ \(x,(y,z)) -> (((x,y),z), lRegroup, const zeroV) - where lRegroup = linear regroup - regroup' = Differentiable $ \((x,y),z) -> ((x,(y,z)), lRegroup, const zeroV) - where lRegroup = linear regroup' + swap = Differentiable $ \(x,y) -> ((y,x), swap, const mempty) + attachUnit = Differentiable $ \x -> ((x, Origin), attachUnit, const mempty) + detachUnit = Differentiable $ \(x, Origin) -> (x, detachUnit, const mempty) + regroup = Differentiable $ \(x,(y,z)) -> (((x,y),z), regroup, const mempty) + regroup' = Differentiable $ \((x,y),z) -> ((x,(y,z)), regroup', const mempty) -instance (MetricScalar s) => Morphism (Differentiable s) where +instance (RealFrac' s) => Morphism (Differentiable s) where Differentiable f *** Differentiable g = Differentiable h - where h (x,y) = ((fx, gy), lPar, devfg) + where h (x,y) = ((fx, gy), f'***g', devfg) where (fx, f', devf) = f x (gy, g', devg) = g y - 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' + devfg δs = transformNorm fst δx + <> transformNorm snd δy + where δx = devf $ transformNorm (id&&&zeroV) δs + δy = devg $ transformNorm (zeroV&&&id) δs AffinDiffable IsDiffableEndo f *** AffinDiffable IsDiffableEndo g = AffinDiffable IsDiffableEndo $ f *** g AffinDiffable _ f *** AffinDiffable _ g = AffinDiffable NotDiffableEndo $ f *** g f *** g = genericiseDifferentiable f *** genericiseDifferentiable g -instance (MetricScalar s) => PreArrow (Differentiable s) where - terminal = Differentiable $ \_ -> (Origin, zeroV, const zeroV) - fst = Differentiable $ \(x,_) -> (x, lfst, const zeroV) - where lfst = linear fst - snd = Differentiable $ \(_,y) -> (y, lsnd, const zeroV) - where lsnd = linear snd +instance (RealFrac' s) => PreArrow (Differentiable s) where + terminal = Differentiable $ \_ -> (Origin, zeroV, const mempty) + fst = Differentiable $ \(x,_) -> (x, fst, const mempty) + snd = Differentiable $ \(_,y) -> (y, snd, const mempty) Differentiable f &&& Differentiable g = Differentiable h - where h x = ((fx, gx), lFanout, devfg) + where h x = ((fx, gx), f'&&&g', devfg) where (fx, f', devf) = f x (gx, g', devg) = g x - devfg δs = (devf $ transformMetric (id&&&zeroV) δs) - ^+^ (devg $ transformMetric (zeroV&&&id) δs) - lFanout = linear $ lapply f'&&&lapply g' + devfg δs = (devf $ transformNorm (id&&&zeroV) δs) + <> (devg $ transformNorm (zeroV&&&id) δs) f &&& g = genericiseDifferentiable f &&& genericiseDifferentiable g -instance (MetricScalar s) => WellPointed (Differentiable s) where +instance (RealFrac' s) => WellPointed (Differentiable s) where unit = Tagged Origin - globalElement x = Differentiable $ \Origin -> (x, zeroV, const zeroV) - const x = Differentiable $ \_ -> (x, zeroV, const zeroV) + globalElement x = Differentiable $ \Origin -> (x, zeroV, const mempty) + const x = Differentiable $ \_ -> (x, zeroV, const mempty) type DfblFuncValue s = GenericAgent (Differentiable s) -instance (MetricScalar s) => HasAgent (Differentiable s) where +instance (RealFrac' s) => HasAgent (Differentiable s) where alg = genericAlg ($~) = genericAgentMap -instance (MetricScalar s) => CartesianAgent (Differentiable s) where +instance (RealFrac' s) => CartesianAgent (Differentiable s) where alg1to2 = genericAlg1to2 alg2to1 = genericAlg2to1 alg2to2 = genericAlg2to2 -instance (MetricScalar s) +instance (RealFrac' s) => PointAgent (DfblFuncValue s) (Differentiable s) a x where point = genericPoint actuallyLinearEndo :: WithField s LinearManifold x - => (x:-*x) -> Differentiable s x x + => (x+>x) -> Differentiable s x x actuallyLinearEndo = AffinDiffable IsDiffableEndo . linearAffine actuallyAffineEndo :: WithField s LinearManifold x - => x -> (x:-*x) -> Differentiable s x 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 + => (x+>y) -> Differentiable s x y actuallyLinear = AffinDiffable NotDiffableEndo . linearAffine actuallyAffine :: ( WithField s LinearManifold x , WithField s AffineManifold y ) - => y -> (x:-*Diff y) -> Differentiable s x y + => y -> (x+>Diff y) -> Differentiable s x y actuallyAffine y₀ f = AffinDiffable NotDiffableEndo $ const y₀ .+^ linearAffine f @@ -419,35 +406,34 @@ dfblFnValsFunc :: ( LocallyScalable s c, LocallyScalable s c', LocallyScalable s d , v ~ Needle c, v' ~ Needle c' - , ε ~ HerMetric v, ε ~ HerMetric v' ) - => (c' -> (c, v':-*v, ε->ε)) -> DfblFuncValue s d c' -> DfblFuncValue s d c + , ε ~ Norm v, ε ~ Norm v' + , RealFrac' s ) + => (c' -> (c, v'+>v, ε->ε)) -> DfblFuncValue s d c' -> DfblFuncValue s d c dfblFnValsFunc f = (Differentiable f $~) dfblFnValsCombine :: forall d c c' c'' v v' v'' ε ε' ε'' s. ( LocallyScalable s c, LocallyScalable s c', LocallyScalable s c'' , LocallyScalable s d , v ~ Needle c, v' ~ Needle c', v'' ~ Needle c'' - , ε ~ HerMetric v , ε' ~ HerMetric v' , ε'' ~ HerMetric v'', ε~ε', ε~ε'' ) - => ( c' -> c'' -> (c, (v',v''):-*v, ε -> (ε',ε'')) ) + , ε ~ Norm v , ε' ~ Norm v' , ε'' ~ Norm v'', ε~ε', ε~ε'' + , RealFrac' s ) + => ( c' -> c'' -> (c, (v',v'')+>v, ε -> (ε',ε'')) ) -> DfblFuncValue s d c' -> DfblFuncValue s d c'' -> DfblFuncValue s d c dfblFnValsCombine cmb (GenericAgent (Differentiable f)) (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'' - jh = convertLinear$->$h' + \d -> let (c', jf, devf) = f d + (c'', jg, devg) = g d + (c, jh, devh) = cmb c' c'' jhl = jh . (id&&&zeroV); jhr = jh . (zeroV&&&id) in ( c - , h' *.* linear (lapply f' &&& lapply g') - , \εc -> let εc' = transformMetric jhl εc - εc'' = transformMetric jhr εc + , jh <<< jf&&&jg + , \εc -> let εc' = transformNorm jhl εc + εc'' = transformNorm jhr εc (δc',δc'') = devh εc - in devf εc' ^+^ devg εc'' - ^+^ transformMetric jf δc' - ^+^ transformMetric jg δc'' + in devf εc' <> devg εc'' + <> transformNorm jf δc' + <> transformNorm jg δc'' ) dfblFnValsCombine cmb (GenericAgent fa) (GenericAgent ga) = dfblFnValsCombine cmb (GenericAgent $ genericiseDifferentiable fa) @@ -457,17 +443,15 @@ -instance (WithField s LinearManifold v, LocallyScalable s a, Floating s) +instance (WithField s LinearManifold v, LocallyScalable s a, RealFloat' s) => AdditiveGroup (DfblFuncValue s a v) where zeroV = point zeroV 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 (^+^) + α^+^β = dfblFnValsCombine (\a b -> (a^+^b, arr addV, const mempty)) α β negateV (GenericAgent (AffinDiffable ef f)) = GenericAgent $ AffinDiffable ef (negateV f) - negateV α = dfblFnValsFunc (\a -> (negateV a, lNegate, const zeroV)) α - where lNegate = linear negateV + negateV α = dfblFnValsFunc (\a -> (negateV a, negateV id, const mempty)) α instance (RealDimension n, LocallyScalable n a) => Num (DfblFuncValue n a n) where @@ -475,8 +459,9 @@ (+) = (^+^) (*) = dfblFnValsCombine $ \a b -> ( a*b - , linear $ \(da,db) -> a*db + b*da - , \d -> let d¹₂ = sqrt d in (d¹₂,d¹₂) + , arr $ addV <<< (scale $ a)***(scale $ b) + , unsafe_dev_ε_δ(show a++"*"++show b) sqrt + >>> \d¹₂ -> (d¹₂,d¹₂) -- ε δa δb = (a+δa)·(b+δb) - (a·b + (a·δa + b·δb)) -- = δa·δb -- so choose δa = δb = √ε @@ -484,14 +469,14 @@ negate = negateV abs = dfblFnValsFunc dfblAbs where dfblAbs a - | a>0 = (a, idL, unsafe_dev_ε_δ("abs "++show a) $ \ε -> a + ε/2) - | a<0 = (-a, negateV idL, unsafe_dev_ε_δ("abs "++show a) $ \ε -> ε/2 - a) - | otherwise = (0, zeroV, (^/ sqrt 2)) + | a>0 = (a, id, unsafe_dev_ε_δ("abs "++show a) $ \ε -> a + ε/2) + | a<0 = (-a, negateV id, unsafe_dev_ε_δ("abs "++show a) $ \ε -> ε/2 - a) + | otherwise = (0, zeroV, scaleNorm (sqrt 0.5)) signum = dfblFnValsFunc dfblSgn where dfblSgn a | a>0 = (1, zeroV, unsafe_dev_ε_δ("signum "++show a) $ const a) | a<0 = (-1, zeroV, unsafe_dev_ε_δ("signum "++show a) $ \_ -> -a) - | otherwise = (0, zeroV, const $ projector 1) + | otherwise = (0, zeroV, const $ spanNorm [1]) @@ -513,19 +498,19 @@ minDblfuncs (Differentiable f) (Differentiable g) = Differentiable h where h x | fx < gx = ( fx, jf - , \d -> devf d ^+^ devg d - ^+^ transformMetric δj - (projector . recip $ recip(metric d 1) + gx - fx) ) + , \d -> devf d <> devg d + <> transformNorm δj + (spanNorm [recip $ recip(d|$|1) + gx - fx]) ) | fx > gx = ( gx, jg - , \d -> devf d ^+^ devg d - ^+^ transformMetric δj - (projector . recip $ recip(metric d 1) + fx - gx) ) + , \d -> devf d <> devg d + <> transformNorm δj + (spanNorm [recip $ recip(d|$|1) + fx - gx]) ) | otherwise = ( fx, (jf^+^jg)^/2 - , \d -> devf d ^+^ devg d - ^+^ transformMetric δj d ) + , \d -> devf d <> devg d + <> transformNorm δj d ) where (fx, jf, devf) = f x (gx, jg, devg) = g x - δj = convertLinear $->$ jf ^-^ jg + δj = jf ^-^ jg postEndo :: ∀ c a b . (HasAgent c, Object c a, Object c b) @@ -580,7 +565,7 @@ positivePreRegion', negativePreRegion' :: (RealDimension s) => PreRegion s s positivePreRegion' = PreRegion $ Differentiable prr where prr x = ( 1 - 1/xp1 - , (1/xp1²) *^ idL + , (1/xp1²) *^ id , unsafe_dev_ε_δ("positivePreRegion@"++show x) δ ) -- ε = (1 − 1/(1+x)) + (-δ · 1/(x+1)²) − (1 − 1/(1+x−δ)) -- = 1/(1+x−δ) − 1/(1+x) − δ · 1/(x+1)² @@ -618,7 +603,7 @@ xp1² = xp1 ^ 2 negativePreRegion' = PreRegion $ ppr . ngt where PreRegion ppr = positivePreRegion' - ngt = actuallyLinearEndo $ linear negate + ngt = actuallyLinearEndo $ negateV id preRegionToInfFrom, preRegionFromMinInfTo :: RealDimension s => s -> PreRegion s s preRegionToInfFrom = RealSubray PositiveHalfSphere @@ -627,16 +612,16 @@ preRegionToInfFrom', preRegionFromMinInfTo' :: RealDimension s => s -> PreRegion s s preRegionToInfFrom' xs = PreRegion $ ppr . trl where PreRegion ppr = positivePreRegion' - trl = actuallyAffineEndo (-xs) idL + trl = actuallyAffineEndo (-xs) id preRegionFromMinInfTo' xe = PreRegion $ ppr . flp where PreRegion ppr = positivePreRegion' - flp = actuallyAffineEndo xe (linear negate) + flp = actuallyAffineEndo xe (negateV id) intervalPreRegion :: RealDimension s => (s,s) -> PreRegion s s intervalPreRegion (lb,rb) = PreRegion $ Differentiable prr where m = lb + radius; radius = (rb - lb)/2 prr x = ( 1 - ((x-m)/radius)^2 - , (2*(m-x)/radius^2) *^ idL + , (2*(m-x)/radius^2) *^ id , unsafe_dev_ε_δ("intervalPreRegion@"++show x) $ (*radius) . sqrt ) @@ -650,7 +635,7 @@ instance (RealDimension s) => Category (RWDiffable s) where - type Object (RWDiffable s) o = LocallyScalable s o + type Object (RWDiffable s) o = (LocallyScalable s o, SimpleSpace (Needle o)) id = RWDiffable $ \x -> (GlobalRegion, pure id) RWDiffable f . RWDiffable g = RWDiffable h where h x₀ = case g x₀ of @@ -663,7 +648,7 @@ -> (rg, fmap (. gr') fhr) (RealSubray diry yl, fhr) -> let hhr = fmap (. gr') fhr - in case lapply ϕg 1 of + in case ϕg $ 1 of y' | y'>0 -> ( unsafePreRegionIntersect rg $ RealSubray diry (x₀ + (yl-y₀)/y') -- y'⋅(xl−x₀) + y₀ ≝ yl @@ -767,7 +752,9 @@ RWDFV_IdVar :: RWDfblFuncValue s c c GenericRWDFV :: RWDiffable s d c -> RWDfblFuncValue s d c -genericiseRWDFV :: (RealDimension s, LocallyScalable s c, LocallyScalable s d) +genericiseRWDFV :: ( RealDimension s + , LocallyScalable s c, SimpleSpace (Needle c) + , LocallyScalable s d, SimpleSpace (Needle d) ) => RWDfblFuncValue s d c -> RWDfblFuncValue s d c genericiseRWDFV (ConstRWDFV c) = GenericRWDFV $ const c genericiseRWDFV RWDFV_IdVar = GenericRWDFV id @@ -795,16 +782,18 @@ :: ( RealDimension s , LocallyScalable s c, LocallyScalable s c', LocallyScalable s d , v ~ Needle c, v' ~ Needle c' - , ε ~ HerMetric v, ε ~ HerMetric v' ) - => (c' -> (c, v':-*v, ε->ε)) -> RWDfblFuncValue s d c' -> RWDfblFuncValue s d c + , SimpleSpace v, SimpleSpace (Needle d) + , ε ~ Norm v, ε ~ Norm v' ) + => (c' -> (c, v'+>v, ε->ε)) -> RWDfblFuncValue s d c' -> RWDfblFuncValue s d c grwDfblFnValsFunc f = (RWDiffable (\_ -> (GlobalRegion, pure (Differentiable f))) $~) grwDfblFnValsCombine :: forall d c c' c'' v v' v'' ε ε' ε'' s. ( LocallyScalable s c, LocallyScalable s c', LocallyScalable s c'' , LocallyScalable s d, RealDimension s , v ~ Needle c, v' ~ Needle c', v'' ~ Needle c'' - , ε ~ HerMetric v , ε' ~ HerMetric v' , ε'' ~ HerMetric v'', ε~ε', ε~ε'' ) - => ( c' -> c'' -> (c, (v',v''):-*v, ε -> (ε',ε'')) ) + , SimpleSpace v, SimpleSpace (Needle d) + , ε ~ Norm v , ε' ~ Norm v' , ε'' ~ Norm v'', ε~ε', ε~ε'' ) + => ( c' -> c'' -> (c, (v',v'')+>v, ε -> (ε',ε'')) ) -> RWDfblFuncValue s d c' -> RWDfblFuncValue s d c'' -> RWDfblFuncValue s d c grwDfblFnValsCombine cmb (GenericRWDFV (RWDiffable fpcs)) (GenericRWDFV (RWDiffable gpcs)) @@ -815,21 +804,18 @@ case (genericiseDifferentiable<$>fmay, genericiseDifferentiable<$>gmay) of (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'' - jh = convertLinear $->$ h' + -> let (c', jf, devf) = f d + (c'',jg, devg) = g d + (c, jh, devh) = cmb c' c'' jhl = jh . (id&&&zeroV); jhr = jh . (zeroV&&&id) in ( c - , h' *.* linear (lapply f' &&& lapply g') - , \εc -> let εc' = transformMetric jhl εc - εc'' = transformMetric jhr εc + , jh <<< jf&&&jg + , \εc -> let εc' = transformNorm jhl εc + εc'' = transformNorm jhr εc (δc',δc'') = devh εc - in devf εc' ^+^ devg εc'' - ^+^ transformMetric jf δc' - ^+^ transformMetric jg δc'' + in devf εc' <> devg εc'' + <> transformNorm jf δc' + <> transformNorm jg δc'' ) _ -> notDefinedHere grwDfblFnValsCombine cmb fv gv @@ -847,7 +833,8 @@ rh = unsafePreRegionIntersect rf rg fgplus :: Differentiable s a v -> Differentiable s a v -> Differentiable s a v fgplus (Differentiable fd) (Differentiable gd) = Differentiable hd - where hd x = (fx^+^gx, jf^+^jg, \ε -> δf(ε^*4) ^+^ δg(ε^*4)) + where hd x = (fx^+^gx, jf^+^jg, \ε -> δf(scaleNorm 2 ε) + <> δg(scaleNorm 2 ε)) where (fx, jf, δf) = fd x (gx, jg, δg) = gd x fgplus (Differentiable fd) (AffinDiffable _ ga) @@ -877,7 +864,8 @@ fneg (AffinDiffable ef af) = AffinDiffable ef $ negateV af postCompRW :: ( RealDimension s - , LocallyScalable s a, LocallyScalable s b, LocallyScalable s c ) + , LocallyScalable s a, LocallyScalable s b, LocallyScalable s c + , SimpleSpace (Needle a), SimpleSpace (Needle b), SimpleSpace (Needle c) ) => RWDiffable s b c -> RWDfblFuncValue s a b -> RWDfblFuncValue s a c postCompRW (RWDiffable f) (ConstRWDFV x) = case f x of (_, Option (Just fd)) -> ConstRWDFV $ fd $ x @@ -885,40 +873,53 @@ postCompRW f (GenericRWDFV g) = GenericRWDFV $ f . g -instance ( WithField s EuclidSpace v, AdditiveGroup v, v ~ Needle (Interior (Needle v)) - , LocallyScalable s a, RealDimension s) +instance ( WithField s EuclidSpace v, SimpleSpace v, v ~ Needle (Interior (Needle v)) + , LocallyScalable s a, SimpleSpace (Needle a), RealDimension s) => AdditiveGroup (RWDfblFuncValue s a v) where zeroV = point zeroV ConstRWDFV c₁ ^+^ ConstRWDFV c₂ = ConstRWDFV (c₁^+^c₂) ConstRWDFV c₁ ^+^ RWDFV_IdVar = GenericRWDFV $ - globalDiffable' (actuallyAffineEndo c₁ idL) + globalDiffable' (actuallyAffineEndo c₁ id) RWDFV_IdVar ^+^ ConstRWDFV c₂ = GenericRWDFV $ - globalDiffable' (actuallyAffineEndo c₂ idL) + globalDiffable' (actuallyAffineEndo c₂ id) ConstRWDFV c₁ ^+^ GenericRWDFV g = GenericRWDFV $ - globalDiffable' (actuallyAffineEndo c₁ idL) . g + globalDiffable' (actuallyAffineEndo c₁ id) . g GenericRWDFV f ^+^ ConstRWDFV c₂ = GenericRWDFV $ - globalDiffable' (actuallyAffineEndo c₂ idL) . f + globalDiffable' (actuallyAffineEndo c₂ id) . 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' (actuallyLinearEndo $ linear negateV) + negateV RWDFV_IdVar = GenericRWDFV $ globalDiffable' (actuallyLinearEndo $ negateV id) negateV (GenericRWDFV f) = GenericRWDFV $ rwDfbl_negateV f -instance (RealDimension n, LocallyScalable n a) +dualCoCoProduct :: ∀ v w s . + ( SimpleSpace v, HilbertSpace v + , SimpleSpace w, Scalar v ~ s, Scalar w ~ s ) + => LinearMap s w v -> LinearMap s w v -> Norm w +dualCoCoProduct s t = Norm $ (tSpread*sSpread) *^ t²Ps²M + where t' = adjoint $ t :: LinearMap s v (DualVector w) + s' = adjoint $ s :: LinearMap s v (DualVector w) + tSpread = sum . map recip_t²PLUSs² $ snd (decomposeLinMap t') [] + sSpread = sum . map recip_t²PLUSs² $ snd (decomposeLinMap s') [] + t²PLUSs²@(Norm t²Ps²M) + = transformNorm t euclideanNorm <> transformNorm s euclideanNorm :: Norm w + recip_t²PLUSs² = normSq (dualNorm t²PLUSs²) :: DualVector w -> s + +instance (RealDimension n, LocallyScalable n a, SimpleSpace (Needle a)) => Num (RWDfblFuncValue n a n) where fromInteger i = point $ fromInteger i (+) = (^+^) ConstRWDFV c₁ * ConstRWDFV c₂ = ConstRWDFV (c₁*c₂) ConstRWDFV c₁ * RWDFV_IdVar = GenericRWDFV $ - globalDiffable' (actuallyLinearEndo $ linear (c₁*)) + globalDiffable' (actuallyLinearEndo . arr $ scale $ c₁) RWDFV_IdVar * ConstRWDFV c₂ = GenericRWDFV $ - globalDiffable' (actuallyLinearEndo $ linear (*c₂)) + globalDiffable' (actuallyLinearEndo . arr $ scale $ c₂) ConstRWDFV c₁ * GenericRWDFV g = GenericRWDFV $ - globalDiffable' (actuallyLinearEndo $ linear (c₁*)) . g + globalDiffable' (actuallyLinearEndo . arr $ scale $ c₁) . g GenericRWDFV f * ConstRWDFV c₂ = GenericRWDFV $ - globalDiffable' (actuallyLinearEndo $ linear (*c₂)) . f + globalDiffable' (actuallyLinearEndo . arr $ scale $ c₂) . f f*g = genericiseRWDFV f ⋅ genericiseRWDFV g - where (⋅) :: ∀ n a . (RealDimension n, LocallyScalable n a) + where (⋅) :: ∀ n a . (RealDimension n, LocallyScalable n a, SimpleSpace (Needle a)) => RWDfblFuncValue n a n -> RWDfblFuncValue n a n -> RWDfblFuncValue n a n GenericRWDFV (RWDiffable fpcs) ⋅ GenericRWDFV (RWDiffable gpcs) = GenericRWDFV . RWDiffable $ @@ -934,28 +935,26 @@ 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' + jf = ϕf $ 1; jg = ϕg $ 1 + invf'g' = recip $ jf*jg in ( fd*gd - , linear.(*)$ fd*g' + gd*f' + , arr $ scale $ fd*jg + gd*jf , 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 + \d -> let (c₁, jf, devf) = f d + (c₂, jg, devg) = g d c = c₁*c₂; c₁² = c₁^2; c₂² = c₂^2 - h' = c₁*^slg ^+^ c₂*^slf + h' = c₁*^jg ^+^ c₂*^jf in ( c , h' - , \εc -> let rε² = metric εc 1 - c₁worst² = c₁² + recip(1 + c₂²*rε²) - c₂worst² = c₂² + recip(1 + c₁²*rε²) - in (4*rε²) *^ dualCoCoProduct jf jg - ^+^ devf (εc^*(4*c₂worst²)) - ^+^ devg (εc^*(4*c₁worst²)) + , \εc -> let rε = εc|$|1 + c₁worst = sqrt $ c₁² + recip(1 + c₂²*rε^2) + c₂worst = sqrt $ c₂² + recip(1 + c₁²*rε^2) + in scaleNorm (2*rε) (dualCoCoProduct jf jg) + <> devf (scaleNorm (2*c₂worst) εc) + <> devg (scaleNorm (2*c₁worst) εc) -- TODO: add formal proof for this (or, if necessary, the correct form) ) mulDi f g = mulDi (genericiseDifferentiable f) (genericiseDifferentiable g) @@ -965,20 +964,20 @@ where absPW a₀ | a₀<0 = (negativePreRegion, pure desc) | otherwise = (positivePreRegion, pure asc) - desc = actuallyLinearEndo $ linear negate - asc = actuallyLinearEndo idL + desc = actuallyLinearEndo $ negateV id + asc = actuallyLinearEndo id signum = (RWDiffable sgnPW $~) where sgnPW a₀ | a₀<0 = (negativePreRegion, pure (const $ -1)) | otherwise = (positivePreRegion, pure (const 1)) -instance (RealDimension n, LocallyScalable n a) +instance (RealDimension n, LocallyScalable n a, SimpleSpace (Needle a)) => Fractional (RWDfblFuncValue n a n) where fromRational i = point $ fromRational i recip = postCompRW . RWDiffable $ \a₀ -> if a₀<0 then (negativePreRegion, pure (Differentiable negp)) else (positivePreRegion, pure (Differentiable posp)) - where negp x = (x'¹, (- x'¹^2) *^ idL, unsafe_dev_ε_δ("1/"++show x) δ) + where negp x = (x'¹, (- x'¹^2) *^ id, unsafe_dev_ε_δ("1/"++show x) δ) -- ε = 1/x − δ/x² − 1/(x+δ) -- ε·x + ε·δ = 1 + δ/x − δ/x − δ²/x² − 1 -- = -δ²/x² @@ -991,7 +990,7 @@ else - x -- numerical underflow of εx³ vs mph -- ≡ ε*x^3 / (2*mph) (Taylor-expansion of the root) x'¹ = recip x - posp x = (x'¹, (- x'¹^2) *^ idL, unsafe_dev_ε_δ("1/"++show x) δ) + posp x = (x'¹, (- x'¹^2) *^ id, unsafe_dev_ε_δ("1/"++show x) δ) where δ ε = let mph = ε*x^2/2 δ₀ = sqrt (mph^2 + ε*x^3) - mph in if δ₀>0 then δ₀ else x @@ -1000,7 +999,7 @@ -instance (RealDimension n, LocallyScalable n a) +instance (RealDimension n, LocallyScalable n a, SimpleSpace (Needle a)) => Floating (RWDfblFuncValue n a n) where pi = point pi @@ -1008,8 +1007,8 @@ $ \x -> let ex = exp x in if ex*2 == ex -- numerical trouble... then if x<0 then ( 0, zeroV, unsafe_dev_ε_δ("exp "++show x) $ \ε -> log ε - x ) - else ( ex, ex*^idL, unsafe_dev_ε_δ("exp "++show x) $ \_ -> 1e-300 ) - else ( ex, ex *^ idL, unsafe_dev_ε_δ("exp "++show x) + else ( ex, ex*^id, unsafe_dev_ε_δ("exp "++show x) $ \_ -> 1e-300 ) + else ( ex, ex *^ id, unsafe_dev_ε_δ("exp "++show x) $ \ε -> case acosh(ε/(2*ex) + 1) of δ | δ==δ -> δ | otherwise -> log ε - x ) @@ -1022,7 +1021,7 @@ log = postCompRW . RWDiffable $ \x -> if x>0 then (positivePreRegion, pure (Differentiable lnPosR)) else (negativePreRegion, notDefinedHere) - where lnPosR x = ( log x, recip x *^ idL, unsafe_dev_ε_δ("log "++show x) $ \ε -> x * sqrt(1 - exp(-ε)) ) + where lnPosR x = ( log x, recip x *^ id, unsafe_dev_ε_δ("log "++show x) $ \ε -> x * sqrt(1 - exp(-ε)) ) -- ε = ln x + (-δ)/x − ln(x−δ) -- = ln (x / ((x−δ) · exp(δ/x))) -- x/e^ε = (x−δ) · exp(δ/x) @@ -1036,13 +1035,13 @@ sqrt = postCompRW . RWDiffable $ \x -> if x>0 then (positivePreRegion, pure (Differentiable sqrtPosR)) else (negativePreRegion, notDefinedHere) - where sqrtPosR x = ( sx, idL ^/ (2*sx), unsafe_dev_ε_δ("sqrt "++show x) $ + where sqrtPosR x = ( sx, id ^/ (2*sx), unsafe_dev_ε_δ("sqrt "++show x) $ \ε -> 2 * (s2 * sqrt sx^3 * sqrt ε + signum (ε*2-sx) * sx * ε) ) where sx = sqrt x; s2 = sqrt 2 -- Exact inverse of O(δ²) remainder. sin = grwDfblFnValsFunc sinDfb - where sinDfb x = ( sx, cx *^ idL, unsafe_dev_ε_δ("sin "++show x) δ ) + where sinDfb x = ( sx, cx *^ id, unsafe_dev_ε_δ("sin "++show x) δ ) where sx = sin x; cx = cos x sx² = sx^2; cx² = cx^2 sx' = abs sx; cx' = abs cx @@ -1057,7 +1056,7 @@ -- Safety margins for overlap between quadratic and cubic model -- (these aren't naturally compatible to be used both together) - cos = sin . (globalDiffable' (actuallyAffineEndo (pi/2) idL) $~) + cos = sin . (globalDiffable' (actuallyAffineEndo (pi/2) id) $~) sinh x = (exp x - exp (-x))/2 {- = grwDfblFnValsFunc sinhDfb @@ -1071,7 +1070,7 @@ cosh x = (exp x + exp (-x))/2 tanh = grwDfblFnValsFunc tanhDfb - where tanhDfb x = ( tnhx, idL ^/ (cosh x^2), unsafe_dev_ε_δ("tan "++show x) δ ) + where tanhDfb x = ( tnhx, id ^/ (cosh x^2), unsafe_dev_ε_δ("tan "++show x) δ ) where tnhx = tanh x c = (tnhx*2/pi)^2 p = 1 + abs x/(2*pi) @@ -1080,7 +1079,7 @@ -- with quite a big margin. TODO: find a tighter definition. atan = grwDfblFnValsFunc atanDfb - where atanDfb x = ( atnx, idL ^/ (1+x^2), unsafe_dev_ε_δ("atan "++show x) δ ) + where atanDfb x = ( atnx, id ^/ (1+x^2), unsafe_dev_ε_δ("atan "++show x) δ ) where atnx = atan x c = (atnx*2/pi)^2 p = 1 + abs x/(2*pi) @@ -1096,7 +1095,7 @@ | x < (-1) -> (preRegionFromMinInfTo (-1), notDefinedHere) | x > 1 -> (preRegionToInfFrom 1, notDefinedHere) | otherwise -> (intervalPreRegion (-1,1), pure (Differentiable asinDefdR)) - where asinDefdR x = ( asinx, asin'x *^ idL, unsafe_dev_ε_δ("asin "++show x) δ ) + where asinDefdR x = ( asinx, asin'x *^ id, unsafe_dev_ε_δ("asin "++show x) δ ) where asinx = asin x; asin'x = recip (sqrt $ 1 - x^2) c = 1 - x^2 δ ε = sqrt ε * c @@ -1106,13 +1105,13 @@ | x < (-1) -> (preRegionFromMinInfTo (-1), notDefinedHere) | x > 1 -> (preRegionToInfFrom 1, notDefinedHere) | otherwise -> (intervalPreRegion (-1,1), pure (Differentiable acosDefdR)) - where acosDefdR x = ( acosx, acos'x *^ idL, unsafe_dev_ε_δ("acos "++show x) δ ) + where acosDefdR x = ( acosx, acos'x *^ id, unsafe_dev_ε_δ("acos "++show x) δ ) where acosx = acos x; acos'x = - recip (sqrt $ 1 - x^2) c = 1 - x^2 δ ε = sqrt ε * c -- Like for asin – it's just a translation/reflection. asinh = grwDfblFnValsFunc asinhDfb - where asinhDfb x = ( asinhx, idL ^/ sqrt(1+x^2), unsafe_dev_ε_δ("asinh "++show x) δ ) + where asinhDfb x = ( asinhx, id ^/ sqrt(1+x^2), unsafe_dev_ε_δ("asinh "++show x) δ ) where asinhx = asinh x δ ε = abs x * sqrt((1 - exp(-ε))*0.8 + ε^2/(3*abs x + 1)) + sqrt(ε/(abs x+0.5)) -- Empirical, modified from log function (the area hyperbolic sine @@ -1121,7 +1120,7 @@ acosh = postCompRW . RWDiffable $ \x -> if x>1 then (preRegionToInfFrom 1, pure (Differentiable acoshDfb)) else (preRegionFromMinInfTo 1, notDefinedHere) - where acoshDfb x = ( acosh x, idL ^/ sqrt(x^2 - 1), unsafe_dev_ε_δ("acosh "++show x) δ ) + where acoshDfb x = ( acosh x, id ^/ sqrt(x^2 - 1), unsafe_dev_ε_δ("acosh "++show x) δ ) where δ ε = (2 - 1/sqrt x) * (s2 * sqrt sx^3 * sqrt(ε/s2) + signum (ε*s2-sx) * sx * ε/s2) sx = sqrt(x-1) s2 = sqrt 2 @@ -1132,7 +1131,7 @@ | x < (-1) -> (preRegionFromMinInfTo (-1), notDefinedHere) | x > 1 -> (preRegionToInfFrom 1, notDefinedHere) | otherwise -> (intervalPreRegion (-1,1), pure (Differentiable atnhDefdR)) - where atnhDefdR x = ( atanh x, recip(1-x^2) *^ idL, unsafe_dev_ε_δ("atanh "++show x) $ \ε -> sqrt(tanh ε)*(1-abs x) ) + where atnhDefdR x = ( atanh x, recip(1-x^2) *^ id, unsafe_dev_ε_δ("atanh "++show x) $ \ε -> sqrt(tanh ε)*(1-abs x) ) -- Empirical, with epsEst upper bound. @@ -1145,21 +1144,27 @@ -- -- However, because this category allows functions to be undefined in some region, -- such decisions can be faked quite well: '?->' restricts a function to --- some region, by simply marking it undefined outside¹, and '?|:' replaces these +-- some region, by simply marking it undefined outside, and '?|:' replaces these -- regions with values from another function. -- -- Example: define a function that is compactly supported on the interval ]-1,1[, -- i.e. exactly zero everywhere outside. -- -- @ --- Graphics.Dynamic.Plot.R2> plotWindow [diffableFnPlot (\\x -> -1 '?<' x '?<' 1 '?->' exp(1/(x^2 - 1)) '?|:' 0)] +-- Graphics.Dynamic.Plot.R2> plotWindow [fnPlot (\\x -> -1 '?<' x '?<' 1 '?->' cos (x*pi/2)^2 '?|:' 0)] -- @ -- --- <> +-- <> -- --- ¹ Note that it may not be necessary to restrict explicitly: for instance if a +-- Note that it may not be necessary to restrict explicitly: for instance if a -- square root appears somewhere in an expression, then the expression is automatically -- restricted so that the root has a positive argument! +-- +-- @ +-- Graphics.Dynamic.Plot.R2> plotWindow [fnPlot (\\x -> sqrt x '?|:' -sqrt (-x))] +-- @ +-- +-- <> infixr 4 ?-> -- | Require the LHS to be defined before considering the RHS as result. @@ -1170,7 +1175,8 @@ -- Just _ 'Control.Applicative.*>' a = a -- _ 'Control.Applicative.*>' a = Nothing -- @ -(?->) :: (RealDimension n, LocallyScalable n a, LocallyScalable n b, LocallyScalable n c) +(?->) :: ( RealDimension n, LocallyScalable n a, LocallyScalable n b, LocallyScalable n c + , SimpleSpace (Needle b), SimpleSpace (Needle c) ) => RWDfblFuncValue n c a -> RWDfblFuncValue n c b -> RWDfblFuncValue n c b ConstRWDFV _ ?-> f = f RWDFV_IdVar ?-> f = f @@ -1196,12 +1202,12 @@ -- allows chaining of comparison operators like in Python.) -- Note that less-than comparison is -- to less-or-equal comparison, because there is no such thing as equality. -(?>) :: (RealDimension n, LocallyScalable n a) +(?>) :: (RealDimension n, LocallyScalable n a, SimpleSpace (Needle a)) => RWDfblFuncValue n a n -> RWDfblFuncValue n a n -> RWDfblFuncValue n a n a ?> b = (positiveRegionalId $~ a-b) ?-> b -- | Return the RHS, if it is greater than the LHS. -(?<) :: (RealDimension n, LocallyScalable n a) +(?<) :: (RealDimension n, LocallyScalable n a, SimpleSpace (Needle 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 @@ -1224,7 +1230,8 @@ -- @ -- -- Basically a weaker and agent-ised version of 'backupRegions'. -(?|:) :: (RealDimension n, LocallyScalable n a, LocallyScalable n b) +(?|:) :: ( RealDimension n, LocallyScalable n a, LocallyScalable n b + , SimpleSpace (Needle a), SimpleSpace (Needle b) ) => RWDfblFuncValue n a b -> RWDfblFuncValue n a b -> RWDfblFuncValue n a b ConstRWDFV c ?|: _ = ConstRWDFV c RWDFV_IdVar ?|: _ = RWDFV_IdVar @@ -1257,16 +1264,12 @@ -- 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)) - +lerp_diffable a b = actuallyAffine a . arr $ flipBilin scale $ b.-.a -isZeroMap :: ∀ v a . (FiniteDimensional v, AdditiveGroup a, Eq a) => (v:-*a) -> Bool -isZeroMap m = all ((==zeroV) . atBasis m) b - where (Tagged b) = completeBasis :: Tagged v [Basis v] 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 @@ -6,8 +6,7 @@ import Data.Semigroup import Data.Function.Affine import Data.VectorSpace -import Data.LinearMap -import Data.LinearMap.HerMetric +import Math.LinearMap.Category import Data.Manifold.Types.Primitive import Data.Manifold.PseudoAffine @@ -55,7 +54,7 @@ -- 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 + , Needle d +> Needle c -- Jacobian , LinDevPropag d c -- Metric showing how far you can go -- from x₀ without deviating from the -- Taylor-1 approximation by more than diff --git a/Data/LinearMap/Category.hs b/Data/LinearMap/Category.hs deleted file mode 100644 --- a/Data/LinearMap/Category.hs +++ /dev/null @@ -1,313 +0,0 @@ --- | --- Module : Data.LinearMap.Category --- Copyright : (c) Justus Sagemüller 2015 --- License : GPL v3 --- --- Maintainer : (@) sagemueller $ geo.uni-koeln.de --- Stability : experimental --- Portability : portable --- -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE UndecidableInstances #-} -{-# LANGUAGE StandaloneDeriving #-} -{-# LANGUAGE DeriveGeneric #-} -{-# LANGUAGE DeriveFunctor #-} -{-# LANGUAGE DeriveFoldable #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE GADTs #-} -{-# LANGUAGE RankNTypes #-} -{-# LANGUAGE TupleSections #-} -{-# LANGUAGE UnicodeSyntax #-} -{-# LANGUAGE CPP #-} -{-# LANGUAGE ConstraintKinds #-} -{-# LANGUAGE PatternGuards #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE DataKinds #-} - -module Data.LinearMap.Category where - -import Data.Tagged - -import Data.VectorSpace -import Data.LinearMap -import Data.VectorSpace.FiniteDimensional -import Data.AffineSpace -import Data.Basis - -import qualified Prelude as Hask hiding(foldl) -import qualified Control.Applicative as Hask -import qualified Control.Monad as Hask -import qualified Data.Foldable as Hask - - -import Control.Category.Constrained.Prelude hiding ((^)) -import Control.Arrow.Constrained - -import Data.Manifold.Types.Primitive -import Data.CoNat -import Data.Embedding - -import qualified Data.Vector as Arr -import qualified Numeric.LinearAlgebra.HMatrix as HMat - - - --- | A linear mapping between finite-dimensional spaces, implemeted as a dense matrix. --- --- 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 - --- | Coerce the matrix representations of two linear mappings; the result makes --- sense iff the spaces are canonically isomorphic (certainly if they --- are good instances of 'Data.Manifold.PseudoAffine.LocallyCoercible'). -unsafeCoerceLinear :: Linear s a b -> Linear s c d -unsafeCoerceLinear (DenseLinear m) = DenseLinear m - -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 - DenseLinear f . DenseLinear g = DenseLinear $ HMat.mul f g - -instance (SmoothScalar s) => Cartesian (Linear s) where - type UnitObject (Linear s) = ZeroDim s - swap = lSwap - where lSwap :: forall v w s - . (FiniteDimensional v, FiniteDimensional w, Scalar v~s, Scalar w~s) - => Linear s (v,w) (w,v) - lSwap = DenseLinear $ HMat.assoc (n,n) 0 l - where l = [ ((i,i+nv), 1) | i<-[0.. nw-1] ] ++ [ ((i+nw,i), 1) | i<-[0.. nv-1] ] - (Tagged nv) = dimension :: Tagged v Int - (Tagged nw) = dimension :: Tagged w Int - n = nv + nw - attachUnit = identMat - detachUnit = identMat - regroup = identMat - regroup' = identMat - -instance (SmoothScalar s) => Morphism (Linear s) where - DenseLinear f *** DenseLinear g = DenseLinear $ HMat.diagBlock [f,g] - -instance (SmoothScalar s) => PreArrow (Linear s) where - DenseLinear f &&& DenseLinear g = DenseLinear $ HMat.fromBlocks [[f], [g]] - fst = lFst - where lFst :: forall v w s - . (FiniteDimensional v, FiniteDimensional w, Scalar v~s, Scalar w~s) - => Linear s (v,w) v - lFst = DenseLinear $ HMat.assoc (nv,n) 0 l - where l = [ ((i,i), 1) | i<-[0.. nv-1] ] - (Tagged nv) = dimension :: Tagged v Int - (Tagged nw) = dimension :: Tagged w Int - n = nv + nw - snd = lSnd - where lSnd :: forall v w s - . (FiniteDimensional v, FiniteDimensional w, Scalar v~s, Scalar w~s) - => Linear s (v,w) w - lSnd = DenseLinear $ HMat.assoc (nw,n) 0 l - where l = [ ((i,i+nv), 1) | i<-[0.. nw-1] ] - (Tagged nv) = dimension :: Tagged v Int - (Tagged nw) = dimension :: Tagged w Int - n = nv + nw - terminal = lTerminal - where lTerminal :: forall v s . (FiniteDimensional v, Scalar v~s) - => Linear s v (ZeroDim s) - lTerminal = DenseLinear $ (0 HMat.>< n) [] - where (Tagged n) = dimension :: Tagged v Int - -instance (SmoothScalar s) => EnhancedCat (->) (Linear s) where - arr (DenseLinear mat) = fromPackedVector . HMat.app mat . asPackedVector - --- | Inverse function application (for isomorphisms), or --- least-square solution of a linear equation. --- Note that least-square is not really well-defined, --- without reference to a norm / scalar product; the operator uses --- the implicit norm induced from the 'FiniteDimensional' representation. -(<\$) :: ( SmoothScalar s, FiniteDimensional v, FiniteDimensional w - , Scalar v ~ s, Scalar w ~ s - ) => Linear s v w -> w -> v -DenseLinear mat <\$ v = fromPackedVector . (mat HMat.<\>) $ asPackedVector v - -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 - - - -canonicalIdentityMatrix :: forall n v s - . (KnownNat n, IsFreeSpace v, FreeDimension v ~ n, Scalar v ~ s) - => Linear s v (FreeVect n s) -canonicalIdentityMatrix = DenseLinear $ HMat.ident n - where (Tagged n) = theNatN :: Tagged n Int - --- | Class of spaces that directly represent a free vector space, i.e. that are simply --- @n@-fold products of the base field. --- This class basically contains 'ℝ', 'ℝ²', 'ℝ³' etc., in future also the complex and --- probably integral versions. -class (FiniteDimensional v, KnownNat (FreeDimension v)) => IsFreeSpace v where - type FreeDimension v :: Nat - identityMatrix :: Isomorphism (Linear (Scalar v)) - v - (FreeVect (FreeDimension v) (Scalar v)) - identityMatrix = fromInversePair emb proj - where emb@(DenseLinear i) = canonicalIdentityMatrix - proj = DenseLinear i - -instance (KnownNat n, Num s, SmoothScalar s) => IsFreeSpace (FreeVect n s) where - type FreeDimension (FreeVect n s) = n - identityMatrix = fromInversePair id id - -instance IsFreeSpace ℝ where - type FreeDimension ℝ = S Z - -instance ( SmoothScalar s, IsFreeSpace v, Scalar v ~ s, FiniteDimensional s, s ~ Scalar s ) - => IsFreeSpace (v,s) where - type FreeDimension (v,s) = S (FreeDimension v) - - - -class VectorSpace v => FreeTuple v where - type Tuplity v :: Nat - freeTuple :: Isomorphism (->) v (FreeVect (Tuplity v) (Scalar v)) - -#define FreeScalar(s) \ -instance FreeTuple (s) where { \ - type Tuplity (s) = S Z; \ - freeTuple = fromInversePair (FreeVect . pure) (\(FreeVect v) -> v Arr.! 0); } - -#define FreePair(s) \ -FreeScalar(s); \ -instance FreeTuple (s,s) where { \ - type Tuplity (s,s) = S(S Z); \ - freeTuple = fromInversePair (\(a,b) -> FreeVect $ Arr.fromList[a,b]) \ - (\(FreeVect v) -> (v Arr.! 0, v Arr.! 1)); } - -#define FreeTriple(s) \ -FreePair(s); \ -instance FreeTuple (s,s,s) where { \ - type Tuplity (s,s,s) = S(S(S Z)); \ - freeTuple = fromInversePair (\(a,b,c) -> FreeVect $ Arr.fromList[a,b,c]) \ - (\(FreeVect v) -> (v Arr.! 0, v Arr.! 1, v Arr.! 2)); };\ -instance FreeTuple (s,(s,s)) where { \ - type Tuplity (s,(s,s)) = S(S(S Z)); \ - freeTuple = fromInversePair (\(a,(b,c)) -> FreeVect $ Arr.fromList[a,b,c]) \ - (\(FreeVect v) -> (v Arr.! 0, (v Arr.! 1, v Arr.! 2))); };\ -instance FreeTuple ((s,s),s) where { \ - type Tuplity ((s,s),s) = S(S(S Z)); \ - freeTuple = fromInversePair (\((a,b),c) -> FreeVect $ Arr.fromList[a,b,c]) \ - (\(FreeVect v) -> ((v Arr.! 0, v Arr.! 1), v Arr.! 2)); } - -FreeTriple(ℝ) -FreeTriple(Int) - - diff --git a/Data/LinearMap/HerMetric.hs b/Data/LinearMap/HerMetric.hs deleted file mode 100644 --- a/Data/LinearMap/HerMetric.hs +++ /dev/null @@ -1,894 +0,0 @@ -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE TupleSections #-} -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE UndecidableInstances #-} -{-# LANGUAGE StandaloneDeriving #-} -{-# LANGUAGE ConstraintKinds #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE UnicodeSyntax #-} -{-# LANGUAGE LambdaCase #-} - - - - -module Data.LinearMap.HerMetric ( - -- * Metric operator types - HerMetric(..), HerMetric'(..) - -- * Evaluating metrics - , toDualWith, fromDualWith - , metricSq, metricSq', metric, metric', metrics, metrics' - -- * Defining metrics - , projector, projector', projectors, projector's - , euclideanMetric' - -- * Metrics induce inner products - , spanHilbertSubspace - , spanSubHilbertSpace - , IsFreeSpace - -- * One-dimensional axes and product spaces - , factoriseMetric, factoriseMetric' - , productMetric, productMetric' - , tryMetricAsLength, metricAsLength, metricFromLength, metric'AsLength - -- * Utility for metrics - , transformMetric, transformMetric', dualCoCoProduct - , dualiseMetric, dualiseMetric' - , recipMetric, recipMetric', safeRecipMetric, safeRecipMetric' - -- ** Eigenvectors - , eigenSpan, eigenSpan' - , eigenCoSpan, eigenCoSpan' - , eigenSystem, HasEigenSystem, EigenVector - -- ** Scaling operations - , metriNormalise, metriNormalise' - , metriScale', metriScale - , volumeRatio, euclideanRelativeMetricVolume - , adjoint - , extendMetric - , applyLinMapMetric, applyLinMapMetric' - , imitateMetricSpanChange - -- * The dual-space class - , HasMetric - , HasMetric'(..) - , (^<.>) --- , riesz, riesz' - -- * Fundamental requirements - , MetricScalar - , FiniteDimensional(..) - -- * Misc - , Stiefel1(..) - , linMapAsTensProd, linMapFromTensProd - , covariance - , outerProducts - , orthogonalComplementSpan - ) where - - - - -import Data.VectorSpace -import Data.LinearMap -import Data.Basis -import Data.Semigroup -import Data.Tagged -import qualified Data.List as List - -import qualified Prelude as Hask -import qualified Control.Applicative as Hask -import qualified Control.Monad as Hask - -import Control.Category.Constrained.Prelude hiding ((^)) -import Control.Arrow.Constrained - -import Data.Manifold.Types.Primitive -import Data.CoNat - -import qualified Data.Vector as Arr -import qualified Numeric.LinearAlgebra.HMatrix as HMat - -import Data.VectorSpace.FiniteDimensional -import Data.LinearMap.Category -import Data.Embedding - - - -infixr 7 <.>^, ^<.> - - --- | 'HerMetric' is a portmanteau of /Hermitian/ and /metric/ (in the sense as --- used in e.g. general relativity – though those particular ones aren't positive --- definite and thus not really metrics). --- --- Mathematically, there are two directly equivalent ways to describe such a metric: --- as a bilinear mapping of two vectors to a scalar, or as a linear mapping from --- a vector space to its dual space. We choose the latter, though you can always --- as well think of metrics as “quadratic dual vectors”. --- --- 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 { - metricMatrix :: Maybe (Linear (Scalar v) v (DualSpace v)) -- @Nothing@ for zero metric. - } - -matrixMetric :: HasMetric v => HMat.Matrix (Scalar v) -> HerMetric v -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 $ negateV <$> m - HerMetric Nothing ^+^ HerMetric n = HerMetric n - HerMetric m ^+^ HerMetric Nothing = HerMetric m - 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 $ (s*^) <$> m - --- | A metric on the dual space; equivalent to a linear mapping from the dual space --- to the original vector space. --- --- 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 (Linear (Scalar v) (DualSpace v) v) - } - -extendMetric :: (HasMetric v, Scalar v~ℝ) => HerMetric v -> v -> HerMetric v -extendMetric (HerMetric Nothing) _ = HerMetric Nothing -extendMetric (HerMetric (Just (DenseLinear m))) v - | isInfinite' detm = HerMetric . Just $ DenseLinear m - | isInfinite' detmninv = singularMetric - | otherwise = HerMetric . Just $ DenseLinear mn - where -- this could probably be done much more efficiently, with only - -- multiplications, no inverses. - (minv, (detm, _)) = HMat.invlndet m - (mn, (detmninv, _)) = HMat.invlndet (minv + HMat.outer vv vv) - vv = asPackedVector v - - -matrixMetric' :: HasMetric v => HMat.Matrix (Scalar v) -> HerMetric' v -matrixMetric' = HerMetric' . Just . DenseLinear - --- | Deprecated -instance (HasMetric v) => AdditiveGroup (HerMetric' v) where - zeroV = HerMetric' Nothing - 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 -instance HasMetric v => VectorSpace (HerMetric' v) where - type Scalar (HerMetric' v) = Scalar v - s *^ (HerMetric' m) = HerMetric' $ (s*^) <$> m - - --- | A metric on @v@ that simply yields the squared overlap of a vector with the --- given dual-space reference. --- --- It will perhaps be the most common way of defining 'HerMetric' values to start --- with such dual-space vectors and superimpose the projectors using the 'VectorSpace' --- instance; e.g. @'projector' (1,0) '^+^' 'projector' (0,2)@ yields a hermitian operator --- describing the ellipsoid span of the vectors /e/₀ and 2⋅/e/₁. --- Metrics generated this way are positive definite if no negative coefficients have --- been introduced with the '*^' scaling operator or with '^-^'. --- --- Note: @projector a ^+^ projector b ^+^ ...@ is more efficiently written as --- @projectors [a, b, ...]@ -projector :: HasMetric v => DualSpace v -> HerMetric v -projector u = HerMetric . pure $ u ⊗ u - -projector' :: HasMetric v => v -> HerMetric' v -projector' v = HerMetric' . pure $ v ⊗ v - --- | Efficient shortcut for the 'sumV' of multiple 'projector's. -projectors :: HasMetric v => [DualSpace v] -> HerMetric v -projectors [] = zeroV -projectors us = HerMetric . pure . outerProducts $ zip us us - -projector's :: HasMetric v => [v] -> HerMetric' v -projector's [] = zeroV -projector's vs = HerMetric' . pure . outerProducts $ zip vs vs - - -singularMetric :: forall v . HasMetric v => HerMetric v -singularMetric = matrixMetric $ HMat.scale (1/0) (HMat.ident dim) - where (Tagged dim) = dimension :: Tagged v Int -singularMetric' :: forall v . HasMetric v => HerMetric' v -singularMetric' = matrixMetric' $ HMat.scale (1/0) (HMat.ident dim) - where (Tagged dim) = dimension :: Tagged v Int - - - --- | Evaluate a vector through a metric. For the canonical metric on a Hilbert space, --- this will be simply 'magnitudeSq'. -metricSq :: HasMetric v => HerMetric v -> v -> Scalar v -metricSq (HerMetric Nothing) _ = 0 -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 (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 --- mathematical metric, i.e. positive definite – otherwise the internally used --- square root may get negative arguments (though it can still produce results if the --- scalars are complex; however, complex spaces aren't supported yet). -metric :: (HasMetric v, Floating (Scalar v)) => HerMetric v -> v -> Scalar v -metric m = sqrt . metricSq m - -metric' :: (HasMetric v, Floating (Scalar v)) => HerMetric' v -> DualSpace v -> Scalar v -metric' m = sqrt . metricSq' m - - -toDualWith :: HasMetric v => HerMetric v -> v -> DualSpace v -toDualWith (HerMetric Nothing) = const zeroV -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 -metriNormalise m v = v ^/ metric m v - -metriNormalise' :: (HasMetric v, Floating (Scalar v)) - => HerMetric' v -> DualSpace v -> DualSpace v -metriNormalise' m v = v ^/ metric' m v - --- | “Anti-normalise” a vector: /multiply/ with its own norm, according to metric. -metriScale :: (HasMetric v, Floating (Scalar v)) => HerMetric v -> v -> v -metriScale m v = metric m v *^ v - -metriScale' :: (HasMetric v, Floating (Scalar v)) - => HerMetric' v -> DualSpace v -> DualSpace v -metriScale' m v = metric' m v *^ v - - --- | Square-sum over the metrics for each dual-space vector. --- --- @ --- metrics m vs ≡ sqrt . sum $ metricSq m '<$>' vs --- @ -metrics :: (HasMetric v, Floating (Scalar v)) => HerMetric v -> [v] -> Scalar v -metrics m vs = sqrt . sum $ metricSq m <$> vs - -metrics' :: (HasMetric v, Floating (Scalar v)) => HerMetric' v -> [DualSpace v] -> Scalar v -metrics' m vs = sqrt . sum $ metricSq' m <$> vs - - -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)) = HerMetric . Just $ adjoint t . m . t - -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)) = 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 ~ 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 tArr = HMat.flatten tmat - sArr = HMat.flatten smat - t²PLUSs² = tmat HMat.<> HMat.tr tmat + smat HMat.<> HMat.tr smat - - --- | This doesn't really do anything at all, since @'HerMetric' v@ is essentially a --- synonym for @'HerMetric' ('DualSpace' v)@. -dualiseMetric :: HasMetric v => HerMetric (DualSpace v) -> HerMetric' v -dualiseMetric (HerMetric m) = HerMetric' m - -dualiseMetric' :: HasMetric v => HerMetric' v -> HerMetric (DualSpace v) -dualiseMetric' (HerMetric' m) = HerMetric m - - --- | 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 m' | Option (Just m) <- safeRecipMetric m' = m -recipMetric _ = singularMetric - -recipMetric' :: HasMetric v => HerMetric v -> HerMetric' v -recipMetric' m | Option (Just m') <- safeRecipMetric' m = m' -recipMetric' _ = singularMetric' - -safeRecipMetric :: HasMetric v => HerMetric' v -> Option (HerMetric v) -safeRecipMetric (HerMetric' Nothing) = empty -safeRecipMetric (HerMetric' (Just (DenseLinear m))) - | isInfinite' detm = empty - | otherwise = return $ matrixMetric minv - where (minv, (detm, _)) = HMat.invlndet m - -safeRecipMetric' :: HasMetric v => HerMetric v -> Option (HerMetric' v) -safeRecipMetric' (HerMetric Nothing) = empty -safeRecipMetric' (HerMetric (Just (DenseLinear m))) - | isInfinite' detm = empty - | otherwise = return $ matrixMetric' minv - where (minv, (detm, _)) = HMat.invlndet m - -isInfinite' :: (Eq a, Num a) => a -> Bool -isInfinite' 0 = False -isInfinite' x = x==x*2 - - - --- | The eigenbasis of a metric, with each eigenvector scaled to the --- square root of the eigenvalue. If the metric is not positive --- definite (i.e. if it has zero eigenvalues), then the 'eigenSpan' --- will contain zero vectors. --- --- This constitutes, in a sense, --- a decomposition of a metric into a set of 'projector'' vectors. If those --- are 'sumV'ed again (use 'projectors's' for this), then the original metric --- is obtained. (This holds even for non-Hilbert/Banach spaces, --- although the concept of eigenbasis and --- “scaled length” doesn't really make sense there.) -eigenSpan :: (HasMetric v, Scalar v ~ ℝ) => HerMetric' v -> [v] -eigenSpan (HerMetric' Nothing) = [] -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 (DenseLinear m))) = map fromPackedVector eigSpan - where (μs,vsm) = HMat.eigSH' m - eigSpan = zipWith (HMat.scale . sqrt) (HMat.toList μs) (HMat.toColumns vsm) - --- | The reciprocal-space counterparts of the nonzero-EV eigenvectors, as can --- be obtained from 'eigenSpan'. The systems of vectors/dual vectors --- behave as orthonormal groups WRT each other, i.e. for each @f@ --- in @'eigenCoSpan' m@ there will be exactly one @v@ in @'eigenSpan' m@ --- such that @f<.>^v ≡ 1@; the other @f<.>^v@ pairings are zero. --- --- Furthermore, @'metric' m f ≡ 1@ for each @f@ in the co-span, which might --- be seen as the actual defining characteristic of these span/co-span systems. -eigenCoSpan :: (HasMetric v, Scalar v ~ ℝ) => HerMetric' v -> [DualSpace v] -eigenCoSpan (HerMetric' Nothing) = [] -eigenCoSpan (HerMetric' (Just (DenseLinear m))) = map fromPackedVector eigSpan - where (μs,vsm) = HMat.eigSH' m - eigSpan = map (uncurry $ HMat.scale . recip . sqrt) - . filter ((>0) . fst) - $ zip (HMat.toList μs) (HMat.toColumns vsm) -eigenCoSpan' :: (HasMetric v, Scalar v ~ ℝ) => HerMetric v -> [v] -eigenCoSpan' (HerMetric Nothing) = [] -eigenCoSpan' (HerMetric (Just (DenseLinear m))) = map fromPackedVector eigSpan - where (μs,vsm) = HMat.eigSH' m - eigSpan = map (uncurry $ HMat.scale . recip . sqrt) - . filter ((>0) . fst) - $ zip (HMat.toList μs) (HMat.toColumns vsm) - - -class HasEigenSystem m where - type EigenVector m :: * - -- | Generalised combination of 'eigenSpan' and 'eigenCoSpan'; this will give a - -- maximum spanning set of vector-covector pairs @(f,v)@ such that @f<.>^v ≡ 1@ - -- and @metric m f ≡ 1@, whereas all @f@ and @v'@ from different tuples - -- are orthogonal. - -- It also yields the /kernel/ of singular metric, spanned by a set of stiefel-manifold - -- points, i.e. vectors of unspecified length that correspond to the eigenvalue 0. - -- - -- You may also consider this as a /factorisation/ of a linear operator - -- @𝐴 : 𝑉 → 𝑉'@ into mappings @𝑅 : 𝑉 → ℝⁿ@ and @𝐿 : ℝⁿ → 𝑉'@ (or, equivalently - -- because ℝⁿ is a Hilbert space, @𝑅' : ℝⁿ → V'@ and @𝐿' : V → ℝⁿ@, which - -- gives you an SVD-style inverse). - eigenSystem :: m -> ( [Stiefel1 (EigenVector m)] - , [(EigenVector m, DualSpace (EigenVector m))] ) - -instance (HasMetric v, Scalar v ~ ℝ) => HasEigenSystem (HerMetric' v) where - type EigenVector (HerMetric' v) = v - eigenSystem (HerMetric' Nothing) = (fmap Stiefel1 completeBasisValues, []) - eigenSystem (HerMetric' (Just (DenseLinear m))) = concat***concat $ unzip eigSpan - where (μs,vsm) = HMat.eigSH' m - eigSpan = zipWith (\μ v - -> if μ>0 - then let sμ = sqrt μ - in ([], [( fromPackedVector $ HMat.scale sμ v - , fromPackedVector $ HMat.scale (recip sμ) v )]) - else ([Stiefel1 $ fromPackedVector v], []) - ) (HMat.toList μs) (HMat.toColumns vsm) - -instance (HasMetric v, Scalar v ~ ℝ) => HasEigenSystem (HerMetric v) where - type EigenVector (HerMetric v) = DualSpace v - eigenSystem (HerMetric Nothing) = (fmap Stiefel1 completeBasisValues, []) - eigenSystem (HerMetric (Just (DenseLinear m))) = concat***concat $ unzip eigSpan - where (μs,vsm) = HMat.eigSH' m - eigSpan = zipWith (\μ v - -> if μ>0 - then let sμ = sqrt μ - in ([], [( fromPackedVector $ HMat.scale sμ v - , fromPackedVector $ HMat.scale (recip sμ) v )]) - else ([Stiefel1 $ fromPackedVector v], []) - ) (HMat.toList μs) (HMat.toColumns vsm) - -instance (HasMetric v, Scalar v ~ ℝ) => HasEigenSystem (HerMetric' v, HerMetric' v) where - type EigenVector (HerMetric' v, HerMetric' v) = v - eigenSystem (n, HerMetric' (Just (DenseLinear m))) | not $ null nSpan - = (++nKernel).concat***concat $ unzip eigSpan - where (μs,vsm) = HMat.eigSH' $ fromv2ℝn HMat.<> m HMat.<> fromℝn2v' - -- m :: v' -> v - eigSpan = zipWith (\μ v - -> if μ>0 - then let sμ = sqrt μ - in ([], [( fromPackedVector $ - fromℝn2v HMat.#> HMat.scale sμ v - , fromPackedVector $ - fromℝn2v' HMat.#> HMat.scale (recip sμ) v ) - ]) - else ([Stiefel1 $ fromPackedVector v], []) - ) (HMat.toList μs) (HMat.toColumns vsm) - fromv2ℝn = HMat.fromRows $ map (asPackedVector . snd) nSpan - fromℝn2v' = HMat.tr fromv2ℝn - fromℝn2v = HMat.fromColumns $ map (asPackedVector . fst) nSpan - (nKernel, nSpan) = eigenSystem n - eigenSystem (_, HerMetric' Nothing) = (fmap Stiefel1 completeBasisValues, []) - -instance (HasMetric v, Scalar v ~ ℝ) => HasEigenSystem (HerMetric v, HerMetric v) where - type EigenVector (HerMetric v, HerMetric v) = DualSpace v - eigenSystem (n, HerMetric (Just (DenseLinear m))) | not $ null nSpan - = (++nKernel).concat***concat $ unzip eigSpan - where (μs,vsm) = HMat.eigSH' $ fromv'2ℝn HMat.<> m HMat.<> fromℝn2v - -- m :: v -> v' - eigSpan = zipWith (\μ v - -> if μ>0 - then let sμ = sqrt μ - in ([], [( fromPackedVector $ - fromℝn2v' HMat.#> HMat.scale sμ v - , fromPackedVector $ - fromℝn2v HMat.#> HMat.scale (recip sμ) v ) - ]) - else ([Stiefel1 $ fromPackedVector v], []) - ) (HMat.toList μs) (HMat.toColumns vsm) - fromv'2ℝn = HMat.fromRows $ map (asPackedVector . snd) nSpan - fromℝn2v = HMat.tr fromv'2ℝn - fromℝn2v' = HMat.fromColumns $ map (asPackedVector . fst) nSpan - (nKernel, nSpan) = eigenSystem n - eigenSystem (_, _) = (fmap Stiefel1 completeBasisValues, []) - - --- | Constraint that a space's scalars need to fulfill so it can be used for 'HerMetric'. -type MetricScalar s = ( SmoothScalar s - , Ord s -- We really rather wouldn't require this... - ) - - -type HasMetric v = (HasMetric' v, HasMetric' (DualSpace v), DualSpace (DualSpace v) ~ v) - - --- | While the main purpose of this class is to express 'HerMetric', it's actually --- all about dual spaces. -class ( FiniteDimensional v, FiniteDimensional (DualSpace v) - , VectorSpace (DualSpace v), HasBasis (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. - -- @v ':-*' 'Scalar' v@. - -- Typically (for all Hilbert- / 'InnerSpace's) this is in turn isomorphic to @v@ - -- itself, which will be rather more efficient (hence the distinction between a - -- vector space and its dual is often neglected or reduced to “column vs row - -- vectors”). - -- Mathematically though, it makes sense to keep the concepts apart, even if ultimately - -- @'DualSpace' v ~ v@ (which needs not /always/ be the case, though!). - type DualSpace v :: * - type DualSpace v = v - - -- | Apply a dual space vector (aka linear functional) to a vector. - (<.>^) :: DualSpace v -> v -> Scalar v - - -- | Interpret a functional as a dual-space vector. Like 'linear', this /assumes/ - -- (completely unchecked) that the supplied function is linear. - functional :: (v -> Scalar v) -> DualSpace v - - -- | While isomorphism between a space and its dual isn't generally canonical, - -- the /double-dual/ space should be canonically isomorphic in pretty much - -- all relevant cases. Indeed, it is recommended that they are the very same type; - -- this condition is enforced by the 'HasMetric' constraint (which is recommended - -- over using 'HasMetric'' itself in signatures). - 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 -ket ^<.> bra = bra <.>^ ket - - -euclideanMetric' :: forall v . (HasMetric v, InnerSpace v) => HerMetric v -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, --- -- as by the Riesz representation theorem. --- -- --- -- Note that usually, Hilbert spaces should just implement @DualSpace v ~ v@, --- -- according to that same correspondence, so 'riesz' is essentially just a more explicit --- -- (and less efficient) way of writing @'id' :: v -> DualSpace v'. --- riesz :: (HasMetric v, InnerSpace v) => v -> DualSpace v --- riesz v = functional (v<.>) --- --- riesz' :: (HasMetric v, InnerSpace v) => DualSpace v -> v --- riesz' f = doubleDual' . functional (f<.>^) - - -instance (MetricScalar k) => HasMetric' (ZeroDim k) where - Origin<.>^Origin = zeroV - functional _ = Origin - doubleDual = id; doubleDual'= id; basisInDual = pure id -instance HasMetric' Double where - (<.>^) = (<.>) - functional f = f 1 - 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 - (<.>^) = (<.>) - functional = fnal - where fnal :: ∀ s n . (SmoothScalar s, KnownNat n) => (s^n -> s) -> s^n - 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; 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 - functional = fnal - where fnal :: ∀ v . HasMetric v => - (FinVecArrRep t v (Scalar v) -> Scalar v) - -> FinVecArrRep t (DualSpace v) (Scalar v) - fnal f = FinVecArrRep . (n HMat.|>) - $ (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, 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_fln m = linear $ \w -> functional $ \v - -> w <.>^lapply m v - - - -metrConst :: forall v. (HasMetric v, v ~ DualSpace v, Num (Scalar v)) - => Scalar v -> HerMetric v -metrConst μ = matrixMetric $ HMat.scale μ (HMat.ident dim) - where (Tagged dim) = dimension :: Tagged v Int - -instance (HasMetric v, v ~ DualSpace v, Num (Scalar v)) => Num (HerMetric v) where - fromInteger = metrConst . fromInteger - (+) = (^+^) - negate = negateV - - -- | This does /not/ work correctly if the metrics don't share an eigenbasis! - 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" - signum = error "signum undefined for HerMetric" - - -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 (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 - fromRational = metrConst . fromRational - recip = metrNumFun recip - -instance (HasMetric v, v ~ Scalar v, v ~ DualSpace v, Floating v) - => Floating (HerMetric v) where - pi = metrConst pi - sqrt = metrNumFun sqrt - exp = metrNumFun exp - log = metrNumFun log - sin = metrNumFun sin - cos = metrNumFun cos - tan = metrNumFun tan - asin = metrNumFun asin - acos = metrNumFun acos - atan = metrNumFun atan - sinh = metrNumFun sinh - cosh = metrNumFun cosh - asinh = metrNumFun asinh - atanh = metrNumFun atanh - acosh = metrNumFun acosh - - - - -normaliseWith :: HasMetric v => HerMetric v -> v -> Option v -normaliseWith m v = case metric m v of - 0 -> empty - μ -> pure (v ^/ μ) - -orthonormalPairsWith :: forall v . HasMetric v => HerMetric v -> [v] -> [(v, DualSpace v)] -orthonormalPairsWith met = mkON - where mkON :: [v] -> [(v, DualSpace v)] -- Generalised Gram-Schmidt process - mkON [] = [] - mkON (v:vs) = let onvs = mkON vs - v' = List.foldl' (\va (vb,pb) -> va ^-^ vb ^* (pb <.>^ va)) v onvs - p' = toDualWith met v' - in case sqrt (p' <.>^ v') of - 0 -> onvs - μ -> (v'^/μ, p'^/μ) : onvs - - - --- | Project a metric on each of the factors of a product space. This works by --- projecting the eigenvectors into both subspaces. -factoriseMetric :: ∀ v w . (HasMetric v, HasMetric w, Scalar v ~ ℝ, Scalar w ~ ℝ) - => HerMetric (v,w) -> (HerMetric v, HerMetric w) -factoriseMetric (HerMetric Nothing) = (HerMetric Nothing, HerMetric Nothing) -factoriseMetric met = (projectors *** projectors) . unzip $ eigenSpan' met - -factoriseMetric' :: ∀ v w . (HasMetric v, HasMetric w, Scalar v ~ ℝ, Scalar w ~ ℝ) - => HerMetric' (v,w) -> (HerMetric' v, HerMetric' w) -factoriseMetric' met = (projector's *** projector's) . unzip $ eigenSpan met - -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 $ 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 $ 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 (Linear ℝ v w) -covariance (HerMetric' Nothing) = pure zeroV -covariance (HerMetric' (Just m)) - | isInfinite' detvnm = empty - | otherwise = return $ snd . m . (id&&&zeroV) . DenseLinear vnorml - where (vnorml, (detvnm, _)) - = HMat.invlndet . getDenseMatrix $ fst . m . (id&&&zeroV) - - -volumeRatio :: HasMetric v => HerMetric v -> HerMetric v -> Scalar v -volumeRatio (HerMetric Nothing) (HerMetric Nothing) = 1 -volumeRatio (HerMetric _) (HerMetric Nothing) = 0 -volumeRatio (HerMetric (Just (DenseLinear m₁))) - (HerMetric (Just (DenseLinear m₂))) - = HMat.det m₂ / HMat.det m₁ -volumeRatio (HerMetric Nothing) (HerMetric _) = 1/0 - -euclideanRelativeMetricVolume :: (HasMetric v, InnerSpace v) => HerMetric v -> Scalar v -euclideanRelativeMetricVolume (HerMetric Nothing) = 1/0 -euclideanRelativeMetricVolume (HerMetric (Just (DenseLinear m))) = recip $ HMat.det m - -tryMetricAsLength :: HerMetric ℝ -> Option ℝ -tryMetricAsLength m = case metricSq m 1 of - o | o > 0 -> pure . sqrt $ recip o - | otherwise -> empty - --- | Unsafe version of 'tryMetricAsLength', only works reliable if the metric --- is strictly positive definite. -metricAsLength :: HerMetric ℝ -> ℝ -metricAsLength m = case metricSq m 1 of - o | o >= 0 -> sqrt $ recip o - | o < 0 -> error "Metric fails to be positive definite!" - | otherwise -> error "Metric yields NaN." - -metricFromLength :: ℝ -> HerMetric ℝ -metricFromLength = projector . recip - -metric'AsLength :: HerMetric' ℝ -> ℝ -metric'AsLength = sqrt . (`metric'`1) - - -spanHilbertSubspace :: ∀ s v w - . (HasMetric v, Scalar v ~ s, IsFreeSpace w, Scalar w ~ s) - => HerMetric v -- ^ Metric to induce the inner product on the Hilbert space. - -> [v] -- ^ @n@ linearly independent vectors, to span the subspace @w@. - -> Option (Embedding (Linear s) w v) - -- ^ An embedding of the @n@-dimensional free subspace @w@ (if the given - -- vectors actually span such a space) into the main space @v@. - -- Regardless of the structure of @v@ (which doesn't need to have an - -- inner product at all!), @w@ will be an 'InnerSpace' with the scalar - -- product defined by the given metric. -spanHilbertSubspace met = emb . orthonormalPairsWith met - where emb onb' - | n'==n = return $ Embedding emb prj . arr identityMatrix - | otherwise = empty - where emb = DenseLinear . HMat.fromColumns $ (asPackedVector . fst) <$> onb - prj = DenseLinear . HMat.fromRows $ (asPackedVector . snd) <$> onb - n' = length onb' - onb = take n onb' - (Tagged n) = theNatN :: Tagged (FreeDimension w) Int - - --- | Same as 'spanHilbertSubspace', but with the standard 'euclideanMetric' (i.e., the --- basis vectors will be orthonormal in the usual sense, in both @w@ and @v@). -spanSubHilbertSpace :: ∀ s v w - . (HasMetric v, InnerSpace v, Scalar v ~ s, IsFreeSpace w, Scalar w ~ s) - => [v] - -> Option (Embedding (Linear s) w v) -spanSubHilbertSpace = spanHilbertSubspace euclideanMetric' - - -orthogonalComplementSpan :: ∀ v . (HasMetric v, Scalar v ~ ℝ) - => [Stiefel1 (DualSpace v)] -> [Stiefel1 v] -orthogonalComplementSpan avoidSpace - = fst ( iterate nextOVect ( [], ( cycle completeBasisValues - , pseudoRieszPair <$> avoidSpace ) ) - !! (d - lav) ) - where Tagged d = dimension :: Tagged v Int - lav = length avoidSpace - nextOVect (result, (v:src, avoid)) - | Option (Just newAvoid@(vfin', _)) <- mkPseudoRieszPair vPurged - = (Stiefel1 vfin':result, (src, newAvoid : avoid)) - where vPurged = foldl (\vp (av', av) -> vp ^-^ av ^* (vp^<.>av')) v avoid - - --- | The /n/-th Stiefel manifold is the space of all possible configurations of --- /n/ orthonormal vectors. In the case /n/ = 1, simply the subspace of normalised --- vectors, i.e. equivalent to the 'UnitSphere'. Even so, it strictly speaking --- requires the containing space to be at least metric (if not Hilbert); we would --- however like to be able to use this concept also in spaces with no inner product, --- therefore we define this space not as normalised vectors, but rather as all --- vectors modulo scaling by positive factors. -newtype Stiefel1 v = Stiefel1 { getStiefel1N :: DualSpace v } - -pseudoRieszPair :: (HasMetric v, Scalar v ~ ℝ) => Stiefel1 v -> (v, DualSpace v) -pseudoRieszPair (Stiefel1 v') - = (fromPackedVector $ HMat.scale (1/HMat.norm_2 vp) vp, v') - where vp = asPackedVector v' - -mkPseudoRieszPair :: (HasMetric v, Scalar v ~ ℝ) => DualSpace v -> Option (v, DualSpace v) -mkPseudoRieszPair v' - | nv' > 0 = pure (fromPackedVector $ HMat.scale (1/nv') vp, v') - | otherwise = empty - where vp = asPackedVector v' - nv' = HMat.norm_2 vp - - - - -instance (HasMetric v, Scalar v ~ Double, Show (DualSpace v)) => Show (HerMetric v) where - showsPrec p m - | null eigSp = showString "zeroV" - | otherwise = showParen (p>5) - . foldr1 ((.) . (.(" ^+^ "++))) - $ ((("projector "++).).showsPrec 10)<$>eigSp - where eigSp = eigenSpan' m - -instance (HasMetric v, Scalar v ~ Double, Show v) => Show (HerMetric' v) where - showsPrec p m - | null eigSp = showString "zeroV" - | otherwise = showParen (p>5) - . foldr1 ((.) . (.(" ^+^ "++))) - $ ((("projector' "++).).showsPrec 10)<$>eigSp - where eigSp = eigenSpan m - - - - - - - - - -linMapAsTensProd :: (FiniteDimensional v, FiniteDimensional w, Scalar v~Scalar w) - => v:-*w -> DualSpace v ⊗ w -linMapAsTensProd f = DensTensProd $ asPackedMatrix f - -linMapFromTensProd :: (FiniteDimensional v, FiniteDimensional w, Scalar v~Scalar w) - => DualSpace v ⊗ w -> v:-*w -linMapFromTensProd (DensTensProd m) = linear $ - asPackedVector >>> HMat.app m >>> fromPackedVector - - - -(⊗) :: (HasMetric v, FiniteDimensional w, Scalar v ~ s, Scalar w ~ s) - => w -> DualSpace v -> Linear s v w -w ⊗ v' = DenseLinear $ HMat.outer wDecomp v'Decomp - where wDecomp = asPackedVector w - v'Decomp = asPackedVector v' - -outerProducts :: (HasMetric v, FiniteDimensional w, Scalar v ~ s, Scalar w ~ s) - => [(w, DualSpace v)] -> Linear s v w -outerProducts [] = zeroV -outerProducts pds = DenseLinear $ HMat.fromColumns (asPackedVector.fst<$>pds) - HMat.<> HMat.fromRows (asPackedVector.snd<$>pds) - -instance ∀ v w s . ( HasMetric v, FiniteDimensional w - , Show (DualSpace v), Show w, Scalar v ~ s, Scalar w ~ s ) - => Show (Linear s v w) where - showsPrec p f = showParen (p>9) $ ("outerProducts "++) - . shows [ (w, v' :: DualSpace v) - | (v,v') <- zip completeBasisValues completeBasisValues - , let w = f $ v ] - - diff --git a/Data/Manifold/Cone.hs b/Data/Manifold/Cone.hs --- a/Data/Manifold/Cone.hs +++ b/Data/Manifold/Cone.hs @@ -36,15 +36,13 @@ import Data.Semigroup import Data.VectorSpace -import Data.LinearMap.HerMetric import Data.Tagged import Data.Manifold.Types.Primitive +import Data.Manifold.Types.Stiefel +import Math.LinearMap.Category import Data.CoNat -import Data.VectorSpace.FiniteDimensional -import qualified Numeric.LinearAlgebra.HMatrix as HMat - import qualified Prelude import qualified Control.Applicative as Hask @@ -57,9 +55,9 @@ -type ConeVecArr m = FinVecArrRep Cℝay (CℝayInterior m) (Scalar (Needle m)) +newtype ConeVecArr m = ConeVecArr {getConeVecArr :: CℝayInterior m} type ConeNeedle m = Needle (ConeVecArr m) -type SConn'dConeVecArr m = FinVecArrRep Cℝay (ℝ, Interior m) ℝ +data SConn'dConeVecArr m = SConn'dConeVecArr ℝ (Interior m) class ( Semimanifold m, Semimanifold (Interior (Interior m)) @@ -84,7 +82,7 @@ -instance (ConeSemimfd m) => Semimanifold (Cℝay m) where +instance ∀ m . (ConeSemimfd m) => Semimanifold (Cℝay m) where type Needle (Cℝay m) = ConeNeedle m type Interior (Cℝay m) = ConeVecArr m fromInterior = fromCℝayInterior @@ -94,6 +92,8 @@ ctp = Tagged ctp' where Tagged ctp' = translateP :: Tagged (ConeVecArr m) (ConeVecArr m -> ConeNeedle m -> ConeVecArr m) + semimanifoldWitness = case semimanifoldWitness :: SemimanifoldWitness (ConeVecArr m) of + SemimanifoldWitness -> SemimanifoldWitness instance (ConeSemimfd m) => Semimanifold (CD¹ m) where type Needle (CD¹ m) = ConeNeedle m @@ -105,137 +105,17 @@ ctp = Tagged ctp' where Tagged ctp' = translateP :: Tagged (ConeVecArr m) (ConeVecArr m -> ConeNeedle m -> ConeVecArr m) + semimanifoldWitness = case semimanifoldWitness :: SemimanifoldWitness (ConeVecArr m) of + SemimanifoldWitness -> SemimanifoldWitness -instance (ConeSemimfd m, SmoothScalar (Scalar (Needle m))) => PseudoAffine (Cℝay m) where - p.-~.i = (.-~.i) =<< toInterior p -instance (ConeSemimfd m, SmoothScalar (Scalar (Needle m))) => PseudoAffine (CD¹ m) where - p.-~.i = (.-~.i) =<< toInterior p -instance ConeSemimfd (ZeroDim ℝ) where - type CℝayInterior (ZeroDim ℝ) = ℝ - fromCℝayInterior (FinVecArrRep qb) | HMat.size qb == 0 = Cℝay 1 Origin - | x <- qb HMat.! 0 = Cℝay (bijectℝtoℝplus x) Origin - toCℝayInterior (Cℝay 0 Origin) = empty - toCℝayInterior (Cℝay y Origin) = pure . FinVecArrRep $ 1 HMat.|>[bijectℝplustoℝ y] -instance ConeSemimfd ℝ where - type CℝayInterior ℝ = ℝ² - fromCℝayInterior (FinVecArrRep qb) = Cℝay (q'+b') (q'-b') - where [q', b'] = HMat.toList $ HMat.cmap ((/2) . bijectℝtoℝplus) qb - toCℝayInterior (Cℝay 0 _) = empty - toCℝayInterior (Cℝay h x) = pure . FinVecArrRep - . HMat.cmap bijectℝplustoℝ $ HMat.fromList [h+x, h-x] - fromCD¹Interior (FinVecArrRep qb) = CD¹ (bijectℝplustoIntv $ q'+b') (q'-b') - where [q', b'] = HMat.toList $ HMat.cmap ((/2) . bijectℝtoℝplus) qb - toCD¹Interior (CD¹ h x) = pure . FinVecArrRep - . HMat.cmap bijectℝplustoℝ $ HMat.fromList [h'+x, h'-x] - where h' = bijectIntvtoℝplus h -instance ConeSemimfd S⁰ where - type CℝayInterior S⁰ = ℝ - fromCℝayInterior xa | x>0 = Cℝay x PositiveHalfSphere - | otherwise = Cℝay (-x) NegativeHalfSphere - where x = getFinVecArrRep xa HMat.! 0 - toCℝayInterior (Cℝay x PositiveHalfSphere) = return . FinVecArrRep $ HMat.scalar x - toCℝayInterior (Cℝay x NegativeHalfSphere) = return . FinVecArrRep . HMat.scalar $ -x - fromCD¹Interior xa | x>0 = CD¹ (bijectℝtoIntv x) PositiveHalfSphere - | otherwise = CD¹ (-bijectℝtoIntv x) NegativeHalfSphere - where x = getFinVecArrRep xa HMat.! 0 - toCD¹Interior (CD¹ 1 _) = empty - toCD¹Interior (CD¹ x PositiveHalfSphere) - = return . FinVecArrRep . HMat.scalar $ bijectIntvtoℝ x - toCD¹Interior (CD¹ x NegativeHalfSphere) - = return . FinVecArrRep . HMat.scalar $ -bijectℝtoIntv x - - -instance ConeSemimfd S¹ where - type CℝayInterior S¹ = ℝ² - fromCℝayInterior (FinVecArrRep xy) = Cℝay r (S¹ $ atan2 y x) - where r = HMat.norm_2 xy - [x,y] = HMat.toList xy - toCℝayInterior (Cℝay r (S¹ φ)) = return . FinVecArrRep - . HMat.scale r $ HMat.fromList [cos φ, sin φ] - fromCD¹Interior (FinVecArrRep xy) = CD¹ (bijectℝtoIntv r) (S¹ $ atan2 y x) - where r = HMat.norm_2 xy - [x,y] = HMat.toList xy - toCD¹Interior (CD¹ 1 _) = empty - toCD¹Interior (CD¹ r (S¹ φ)) = return . FinVecArrRep - . HMat.scale r' $ HMat.fromList [cos φ, sin φ] - where r' = bijectIntvtoℝ r - - -instance ConeSemimfd S² where - type CℝayInterior S² = ℝ³ - fromCℝayInterior (FinVecArrRep xyz) = Cℝay r (S² (acos $ z/r) (atan2 y x)) - where r = HMat.norm_2 xyz - [x,y,z] = HMat.toList xyz - toCℝayInterior (Cℝay r (S² ϑ φ)) = return . FinVecArrRep - . HMat.scale r $ HMat.fromList [w*x₀, w*y₀, z₀] - where x₀ = cos φ; y₀ = sin φ; z₀ = cos ϑ; w = sin ϑ - --- | Products of simply connected spaces. -instance ( PseudoAffine x, PseudoAffine y - , WithField ℝ HilbertSpace (Interior x), WithField ℝ HilbertSpace (Interior y) - , LinearManifold (FinVecArrRep Cℝay (ℝ, (Interior x, Interior y)) ℝ) - ) => ConeSemimfd (x,y) where - type CℝayInterior (x,y) = (ℝ, (Interior x, Interior y)) - fromCℝayInterior = simplyCncted_fromCℝayInterior - toCℝayInterior = simplyCncted_toCℝayInterior -instance ( KnownNat n ) => ConeSemimfd (ℝ^n) where - type CℝayInterior (ℝ^n) = (ℝ, ℝ^n) - fromCℝayInterior = simplyCncted_fromCℝayInterior - toCℝayInterior = simplyCncted_toCℝayInterior -instance ( HilbertSpace (FinVecArrRep t v ℝ) ) => ConeSemimfd (FinVecArrRep t v ℝ) where - type CℝayInterior (FinVecArrRep t v ℝ) = (ℝ, FinVecArrRep t v ℝ) - fromCℝayInterior = simplyCncted_fromCℝayInterior - toCℝayInterior = simplyCncted_toCℝayInterior - - - -instance ( WithField ℝ ConeSemimfd x, PseudoAffine (Cℝay x) - , HilbertSpace (CℝayInterior x) - , HilbertSpace (FinVecArrRep Cℝay (CℝayInterior x) ℝ) - ) => ConeSemimfd (CD¹ x) where - type CℝayInterior (CD¹ x) = (ℝ, ConeVecArr x) - fromCℝayInterior i = Cℝay h (embCℝayToCD¹ o) - where (Cℝay h o) = simplyCncted_fromCℝayInterior i - toCℝayInterior (Cℝay _ (CD¹ 1 _)) = empty - toCℝayInterior (Cℝay h p) = simplyCncted_toCℝayInterior $ Cℝay h (projCD¹ToCℝay p) - - -instance ( WithField ℝ ConeSemimfd x, PseudoAffine (Cℝay x) - , HilbertSpace (CℝayInterior x) - , HilbertSpace (FinVecArrRep Cℝay (CℝayInterior x) ℝ) - ) => ConeSemimfd (Cℝay x) where - type CℝayInterior (Cℝay x) = (ℝ, ConeVecArr x) - fromCℝayInterior = simplyCncted_fromCℝayInterior - toCℝayInterior = simplyCncted_toCℝayInterior - - -simplyCncted_fromCℝayInterior :: (PseudoAffine x, WithField ℝ HilbertSpace (Interior x)) - => SConn'dConeVecArr x -> Cℝay x -simplyCncted_fromCℝayInterior (FinVecArrRep ri) = Cℝay h . fromInterior . fromPackedVector - $ subtract (h/n) `Arr.map` Arr.tail cmps - where h = Arr.sum cmps - cmps = bijectℝtoℝplus `HMat.cmap` ri - n = fromIntegral $ Arr.length cmps - -simplyCncted_toCℝayInterior :: (PseudoAffine x, WithField ℝ HilbertSpace (Interior x)) - => Cℝay x -> Option (SConn'dConeVecArr x) -simplyCncted_toCℝayInterior (Cℝay h v) | h/=0, Option (Just vi) <- toInterior v - = let cmps'' = asPackedVector vi - cmps' = (+ h/n) `HMat.cmap` cmps'' - cmps = (h - Arr.sum cmps') `Arr.cons` cmps - n = fromIntegral $ Arr.length cmps - in return $ FinVecArrRep (bijectℝplustoℝ `Arr.map` cmps) -simplyCncted_toCℝayInterior (Cℝay _ _) = empty - - -- Some essential homeomorphisms bijectℝtoℝplus , bijectℝplustoℝ , bijectIntvtoℝplus, bijectℝplustoIntv @@ -265,15 +145,15 @@ stiefel1Project :: LinearManifold v => - DualSpace v -- ^ Must be nonzero. + DualVector v -- ^ Must be nonzero. -> Stiefel1 v stiefel1Project = Stiefel1 -stiefel1Embed :: HilbertSpace v => Stiefel1 v -> v +stiefel1Embed :: (HilbertSpace v, RealFloat (Scalar v)) => Stiefel1 v -> v stiefel1Embed (Stiefel1 n) = normalized n -class (PseudoAffine v, InnerSpace v, NaturallyEmbedded (UnitSphere v) (DualSpace v)) +class (PseudoAffine v, InnerSpace v, NaturallyEmbedded (UnitSphere v) (DualVector v)) => HasUnitSphere v where type UnitSphere v :: * stiefel :: UnitSphere v -> Stiefel1 v @@ -282,13 +162,10 @@ unstiefel = coEmbed . getStiefel1N instance HasUnitSphere ℝ where type UnitSphere ℝ = S⁰ -instance HasUnitSphere (FinVecArrRep t ℝ ℝ) where type UnitSphere (FinVecArrRep t ℝ ℝ) = S⁰ instance HasUnitSphere ℝ² where type UnitSphere ℝ² = S¹ -instance HasUnitSphere (FinVecArrRep t ℝ² ℝ) where type UnitSphere (FinVecArrRep t ℝ² ℝ) = S¹ instance HasUnitSphere ℝ³ where type UnitSphere ℝ³ = S² -instance HasUnitSphere (FinVecArrRep t ℝ³ ℝ) where type UnitSphere (FinVecArrRep t ℝ³ ℝ) = S² diff --git a/Data/Manifold/DifferentialEquation.hs b/Data/Manifold/DifferentialEquation.hs --- a/Data/Manifold/DifferentialEquation.hs +++ b/Data/Manifold/DifferentialEquation.hs @@ -49,8 +49,7 @@ import Data.Semigroup import Data.VectorSpace -import Data.LinearMap.HerMetric -import Data.LinearMap.Category +import Math.LinearMap.Category import Data.AffineSpace import Data.Basis @@ -61,7 +60,6 @@ import Data.Manifold.TreeCover import Data.Manifold.Web -import qualified Numeric.LinearAlgebra.HMatrix as HMat import qualified Data.List as List import qualified Prelude as Hask hiding(foldl, sum, sequence) @@ -78,13 +76,15 @@ import Data.Traversable.Constrained (Traversable, traverse) -constLinearDEqn :: (WithField ℝ LinearManifold x, WithField ℝ LinearManifold y) - => Linear ℝ (DualSpace y) (Linear ℝ y x) -> DifferentialEqn x y +constLinearDEqn :: ( WithField ℝ LinearManifold x, SimpleSpace x + , WithField ℝ LinearManifold y, SimpleSpace y ) + => (DualVector y +> (y +> x)) -> DifferentialEqn x y constLinearDEqn bwt = factoriseShade - >>> \(_x, Shade y δy) -> let j = bwt'm HMat.<\> (asPackedVector y) - δj = bwt' `transformMetric` recipMetric δy - in Shade' (fromPackedVector j) δj - where bwt'@(DenseLinear bwt'm) = adjoint bwt + >>> \(_x, Shade y δy) -> let j = bwt'inv y + δj = bwt' `transformNorm` dualNorm δy + in Shade' j δj + where bwt' = adjoint $ bwt + bwt'inv = (bwt'\$) -- | A function that variates, relatively speaking, most strongly @@ -93,25 +93,33 @@ -- approaches 0. -- -- The idea is that if you consider the ratio of two function values, --- it will be close to 1 if both arguments on the same side of 1, --- even if their ratio is large. +-- it will be close to 1 if either both arguments are much smaller or both +-- much larger than 1, even if the ratio of these arguments is large. -- Only if both arguments are close to 1, or lie on opposite sides -- of it, will the ratio of the function values will be significant. goalSensitive :: ℝ -> ℝ goalSensitive η = 0.3 + sqrt (η * (1 + η/(1+η)) / (3 + η)) -euclideanVolGoal :: WithField ℝ EuclidSpace y => ℝ -> x -> Shade' y -> ℝ +euclideanVolGoal :: (WithField ℝ EuclidSpace y, SimpleSpace (Needle y)) + => ℝ -> x -> Shade' y -> ℝ euclideanVolGoal vTgt _ (Shade' _ shy) = goalSensitive η where η = euclideanRelativeMetricVolume shy / vTgt -maxDeviationsGoal :: WithField ℝ EuclidSpace y => [Needle y] -> x -> Shade' y -> ℝ -maxDeviationsGoal = uncertaintyGoal . projector's +euclideanRelativeMetricVolume :: (SimpleSpace y, HilbertSpace y) => Norm y -> Scalar y +euclideanRelativeMetricVolume (Norm m) = recip . roughDet . arr $ ue . m + where Norm ue = euclideanNorm -uncertaintyGoal :: WithField ℝ EuclidSpace y => Metric' y -> x -> Shade' y -> ℝ +maxDeviationsGoal :: (WithField ℝ EuclidSpace y, SimpleSpace (Needle y)) + => [Needle y] -> x -> Shade' y -> ℝ +maxDeviationsGoal = uncertaintyGoal . spanNorm + +uncertaintyGoal :: (WithField ℝ EuclidSpace y, SimpleSpace (Needle y)) + => Metric' y -> x -> Shade' y -> ℝ uncertaintyGoal = uncertaintyGoal' . const -uncertaintyGoal' :: WithField ℝ EuclidSpace y => (x -> Metric' y) -> x -> Shade' y -> ℝ +uncertaintyGoal' :: (WithField ℝ EuclidSpace y, SimpleSpace (Needle y)) + => (x -> Metric' y) -> x -> Shade' y -> ℝ uncertaintyGoal' f x (Shade' _ shy) - = List.sum [goalSensitive $ 1 / metricSq' m q | q <- shySpan] - where shySpan = eigenSpan' shy + = List.sum [goalSensitive $ 1 / normSq m q | q <- shySpan] + where shySpan = normSpanningSystem shy m = f x diff --git a/Data/Manifold/Griddable.hs b/Data/Manifold/Griddable.hs --- a/Data/Manifold/Griddable.hs +++ b/Data/Manifold/Griddable.hs @@ -38,7 +38,7 @@ import Data.List hiding (filter, all, elem, sum) import Data.Maybe -import Data.LinearMap.HerMetric +import Math.LinearMap.Category import Data.Manifold.Types import Data.Manifold.Types.Primitive ((^), (^.)) @@ -97,7 +97,7 @@ where l = c - expa r = c + expa - expa = metric'AsLength expa' + expa = normalLength expa' (Just ax) = find ((>=n) . axisGrLength) $ [ let qe = 10^^lqe' * nb @@ -110,20 +110,21 @@ | n < 0 = floor $ lg (-n) -instance (Griddable m a, Griddable n a) => Griddable (m,n) a where +instance ( SimpleSpace (Needle m), SimpleSpace (Needle n), SimpleSpace (Needle a) + , Griddable m a, Griddable n a ) => Griddable (m,n) a where data GriddingParameters (m,n) a = PairGriddingParameters { fstGriddingParams :: GriddingParameters m a , sndGriddingParams :: GriddingParameters n a } mkGridding (PairGriddingParameters p₁ p₂) n (Shade (c₁,c₂) e₁e₂) = ( gshmap ( uncurry fullShade . ( (,c₂).(^.shadeCtr) - &&& (`productMetric'`e₂).(^.shadeExpanse)) ) + &&& (`sumSubspaceNorms`e₂).(^.shadeExpanse)) ) <$> g₁s ) ++ ( gshmap ( uncurry fullShade . ( (c₁,).(^.shadeCtr) - &&& ( productMetric' e₁).(^.shadeExpanse)) ) + &&& ( sumSubspaceNorms e₁).(^.shadeExpanse)) ) <$> g₂s ) where g₁s = mkGridding p₁ n $ fullShade c₁ e₁ g₂s = mkGridding p₂ n $ fullShade c₂ e₂ - (e₁,e₂) = factoriseMetric' e₁e₂ + (e₁,e₂) = summandSpaceNorms e₁e₂ prettyFloatShow :: Int -> Double -> String prettyFloatShow _ 0 = "0" @@ -144,11 +145,11 @@ shade2Intvl :: Shade ℝ -> Interval shade2Intvl sh = Interval l r where c = sh ^. shadeCtr - expa = metric'AsLength $ sh ^. shadeExpanse + expa = normalLength $ sh ^. shadeExpanse l = c - expa; r = c + expa intvl2Shade :: Interval -> Shade ℝ -intvl2Shade (Interval l r) = fullShade c (projector' expa) +intvl2Shade (Interval l r) = fullShade c (spanNorm [expa]) where c = (l+r) / 2 expa = (r-l) / 2 diff --git a/Data/Manifold/PseudoAffine.hs b/Data/Manifold/PseudoAffine.hs --- a/Data/Manifold/PseudoAffine.hs +++ b/Data/Manifold/PseudoAffine.hs @@ -37,6 +37,7 @@ {-# LANGUAGE RankNTypes #-} {-# LANGUAGE TupleSections #-} {-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE DefaultSignatures #-} {-# LANGUAGE PatternGuards #-} {-# LANGUAGE TypeOperators #-} {-# LANGUAGE UnicodeSyntax #-} @@ -56,17 +57,18 @@ , Metric, Metric', euclideanMetric , RieMetric, RieMetric' -- ** Constraints + , SemimanifoldWitness(..) , RealDimension, AffineManifold , LinearManifold , WithField - , HilbertSpace + , HilbertManifold , EuclidSpace , LocallyScalable -- ** Local functions , LocalLinear, LocalAffine -- * Misc - , alerpB, palerp, palerpB, LocallyCoercible(..) - , ImpliesMetric(..) + , alerpB, palerp, palerpB, LocallyCoercible(..), CanonicalDiffeomorphism(..) + , ImpliesMetric(..), coerceMetric, coerceMetric' ) where @@ -76,16 +78,20 @@ import Data.Fixed import Data.VectorSpace +import Linear.V0 +import Linear.V1 +import Linear.V2 +import Linear.V3 +import Linear.V4 +import qualified Linear.Affine as LinAff import Data.Embedding import Data.LinearMap -import Data.LinearMap.HerMetric -import Data.LinearMap.Category +import Math.LinearMap.Category import Data.AffineSpace import Data.Tagged import Data.Manifold.Types.Primitive import Data.CoNat -import Data.VectorSpace.FiniteDimensional import qualified Prelude import qualified Control.Applicative as Hask @@ -99,11 +105,20 @@ +-- | This is the reified form of the property that the interior of a semimanifold +-- is a manifold. +data SemimanifoldWitness x where + SemimanifoldWitness :: + ( Semimanifold (Interior x), Semimanifold (Needle x) + , Interior (Interior x) ~ Interior x, Needle (Interior x) ~ Needle x + , Interior (Needle x) ~ Needle x ) + => SemimanifoldWitness x + + infix 6 .-~. infixl 6 .+~^, .-~^ -class ( AdditiveGroup (Needle x), Interior (Interior x) ~ Interior x ) - => Semimanifold x where +class AdditiveGroup (Needle x) => Semimanifold x where {-# MINIMAL ((.+~^) | fromInterior), toInterior, translateP #-} -- | The space of “natural” ways starting from some reference point -- and going to some particular target point. Hence, @@ -160,6 +175,14 @@ -- instance). (.-~^) :: Interior x -> Needle x -> x p .-~^ v = p .+~^ negateV v + + semimanifoldWitness :: SemimanifoldWitness x + default semimanifoldWitness :: + ( Semimanifold (Interior x), Semimanifold (Needle x) + , Interior (Interior x) ~ Interior x, Needle (Interior x) ~ Needle x + , Interior (Needle x) ~ Needle x ) + => SemimanifoldWitness x + semimanifoldWitness = SemimanifoldWitness -- | This is the class underlying manifolds. ('Manifold' only precludes boundaries @@ -228,42 +251,66 @@ -- /canonically isomorphic/ tangent spaces, so that -- @'fromPackedVector' . 'asPackedVector' :: 'Needle' x -> 'Needle' ξ@ -- defines a meaningful “representational identity“ between these spaces. -class (PseudoAffine x, PseudoAffine ξ, Scalar (Needle x) ~ Scalar (Needle ξ)) +class ( Semimanifold x, Semimanifold ξ, LSpace (Needle x), LSpace (Needle ξ) + , Scalar (Needle x) ~ Scalar (Needle ξ) ) => LocallyCoercible x ξ where - -- | Must be compatible with the canonical isomorphism on the tangent spaces, - -- i.e. + -- | Must be compatible with the isomorphism on the tangent spaces, i.e. -- @ - -- locallyTrivialDiffeomorphism (p .+~^ 'fromPackedVector' v) - -- ≡ locallyTrivialDiffeomorphism p .+~^ 'fromPackedVector' v + -- locallyTrivialDiffeomorphism (p .+~^ v) + -- ≡ locallyTrivialDiffeomorphism p .+~^ 'coerceNeedle' v -- @ locallyTrivialDiffeomorphism :: x -> ξ - -instance LocallyCoercible ℝ ℝ where locallyTrivialDiffeomorphism = id -instance LocallyCoercible (ℝ,ℝ) (ℝ,ℝ) where locallyTrivialDiffeomorphism = id -instance LocallyCoercible (ℝ,(ℝ,ℝ)) (ℝ,(ℝ,ℝ)) where locallyTrivialDiffeomorphism = id -instance LocallyCoercible ((ℝ,ℝ),ℝ) ((ℝ,ℝ),ℝ) where locallyTrivialDiffeomorphism = id + coerceNeedle :: Functor p (->) (->) => p (x,ξ) -> (Needle x -+> Needle ξ) + coerceNeedle' :: Functor p (->) (->) => p (x,ξ) -> (Needle' x -+> Needle' ξ) + oppositeLocalCoercion :: CanonicalDiffeomorphism ξ x + default oppositeLocalCoercion :: LocallyCoercible ξ x => CanonicalDiffeomorphism ξ x + oppositeLocalCoercion = CanonicalDiffeomorphism + interiorLocalCoercion :: Functor p (->) (->) + => p (x,ξ) -> CanonicalDiffeomorphism (Interior x) (Interior ξ) + default interiorLocalCoercion :: LocallyCoercible (Interior x) (Interior ξ) + => p (x,ξ) -> CanonicalDiffeomorphism (Interior x) (Interior ξ) + interiorLocalCoercion _ = CanonicalDiffeomorphism +#define identityCoercion(c,t) \ +instance (c) => LocallyCoercible (t) (t) where { \ + locallyTrivialDiffeomorphism = id; \ + coerceNeedle _ = id; \ + coerceNeedle' _ = id; \ + oppositeLocalCoercion = CanonicalDiffeomorphism; \ + interiorLocalCoercion _ = CanonicalDiffeomorphism } +identityCoercion(NumberManifold s, ZeroDim s) +identityCoercion(NumberManifold s, V0 s) +identityCoercion((), ℝ) +identityCoercion(NumberManifold s, V1 s) +identityCoercion((), (ℝ,ℝ)) +identityCoercion(NumberManifold s, V2 s) +identityCoercion((), (ℝ,(ℝ,ℝ))) +identityCoercion((), ((ℝ,ℝ),ℝ)) +identityCoercion(NumberManifold s, V3 s) +identityCoercion(NumberManifold s, V4 s) +data CanonicalDiffeomorphism a b where + CanonicalDiffeomorphism :: LocallyCoercible a b => CanonicalDiffeomorphism a b + + type LocallyScalable s x = ( PseudoAffine x - , HasMetric (Needle x) - , s ~ Scalar (Needle x) ) + , LSpace (Needle x) + , s ~ Scalar (Needle x) + , Num''' s ) -type LocalLinear x y = Linear (Scalar (Needle x)) (Needle x) (Needle y) +type LocalLinear x y = LinearMap (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. --- --- (Actually, 'LinearManifold' is stronger than 'VectorSpace' at the moment, since --- 'HasMetric' requires 'FiniteDimensional'. This might be lifted in the future.) -type LinearManifold x = ( AffineManifold x, Needle x ~ x, HasMetric x ) +type LinearManifold x = ( AffineManifold x, Needle x ~ x, LSpace x ) type LinearManifold' x = ( PseudoAffine x, AffineSpace x, Diff x ~ x - , Interior x ~ x, Needle x ~ x, HasMetric x ) + , Interior x ~ x, Needle x ~ x, LSpace x ) -- | Require some constraint on a manifold, and also fix the type of the manifold's --- underlying field. For example, @WithField ℝ 'HilbertSpace' v@ constrains +-- underlying field. For example, @WithField ℝ 'HilbertManifold' v@ constrains -- @v@ to be a real (i.e., 'Double'-) Hilbert space. -- Note that for this to compile, you will in -- general need the @-XLiberalTypeSynonyms@ extension (except if the constraint @@ -272,9 +319,7 @@ type WithField s c x = ( c x, s ~ Scalar (Needle x) ) -- | The 'RealFloat' class plus manifold constraints. -type RealDimension r = ( PseudoAffine r, Interior r ~ r, Needle r ~ r - , HasMetric r, DualSpace r ~ r, Scalar r ~ r - , RealFloat r, r ~ ℝ) +type RealDimension r = ( PseudoAffine r, Interior r ~ r, Needle r ~ r, r ~ ℝ) -- | The 'AffineSpace' class plus manifold constraints. type AffineManifold m = ( PseudoAffine m, Interior m ~ m, AffineSpace m @@ -286,27 +331,36 @@ -- (Stricly speaking, that doesn't have much to do with the completeness criterion; -- but since 'Manifold's are at the moment confined to finite dimension, they are in -- fact (trivially) complete.) -type HilbertSpace x = ( LinearManifold x, InnerSpace x - , Interior x ~ x, Needle x ~ x, DualSpace x ~ x - , Floating (Scalar x) ) +type HilbertManifold x = ( LinearManifold x, InnerSpace x + , Interior x ~ x, Needle x ~ x, DualVector x ~ x + , Floating (Scalar x) ) -- | An euclidean space is a real affine space whose tangent space is a Hilbert space. type EuclidSpace x = ( AffineManifold x, InnerSpace (Diff x) - , DualSpace (Diff x) ~ Diff x, Floating (Scalar (Diff x)) ) + , DualVector (Diff x) ~ Diff x, Floating (Scalar (Diff x)) ) +type NumberManifold n = ( Num''' n, Manifold n, Interior n ~ n, Needle n ~ n + , LSpace n, DualVector n ~ n, Scalar n ~ n ) + euclideanMetric :: EuclidSpace x => proxy x -> Metric x -euclideanMetric _ = euclideanMetric' +euclideanMetric _ = euclideanNorm -- | A co-needle can be understood as a “paper stack”, with which you can measure -- the length that a needle reaches in a given direction by counting the number -- of holes punched through them. -type Needle' x = DualSpace (Needle x) +type Needle' x = DualVector (Needle x) --- | The word “metric” is used in the sense as in general relativity. Cf. 'HerMetric'. -type Metric x = HerMetric (Needle x) -type Metric' x = HerMetric' (Needle x) +-- | The word “metric” is used in the sense as in general relativity. +-- Actually this is just the type of scalar products on the tangent space. +-- The actual metric is the function @x -> x -> Scalar (Needle x)@ defined by +-- +-- @ +-- \\p q -> m '|$|' (p.-~!q) +-- @ +type Metric x = Norm (Needle x) +type Metric' x = Variance (Needle x) -- | A Riemannian metric assigns each point on a manifold a scalar product on the tangent space. -- Note that this association is /not/ continuous, because the charts/tangent spaces in the bundle @@ -315,6 +369,23 @@ type RieMetric x = x -> Metric x type RieMetric' x = x -> Metric' x + +coerceMetric :: ∀ x ξ . (LocallyCoercible x ξ, LSpace (Needle ξ)) + => RieMetric ξ -> RieMetric x +coerceMetric m x = case m $ locallyTrivialDiffeomorphism x of + Norm sc -> Norm $ bw . sc . fw + where fw = coerceNeedle ([]::[(x,ξ)]) + bw = case oppositeLocalCoercion :: CanonicalDiffeomorphism ξ x of + CanonicalDiffeomorphism -> coerceNeedle' ([]::[(ξ,x)]) +coerceMetric' :: ∀ x ξ . (LocallyCoercible x ξ, LSpace (Needle ξ)) + => RieMetric' ξ -> RieMetric' x +coerceMetric' m x = case m $ locallyTrivialDiffeomorphism x of + Norm sc -> Norm $ bw . sc . fw + where fw = coerceNeedle' ([]::[(x,ξ)]) + bw = case oppositeLocalCoercion :: CanonicalDiffeomorphism ξ x of + CanonicalDiffeomorphism -> coerceNeedle ([]::[(ξ,x)]) + + -- | Interpolate between points, approximately linearly. For -- points that aren't close neighbours (i.e. lie in an almost -- flat region), the pathway is basically undefined – save for @@ -347,37 +418,71 @@ hugeℝVal :: ℝ hugeℝVal = 1e+100 -#define deriveAffine(t) \ -instance Semimanifold (t) where { \ - type Needle (t) = Diff (t); \ - fromInterior = id; \ - toInterior = pure; \ - translateP = Tagged (.+^); \ - (.+~^) = (.+^) }; \ -instance PseudoAffine (t) where { \ +#define deriveAffine(c,t) \ +instance (c) => Semimanifold (t) where { \ + type Needle (t) = Diff (t); \ + fromInterior = id; \ + toInterior = pure; \ + translateP = Tagged (.+^); \ + (.+~^) = (.+^) }; \ +instance (c) => PseudoAffine (t) where { \ a.-~.b = pure (a.-.b); } -deriveAffine(Double) -deriveAffine(Rational) +deriveAffine((),Double) +deriveAffine((),Rational) +deriveAffine(NumberManifold s, V1 s) +deriveAffine(NumberManifold s, V2 s) +deriveAffine(NumberManifold s, V3 s) +deriveAffine(NumberManifold s, V4 s) -instance SmoothScalar s => Semimanifold (FinVecArrRep t b s) where - type Needle (FinVecArrRep t b s) = FinVecArrRep t b s - type Interior (FinVecArrRep t b s) = FinVecArrRep t b s - fromInterior = id - toInterior = pure - translateP = Tagged (.+^) - (.+~^) = (.+^) -instance SmoothScalar s => PseudoAffine (FinVecArrRep t b s) where - a.-~.b = pure (a.-.b) -instance SmoothScalar s => LocallyCoercible (FinVecArrRep t b s) (FinVecArrRep t b s) where - locallyTrivialDiffeomorphism = id -instance (SmoothScalar s, LinearManifold b, Scalar b ~ s) - => LocallyCoercible (FinVecArrRep t b s) b where - locallyTrivialDiffeomorphism = (concreteArrRep$<-$) -instance (SmoothScalar s, LinearManifold b, Scalar b ~ s) - => LocallyCoercible b (FinVecArrRep t b s) where - locallyTrivialDiffeomorphism = (concreteArrRep$->$) - +instance (NumberManifold s) => LocallyCoercible (ZeroDim s) (V0 s) where + locallyTrivialDiffeomorphism Origin = V0 + coerceNeedle _ = LinearFunction $ \Origin -> V0 + coerceNeedle' _ = LinearFunction $ \Origin -> V0 +instance (NumberManifold s) => LocallyCoercible (V0 s) (ZeroDim s) where + locallyTrivialDiffeomorphism V0 = Origin + coerceNeedle _ = LinearFunction $ \V0 -> Origin + coerceNeedle' _ = LinearFunction $ \V0 -> Origin +instance LocallyCoercible ℝ (V1 ℝ) where + locallyTrivialDiffeomorphism = V1 + coerceNeedle _ = LinearFunction V1 + coerceNeedle' _ = LinearFunction V1 +instance LocallyCoercible (V1 ℝ) ℝ where + locallyTrivialDiffeomorphism (V1 n) = n + coerceNeedle _ = LinearFunction $ \(V1 n) -> n + coerceNeedle' _ = LinearFunction $ \(V1 n) -> n +instance LocallyCoercible (ℝ,ℝ) (V2 ℝ) where + locallyTrivialDiffeomorphism = uncurry V2 + coerceNeedle _ = LinearFunction $ uncurry V2 + coerceNeedle' _ = LinearFunction $ uncurry V2 +instance LocallyCoercible (V2 ℝ) (ℝ,ℝ) where + locallyTrivialDiffeomorphism (V2 x y) = (x,y) + coerceNeedle _ = LinearFunction $ \(V2 x y) -> (x,y) + coerceNeedle' _ = LinearFunction $ \(V2 x y) -> (x,y) +instance LocallyCoercible ((ℝ,ℝ),ℝ) (V3 ℝ) where + locallyTrivialDiffeomorphism ((x,y),z) = V3 x y z + coerceNeedle _ = LinearFunction $ \((x,y),z) -> V3 x y z + coerceNeedle' _ = LinearFunction $ \((x,y),z) -> V3 x y z +instance LocallyCoercible (ℝ,(ℝ,ℝ)) (V3 ℝ) where + locallyTrivialDiffeomorphism (x,(y,z)) = V3 x y z + coerceNeedle _ = LinearFunction $ \(x,(y,z)) -> V3 x y z + coerceNeedle' _ = LinearFunction $ \(x,(y,z)) -> V3 x y z +instance LocallyCoercible (V3 ℝ) ((ℝ,ℝ),ℝ) where + locallyTrivialDiffeomorphism (V3 x y z) = ((x,y),z) + coerceNeedle _ = LinearFunction $ \(V3 x y z) -> ((x,y),z) + coerceNeedle' _ = LinearFunction $ \(V3 x y z) -> ((x,y),z) +instance LocallyCoercible (V3 ℝ) (ℝ,(ℝ,ℝ)) where + locallyTrivialDiffeomorphism (V3 x y z) = (x,(y,z)) + coerceNeedle _ = LinearFunction $ \(V3 x y z) -> (x,(y,z)) + coerceNeedle' _ = LinearFunction $ \(V3 x y z) -> (x,(y,z)) +instance LocallyCoercible ((ℝ,ℝ),(ℝ,ℝ)) (V4 ℝ) where + locallyTrivialDiffeomorphism ((x,y),(z,w)) = V4 x y z w + coerceNeedle _ = LinearFunction $ \((x,y),(z,w)) -> V4 x y z w + coerceNeedle' _ = LinearFunction $ \((x,y),(z,w)) -> V4 x y z w +instance LocallyCoercible (V4 ℝ) ((ℝ,ℝ),(ℝ,ℝ)) where + locallyTrivialDiffeomorphism (V4 x y z w) = ((x,y),(z,w)) + coerceNeedle _ = LinearFunction $ \(V4 x y z w) -> ((x,y),(z,w)) + coerceNeedle' _ = LinearFunction $ \(V4 x y z w) -> ((x,y),(z,w)) instance Semimanifold (ZeroDim k) where type Needle (ZeroDim k) = ZeroDim k @@ -388,105 +493,114 @@ translateP = Tagged (.+~^) instance PseudoAffine (ZeroDim k) where Origin .-~. Origin = pure Origin +instance Num k => Semimanifold (V0 k) where + type Needle (V0 k) = V0 k + fromInterior = id + toInterior = pure + V0 .+~^ V0 = V0 + V0 .-~^ V0 = V0 + translateP = Tagged (.+~^) +instance Num k => PseudoAffine (V0 k) where + V0 .-~. V0 = pure V0 -instance (Semimanifold a, Semimanifold b) => Semimanifold (a,b) where +instance ∀ a b . (Semimanifold a, Semimanifold b) => Semimanifold (a,b) where type Needle (a,b) = (Needle a, Needle b) type Interior (a,b) = (Interior a, Interior b) (a,b).+~^(v,w) = (a.+~^v, b.+~^w) (a,b).-~^(v,w) = (a.-~^v, b.-~^w) fromInterior (i,j) = (fromInterior i, fromInterior j) toInterior (a,b) = fzip (toInterior a, toInterior b) - translateP = tp - where tp :: ∀ a b . (Semimanifold a, Semimanifold b) - => Tagged (a,b) ( (Interior a, Interior b) - -> (Needle a, Needle b) - -> (Interior a, Interior b) ) - tp = Tagged $ \(a,b) (v,w) -> (ta a v, tb b w) - where Tagged ta = translateP :: Tagged a (Interior a -> Needle a -> Interior a) - Tagged tb = translateP :: Tagged b (Interior b -> Needle b -> Interior b) + translateP = Tagged $ \(a,b) (v,w) -> (ta a v, tb b w) + where Tagged ta = translateP :: Tagged a (Interior a -> Needle a -> Interior a) + Tagged tb = translateP :: Tagged b (Interior b -> Needle b -> Interior b) + semimanifoldWitness = case ( semimanifoldWitness :: SemimanifoldWitness a + , semimanifoldWitness :: SemimanifoldWitness b ) of + (SemimanifoldWitness, SemimanifoldWitness) -> SemimanifoldWitness instance (PseudoAffine a, PseudoAffine b) => PseudoAffine (a,b) where (a,b).-~.(c,d) = liftA2 (,) (a.-~.c) (b.-~.d) -instance (PseudoAffine a, PseudoAffine b, PseudoAffine c) - => LocallyCoercible (a,(b,c)) ((a,b),c) where locallyTrivialDiffeomorphism = regroup -instance (PseudoAffine a, PseudoAffine b, PseudoAffine c) - => LocallyCoercible ((a,b),c) (a,(b,c)) where locallyTrivialDiffeomorphism = regroup' +instance ( Semimanifold a, Semimanifold b, Semimanifold c + , LSpace (Needle a), LSpace (Needle b), LSpace (Needle c) + , Scalar (Needle a) ~ Scalar (Needle b), Scalar (Needle b) ~ Scalar (Needle c) ) + => LocallyCoercible (a,(b,c)) ((a,b),c) where + locallyTrivialDiffeomorphism = regroup + coerceNeedle _ = regroup + coerceNeedle' _ = regroup + oppositeLocalCoercion = CanonicalDiffeomorphism + interiorLocalCoercion _ = case ( semimanifoldWitness :: SemimanifoldWitness a + , semimanifoldWitness :: SemimanifoldWitness b + , semimanifoldWitness :: SemimanifoldWitness c ) of + (SemimanifoldWitness, SemimanifoldWitness, SemimanifoldWitness) + -> CanonicalDiffeomorphism +instance ∀ a b c . + ( Semimanifold a, Semimanifold b, Semimanifold c + , LSpace (Needle a), LSpace (Needle b), LSpace (Needle c) + , Scalar (Needle a) ~ Scalar (Needle b), Scalar (Needle b) ~ Scalar (Needle c) ) + => LocallyCoercible ((a,b),c) (a,(b,c)) where + locallyTrivialDiffeomorphism = regroup' + coerceNeedle _ = regroup' + coerceNeedle' _ = regroup' + oppositeLocalCoercion = CanonicalDiffeomorphism + interiorLocalCoercion _ = case ( semimanifoldWitness :: SemimanifoldWitness a + , semimanifoldWitness :: SemimanifoldWitness b + , semimanifoldWitness :: SemimanifoldWitness c ) of + (SemimanifoldWitness, SemimanifoldWitness, SemimanifoldWitness) + -> CanonicalDiffeomorphism -instance (Semimanifold a, Semimanifold b, Semimanifold c) => Semimanifold (a,b,c) where +instance ∀ a b c . (Semimanifold a, Semimanifold b, Semimanifold c) + => Semimanifold (a,b,c) where type Needle (a,b,c) = (Needle a, Needle b, Needle c) type Interior (a,b,c) = (Interior a, Interior b, Interior c) (a,b,c).+~^(v,w,x) = (a.+~^v, b.+~^w, c.+~^x) (a,b,c).-~^(v,w,x) = (a.-~^v, b.-~^w, c.-~^x) fromInterior (i,j,k) = (fromInterior i, fromInterior j, fromInterior k) toInterior (a,b,c) = liftA3 (,,) (toInterior a) (toInterior b) (toInterior c) - translateP = tp - where tp :: ∀ a b v . (Semimanifold a, Semimanifold b, Semimanifold c) - => Tagged (a,b,c) ( (Interior a, Interior b, Interior c) - -> (Needle a, Needle b, Needle c) - -> (Interior a, Interior b, Interior c) ) - tp = Tagged $ \(a,b,c) (v,w,x) -> (ta a v, tb b w, tc c x) - where Tagged ta = translateP :: Tagged a (Interior a -> Needle a -> Interior a) - Tagged tb = translateP :: Tagged b (Interior b -> Needle b -> Interior b) - Tagged tc = translateP :: Tagged c (Interior c -> Needle c -> Interior c) + translateP = Tagged $ \(a,b,c) (v,w,x) -> (ta a v, tb b w, tc c x) + where Tagged ta = translateP :: Tagged a (Interior a -> Needle a -> Interior a) + Tagged tb = translateP :: Tagged b (Interior b -> Needle b -> Interior b) + Tagged tc = translateP :: Tagged c (Interior c -> Needle c -> Interior c) + semimanifoldWitness = case ( semimanifoldWitness :: SemimanifoldWitness a + , semimanifoldWitness :: SemimanifoldWitness b + , semimanifoldWitness :: SemimanifoldWitness c ) of + (SemimanifoldWitness, SemimanifoldWitness, SemimanifoldWitness) + -> SemimanifoldWitness instance (PseudoAffine a, PseudoAffine b, PseudoAffine c) => PseudoAffine (a,b,c) where (a,b,c).-~.(d,e,f) = liftA3 (,,) (a.-~.d) (b.-~.e) (c.-~.f) -instance (PseudoAffine a, PseudoAffine b, PseudoAffine c) - => LocallyCoercible (a,b,c) ((a,b),c) where - locallyTrivialDiffeomorphism (a,b,c) = ((a,b),c) -instance (PseudoAffine a, PseudoAffine b, PseudoAffine c) - => LocallyCoercible (a,b,c) (a,(b,c)) where - locallyTrivialDiffeomorphism (a,b,c) = (a,(b,c)) -instance (PseudoAffine a, PseudoAffine b, PseudoAffine c) - => LocallyCoercible ((a,b),c) (a,b,c) where - locallyTrivialDiffeomorphism ((a,b),c) = (a,b,c) -instance (PseudoAffine a, PseudoAffine b, PseudoAffine c) - => LocallyCoercible (a,(b,c)) (a,b,c) where - locallyTrivialDiffeomorphism (a,(b,c)) = (a,b,c) -instance (MetricScalar a, KnownNat n) => Semimanifold (FreeVect n a) where - type Needle (FreeVect n a) = FreeVect n a + +instance LinearManifold (a n) => Semimanifold (LinAff.Point a n) where + type Needle (LinAff.Point a n) = a n fromInterior = id toInterior = pure - translateP = Tagged (.+~^) - (.+~^) = (.+^) -instance (MetricScalar a, KnownNat n) => PseudoAffine (FreeVect n a) where - a.-~.b = pure (a.-.b) -instance LocallyCoercible ℝ (ℝ ^ S Z) where - locallyTrivialDiffeomorphism = replicVector -instance LocallyCoercible (ℝ ^ S Z) ℝ where - locallyTrivialDiffeomorphism = (<.>^replicVector 1) + LinAff.P v .+~^ w = LinAff.P $ v ^+^ w + translateP = Tagged $ \(LinAff.P v) w -> LinAff.P $ v ^+^ w +instance LinearManifold (a n) => PseudoAffine (LinAff.Point a n) where + LinAff.P v .-~. LinAff.P w = return $ v ^-^ w -instance (HasMetric a, FiniteDimensional b, Scalar a~Scalar b) => Semimanifold (a⊗b) where - type Needle (a⊗b) = a ⊗ b +instance (LSpace a, LSpace b, s~Scalar a, s~Scalar b) + => Semimanifold (Tensor s a b) where + type Needle (Tensor s a b) = Tensor s a b fromInterior = id toInterior = pure translateP = Tagged (.+~^) (.+~^) = (^+^) -instance (HasMetric a, FiniteDimensional b, Scalar a~Scalar b) => PseudoAffine (a⊗b) where +instance (LSpace a, LSpace b, s~Scalar a, s~Scalar b) + => PseudoAffine (Tensor s a b) where a.-~.b = pure (a^-^b) -instance (HasMetric a, FiniteDimensional b, Scalar a~Scalar b) => Semimanifold (a:-*b) where - type Needle (a:-*b) = DualSpace a ⊗ b - fromInterior = id - toInterior = pure - translateP = Tagged (.+~^) - 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 +instance (LSpace a, LSpace b, Scalar a~s, Scalar b~s) + => Semimanifold (LinearMap s a b) where + type Needle (LinearMap s a b) = LinearMap 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 +instance (LSpace a, LSpace b, Scalar a~s, Scalar b~s) + => PseudoAffine (LinearMap s a b) where a.-~.b = pure (a^-^b) instance Semimanifold S⁰ where - type Needle S⁰ = ℝ⁰ + type Needle S⁰ = ZeroDim ℝ fromInterior = id toInterior = pure translateP = Tagged (.+~^) @@ -555,7 +669,7 @@ fromInterior = id toInterior = pure translateP = Tagged (.+~^) - ℝP² r₀ φ₀ .+~^ (δr, δφ) + ℝP² r₀ φ₀ .+~^ V2 δr δφ | r₀ > 1/2 = case r₀ + δr of r₁ | r₁ > 1 -> ℝP² (2-r₁) (toS¹range $ φ₀+δφ+pi) | otherwise -> ℝP² r₁ (toS¹range $ φ₀+δφ) @@ -566,11 +680,11 @@ instance PseudoAffine ℝP² where ℝP² r₁ φ₁ .-~. ℝP² r₀ φ₀ | r₀ > 1/2 = pure `id` case φ₁-φ₀ of - δφ | δφ > 3*pi/2 -> ( r₁ - r₀, δφ - 2*pi) - | δφ < -3*pi/2 -> ( r₁ - r₀, δφ + 2*pi) - | δφ > pi/2 -> (2-r₁ - r₀, δφ - pi ) - | δφ < -pi/2 -> (2-r₁ - r₀, δφ + pi ) - | otherwise -> ( r₁ - r₀, δφ ) + δφ | δφ > 3*pi/2 -> V2 ( r₁ - r₀) (δφ - 2*pi) + | δφ < -3*pi/2 -> V2 ( r₁ - r₀) (δφ + 2*pi) + | δφ > pi/2 -> V2 (2-r₁ - r₀) (δφ - pi ) + | δφ < -pi/2 -> V2 (2-r₁ - r₀) (δφ + pi ) + | otherwise -> V2 ( r₁ - r₀) (δφ ) | otherwise = pure ( r₁*^embed(S¹ φ₁) ^-^ r₀*^embed(S¹ φ₀) ) @@ -597,22 +711,16 @@ class ImpliesMetric s where - {-# MINIMAL inferMetric | inferMetric' #-} type MetricRequirement s x :: Constraint type MetricRequirement s x = Semimanifold x - inferMetric :: (MetricRequirement s x, HasMetric (Needle x)) - => s x -> Option (Metric x) - inferMetric = safeRecipMetric <=< inferMetric' - inferMetric' :: (MetricRequirement s x, HasMetric (Needle x)) - => s x -> Option (Metric' x) - inferMetric' = safeRecipMetric' <=< inferMetric - -instance ImpliesMetric HerMetric where - type MetricRequirement HerMetric x = x ~ Needle x - inferMetric = pure + inferMetric :: (MetricRequirement s x, LSpace (Needle x)) + => s x -> Metric x + inferMetric' :: (MetricRequirement s x, LSpace (Needle x)) + => s x -> Metric' x -instance ImpliesMetric HerMetric' where - type MetricRequirement HerMetric' x = x ~ Needle x - inferMetric' = pure +instance ImpliesMetric Norm where + type MetricRequirement Norm x = (SimpleSpace x, x ~ Needle x) + inferMetric = id + inferMetric' = dualNorm diff --git a/Data/Manifold/Riemannian.hs b/Data/Manifold/Riemannian.hs --- a/Data/Manifold/Riemannian.hs +++ b/Data/Manifold/Riemannian.hs @@ -51,13 +51,14 @@ import Data.Semigroup import Data.VectorSpace -import Data.LinearMap.HerMetric +import Data.VectorSpace.Free import Data.AffineSpace +import Math.LinearMap.Category import Data.Manifold.Types import Data.Manifold.Types.Primitive ((^), empty, embed, coEmbed) +import Data.Manifold.Types.Stiefel import Data.Manifold.PseudoAffine -import Data.VectorSpace.FiniteDimensional import Data.CoNat @@ -68,8 +69,6 @@ import qualified Data.Foldable as Hask import qualified Data.Traversable as Hask -import qualified Numeric.LinearAlgebra.HMatrix as HMat - import Control.Category.Constrained.Prelude hiding ((^), all, elem, sum, forM, Foldable(..), Traversable) import Control.Arrow.Constrained @@ -111,18 +110,12 @@ = liftA3 (\ia ib ic t -> (ia t, ib t, ic t)) (geodesicBetween a α) (geodesicBetween b β) (geodesicBetween c γ) -instance (KnownNat n) => Geodesic (FreeVect n ℝ) where - geodesicBetween (FreeVect v) (FreeVect w) - = return $ \(D¹ t) -> let μv = (1-t)/2; μw = (t+1)/2 - in FreeVect $ Arr.zipWith (\vi wi -> μv*vi + μw*wi) v w - -instance (PseudoAffine v) => Geodesic (FinVecArrRep t v ℝ) where - geodesicBetween (FinVecArrRep v) (FinVecArrRep w) - | HMat.size v>0 && HMat.size w>0 - = return $ \(D¹ t) -> let μv = (1-t)/2; μw = (t+1)/2 - in FinVecArrRep $ HMat.scale μv v + HMat.scale μw w +-- instance (KnownNat n) => Geodesic (FreeVect n ℝ) where +-- geodesicBetween (FreeVect v) (FreeVect w) +-- = return $ \(D¹ t) -> let μv = (1-t)/2; μw = (t+1)/2 +-- in FreeVect $ Arr.zipWith (\vi wi -> μv*vi + μw*wi) v w -instance (Geodesic v, WithField ℝ HilbertSpace v) +instance (Geodesic v, FiniteFreeSpace v, WithField ℝ HilbertManifold v) => Geodesic (Stiefel1 v) where geodesicBetween (Stiefel1 p') (Stiefel1 q') = (\f -> \(D¹ t) -> Stiefel1 . f . D¹ $ g * tan (ϑ*t)) @@ -147,39 +140,39 @@ <$> geodesicBetween (-pi-φ) (pi-ϕ) -instance Geodesic (Cℝay S⁰) where - geodesicBetween p q = (>>> fromℝ) <$> geodesicBetween (toℝ p) (toℝ q) - where toℝ (Cℝay h PositiveHalfSphere) = h - toℝ (Cℝay h NegativeHalfSphere) = -h - fromℝ x | x>0 = Cℝay x PositiveHalfSphere - | otherwise = Cℝay (-x) NegativeHalfSphere - -instance Geodesic (CD¹ S⁰) where - geodesicBetween p q = (>>> fromI) <$> geodesicBetween (toI p) (toI q) - where toI (CD¹ h PositiveHalfSphere) = h - toI (CD¹ h NegativeHalfSphere) = -h - fromI x | x>0 = CD¹ x PositiveHalfSphere - | otherwise = CD¹ (-x) NegativeHalfSphere - -instance Geodesic (Cℝay S¹) where - geodesicBetween p q = (>>> fromP) <$> geodesicBetween (toP p) (toP q) - where fromP = fromInterior - toP w = case toInterior w of {Option (Just i) -> i} - -instance Geodesic (CD¹ S¹) where - geodesicBetween p q = (>>> fromI) <$> geodesicBetween (toI p) (toI q) - where toI (CD¹ h (S¹ φ)) = (h*cos φ, h*sin φ) - fromI (x,y) = CD¹ (sqrt $ x^2+y^2) (S¹ $ atan2 y x) - -instance Geodesic (Cℝay S²) where - geodesicBetween p q = (>>> fromP) <$> geodesicBetween (toP p) (toP q) - where fromP = fromInterior - toP w = case toInterior w of {Option (Just i) -> i} - -instance Geodesic (CD¹ S²) where - geodesicBetween p q = (>>> fromI) <$> geodesicBetween (toI p) (toI q :: ℝ³) - where toI (CD¹ h sph) = h *^ embed sph - fromI v = CD¹ (magnitude v) (coEmbed v) +-- instance Geodesic (Cℝay S⁰) where +-- geodesicBetween p q = (>>> fromℝ) <$> geodesicBetween (toℝ p) (toℝ q) +-- where toℝ (Cℝay h PositiveHalfSphere) = h +-- toℝ (Cℝay h NegativeHalfSphere) = -h +-- fromℝ x | x>0 = Cℝay x PositiveHalfSphere +-- | otherwise = Cℝay (-x) NegativeHalfSphere +-- +-- instance Geodesic (CD¹ S⁰) where +-- geodesicBetween p q = (>>> fromI) <$> geodesicBetween (toI p) (toI q) +-- where toI (CD¹ h PositiveHalfSphere) = h +-- toI (CD¹ h NegativeHalfSphere) = -h +-- fromI x | x>0 = CD¹ x PositiveHalfSphere +-- | otherwise = CD¹ (-x) NegativeHalfSphere +-- +-- instance Geodesic (Cℝay S¹) where +-- geodesicBetween p q = (>>> fromP) <$> geodesicBetween (toP p) (toP q) +-- where fromP = fromInterior +-- toP w = case toInterior w of {Option (Just i) -> i} +-- +-- instance Geodesic (CD¹ S¹) where +-- geodesicBetween p q = (>>> fromI) <$> geodesicBetween (toI p) (toI q) +-- where toI (CD¹ h (S¹ φ)) = (h*cos φ, h*sin φ) +-- fromI (x,y) = CD¹ (sqrt $ x^2+y^2) (S¹ $ atan2 y x) +-- +-- instance Geodesic (Cℝay S²) where +-- geodesicBetween p q = (>>> fromP) <$> geodesicBetween (toP p) (toP q) +-- where fromP = fromInterior +-- toP w = case toInterior w of {Option (Just i) -> i} +-- +-- instance Geodesic (CD¹ S²) where +-- geodesicBetween p q = (>>> fromI) <$> geodesicBetween (toI p) (toI q :: ℝ³) +-- where toI (CD¹ h sph) = h *^ embed sph +-- fromI v = CD¹ (magnitude v) (coEmbed v) #define geoVSpCone(c,t) \ instance (c) => Geodesic (Cℝay (t)) where { \ @@ -193,11 +186,11 @@ ; fromP (x,h) = CD¹ h (x^/h) \ ; toP (CD¹ h w) = ( h*^w, h ) } } -geoVSpCone ((), ℝ) -geoVSpCone ((), ℝ⁰) -geoVSpCone ((WithField ℝ HilbertSpace a, WithField ℝ HilbertSpace b, Geodesic (a,b)), (a,b)) -geoVSpCone (KnownNat n, FreeVect n ℝ) -geoVSpCone ((Geodesic v, WithField ℝ HilbertSpace v), FinVecArrRep t v ℝ) +-- geoVSpCone ((), ℝ) +-- geoVSpCone ((), ℝ⁰) +-- geoVSpCone ((WithField ℝ HilbertManifold a, WithField ℝ HilbertManifold b +-- , Geodesic (a,b)), (a,b)) +-- geoVSpCone (KnownNat n, FreeVect n ℝ) @@ -208,16 +201,16 @@ instance IntervalLike D¹ where toClosedInterval = id -instance IntervalLike (CD¹ S⁰) where - toClosedInterval (CD¹ h PositiveHalfSphere) = D¹ h - toClosedInterval (CD¹ h NegativeHalfSphere) = D¹ (-h) -instance IntervalLike (Cℝay S⁰) where - toClosedInterval (Cℝay h PositiveHalfSphere) = D¹ $ tanh h - toClosedInterval (Cℝay h NegativeHalfSphere) = D¹ $ -tanh h -instance IntervalLike (CD¹ ℝ⁰) where - toClosedInterval (CD¹ h Origin) = D¹ $ h*2 - 1 -instance IntervalLike (Cℝay ℝ⁰) where - toClosedInterval (Cℝay h Origin) = D¹ $ 1 - 2/(h+1) +-- instance IntervalLike (CD¹ S⁰) where +-- toClosedInterval (CD¹ h PositiveHalfSphere) = D¹ h +-- toClosedInterval (CD¹ h NegativeHalfSphere) = D¹ (-h) +-- instance IntervalLike (Cℝay S⁰) where +-- toClosedInterval (Cℝay h PositiveHalfSphere) = D¹ $ tanh h +-- toClosedInterval (Cℝay h NegativeHalfSphere) = D¹ $ -tanh h +-- instance IntervalLike (CD¹ ℝ⁰) where +-- toClosedInterval (CD¹ h Origin) = D¹ $ h*2 - 1 +-- instance IntervalLike (Cℝay ℝ⁰) where +-- toClosedInterval (Cℝay h Origin) = D¹ $ 1 - 2/(h+1) instance IntervalLike ℝ where toClosedInterval x = D¹ $ tanh x @@ -229,4 +222,4 @@ rieMetric :: RieMetric m instance Riemannian ℝ where - rieMetric = const m where m = projector 1 + rieMetric = const euclideanNorm diff --git a/Data/Manifold/TreeCover.hs b/Data/Manifold/TreeCover.hs --- a/Data/Manifold/TreeCover.hs +++ b/Data/Manifold/TreeCover.hs @@ -42,7 +42,7 @@ -- ** Lenses , shadeCtr, shadeExpanse, shadeNarrowness -- ** Construction - , fullShade, fullShade', pointsShades, pointsCovers + , fullShade, fullShade', pointsShades, pointsShade's, pointsCovers, pointsCover's -- ** Evaluation , occlusion -- ** Misc @@ -56,13 +56,14 @@ , SimpleTree, Trees, NonEmptyTree, GenericTree(..) -- * Misc , sShSaw, chainsaw, HasFlatView(..), shadesMerge, smoothInterpolate - , twigsWithEnvirons, completeTopShading, flexTwigsShading + , twigsWithEnvirons, Twig, TwigEnviron + , completeTopShading, flexTwigsShading , WithAny(..), Shaded, fmapShaded, stiAsIntervalMapping, spanShading , constShaded, stripShadedUntopological , DifferentialEqn, propagateDEqnSolution_loc -- ** Triangulation-builders - , TriangBuild, doTriangBuild, singleFullSimplex, autoglueTriangulation - , AutoTriang, elementaryTriang, breakdownAutoTriang + , TriangBuild, doTriangBuild + , AutoTriang, breakdownAutoTriang ) where @@ -79,8 +80,7 @@ import Data.VectorSpace import Data.AffineSpace -import Data.LinearMap.HerMetric -import Data.LinearMap.Category +import Math.LinearMap.Category import Data.Tagged import Data.SimplicialComplex @@ -107,8 +107,6 @@ import qualified Data.Traversable as Hask import Data.Traversable (forM) -import qualified Numeric.LinearAlgebra.HMatrix as HMat - import Control.Category.Constrained.Prelude hiding ((^), all, elem, sum, forM, Foldable(..), foldr1, Traversable, traverse) import Control.Arrow.Constrained @@ -117,6 +115,7 @@ import Data.Traversable.Constrained (traverse) import GHC.Generics (Generic) +import Data.Type.Coercion -- | Possibly / Partially / asymPtotically singular metric. @@ -135,15 +134,14 @@ -- 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) +deriving instance (Show x, Show (Metric' 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) +deriving instance (Show x, Show (Metric x), WithField ℝ Manifold x) => Show (Shade' x) class IsShade shade where -- type (*) shade :: *->* @@ -153,10 +151,12 @@ -- 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 ) + occlusion :: ( Manifold x, SimpleSpace (Needle x) + , s ~ (Scalar (Needle x)), RealDimension s ) => shade x -> x -> s - factoriseShade :: ( Manifold x, RealDimension (Scalar (Needle x)) - , Manifold y, RealDimension (Scalar (Needle y)) ) + factoriseShade :: ( Manifold x, SimpleSpace (Needle x) + , Manifold y, SimpleSpace (Needle y) + , Scalar (Needle x) ~ Scalar (Needle y) ) => shade (x,y) -> (shade x, shade y) coerceShade :: (Manifold x, Manifold y, LocallyCoercible x y) => shade x -> shade y @@ -164,23 +164,31 @@ 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 + Option(Just vd) | mSq <- normSq δinv vd , mSq == mSq -- avoid NaN -> exp (negate mSq) _ -> zeroV - δinv = recipMetric δ + δinv = dualNorm δ factoriseShade (Shade (x₀,y₀) δxy) = (Shade x₀ δx, Shade y₀ δy) - where (δx,δy) = factoriseMetric' δxy - coerceShade (Shade x (HerMetric' δxym)) - = Shade (locallyTrivialDiffeomorphism x) (HerMetric' $ unsafeCoerceLinear<$>δxym) + where (δx,δy) = summandSpaceNorms δxy + coerceShade = cS + where cS :: ∀ x y . (LocallyCoercible x y) => Shade x -> Shade y + cS = \(Shade x δxym) -> Shade (internCoerce x) (tN δxym) + where tN = case oppositeLocalCoercion :: CanonicalDiffeomorphism y x of + CanonicalDiffeomorphism -> + transformNorm . arr $ coerceNeedle' ([]::[(y,x)]) + internCoerce = case interiorLocalCoercion ([]::[(x,y)]) of + CanonicalDiffeomorphism -> locallyTrivialDiffeomorphism instance ImpliesMetric Shade where - type MetricRequirement Shade x = Manifold x - inferMetric' (Shade _ e) = pure e + type MetricRequirement Shade x = (Manifold x, SimpleSpace (Needle x)) + inferMetric' (Shade _ e) = e + inferMetric (Shade _ e) = dualNorm e instance ImpliesMetric Shade' where - type MetricRequirement Shade' x = Manifold x - inferMetric (Shade' _ e) = pure e + type MetricRequirement Shade' x = (Manifold x, SimpleSpace (Needle x)) + inferMetric (Shade' _ e) = e + inferMetric' (Shade' _ e) = dualNorm e shadeExpanse :: Lens' (Shade x) (Metric' x) shadeExpanse f (Shade c e) = fmap (Shade c) $ f e @@ -189,14 +197,20 @@ 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 + Option(Just vd) | mSq <- normSq δ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 - coerceShade (Shade' x (HerMetric δxym)) - = Shade' (locallyTrivialDiffeomorphism x) (HerMetric $ unsafeCoerceLinear<$>δxym) + where (δx,δy) = summandSpaceNorms δxy + coerceShade = cS + where cS :: ∀ x y . (LocallyCoercible x y) => Shade' x -> Shade' y + cS = \(Shade' x δxym) -> Shade' (internCoerce x) (tN δxym) + where tN = case oppositeLocalCoercion :: CanonicalDiffeomorphism y x of + CanonicalDiffeomorphism -> + transformNorm . arr $ coerceNeedle ([]::[(y,x)]) + internCoerce = case interiorLocalCoercion ([]::[(x,y)]) of + CanonicalDiffeomorphism -> locallyTrivialDiffeomorphism shadeNarrowness :: Lens' (Shade' x) (Metric x) shadeNarrowness f (Shade' c e) = fmap (Shade' c) $ f e @@ -209,13 +223,13 @@ Shade c e .+~^ v = Shade (c.+^v) e Shade c e .-~^ v = Shade (c.-^v) e -instance (WithField ℝ AffineManifold x, Geodesic x) => Geodesic (Shade x) where +instance (WithField ℝ AffineManifold x, Geodesic x, SimpleSpace (Needle x)) + => Geodesic (Shade x) where geodesicBetween (Shade c e) (Shade ζ η) = pure interp - where ([], sharedSpan) = eigenSystem (e,η) + where sharedSpan = sharedNormSpanningSystem e η interp t = Shade (pinterp t) - (projector's [ v ^* (alerpB qe qη t) - | ([qe,qη], (v,_)) <- zip coeffs sharedSpan ]) - coeffs = [ [metric' m v' | m <- [e,η]] | (_,v') <- sharedSpan ] + (spanNorm [ v ^* (alerpB 1 qη t) + | (v,qη) <- sharedSpan ]) Option (Just pinterp) = geodesicBetween c ζ instance (AffineManifold x) => Semimanifold (Shade' x) where @@ -226,13 +240,13 @@ Shade' c e .+~^ v = Shade' (c.+^v) e Shade' c e .-~^ v = Shade' (c.-^v) e -instance (WithField ℝ AffineManifold x, Geodesic x) => Geodesic (Shade' x) where +instance (WithField ℝ AffineManifold x, Geodesic x, SimpleSpace (Needle x)) + => Geodesic (Shade' x) where geodesicBetween (Shade' c e) (Shade' ζ η) = pure interp - where ([], sharedSpan) = eigenSystem (e,η) + where sharedSpan = sharedNormSpanningSystem e η interp t = Shade' (pinterp t) - (projectors [ v' ^/ (alerpB qe qη t) - | ([qe,qη], (v',_)) <- zip coeffs sharedSpan ]) - coeffs = [ [recip $ metric m v | m <- [e,η]] | (_,v) <- sharedSpan ] + (spanNorm [ v ^/ (alerpB 1 (recip qη) t) + | (v,qη) <- sharedSpan ]) Option (Just pinterp) = geodesicBetween c ζ fullShade :: WithField ℝ Manifold x => x -> Metric' x -> Shade x @@ -243,9 +257,10 @@ -- | Span a 'Shade' from a center point and multiple deviation-vectors. -pattern (:±) :: () => WithField ℝ Manifold x => x -> [Needle x] -> Shade x -pattern x :± shs <- Shade x (eigenSpan -> shs) - where x :± shs = fullShade x $ projector's shs +pattern (:±) :: () => (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => x -> [Needle x] -> Shade x +pattern x :± shs <- Shade x (normSpanningSystem -> shs) + where x :± shs = fullShade x $ spanVariance shs -- | Similar to ':±', but instead of expanding the shade, each vector /restricts/ it. @@ -255,7 +270,7 @@ -- Note that '|±|' is only possible, as such, in an inner-product space; in -- general you need reciprocal vectors ('Needle'') to define a 'Shade''. (|±|) :: WithField ℝ EuclidSpace x => x -> [Needle x] -> Shade' x -x |±| shs = Shade' x $ projectors [v^/(v<.>v) | v<-shs] +x |±| shs = Shade' x $ spanNorm [v^/(v<.>v) | v<-shs] @@ -267,8 +282,9 @@ in (iu, if vl>0 then UpperBulb else LowerBulb) _ -> (-1, error "Trying to obtain the subshadeId of a point not actually included in the shade.") -subshadeId :: WithField ℝ Manifold x => Shade x -> x -> (Int, HourglassBulb) -subshadeId (Shade c expa) = subshadeId' c . NE.fromList $ eigenCoSpan expa +subshadeId :: (WithField ℝ Manifold x, FiniteDimensional (Needle' x)) + => Shade x -> x -> (Int, HourglassBulb) +subshadeId (Shade c expa) = subshadeId' c . NE.fromList $ normSpanningSystem' expa @@ -282,31 +298,37 @@ -- For /nonconnected/ manifolds it will be necessary to yield separate shades -- for each connected component. And for an empty input list, there is no shade! -- Hence the result type is a list. -pointsShades :: WithField ℝ Manifold x => [x] -> [Shade x] -pointsShades = map snd . pointsShades' zeroV +pointsShades :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => [x] -> [Shade x] +pointsShades = map snd . pointsShades' mempty -- | Like 'pointsShades', but ensure that all points are actually in -- the shade, i.e. if @['Shade' x₀ ex]@ is the result then -- @'metric' (recipMetric ex) (p-x₀) ≤ 1@ for all @p@ in the list. -pointsCovers :: ∀ x . WithField ℝ Manifold x => [x] -> [Shade x] -pointsCovers = map guaranteeIn . pointsShades' zeroV +pointsCovers :: ∀ x . (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => [x] -> [Shade x] +pointsCovers = map guaranteeIn . pointsShades' mempty where guaranteeIn (ps, Shade x₀ ex) = case ps >>= \p -> let Option (Just v) = p.-~.x₀ - in guard (metric ex' v > 1) >> [(p,projector' v)] + in guard ((ex'|$|v) > 1) >> [(p, spanVariance [v])] of [] -> Shade x₀ ex outs -> guaranteeIn ( fst<$>outs , Shade x₀ - $ ex ^+^ sumV (snd<$>outs) - ^/ fromIntegral (2 * length outs) ) - where ex' = recipMetric ex + $ ex <> scaleNorm + (sqrt . recip . fromIntegral + $ 2 * length outs) + (mconcat $ snd<$>outs) + ) + where ex' = dualNorm ex -pointsShade's :: WithField ℝ Manifold x => [x] -> [Shade' x] -pointsShade's = map (\(Shade c e) -> Shade' c $ recipMetric e) . pointsShades +pointsShade's :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) => [x] -> [Shade' x] +pointsShade's = map (\(Shade c e) -> Shade' c $ dualNorm e) . pointsShades -pointsCover's :: WithField ℝ Manifold x => [x] -> [Shade' x] -pointsCover's = map (\(Shade c e) -> Shade' c $ recipMetric e) . pointsCovers +pointsCover's :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) => [x] -> [Shade' x] +pointsCover's = map (\(Shade c e) -> Shade' c $ dualNorm e) . pointsCovers -pseudoECM :: WithField ℝ Manifold x => NonEmpty x -> (x, ([x],[x])) +pseudoECM :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => NonEmpty x -> (x, ([x],[x])) pseudoECM (p₀ NE.:| psr) = foldl' ( \(acc, (rb,nr)) (i,p) -> case p.-~.acc of Option (Just δ) -> (acc .+~^ δ^/i, (p:rb, nr)) @@ -314,7 +336,8 @@ (p₀, mempty) ( zip [1..] $ p₀:psr ) -pointsShades' :: WithField ℝ Manifold x => Metric' x -> [x] -> [([x], Shade x)] +pointsShades' :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => Metric' x -> [x] -> [([x], Shade x)] pointsShades' _ [] = [] pointsShades' minExt ps = case expa of Option (Just e) -> (ps, fullShade ctr e) @@ -322,14 +345,14 @@ _ -> pointsShades' minExt inc'd ++ pointsShades' minExt unreachable where (ctr,(inc'd,unreachable)) = pseudoECM $ NE.fromList ps - expa = ( (^+^minExt) . (^/ fromIntegral(length ps)) . projector's ) + expa = ( (<>minExt) . spanVariance . map (^/ fromIntegral (length ps)) ) <$> mapM (.-~.ctr) ps -- | Attempt to reduce the number of shades to fewer (ideally, a single one). -- In the simplest cases these should guaranteed cover the same area; -- for non-flat manifolds it only works in a heuristic sense. -shadesMerge :: WithField ℝ Manifold x +shadesMerge :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) => ℝ -- ^ How near (inverse normalised distance, relative to shade expanse) -- two shades must be to be merged. If this is zero, any shades -- in the same connected region of a manifold are merged. @@ -343,14 +366,14 @@ where tryMerge (Shade c₂ e₂) | Option (Just v) <- c₁.-~.c₂ , Option (Just v') <- c₂.-~.c₁ - , [e₁',e₂'] <- recipMetric<$>[e₁, e₂] - , b₁ <- metric e₂' v - , b₂ <- metric e₁' v + , [e₁',e₂'] <- dualNorm<$>[e₁, e₂] + , b₁ <- e₂'|$|v + , b₂ <- e₁'|$|v , fuzz*b₁*b₂ <= b₁ + b₂ = Just $ let cc = c₂ .+~^ v ^/ 2 Option (Just cv₁) = c₁.-~.cc Option (Just cv₂) = c₂.-~.cc - in Shade cc $ e₁ ^+^ e₂ ^+^ projector's [cv₁, cv₂] + in Shade cc $ e₁ <> e₂ <> spanVariance [cv₁, cv₂] | otherwise = Nothing shadesMerge _ shs = shs @@ -364,19 +387,20 @@ => Shade' x -> x -> s minusLogOcclusion' (Shade' p₀ δinv) = occ where occ p = case p .-~. p₀ of - Option(Just vd) | mSq <- metricSq δinv vd + Option(Just vd) | mSq <- normSq δinv vd , mSq == mSq -- avoid NaN -> mSq _ -> 1/0 -minusLogOcclusion :: ( Manifold x, s ~ (Scalar (Needle x)), RealDimension s ) +minusLogOcclusion :: ( Manifold x, SimpleSpace (Needle 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 + Option(Just vd) | mSq <- normSq δinv vd , mSq == mSq -- avoid NaN -> mSq _ -> 1/0 - δinv = recipMetric δ + δinv = dualNorm δ @@ -490,13 +514,13 @@ DisjointBranches n br .+~^ v = DisjointBranches n $ (.+~^v)<$>br -- | WRT union. -instance WithField ℝ Manifold x => Semigroup (ShadeTree x) where +instance (WithField ℝ Manifold x, SimpleSpace (Needle x)) => Semigroup (ShadeTree x) where PlainLeaves [] <> t = t t <> PlainLeaves [] = t t <> s = fromLeafPoints $ onlyLeaves t ++ onlyLeaves s -- Could probably be done more efficiently sconcat = mconcat . NE.toList -instance WithField ℝ Manifold x => Monoid (ShadeTree x) where +instance (WithField ℝ Manifold x, SimpleSpace (Needle x)) => Monoid (ShadeTree x) where mempty = PlainLeaves [] mappend = (<>) mconcat l = case filter ne l of @@ -511,7 +535,8 @@ -- Example: https://nbviewer.jupyter.org/github/leftaroundabout/manifolds/blob/master/test/Trees-and-Webs.ipynb#pseudorandomCloudTree -- -- <> -fromLeafPoints :: ∀ x. WithField ℝ Manifold x => [x] -> ShadeTree x +fromLeafPoints :: ∀ x. (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => [x] -> ShadeTree x fromLeafPoints = fromLeafPoints' sShIdPartition @@ -541,7 +566,7 @@ -- | “Inverse indexing” of a tree. This is roughly a nearest-neighbour search, -- but not guaranteed to give the correct result unless evaluated at the -- precise position of a tree leaf. -positionIndex :: ∀ x . WithField ℝ Manifold x +positionIndex :: ∀ x . (WithField ℝ Manifold x, SimpleSpace (Needle x)) => Option (Metric x) -- ^ For deciding (at the lowest level) what “close” means; -- this is optional for any tree of depth >1. -> ShadeTree x -- ^ The tree to index into @@ -551,7 +576,7 @@ -- environment trees leading down to its position (in decreasing -- order of size), and actual position of the found node. positionIndex (Option (Just m)) sh@(PlainLeaves lvs) x - = case catMaybes [ ((i,p),) . metricSq m <$> getOption (p.-~.x) + = case catMaybes [ ((i,p),) . normSq m <$> getOption (p.-~.x) | (i,p) <- zip [0..] lvs] of [] -> empty l | ((i,p),_) <- minimumBy (comparing snd) l @@ -571,26 +596,27 @@ , let ω = d<.>^vx , (t',σ) <- [(t'u, 1), (t'd, -1)] ] in ((+i₀) *** first (sh:)) - <$> positionIndex (return $ recipMetric ce) t' x + <$> positionIndex (return $ dualNorm ce) t' x positionIndex _ _ _ = empty -fromFnGraphPoints :: ∀ x y . (WithField ℝ Manifold x, WithField ℝ Manifold y) +fromFnGraphPoints :: ∀ x y . ( WithField ℝ Manifold x, WithField ℝ Manifold y + , SimpleSpace (Needle x), SimpleSpace (Needle y) ) => [(x,y)] -> ShadeTree (x,y) fromFnGraphPoints = fromLeafPoints' fg_sShIdPart where fg_sShIdPart :: Shade (x,y) -> [(x,y)] -> NonEmpty (DBranch' (x,y) [(x,y)]) fg_sShIdPart (Shade c expa) xs | b:bs <- [DBranch (v, zeroV) mempty - | v <- eigenCoSpan - (transformMetric' fst expa :: Metric' x) ] + | v <- normSpanningSystem' + (transformNorm (id&&&zeroV) expa :: Metric' x) ] = sShIdPartition' c xs $ b:|bs -fromLeafPoints' :: ∀ x. WithField ℝ Manifold x => +fromLeafPoints' :: ∀ x. (WithField ℝ Manifold x, SimpleSpace (Needle x)) => (Shade x -> [x] -> NonEmpty (DBranch' x [x])) -> [x] -> ShadeTree x -fromLeafPoints' sShIdPart = go zeroV +fromLeafPoints' sShIdPart = go mempty where go :: Metric' x -> [x] -> ShadeTree x - go preShExpa = \xs -> case pointsShades' (preShExpa^/10) xs of + go preShExpa = \xs -> case pointsShades' (scaleNorm (1/3) preShExpa) xs of [] -> mempty [(_,rShade)] -> let trials = sShIdPart rShade xs in case reduce rShade trials of @@ -601,7 +627,7 @@ _ -> PlainLeaves xs partitions -> DisjointBranches (length xs) . NE.fromList - $ map (\(xs',pShade) -> go zeroV xs') partitions + $ map (\(xs',pShade) -> go mempty xs') partitions where branchProc redSh = fmap (fmap $ go redSh) @@ -630,9 +656,10 @@ i ) st xs where ssi = subshadeId' c (boughDirection<$>st) -sShIdPartition :: WithField ℝ Manifold x => Shade x -> [x] -> NonEmpty (DBranch' x [x]) +sShIdPartition :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => Shade x -> [x] -> NonEmpty (DBranch' x [x]) sShIdPartition (Shade c expa) xs - | b:bs <- [DBranch v mempty | v <- eigenCoSpan expa] + | b:bs <- [DBranch v mempty | v <- normSpanningSystem' expa] = sShIdPartition' c xs $ b:|bs @@ -669,7 +696,8 @@ sortByKey = map snd . sortBy (comparing fst) -trunks :: ∀ x. WithField ℝ Manifold x => ShadeTree x -> [Shade x] +trunks :: ∀ x. (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => ShadeTree x -> [Shade x] trunks (PlainLeaves lvs) = pointsCovers lvs trunks (DisjointBranches _ brs) = Hask.foldMap trunks brs trunks (OverlappingBranches _ sh _) = [sh] @@ -682,11 +710,16 @@ instance ImpliesMetric ShadeTree where - type MetricRequirement ShadeTree x = WithField ℝ Manifold x - inferMetric' (OverlappingBranches _ (Shade _ e) _) = pure e + type MetricRequirement ShadeTree x = (WithField ℝ Manifold x, SimpleSpace (Needle x)) + inferMetric (OverlappingBranches _ (Shade _ e) _) = dualNorm e + inferMetric (PlainLeaves lvs) = case pointsShades lvs of + (Shade _ sh:_) -> dualNorm sh + _ -> mempty + inferMetric (DisjointBranches _ (br:|_)) = inferMetric br + inferMetric' (OverlappingBranches _ (Shade _ e) _) = e inferMetric' (PlainLeaves lvs) = case pointsShades lvs of - (Shade _ sh:_) -> pure sh - _ -> empty + (Shade _ sh:_) -> sh + _ -> mempty inferMetric' (DisjointBranches _ (br:|_)) = inferMetric' br @@ -720,36 +753,36 @@ -- | Class of manifolds which can use 'Shade'' as a basic set type. -- This is easily possible for vector spaces with the default implementations. -class (WithField ℝ Manifold y) => Refinable y where +class (WithField ℝ Manifold y, SimpleSpace (Needle y)) => Refinable y where -- | @a `subShade'` b ≡ True@ means @a@ is fully contained in @b@, i.e. from -- @'minusLogOcclusion'' a p < 1@ follows also @minusLogOcclusion' b p < 1@. subShade' :: Shade' y -> Shade' y -> Bool subShade' (Shade' ac ae) tsh = all ((<1) . minusLogOcclusion' tsh) - [ ac.+~^σ*^v | σ<-[-1,1], v<-eigenCoSpan' ae ] + [ ac.+~^σ*^v | σ<-[-1,1], v<-normSpanningSystem' ae ] refineShade' :: Shade' y -> Shade' y -> Option (Shade' y) - refineShade' (Shade' c₀ (HerMetric (Just e₁))) - (Shade' c₀₂ (HerMetric (Just e₂))) + refineShade' (Shade' c₀ (Norm e₁)) + (Shade' c₀₂ (Norm e₂)) | Option (Just c₂) <- c₀₂.-~.c₀ , e₁c₂ <- e₁ $ c₂ , e₂c₂ <- e₂ $ c₂ - , cc <- σe <\$ e₂c₂ + , cc <- σe \$ e₂c₂ , cc₂ <- cc ^-^ c₂ , e₁cc <- e₁ $ cc , e₂cc <- e₂ $ cc , α <- 2 + cc₂<.>^e₂c₂ , α > 0 , ee <- σe ^/ α - , c₂e₁c₂ <- c₂^<.>e₁c₂ - , c₂e₂c₂ <- c₂^<.>e₂c₂ + , c₂e₁c₂ <- c₂<.>^e₁c₂ + , c₂e₂c₂ <- c₂<.>^e₂c₂ , c₂eec₂ <- (c₂e₁c₂ + c₂e₂c₂) / α , [γ₁,γ₂] <- middle . sort $ quadraticEqnSol c₂e₁c₂ - (2 * (c₂^<.>e₁cc)) - (cc^<.>e₁cc - 1) + (2 * (c₂<.>^e₁cc)) + (cc<.>^e₁cc - 1) ++ quadraticEqnSol c₂e₂c₂ - (2 * (c₂^<.>e₂cc - c₂e₂c₂)) - (cc^<.>e₂cc - 2 * (cc^<.>e₂c₂) + c₂e₂c₂ - 1) + (2 * (c₂<.>^e₂cc - c₂e₂c₂)) + (cc<.>^e₂cc - 2 * (cc<.>^e₂c₂) + c₂e₂c₂ - 1) , cc' <- cc ^+^ ((γ₁+γ₂)/2)*^c₂ , rγ <- abs (γ₁ - γ₂) / 2 , η <- if rγ * c₂eec₂ /= 0 && 1 - rγ^2 * c₂eec₂ > 0 @@ -757,10 +790,9 @@ else 0 = return $ Shade' (c₀.+~^cc') - (HerMetric (Just ee) ^+^ projector (ee $ c₂^*η)) - + (Norm (arr ee) <> spanNorm [ee $ c₂^*η]) | otherwise = empty - where σe = e₁^+^e₂ + where σe = arr $ e₁^+^e₂ quadraticEqnSol a b c | a /= 0 && disc > 0 = [ (σ * sqrt disc - b) / (2*a) | σ <- [-1, 1] ] @@ -768,8 +800,6 @@ where disc = b^2 - 4*a*c middle (_:x:y:_) = [x,y] middle l = l - refineShade' (Shade' _ (HerMetric Nothing)) s₂ = pure s₂ - refineShade' s₁ (Shade' _ (HerMetric Nothing)) = pure s₁ -- ⟨x−c₁|e₁|x−c₁⟩ < 1 ∧ ⟨x−c₂|e₂|x−c₂⟩ < 1 -- We search (cc,ee) such that this implies -- ⟨x−cc|ee|x−cc⟩ < 1. @@ -854,11 +884,11 @@ convolveShade' :: Shade' y -> Shade' (Needle y) -> Shade' y convolveShade' (Shade' y₀ ey) (Shade' δ₀ eδ) = Shade' (y₀.+~^δ₀) - ( projectors [ f ^* ζ crl - | (f,_) <- eδsp - | crl <- corelap ] ) - where (_,eδsp) = eigenSystem (ey,eδ) - corelap = map (metric ey . snd) eδsp + ( spanNorm [ f ^* ζ crl + | (f,_) <- eδsp + | crl <- corelap ] ) + where eδsp = sharedNormSpanningSystem ey eδ + corelap = map snd eδsp ζ = case filter (>0) corelap of [] -> const 0 nzrelap @@ -873,7 +903,7 @@ instance Refinable ℝ where refineShade' (Shade' cl el) (Shade' cr er) - = case (metricSq el 1, metricSq er 1) of + = case (normSq el 1, normSq er 1) of (0, _) -> return $ Shade' cr er (_, 0) -> return $ Shade' cl el (ql,qr) | ql>0, qr>0 @@ -883,7 +913,7 @@ in guard (b> let cm = (b+t)/2 rm = (t-b)/2 - in return $ Shade' cm (projector $ recip rm) + in return $ Shade' cm (spanNorm [recip rm]) -- convolveShade' (Shade' y₀ ey) (Shade' δ₀ eδ) -- = case (metricSq ey 1, metricSq eδ 1) of -- (wy,wδ) | wy>0, wδ>0 @@ -894,6 +924,11 @@ instance (Refinable a, Refinable b) => Refinable (a,b) +instance Refinable ℝ⁰ +instance Refinable ℝ¹ +instance Refinable ℝ² +instance Refinable ℝ³ +instance Refinable ℝ⁴ intersectShade's :: ∀ y . Refinable y => NonEmpty (Shade' y) -> Option (Shade' y) @@ -905,7 +940,8 @@ type DifferentialEqn x y = Shade (x,y) -> Shade' (LocalLinear x y) -propagateDEqnSolution_loc :: ∀ x y . (WithField ℝ Manifold x, Refinable y) +propagateDEqnSolution_loc :: ∀ x y . ( WithField ℝ Manifold x, Refinable y + , SimpleSpace (Needle x) ) => DifferentialEqn x y -> ((x, Shade' y), NonEmpty (Needle x, Shade' y)) -> NonEmpty (Shade' y) propagateDEqnSolution_loc f ((x, shy@(Shade' y _)), neighbours) = ycs @@ -913,7 +949,7 @@ [shxy] = pointsCovers [ (xs, ys') | (xs, Shade' ys yse) <- (x,shy):(first (x.+~^)<$>NE.toList neighbours) - , δy <- eigenCoSpan' yse + , δy <- normSpanningSystem' yse , ys' <- [ys.+~^δy, ys.-~^δy] ] [Shade' _ expax] = pointsCover's $ x : ((x.+~^).fst<$>NE.toList neighbours) marginδs :: NonEmpty (Needle x, (Needle y, Metric y)) @@ -925,28 +961,38 @@ back2Centre (δx, (δym, expany)) = convolveShade' (Shade' y expany) - (Shade' δyb $ applyLinMapMetric jExpa (δx'^/(δx'<.>^δx))) + (Shade' δyb $ applyLinMapNorm jExpa (δx'^/(δx'<.>^δx))) where δyb = δym ^-^ (j₀ $ δx) - δx' = toDualWith expax δx + δx' = expax<$|δx ycs :: NonEmpty (Shade' y) ycs = back2Centre <$> marginδs - xSpan = eigenCoSpan' expax + xSpan = normSpanningSystem expax +applyLinMapNorm :: (LSpace x, LSpace y, Scalar x ~ Scalar y) + => Norm (x+>y) -> DualVector x -> Norm y +applyLinMapNorm n dx + = transformNorm (fmap (arr Coercion . transposeTensor) . blockVectSpan' $ dx) n --- 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 -> [((Int, ShadeTree x), [(Int, ShadeTree x)])] + +type Twig x = (Int, ShadeTree x) +type TwigEnviron x = [Twig x] + +-- Formerly, 'twigsWithEnvirons' what has now become 'traverseTwigsWithEnvirons'. +-- The simple list-yielding version (see rev. b4a427d59ec82889bab2fde39225b14a57b694df) +-- may well be more efficient than the current traversal-derived version. + +-- | Example: https://nbviewer.jupyter.org/github/leftaroundabout/manifolds/blob/master/test/Trees-and-Webs.ipynb#pseudorandomCloudTree +-- +-- <> +twigsWithEnvirons :: ∀ x. (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => ShadeTree x -> [(Twig x, TwigEnviron x)] twigsWithEnvirons = execWriter . traverseTwigsWithEnvirons (writer . (snd.fst&&&pure)) traverseTwigsWithEnvirons :: ∀ x f . - (WithField ℝ Manifold x, Hask.Applicative f) - => ( ((Int, ShadeTree x), [(Int, ShadeTree x)]) -> f (ShadeTree x)) - -> ShadeTree x -> f (ShadeTree x) + (WithField ℝ Manifold x, SimpleSpace (Needle x), Hask.Applicative f) + => ( (Twig x, TwigEnviron x) -> f (ShadeTree x) ) -> ShadeTree x -> f (ShadeTree x) traverseTwigsWithEnvirons f = fst . go [] . (0,) - where go :: [(Int, ShadeTree x)] -> (Int, ShadeTree x) - -> (f (ShadeTree x), Bool) + where go :: TwigEnviron x -> Twig x -> (f (ShadeTree x), Bool) go _ (i₀, DisjointBranches nlvs djbs) = ( fmap (DisjointBranches nlvs) . Hask.traverse (fst . go []) $ NE.zip ioffs djbs @@ -967,7 +1013,7 @@ where envi'' = filter (snd >>> trunks >>> \(Shade ce _:_) -> let Option (Just δyenv) = ce.-~.robc qq = vy<.>^δyenv - in qq > -1 && qq < 5 + in qq > -1 ) envi' ++ map ((+i₀)***snd) alts envi' = approach =<< envi @@ -975,13 +1021,14 @@ = first (+i₀e) <$> twigsaveTrim hither apt where Option (Just δxenv) = robc .-~. envc hither (DBranch bdir (Hourglass bdc₁ bdc₂)) - | bdir<.>^δxenv > 0 = [(0 , bdc₁)] - | otherwise = [(nLeaves bdc₁, bdc₂)] + = [(0 , bdc₁) | overlap > -1] + ++ [(nLeaves bdc₁, bdc₂) | overlap < 1] + where overlap = bdir<.>^δxenv approach q = [q] go envi plvs@(i₀, (PlainLeaves _)) = (f $ purgeRemotes (plvs, envi), True) - twigProximæ :: x -> ShadeTree x -> [(Int, ShadeTree x)] + twigProximæ :: x -> ShadeTree x -> TwigEnviron x twigProximæ x₀ (DisjointBranches _ djbs) = Hask.foldMap (\(i₀,st) -> first (+i₀) <$> twigProximæ x₀ st) $ NE.zip ioffs djbs @@ -990,13 +1037,12 @@ = twigsaveTrim hither ct where Option (Just δxb) = x₀ .-~. xb hither (DBranch bdir (Hourglass bdc₁ bdc₂)) - | bdir<.>^δxb > 0 = twigProximæ x₀ bdc₁ - | otherwise = first (+nLeaves bdc₁) - <$> twigProximæ x₀ bdc₂ + = ((guard (overlap > -1)) >> twigProximæ x₀ bdc₁) + ++ ((guard (overlap < 1)) >> first (+nLeaves bdc₁)<$>twigProximæ x₀ bdc₂) + where overlap = bdir<.>^δxb twigProximæ _ plainLeaves = [(0, plainLeaves)] - twigsaveTrim :: (DBranch x -> [(Int,ShadeTree x)]) - -> ShadeTree x -> [(Int,ShadeTree x)] + twigsaveTrim :: (DBranch x -> TwigEnviron x) -> ShadeTree x -> TwigEnviron x twigsaveTrim f ct@(OverlappingBranches _ _ dbs) = case Hask.mapM (\(i₀,dbr) -> noLeaf $ first(+i₀)<$>f dbr) $ NE.zip ioffs dbs of @@ -1006,11 +1052,11 @@ noLeaf bqs = pure bqs ioffs = NE.scanl (\i -> (+i) . sum . fmap nLeaves . toList) 0 dbs - purgeRemotes :: ((Int,ShadeTree x), [(Int,ShadeTree x)]) - -> ((Int,ShadeTree x), [(Int,ShadeTree x)]) + purgeRemotes :: (Twig x, TwigEnviron x) -> (Twig x, TwigEnviron x) purgeRemotes = id -- See 7d1f3a4 for the implementation; this didn't work reliable. -completeTopShading :: (WithField ℝ Manifold x, WithField ℝ Manifold y) +completeTopShading :: ( WithField ℝ Manifold x, WithField ℝ Manifold y + , SimpleSpace (Needle x), SimpleSpace (Needle y) ) => x`Shaded`y -> [Shade' (x,y)] completeTopShading (PlainLeaves plvs) = pointsShade's $ (_topological &&& _untopological) <$> plvs @@ -1018,7 +1064,12 @@ = take 1 . completeTopShading =<< NE.toList bqs completeTopShading t = pointsCover's . map (_topological &&& _untopological) $ onlyLeaves t + +transferAsNormsDo :: LSpace v => Norm v -> Variance v -> v-+>v +transferAsNormsDo (Norm m) (Norm n) = n . m + flexTopShading :: ∀ x y f . ( WithField ℝ Manifold x, WithField ℝ Manifold y + , SimpleSpace (Needle x), SimpleSpace (Needle y) , Applicative f (->) (->) ) => (Shade' (x,y) -> f (x, (Shade' y, LocalLinear x y))) -> x`Shaded`y -> f (x`Shaded`y) @@ -1027,13 +1078,13 @@ 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₀ + where expa'₀ = dualNorm expa₀ j₀ :: LocalLinear x y - Option (Just j₀) = covariance expa'₀ - (_,expay₀) = factoriseMetric expa₀ + j₀ = dependence expa'₀ + (_,expay₀) = summandSpaceNorms expa₀ fts (xc, (Shade' yc expay, jtg)) = unsafeFmapLeaves applδj t where Option (Just δyc) = yc.-~.yc₀ - tfm = imitateMetricSpanChange expay₀ (recipMetric' expay) + tfm = transferAsNormsDo expay₀ (dualNorm expay) applδj (WithAny y x) = WithAny (yc₀ .+~^ ((tfm$δy) ^+^ (jtg$δx) ^+^ δyc)) x where Option (Just δx) = x.-~.xc @@ -1047,6 +1098,7 @@ assert_connected (PlainLeaves _) = () flexTwigsShading :: ∀ x y f . ( WithField ℝ Manifold x, WithField ℝ Manifold y + , SimpleSpace (Needle x), SimpleSpace (Needle y) , Hask.Applicative f ) => (Shade' (x,y) -> f (x, (Shade' y, LocalLinear x y))) -> x`Shaded`y -> f (x`Shaded`y) @@ -1083,6 +1135,7 @@ toInterior = pure translateP = Tagged (.+~^) (.+~^) = (.+^) + semimanifoldWitness = undefined instance (KnownNat n) => PseudoAffine (BaryCoords n) where (.-~.) = pure .: (.-.) @@ -1118,30 +1171,6 @@ --- startTriangulation :: forall n x . (KnownNat n, WithField ℝ Manifold x) --- => ISimplex n x -> TriangBuilder n x --- startTriangulation ispl@(ISimplex emb) = startWith $ fromISimplex ispl --- where startWith (ZeroSimplex p) = TriangVerticesSt [p] --- startWith s@(Simplex _ _) --- = TriangBuilder (Triangulation [s]) --- (splxVertices s) --- [ (s', expandInDir j) --- | j<-[0..n] --- | s' <- getTriangulation $ simplexFaces s ] --- where expandInDir j xs = case sortBy (comparing snd) $ filter ((> -1) . snd) xs_bc of --- ((x, q) : _) | q<0 -> pure x --- _ -> empty --- where xs_bc = map (\x -> (x, getBaryCoord (emb >-$ x) j)) xs --- (Tagged n) = theNatN :: Tagged n Int - --- extendTriangulation :: forall n x . (KnownNat n, WithField ℝ Manifold x) --- => [x] -> TriangBuilder n x -> TriangBuilder n x --- extendTriangulation xs (TriangBuilder tr tb te) = foldr tryex (TriangBuilder tr tb []) te --- where tryex (bspl, expd) (TriangBuilder (Triangulation tr') tb' te') --- | Option (Just fav) <- expd xs --- = let snew = Simplex fav bspl --- in TriangBuilder (Triangulation $ snew:tr') (fav:tb') undefined - bottomExtendSuitability :: (KnownNat n, WithField ℝ Manifold x) => ISimplex (S n) x -> x -> ℝ @@ -1158,17 +1187,7 @@ qs -> pure . fst . maximumBy (comparing snd) $ qs -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 = simplexBarycenter s - spread = init . map ((.-~.bc) >>> \(Option (Just v)) -> v) $ splxVertices s - embedding = case spanHilbertSubspace m spread of - (Option (Just e)) -> e - _ -> error "Trying to obtain simplexPlane from zero-volume\ - \ 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 @@ -1183,23 +1202,6 @@ Tagged n = theNatN :: Tagged n ℝ x' – x = case x'.-~.x of {Option(Just v)->v} -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 = simplexBarycenter s - (Embedding emb (DenseLinear prj)) - = simplexPlane m s - (r₀:rs) = [ prj HMat.#> asPackedVector v - | x <- splxVertices s, let (Option (Just v)) = x.-~.bc ] - tmat = HMat.inv $ HMat.fromColumns [ r - r₀ | r<-rs ] - toBrc x = case x.-~.bc of - Option (Just v) -> let rx = prj HMat.#> asPackedVector v - r₀ - in finalise $ tmat HMat.#> rx - finalise v = case freeVector $ HMat.toList v of - Option (Just bv) -> BaryCoords bv - fromBrc bccs = bc .+~^ (emb $ v) - where v = linearCombo $ (fromPackedVector r₀, b₀) : zip (fromPackedVector<$>rs) bs - (b₀:bs) = getBaryCoords' bccs fromISimplex :: forall x n . (KnownNat n, WithField ℝ Manifold x) => ISimplex n x -> Simplex n x @@ -1231,22 +1233,6 @@ doTriangBuild t = runIdentity (fst <$> doTriangT (unliftInTriangT (`evalStateT`mempty) t >> simplexITList >>= mapM lookSimplex)) -singleFullSimplex :: ∀ t n x . (KnownNat n, WithField ℝ Manifold x) - => ISimplex n x -> FullTriang t n x (SimplexIT t n x) -singleFullSimplex is = do - frame <- disjointSimplex (fromISimplex is) - lift . modify' $ Map.insert frame is - return frame - -fullOpenSimplex :: ∀ t n x . (KnownNat n, WithField ℝ Manifold x) - => Metric x -> Simplex (S n) x -> TriangBuild t n x [SimplexIT t n x] -fullOpenSimplex m s = do - let is = toISimplex m s - frame <- disjointSimplex (fromISimplex is) - fsides <- toList <$> lookSplxFacesIT frame - lift . forM (zip fsides $ iSimplexSideViews is) - $ \(fside,is') -> modify' $ Map.insert fside (m,is') - return fsides hypotheticalSimplexScore :: ∀ t n n' x . (KnownNat n', WithField ℝ Manifold x, n~S n') @@ -1268,90 +1254,15 @@ _ -> empty return . fmap sum $ Hask.sequence scores -spanSemiOpenSimplex :: ∀ t n n' x . (KnownNat n', WithField ℝ Manifold x, n~S n') - => SimplexIT t Z x -- ^ Tip of the desired simplex. - -> SimplexIT t n x -- ^ Base of the desired simplex. - -> TriangBuild t n x [SimplexIT t n x] - -- ^ Return the exposed faces of the new simplices. -spanSemiOpenSimplex p b = do - m <- lift $ fst <$> (Map.!b) <$> get - neighbours <- filterM isAdjacent =<< lookSupersimplicesIT p - let bs = b:|neighbours - frame <- webinateTriang p b - backSplx <- lookSimplex frame - let iSplx = toISimplex m backSplx - fsides <- toList <$> lookSplxFacesIT frame - let sviews = filter (not . (`elem`bs) . fst) $ zip fsides (iSimplexSideViews iSplx) - lift . forM sviews $ \(fside,is') -> modify' $ Map.insert fside (m,is') - lift . Hask.forM_ bs $ \fside -> modify' $ Map.delete fside - return $ fst <$> sviews - where isAdjacent = fmap (isJust . getOption) . sharedBoundary b -multiextendTriang :: ∀ t n n' x . (KnownNat n', WithField ℝ Manifold x, n~S n') - => [SimplexIT t Z x] -> TriangBuild t n x () -multiextendTriang vs = do - ps <- mapM lookVertexIT vs - sides <- lift $ Map.toList <$> get - forM_ sides $ \(f,(m,s)) -> - case optimalBottomExtension s ps of - Option (Just c) -> spanSemiOpenSimplex (vs !! c) f - _ -> return [] --- | BUGGY: this does connect the supplied triangulations, but it doesn't choose --- the right boundary simplices yet. Probable cause: inconsistent internal --- numbering of the subsimplices. -autoglueTriangulation :: ∀ t n n' n'' x - . (KnownNat n'', WithField ℝ Manifold x, n~S n', n'~S n'') - => (∀ t' . TriangBuild t' n' x ()) -> TriangBuild t n' x () -autoglueTriangulation tb = do - mbBounds <- Map.toList <$> lift get - mps <- pointsOfSurf mbBounds - - WriterT gbBounds <- liftInTriangT $ mixinTriangulation tb' - lift . forM_ gbBounds $ \(i,ms) -> do - modify' $ Map.insert i ms - gps <- pointsOfSurf gbBounds - - autoglue mps gbBounds - autoglue gps mbBounds - - where tb' :: ∀ s . TriangT s n x Identity - (WriterT (Metric x, ISimplex n x) [] (SimplexIT s n' x)) - tb' = unliftInTriangT (`evalStateT`mempty) $ - tb >> (WriterT . Map.toList) <$> lift get - - pointsOfSurf s = fnubConcatMap Hask.toList <$> forM s (lookSplxVerticesIT . fst) - - autoglue :: [SimplexIT t Z x] -> [(SimplexIT t n' x, (Metric x, ISimplex n x))] - -> TriangBuild t n' x () - autoglue vs sides = do - forM_ sides $ \(f,_) -> do - possibs <- forM vs $ \p -> fmap(p,) <$> hypotheticalSimplexScore p f - case catOptions possibs of - [] -> return () - qs -> do - spanSemiOpenSimplex (fst `id` maximumBy (comparing $ snd) qs) f - return () data AutoTriang n x where AutoTriang :: { getAutoTriang :: ∀ t . TriangBuild t n x () } -> AutoTriang (S n) x -instance (KnownNat n, WithField ℝ Manifold x) => Semigroup (AutoTriang (S (S n)) x) where - (<>) = autoTriangMappend -autoTriangMappend :: ∀ n n' n'' x . ( KnownNat n'', n ~ S n', n' ~ S n'' - , WithField ℝ Manifold x ) - => AutoTriang n x -> AutoTriang n x -> AutoTriang n x -AutoTriang a `autoTriangMappend` AutoTriang b = AutoTriang c - where c :: ∀ t . TriangBuild t n' x () - c = a >> autoglueTriangulation b -elementaryTriang :: ∀ n n' x . (KnownNat n', n~S n', WithField ℝ EuclidSpace x) - => Simplex n x -> AutoTriang n x -elementaryTriang t = AutoTriang (fullOpenSimplex m t >> return ()) - where m = euclideanMetric t - breakdownAutoTriang :: ∀ n n' x . (KnownNat n', n ~ S n') => AutoTriang n x -> [Simplex n x] breakdownAutoTriang (AutoTriang t) = doTriangBuild t @@ -1391,28 +1302,10 @@ --- triangulate :: forall x n . (KnownNat n, WithField ℝ Manifold x) --- => ShadeTree x -> Triangulation n x --- triangulate (DisjointBranches _ brs) --- = Triangulation $ Hask.foldMap (getTriangulation . triangulate) brs --- triangulate (PlainLeaves xs) = primitiveTriangulation xs --- triangBranches :: WithField ℝ Manifold x --- => ShadeTree x -> Branchwise x (Triangulation x) n --- triangBranches _ = undefined --- --- tringComplete :: WithField ℝ Manifold x --- => Triangulation x (n-1) -> Triangulation x n -> Triangulation x n --- tringComplete (Triangulation trr) (Triangulation tr) = undefined --- where --- bbSimplices = Map.fromList [(i, Left s) | s <- tr | i <- [0::Int ..] ] --- bbVertices = [(i, splxVertices s) | s <- tr | i <- [0::Int ..] ] --- - - -- | -- @ -- 'SimpleTree' x ≅ Maybe (x, 'Trees' x) @@ -1441,7 +1334,7 @@ deriving instance Show (c (x, GenericTree b b x)) => Show (GenericTree c b x) -- | Imitate the specialised 'ShadeTree' structure with a simpler, generic tree. -onlyNodes :: WithField ℝ Manifold x => ShadeTree x -> Trees x +onlyNodes :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) => ShadeTree x -> Trees x onlyNodes (PlainLeaves []) = GenericTree [] onlyNodes (PlainLeaves ps) = let (ctr,_) = pseudoECM $ NE.fromList ps in GenericTree [ (ctr, GenericTree $ (,mempty) <$> ps) ] @@ -1475,7 +1368,8 @@ mappend = (<>) -chainsaw :: WithField ℝ Manifold x => Cutplane x -> ShadeTree x -> Sawbones x +chainsaw :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => Cutplane x -> ShadeTree x -> Sawbones x chainsaw cpln (PlainLeaves xs) = Sawbones (sd1++) (sd2++) sd2 sd1 where (sd1,sd2) = partition (\x -> sideOfCut cpln x == Option(Just PositiveHalfSphere)) xs chainsaw cpln (DisjointBranches _ brs) = Hask.foldMap (chainsaw cpln) brs @@ -1490,7 +1384,7 @@ where shelter dpCutDist dq = case ptsDist dp dq of Option (Just d) -> d < abs dpCutDist _ -> False - ptsDist = fmap (metric $ recipMetric bexpa) .: (.-~.) + ptsDist = fmap (dualNorm bexpa|$|) .: (.-~.) fathomCD = fathomCutDistance cpln bexpa @@ -1510,7 +1404,7 @@ -- | Saw a tree into the domains covered by the respective branches of another tree. -sShSaw :: WithField ℝ Manifold x +sShSaw :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) => ShadeTree x -- ^ “Reference tree”, defines the cut regions. -- Must be at least one level of 'OverlappingBranches' deep. -> ShadeTree x -- ^ Tree to take the actual contents from. @@ -1551,7 +1445,7 @@ where shelter dpCutDist dq = case ptsDist dp dq of Option (Just d) -> d < abs dpCutDist _ -> False - ptsDist = fmap (metric $ recipMetric bexpa) .: (.-~.) + ptsDist = fmap (dualNorm bexpa|$|) .: (.-~.) fathomCD = fathomCutDistance cpl bexpa sShSaw _ _ = error "`sShSaw` is not supposed to cut anything else but `OverlappingBranches`" @@ -1566,7 +1460,7 @@ instance (NFData x, NFData y) => NFData (WithAny x y) -instance (Semimanifold x) => Semimanifold (x`WithAny`y) where +instance ∀ x y . (Semimanifold x) => Semimanifold (x`WithAny`y) where type Needle (WithAny x y) = Needle x type Interior (WithAny x y) = Interior x `WithAny` y WithAny y x .+~^ δx = WithAny y $ x.+~^δx @@ -1577,6 +1471,8 @@ (Interior x`WithAny`y -> Needle x -> Interior x`WithAny`y) tpWD = Tagged `id` \(WithAny y x) δx -> WithAny y $ tpx x δx where Tagged tpx = translateP :: Tagged x (Interior x -> Needle x -> Interior x) + semimanifoldWitness = case semimanifoldWitness :: SemimanifoldWitness x of + SemimanifoldWitness -> SemimanifoldWitness instance (PseudoAffine x) => PseudoAffine (x`WithAny`y) where WithAny _ x .-~. WithAny _ ξ = x.-~.ξ @@ -1624,7 +1520,8 @@ -- | This is to 'ShadeTree' as 'Data.Map.Map' is to 'Data.Set.Set'. type x`Shaded`y = ShadeTree (x`WithAny`y) -stiWithDensity :: (WithField ℝ Manifold x, WithField ℝ LinearManifold y) +stiWithDensity :: ( WithField ℝ Manifold x, WithField ℝ LinearManifold y + , SimpleSpace (Needle x) ) => x`Shaded`y -> x -> Cℝay y stiWithDensity (PlainLeaves lvs) | [locShape@(Shade baryc expa)] <- pointsShades $ _topological <$> lvs @@ -1641,12 +1538,12 @@ stiWithDensity (OverlappingBranches n (Shade (WithAny _ bc) extend) brs) = ovbSWD where ovbSWD x = case x .-~. bc of Option (Just v) - | dist² <- metricSq ε v + | dist² <- normSq ε v , dist² < 9 , att <- exp(1/(dist²-9)+1/9) -> qGather att $ fmap ($x) downPrepared _ -> coneTip - ε = recipMetric extend + ε = dualNorm extend downPrepared = dp =<< brs where dp (DBranch _ (Hourglass up dn)) = fmap stiWithDensity $ up:|[dn] @@ -1655,14 +1552,14 @@ where dens = sum (hParamCℝay <$> contribs) stiAsIntervalMapping :: (x ~ ℝ, y ~ ℝ) - => x`Shaded`y -> [(x, ((y, Diff y), Linear ℝ x y))] + => x`Shaded`y -> [(x, ((y, Diff y), LinearMap ℝ x y))] stiAsIntervalMapping = twigsWithEnvirons >=> pure.snd.fst >=> completeTopShading >=> pure. \(Shade' (xloc, yloc) shd) - -> ( xloc, ( (yloc, recip $ metric shd (0,1)) - , case covariance (recipMetric' shd) of - {Option(Just j)->j} ) ) + -> ( xloc, ( (yloc, recip $ shd|$|(0,1)) + , dependence (dualNorm shd) ) ) -smoothInterpolate :: (WithField ℝ Manifold x, WithField ℝ LinearManifold y) +smoothInterpolate :: ( WithField ℝ Manifold x, WithField ℝ LinearManifold y + , SimpleSpace (Needle x) ) => NonEmpty (x,y) -> x -> y smoothInterpolate l = \x -> case ltr x of @@ -1674,7 +1571,8 @@ ltr = stiWithDensity $ fromLeafPoints l' -spanShading :: ∀ x y . (WithField ℝ Manifold x, WithField ℝ Manifold y) +spanShading :: ∀ x y . ( WithField ℝ Manifold x, WithField ℝ Manifold y + , SimpleSpace (Needle x), SimpleSpace (Needle y) ) => (Shade x -> Shade y) -> ShadeTree x -> x`Shaded`y spanShading f = unsafeFmapTree addYs id addYSh where addYs :: NonEmpty x -> NonEmpty (x`WithAny`y) @@ -1683,7 +1581,7 @@ where [xsh@(Shade xmid _)] = pointsCovers $ toList l Shade ymid yexpa = f xsh yexamp = [ ymid .+~^ σ*^δy - | δy <- eigenSpan yexpa, σ <- [-1,1] ] + | δy <- normSpanningSystem yexpa, σ <- [-1,1] ] addYSh :: Shade x -> Shade (x`WithAny`y) addYSh xsh = shadeWithAny (_shadeCtr $ f xsh) xsh diff --git a/Data/Manifold/Types.hs b/Data/Manifold/Types.hs --- a/Data/Manifold/Types.hs +++ b/Data/Manifold/Types.hs @@ -26,7 +26,7 @@ {-# LANGUAGE PatternGuards #-} {-# LANGUAGE TypeOperators #-} {-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE RecordWildCards #-} +{-# LANGUAGE UnicodeSyntax #-} module Data.Manifold.Types ( @@ -37,10 +37,10 @@ , Disk1, Disk2, Cone, OpenCone -- * Linear manifolds , ZeroDim(..) - , ℝ⁰, ℝ, ℝ², ℝ³ + , ℝ, ℝ⁰, ℝ¹, ℝ², ℝ³, ℝ⁴ -- * Hyperspheres -- ** General form: Stiefel manifolds - , Stiefel1, stiefel1Project, stiefel1Embed + , Stiefel1(..), stiefel1Project, stiefel1Embed -- ** Specific examples , HasUnitSphere(..) , S⁰(..), S¹(..), S²(..) @@ -57,27 +57,29 @@ , Cutplane(..) , fathomCutDistance, sideOfCut, cutPosBetween -- * Linear mappings - , Linear, LocalLinear, denseLinear + , LinearMap, LocalLinear ) where import Data.VectorSpace +import Data.VectorSpace.Free import Data.AffineSpace import Data.MemoTrie (HasTrie(..)) import Data.Basis import Data.Fixed import Data.Tagged import Data.Semigroup -import qualified Numeric.LinearAlgebra.HMatrix as HMat import qualified Data.Vector.Generic as Arr import qualified Data.Vector +import qualified Data.Vector.Unboxed as UArr +import Data.List (maximumBy) +import Data.Ord (comparing) import Data.Manifold.Types.Primitive +import Data.Manifold.Types.Stiefel import Data.Manifold.PseudoAffine import Data.Manifold.Cone -import Data.LinearMap.HerMetric -import Data.VectorSpace.FiniteDimensional -import Data.LinearMap.Category (Linear, denseLinear) +import Math.LinearMap.Category import qualified Prelude @@ -86,6 +88,8 @@ import Control.Monad.Constrained import Data.Foldable.Constrained +import Data.Type.Coercion + #define deriveAffine(c,t) \ instance (c) => Semimanifold (t) where { \ type Needle (t) = Diff (t); \ @@ -97,127 +101,173 @@ a.-~.b = pure (a.-.b); } -newtype Stiefel1Needle v = Stiefel1Needle { getStiefel1Tangent :: HMat.Vector (Scalar v) } +newtype Stiefel1Needle v = Stiefel1Needle { getStiefel1Tangent :: UArr.Vector (Scalar v) } newtype Stiefel1Basis v = Stiefel1Basis { getStiefel1Basis :: Int } -s1bTrie :: forall v b. FiniteDimensional v => (Stiefel1Basis v->b) -> Stiefel1Basis v:->:b +s1bTrie :: ∀ v b. FiniteFreeSpace v => (Stiefel1Basis v->b) -> Stiefel1Basis v:->:b s1bTrie = \f -> St1BTrie $ fmap (f . Stiefel1Basis) allIs - where (Tagged d) = dimension :: Tagged v Int + where d = freeDimension ([]::[v]) allIs = Arr.fromList [0 .. d-2] -instance FiniteDimensional v => HasTrie (Stiefel1Basis v) where +instance FiniteFreeSpace v => HasTrie (Stiefel1Basis v) where data (Stiefel1Basis v :->: a) = St1BTrie ( Array a ) trie = s1bTrie; untrie (St1BTrie a) (Stiefel1Basis i) = a Arr.! i enumerate (St1BTrie a) = Arr.ifoldr (\i x l -> (Stiefel1Basis i,x):l) [] a type Array = Data.Vector.Vector -instance(SmoothScalar(Scalar v),FiniteDimensional v)=>AdditiveGroup(Stiefel1Needle v) where - Stiefel1Needle v ^+^ Stiefel1Needle w = Stiefel1Needle $ v + w - zeroV = s1nZ; negateV (Stiefel1Needle v) = Stiefel1Needle $ negate v -s1nZ :: forall v. FiniteDimensional v => Stiefel1Needle v -s1nZ=Stiefel1Needle .HMat.fromList$replicate(d-1)0 where(Tagged d)=dimension::Tagged v Int +instance (FiniteFreeSpace v, UArr.Unbox (Scalar v)) + => AdditiveGroup(Stiefel1Needle v) where + Stiefel1Needle v ^+^ Stiefel1Needle w = Stiefel1Needle $ uarrAdd v w + Stiefel1Needle v ^-^ Stiefel1Needle w = Stiefel1Needle $ uarrSubtract v w + zeroV = s1nZ; negateV (Stiefel1Needle v) = Stiefel1Needle $ UArr.map negate v -instance (SmoothScalar(Scalar v),FiniteDimensional v)=>VectorSpace(Stiefel1Needle v) where +uarrAdd :: (Num n, UArr.Unbox n) => UArr.Vector n -> UArr.Vector n -> UArr.Vector n +uarrAdd = UArr.zipWith (+) +uarrSubtract :: (Num n, UArr.Unbox n) => UArr.Vector n -> UArr.Vector n -> UArr.Vector n +uarrSubtract = UArr.zipWith (-) + +s1nZ :: ∀ v. (FiniteFreeSpace v, UArr.Unbox (Scalar v)) => Stiefel1Needle v +s1nZ = Stiefel1Needle . UArr.fromList $ replicate (d-1) 0 + where d = freeDimension ([]::[v]) + +instance (FiniteFreeSpace v, UArr.Unbox (Scalar v)) => VectorSpace (Stiefel1Needle v) where type Scalar (Stiefel1Needle v) = Scalar v - μ *^ Stiefel1Needle v = Stiefel1Needle $ HMat.scale μ v + μ *^ Stiefel1Needle v = Stiefel1Needle $ uarrScale μ v -instance (SmoothScalar (Scalar v), FiniteDimensional v)=>HasBasis (Stiefel1Needle v) where +uarrScale :: (Num n, UArr.Unbox n) => n -> UArr.Vector n -> UArr.Vector n +uarrScale μ = UArr.map (*μ) + +instance (FiniteFreeSpace v, UArr.Unbox (Scalar v)) => HasBasis (Stiefel1Needle v) where type Basis (Stiefel1Needle v) = Stiefel1Basis v basisValue = s1bV - decompose (Stiefel1Needle v) = zipWith ((,).Stiefel1Basis) [0..] $ HMat.toList v - decompose' (Stiefel1Needle v) (Stiefel1Basis i) = v HMat.! i -s1bV :: forall v b. FiniteDimensional v => Stiefel1Basis v -> Stiefel1Needle v + decompose (Stiefel1Needle v) = zipWith ((,).Stiefel1Basis) [0..] $ UArr.toList v + decompose' (Stiefel1Needle v) (Stiefel1Basis i) = v UArr.! i + +s1bV :: ∀ v b. (FiniteFreeSpace v, UArr.Unbox (Scalar v)) + => Stiefel1Basis v -> Stiefel1Needle v s1bV = \(Stiefel1Basis i) -> Stiefel1Needle - $ HMat.fromList [ if k==i then 1 else 0 | k<-[0..d-2] ] - where (Tagged d) = dimension :: Tagged v Int + $ UArr.fromList [ if k==i then 1 else 0 | k<-[0..d-2] ] + where d = freeDimension ([]::[v]) -instance (SmoothScalar (Scalar v), FiniteDimensional v) - => FiniteDimensional (Stiefel1Needle v) where - dimension = s1nD - basisIndex = Tagged $ \(Stiefel1Basis i) -> i - indexBasis = Tagged Stiefel1Basis - fromPackedVector = Stiefel1Needle - asPackedVector = getStiefel1Tangent -s1nD :: forall v. FiniteDimensional v => Tagged (Stiefel1Needle v) Int -s1nD = Tagged (d - 1) where (Tagged d) = dimension :: Tagged v Int +instance (FiniteFreeSpace v, UArr.Unbox (Scalar v)) + => FiniteFreeSpace (Stiefel1Needle v) where + freeDimension = s1nD + toFullUnboxVect = getStiefel1Tangent + unsafeFromFullUnboxVect = Stiefel1Needle +s1nD :: ∀ v p . FiniteFreeSpace v => p (Stiefel1Needle v) -> Int +s1nD _ = freeDimension ([]::[v]) - 1 -instance (SmoothScalar (Scalar v), FiniteDimensional v) - => AffineSpace (Stiefel1Needle v) where +instance (FiniteFreeSpace v, UArr.Unbox (Scalar v)) => AffineSpace (Stiefel1Needle v) where type Diff (Stiefel1Needle v) = Stiefel1Needle v (.+^) = (^+^) (.-.) = (^-^) -deriveAffine((SmoothScalar (Scalar v), FiniteDimensional v), Stiefel1Needle v) +deriveAffine((FiniteFreeSpace v, UArr.Unbox (Scalar v)), Stiefel1Needle v) -instance (MetricScalar (Scalar v), FiniteDimensional v) - => HasMetric' (Stiefel1Needle v) where - type DualSpace (Stiefel1Needle v) = Stiefel1Needle v - Stiefel1Needle v <.>^ Stiefel1Needle w = HMat.dot v w - functional = s1nF - doubleDual = id; doubleDual' = id -s1nF :: forall v. FiniteDimensional v => (Stiefel1Needle v->Scalar v)->Stiefel1Needle v -s1nF = \f -> Stiefel1Needle $ HMat.fromList [f $ basisValue b | b <- cb] - where (Tagged cb) = completeBasis :: Tagged (Stiefel1Needle v) [Stiefel1Basis v] +instance ∀ v . (FiniteFreeSpace v, UArr.Unbox (Scalar v)) + => TensorSpace (Stiefel1Needle v) where + type TensorProduct (Stiefel1Needle v) w = Array w + zeroTensor = Tensor $ Arr.replicate (freeDimension ([]::[v]) - 1) zeroV + toFlatTensor = LinearFunction $ Tensor . Arr.convert . getStiefel1Tangent + fromFlatTensor = LinearFunction $ Stiefel1Needle . Arr.convert . getTensorProduct + addTensors (Tensor a) (Tensor b) = Tensor $ Arr.zipWith (^+^) a b + scaleTensor = bilinearFunction $ \μ (Tensor a) -> Tensor $ Arr.map (μ*^) a + negateTensor = LinearFunction $ \(Tensor a) -> Tensor $ Arr.map negateV a + tensorProduct = bilinearFunction $ \(Stiefel1Needle n) w + -> Tensor $ Arr.map (*^w) $ Arr.convert n + transposeTensor = LinearFunction $ \(Tensor a) -> Arr.foldl' (^+^) zeroV + $ Arr.imap ( \i w -> (tensorProduct $ w) $ Stiefel1Needle + $ UArr.generate d (\j -> if i==j then 1 else 0) ) a + where d = freeDimension ([]::[v]) - 1 + fmapTensor = bilinearFunction $ \f (Tensor a) -> Tensor $ Arr.map (f$) a + fzipTensorWith = bilinearFunction $ \f (Tensor a, Tensor b) + -> Tensor $ Arr.zipWith (curry $ arr f) a b + coerceFmapTensorProduct _ Coercion = Coercion + +instance ∀ v . (FiniteFreeSpace v, UArr.Unbox (Scalar v), Num''' (Scalar v)) + => LinearSpace (Stiefel1Needle v) where + type DualVector (Stiefel1Needle v) = Stiefel1Needle v + linearId = LinearMap . Arr.generate d $ \i -> Stiefel1Needle . Arr.generate d $ + \j -> if i==j then 1 else 0 + where d = freeDimension ([]::[v]) - 1 + coerceDoubleDual = Coercion + blockVectSpan = LinearFunction $ \w -> Tensor . Arr.generate d + $ \i -> LinearMap . Arr.generate d + $ \j -> if i==j then w else zeroV + where d = freeDimension ([]::[v]) - 1 + blockVectSpan'= LinearFunction $ \w -> LinearMap . Arr.generate d + $ \i -> Tensor . Arr.generate d + $ \j -> if i==j then w else zeroV + where d = freeDimension ([]::[v]) - 1 + contractTensorMap = LinearFunction $ \(LinearMap m) + -> Arr.ifoldl' (\acc i (Tensor t) -> acc ^+^ t Arr.! i) zeroV m + contractMapTensor = LinearFunction $ \(Tensor m) + -> Arr.ifoldl' (\acc i (LinearMap t) -> acc ^+^ t Arr.! i) zeroV m + contractLinearMapAgainst = bilinearFunction $ \(LinearMap m) f + -> Arr.ifoldl' (\acc i w -> case f $ w of + Stiefel1Needle n -> n UArr.! i ) 0 m + applyDualVector = bilinearFunction $ \(Stiefel1Needle v) (Stiefel1Needle w) + -> UArr.sum $ UArr.zipWith (*) v w + applyLinear = bilinearFunction $ \(LinearMap m) (Stiefel1Needle v) + -> Arr.ifoldl' (\acc i w -> acc ^+^ v UArr.! i *^ w) zeroV m + composeLinear = bilinearFunction $ \f (LinearMap g) -> LinearMap $ Arr.map (f$) g -instance (WithField k LinearManifold v, Real k) => Semimanifold (Stiefel1 v) where +instance ( WithField k LinearManifold v, FiniteFreeSpace v, FiniteFreeSpace (DualVector v) + , RealFloat k, UArr.Unbox k + ) => Semimanifold (Stiefel1 v) where type Needle (Stiefel1 v) = Stiefel1Needle v fromInterior = id toInterior = pure translateP = Tagged (.+~^) - Stiefel1 s .+~^ Stiefel1Needle n = Stiefel1 . fromPackedVector . HMat.scale (signum s'i) + Stiefel1 s .+~^ Stiefel1Needle n = Stiefel1 . unsafeFromFullUnboxVect . uarrScale (signum s'i) $ if| ν==0 -> s' -- ν'≡0 is a special case of this, so we can otherwise assume ν'>0. - | ν<=2 -> let m = HMat.scale ιmν spro + HMat.scale ((1-abs ιmν)/ν') n + | ν<=2 -> let m = uarrScale ιmν spro `uarrAdd` uarrScale ((1-abs ιmν)/ν') n ιmν = 1-ν in insi ιmν m - | otherwise -> let m = HMat.scale ιmν spro + HMat.scale ((abs ιmν-1)/ν') n + | otherwise -> let m = uarrScale ιmν spro `uarrAdd` uarrScale ((abs ιmν-1)/ν') n ιmν = ν-3 in insi ιmν m - where d = HMat.size s' - s'= asPackedVector s + where d = UArr.length s' + s'= toFullUnboxVect s ν' = l2norm n quop = signum s'i / ν' ν = ν' `mod'` 4 - im = HMat.maxIndex $ HMat.cmap abs s' - s'i = s' HMat.! im - spro = let v = deli s' in HMat.scale (recip s'i) v + im = UArr.maxIndex $ UArr.map abs s' + s'i = s' UArr.! im + spro = let v = deli s' in uarrScale (recip s'i) v deli v = Arr.take im v Arr.++ Arr.drop (im+1) v insi ti v = Arr.generate d $ \i -> if | i v Arr.! i | i>im -> v Arr.! (i-1) | otherwise -> ti -instance (WithField k LinearManifold v, Real k) => PseudoAffine (Stiefel1 v) where - Stiefel1 s .-~. Stiefel1 t = pure . Stiefel1Needle $ case s' HMat.! im of - 0 -> HMat.scale (recip $ l2norm delis) delis - s'i | v <- HMat.scale (recip s'i) delis - tpro +instance ( WithField k LinearManifold v, FiniteFreeSpace v, FiniteFreeSpace (DualVector v) + , RealFloat k, UArr.Unbox k + ) => PseudoAffine (Stiefel1 v) where + Stiefel1 s .-~. Stiefel1 t = pure . Stiefel1Needle $ case s' UArr.! im of + 0 -> uarrScale (recip $ l2norm delis) delis + s'i | v <- uarrScale (recip s'i) delis `uarrSubtract` tpro , absv <- l2norm v , absv > 0 -> let μ = (signum (t'i/s'i) - recip(absv + 1)) / absv - in HMat.scale μ v + in uarrScale μ v | t'i/s'i > 0 -> samePoint | otherwise -> antipode - where d = HMat.size t' - s'= asPackedVector s; t' = asPackedVector t - im = HMat.maxIndex $ HMat.cmap abs t' - t'i = t' HMat.! im - tpro = let v = deli t' in HMat.scale (recip t'i) v + where d = UArr.length t' + s'= toFullUnboxVect s; t' = toFullUnboxVect t + im = UArr.maxIndex $ UArr.map abs t' + t'i = t' UArr.! im + tpro = let v = deli t' in uarrScale (recip t'i) v delis = deli s' deli v = Arr.take im v Arr.++ Arr.drop (im+1) v - samePoint = (d-1) HMat.|> repeat 0 - antipode = (d-1) HMat.|> (2 : repeat 0) + samePoint = UArr.replicate (d-1) 0 + antipode = (d-1) `UArr.fromListN` (2 : repeat 0) -instance ( WithField ℝ HilbertSpace x ) => ConeSemimfd (Stiefel1 x) where - type CℝayInterior (Stiefel1 x) = x - fromCℝayInterior (FinVecArrRep v) = case HMat.size v of - 0 -> Cℝay 0 $ Stiefel1 zeroV - _ -> Cℝay (HMat.norm_2 v) $ Stiefel1 (fromPackedVector v) - toCℝayInterior (Cℝay 0 _) = pure zeroV - toCℝayInterior (Cℝay l (Stiefel1 v)) - = pure.FinVecArrRep $ HMat.scale (l/HMat.norm_2 v') v' - where v' = asPackedVector v +-- instance ( WithField ℝ HilbertManifold x ) => ConeSemimfd (Stiefel1 x) where +-- type CℝayInterior (Stiefel1 x) = x -l2norm :: MetricScalar s => HMat.Vector s -> s -l2norm = realToFrac . HMat.norm_2 +l2norm :: (Floating s, UArr.Unbox s) => UArr.Vector s -> s +l2norm = sqrt . UArr.sum . UArr.map (^2) @@ -247,7 +297,7 @@ fathomCutDistance :: WithField ℝ Manifold x => Cutplane x -- ^ Hyperplane to measure the distance from. - -> HerMetric'(Needle x) -- ^ Metric to use for measuring that distance. + -> Metric' x -- ^ Metric to use for measuring that distance. -- This can only be accurate if the metric -- is valid both around the cut-plane's 'sawHandle', and -- around the points you measure. @@ -259,7 +309,7 @@ -- 'Nothing' if the point isn't reachable from the plane. fathomCutDistance (Cutplane sh (Stiefel1 cn)) met = \x -> fmap fathom $ x .-~. sh where fathom v = (cn <.>^ v) / scaleDist - scaleDist = metric' met cn + scaleDist = met|$|cn cutPosBetween :: WithField ℝ Manifold x => Cutplane x -> (x,x) -> Option D¹ @@ -270,7 +320,16 @@ | otherwise = empty -lineAsPlaneIntersection :: WithField ℝ Manifold x => Line x -> [Cutplane x] -lineAsPlaneIntersection (Line h dir) - = [Cutplane h nrml | nrml <- orthogonalComplementSpan [dir]] +lineAsPlaneIntersection :: + (WithField ℝ Manifold x, FiniteDimensional (Needle' x)) + => Line x -> [Cutplane x] +lineAsPlaneIntersection (Line h (Stiefel1 dir)) + = [ Cutplane h . Stiefel1 + $ candidate ^-^ worstCandidate ^* (overlap/worstOvlp) + | (i, (candidate, overlap)) <- zip [0..] $ zip candidates overlaps + , i /= worstId ] + where candidates = enumerateSubBasis entireBasis + overlaps = (<.>^dir) <$> candidates + (worstId, worstOvlp) = maximumBy (comparing $ abs . snd) $ zip [0..] overlaps + worstCandidate = candidates !! worstId diff --git a/Data/Manifold/Types/Primitive.hs b/Data/Manifold/Types/Primitive.hs --- a/Data/Manifold/Types/Primitive.hs +++ b/Data/Manifold/Types/Primitive.hs @@ -37,8 +37,8 @@ , Projective1, Projective2 , Disk1, Disk2, Cone, OpenCone -- * Linear manifolds - , ZeroDim(..), isoAttachZeroDim - , ℝ⁰, ℝ, ℝ², ℝ³ + , ZeroDim(..) + , ℝ, ℝ⁰, ℝ¹, ℝ², ℝ³, ℝ⁴ -- * Hyperspheres , S⁰(..), otherHalfSphere, S¹(..), S²(..) -- * Projective spaces @@ -57,12 +57,15 @@ import Data.VectorSpace +import Data.VectorSpace.Free +import Linear.V2 +import Linear.V3 +import Math.VectorSpace.ZeroDimensional import Data.AffineSpace import Data.Basis import Data.Void import Data.Monoid - -import qualified Numeric.LinearAlgebra.HMatrix as HMat +import Math.LinearMap.Category ((⊗)()) import Control.Applicative (Const(..), Alternative(..)) @@ -91,35 +94,7 @@ --- | A single point. Can be considered a zero-dimensional vector space, WRT any scalar. -data ZeroDim k = Origin deriving(Eq, Show) -instance Monoid (ZeroDim k) where - mempty = Origin - mappend Origin Origin = Origin -instance AffineSpace (ZeroDim k) where - type Diff (ZeroDim k) = ZeroDim k - Origin .+^ Origin = Origin - Origin .-. Origin = Origin -instance AdditiveGroup (ZeroDim k) where - zeroV = Origin - Origin ^+^ Origin = Origin - negateV Origin = Origin -instance VectorSpace (ZeroDim k) where - type Scalar (ZeroDim k) = k - _ *^ Origin = Origin -instance HasBasis (ZeroDim k) where - type Basis (ZeroDim k) = Void - basisValue = absurd - decompose Origin = [] - decompose' Origin = absurd -{-# INLINE isoAttachZeroDim #-} -isoAttachZeroDim :: ( WellPointed c, UnitObject c ~ (), ObjectPair c a () - , Object c (ZeroDim k), ObjectPair c a (ZeroDim k) - , PointObject c (ZeroDim k) ) - => Isomorphism c a (a, ZeroDim k) -isoAttachZeroDim = second (Isomorphism (const Origin) terminal) . attachUnit - -- | The zero-dimensional sphere is actually just two points. Implementation might -- therefore change to @ℝ⁰ 'Control.Category.Constrained.+' ℝ⁰@: the disjoint sum of two -- single-point spaces. @@ -188,11 +163,8 @@ --- | Dense tensor product of two vector spaces. -newtype x⊗y = DensTensProd { getDensTensProd :: HMat.Matrix (Scalar y) } - class NaturallyEmbedded m v where embed :: m -> v coEmbed :: v -> m @@ -211,16 +183,16 @@ coEmbed x | x>=0 = PositiveHalfSphere | otherwise = NegativeHalfSphere instance NaturallyEmbedded S¹ ℝ² where - embed (S¹ φ) = (cos φ, sin φ) - coEmbed (x,y) = S¹ $ atan2 y x + embed (S¹ φ) = V2 (cos φ) (sin φ) + coEmbed (V2 x y) = S¹ $ atan2 y x instance NaturallyEmbedded S² ℝ³ where - embed (S² ϑ φ) = ((cos φ * sin ϑ, sin φ * sin ϑ), cos ϑ) - coEmbed ((x,y),z) = S² (acos $ z/r) (atan2 y x) + embed (S² ϑ φ) = V3 (cos φ * sin ϑ) (sin φ * sin ϑ) (cos ϑ) + coEmbed (V3 x y z) = S² (acos $ z/r) (atan2 y x) where r = sqrt $ x^2 + y^2 + z^2 instance NaturallyEmbedded ℝP² ℝ³ where - embed (ℝP² r φ) = ((r * cos φ, r * sin φ), sqrt $ 1-r^2) - coEmbed ((x,y),z) = ℝP² (sqrt $ 1-(z/r)^2) (atan2 (y/r) (x/r)) + embed (ℝP² r φ) = V3 (r * cos φ) (r * sin φ) (sqrt $ 1-r^2) + coEmbed (V3 x y z) = ℝP² (sqrt $ 1-(z/r)^2) (atan2 (y/r) (x/r)) where r = sqrt $ x^2 + y^2 + z^2 instance NaturallyEmbedded D¹ ℝ where @@ -236,10 +208,12 @@ type Endomorphism a = a->a -type ℝ⁰ = ZeroDim ℝ type ℝ = Double -type ℝ² = (ℝ,ℝ) -type ℝ³ = (ℝ²,ℝ) +type ℝ⁰ = ZeroDim ℝ +type ℝ¹ = V1 ℝ +type ℝ² = V2 ℝ +type ℝ³ = V3 ℝ +type ℝ⁴ = V4 ℝ -- | Better known as ℝ⁺ (which is not a legal Haskell name), the ray diff --git a/Data/Manifold/Types/Stiefel.hs b/Data/Manifold/Types/Stiefel.hs new file mode 100644 --- /dev/null +++ b/Data/Manifold/Types/Stiefel.hs @@ -0,0 +1,48 @@ +-- | +-- Module : Data.Manifold.Types.Stiefel +-- Copyright : (c) Justus Sagemüller 2015 +-- License : GPL v3 +-- +-- Maintainer : (@) sagemueller $ geo.uni-koeln.de +-- Stability : experimental +-- Portability : portable +-- +-- Stiefel manifolds are a generalisation of the concept of the 'UnitSphere' +-- in real vector spaces. +-- The /n/-th Stiefel manifold is the space of all possible configurations of +-- /n/ orthonormal vectors. In the case /n/ = 1, simply a single normalised vector, +-- i.e. a vector on the unit sphere. +-- +-- Alternatively, the stiefel manifolds can be defined as quotient spaces under +-- scalings, and we prefer that definition since it doesn't require a notion of +-- unit length (which is only defined in inner-product spaces). + + + + +module Data.Manifold.Types.Stiefel where + + +import Data.Maybe +import qualified Data.Vector as Arr +import Data.Semigroup + +import Data.VectorSpace +import Data.AffineSpace +import Math.LinearMap.Category + +import Data.Manifold.Types.Primitive ((^), empty, embed, coEmbed) +import Data.Manifold.PseudoAffine + +import qualified Prelude as Hask hiding(foldl, sum, sequence) +import qualified Control.Applicative as Hask +import qualified Control.Monad as Hask hiding(forM_, sequence) + +import Control.Category.Constrained.Prelude hiding + ((^), all, elem, sum, forM, Foldable(..), Traversable) +import Control.Arrow.Constrained +import Control.Monad.Constrained hiding (forM) +import Data.Foldable.Constrained + + +newtype Stiefel1 v = Stiefel1 { getStiefel1N :: DualVector v } diff --git a/Data/Manifold/Web.hs b/Data/Manifold/Web.hs --- a/Data/Manifold/Web.hs +++ b/Data/Manifold/Web.hs @@ -66,7 +66,8 @@ import Control.DeepSeq import Data.VectorSpace -import Data.LinearMap.HerMetric +import Math.LinearMap.Category + import Data.Tagged import Data.Function (on) import Data.Fixed (mod') @@ -111,7 +112,7 @@ } deriving (Generic) -instance (NFData x, NFData (HerMetric (Needle x))) => NFData (Neighbourhood x) +instance (NFData x, NFData (Metric x)) => NFData (Neighbourhood x) -- | A 'PointsWeb' is almost, but not quite a mesh. It is a stongly connected† -- directed graph, backed by a tree for fast nearest-neighbour lookup of points. @@ -125,7 +126,7 @@ } -> PointsWeb x y deriving (Generic, Hask.Functor, Hask.Foldable, Hask.Traversable) -instance (NFData x, NFData (HerMetric (Needle x)), NFData (Needle' x), NFData y) => NFData (PointsWeb x y) +instance (NFData x, NFData (Metric x), NFData (Needle' x), NFData y) => NFData (PointsWeb x y) instance Foldable (PointsWeb x) (->) (->) where ffoldl = uncurry . Hask.foldl' . curry @@ -141,23 +142,23 @@ type MetricChoice x = Shade x -> Metric x -fromWebNodes :: ∀ x y . WithField ℝ Manifold x +fromWebNodes :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x)) => (MetricChoice x) -> [(x,y)] -> PointsWeb x y fromWebNodes mf = fromShaded mf . fromLeafPoints . map (uncurry WithAny . swap) -fromTopWebNodes :: ∀ x y . WithField ℝ Manifold x +fromTopWebNodes :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x)) => (MetricChoice x) -> [((x,[Needle x]),y)] -> PointsWeb x y fromTopWebNodes mf = fromTopShaded mf . fromLeafPoints . map (uncurry WithAny . swap . regroup') -fromShadeTree_auto :: ∀ x . WithField ℝ Manifold x => ShadeTree x -> PointsWeb x () -fromShadeTree_auto = fromShaded (recipMetric . _shadeExpanse) . constShaded () +fromShadeTree_auto :: ∀ x . (WithField ℝ Manifold x, SimpleSpace (Needle x)) => ShadeTree x -> PointsWeb x () +fromShadeTree_auto = fromShaded (dualNorm . _shadeExpanse) . constShaded () -fromShadeTree :: ∀ x . WithField ℝ Manifold x +fromShadeTree :: ∀ x . (WithField ℝ Manifold x, SimpleSpace (Needle x)) => (Shade x -> Metric x) -> ShadeTree x -> PointsWeb x () fromShadeTree mf = fromShaded mf . constShaded () -fromShaded :: ∀ x y . WithField ℝ Manifold x +fromShaded :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x)) => (MetricChoice x) -- ^ Local scalar-product generator. You can always -- use @'recipMetric' . '_shadeExpanse'@ (but this -- may give distortions compared to an actual @@ -166,7 +167,7 @@ -> PointsWeb x y fromShaded metricf = fromTopShaded metricf . fmapShaded ([],) -fromTopShaded :: ∀ x y . WithField ℝ Manifold x +fromTopShaded :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x)) => (MetricChoice x) -> (x`Shaded`([Needle x], y)) -- ^ Source tree, with a priori topology information -- (needles pointing to already-known neighbour candidates) @@ -203,7 +204,7 @@ oldNgbs <- get when (all (\(_,(_,nw)) -> visibleOverlap nw v) oldNgbs) `id`do let w = w₀ ^/ (w₀<.>^v) - where w₀ = toDualWith locRieM v + where w₀ = locRieM<$|v put $ (iNgb, (v,w)) : [ neighbour | neighbour@(_,(nv,_))<-oldNgbs @@ -226,7 +227,8 @@ ++ Hask.foldMap (onlyLeaves . snd) neighRegions of [sh₀] -> metricf sh₀ -indexWeb :: WithField ℝ Manifold x => PointsWeb x y -> WebNodeId -> Option (x,y) +indexWeb :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => PointsWeb x y -> WebNodeId -> Option (x,y) indexWeb (PointsWeb rsc assocD) i | i>=0, i WebNodeId -> y unsafeIndexWebData (PointsWeb _ asd) i = fst (asd Arr.! i) -webEdges :: ∀ x y . WithField ℝ Manifold x +webEdges :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x)) => PointsWeb x y -> [((x,y), (x,y))] webEdges web@(PointsWeb rsc assoc) = (lookId***lookId) <$> toList allEdges where allEdges :: Set.Set (WebNodeId,WebNodeId) @@ -269,7 +271,8 @@ -- | Fetch a point between any two neighbouring web nodes on opposite -- sides of the plane, and linearly interpolate the values onto the -- cut plane. -sliceWeb_lin :: ∀ x y . (WithField ℝ Manifold x, Geodesic x, Geodesic y) +sliceWeb_lin :: ∀ x y . ( WithField ℝ Manifold x, SimpleSpace (Needle x) + , Geodesic x, Geodesic y ) => PointsWeb x y -> Cutplane x -> [(x,y)] sliceWeb_lin web = sliceEdgs where edgs = webEdges web @@ -300,14 +303,16 @@ = GridSetup (x₀,y₀) [ GridPlanes (0,1) (0, (y₁-y₀)/fromIntegral ny) ny , GridPlanes (1,0) ((x₁-x₀)/fromIntegral nx, 0) ny ] -splitToGridLines :: (WithField ℝ Manifold x, Geodesic x, Geodesic y) +splitToGridLines :: ( WithField ℝ Manifold x, SimpleSpace (Needle x) + , Geodesic x, Geodesic y ) => PointsWeb x y -> GridSetup x -> [((x, GridPlanes x), [(x,y)])] splitToGridLines web (GridSetup x₀ [GridPlanes dirΩ spcΩ nΩ, linePln]) = [ ((x₀', linePln), sliceWeb_lin web $ Cutplane x₀' (Stiefel1 dirΩ)) | k <- [0 .. nΩ-1] , let x₀' = x₀.+~^(fromIntegral k *^ spcΩ) ] -sampleWebAlongGrid_lin :: ∀ x y . (WithField ℝ Manifold x, Geodesic x, Geodesic y) +sampleWebAlongGrid_lin :: ∀ x y . ( WithField ℝ Manifold x, SimpleSpace (Needle x) + , Geodesic x, Geodesic y ) => PointsWeb x y -> GridSetup x -> [(x,Option y)] sampleWebAlongGrid_lin web grid = finalLine =<< splitToGridLines web grid where finalLine :: ((x, GridPlanes x), [(x,y)]) -> [(x,Option y)] @@ -315,13 +320,13 @@ | length verts < 2 = take nSpl $ (,empty)<$>iterate (.+~^dir) x₀ finalLine ((x₀, GridPlanes _ dir nSpl), verts) = take nSpl $ go (x₀,0) intpseq where intpseq = mkInterpolationSeq_lin - [ (metric metr $ x.-~!x₀, y) | (x,y) <- verts ] + [ (metr |$| x.-~!x₀, y) | (x,y) <- verts ] go (x,_) [] = (,empty)<$>iterate (.+~^dir) x go xt (InterpolationIv (_,te) f:fs) = case break (( ((id***pure.f)<$>thisRange) ++ go xtn fs - Option (Just metr) = inferMetric $ webNodeRsc web + metr = inferMetric $ webNodeRsc web sampleWeb_2Dcartesian_lin :: (x~ℝ, y~ℝ, Geodesic z) => PointsWeb (x,y) z -> ((x,x),Int) -> ((y,y),Int) -> [(y,[(x,Option z)])] @@ -359,7 +364,7 @@ }, ngbH ) anyUnopposed rieM ngbCo = (`any`ngbCo) $ \(v,_) -> not $ (`any`ngbCo) $ \(v',_) - -> toDualWith rieM v <.>^ v' < 0 + -> (rieM<$|v) <.>^ v' < 0 localFocusWeb :: WithField ℝ Manifold x => PointsWeb x y -> PointsWeb x ((x,y), [(Needle x, y)]) @@ -376,16 +381,16 @@ ) asd' -nearestNeighbour :: WithField ℝ Manifold x +nearestNeighbour :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) => PointsWeb x y -> x -> Option (x,y) nearestNeighbour (PointsWeb rsc asd) x = fmap lkBest $ positionIndex empty rsc x where lkBest (iEst, (_, xEst)) = (xProx, yProx) where (iProx, (xProx, _)) = minimumBy (comparing $ snd . snd) - $ (iEst, (xEst, metricSq locMetr vEst)) + $ (iEst, (xEst, normSq locMetr vEst)) : neighbours (yProx, _) = asd Arr.! iProx (_, Neighbourhood neighbourIds locMetr) = asd Arr.! iEst - neighbours = [ (i, (xNgb, metricSq locMetr v)) + neighbours = [ (i, (xNgb, normSq locMetr v)) | i <- UArr.toList neighbourIds , let Right (_, xNgb) = indexShadeTree rsc i Option (Just v) = xNgb.-~.x @@ -417,7 +422,8 @@ -toGraph :: WithField ℝ Manifold x => PointsWeb x y -> (Graph, Vertex -> (x, y)) +toGraph :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => PointsWeb x y -> (Graph, Vertex -> (x, y)) toGraph wb = second (>>> \(i,_,_) -> case indexWeb wb i of {Option (Just xy) -> xy}) (graphFromEdges' edgs) where edgs :: [(Int, Int, [Int])] @@ -468,13 +474,15 @@ dupHead (x:|xs) = x:|x:xs -iterateFilterDEqn_static :: (WithField ℝ Manifold x, Refinable y) +iterateFilterDEqn_static :: ( WithField ℝ Manifold x, SimpleSpace (Needle x) + , Refinable y ) => DifferentialEqn x y -> PointsWeb x (Shade' y) -> [PointsWeb x (Shade' y)] iterateFilterDEqn_static f = map (fmap convexSetHull) . itWhileJust (filterDEqnSolutions_static f) . fmap (`ConvexSet`[]) -filterDEqnSolution_static :: (WithField ℝ Manifold x, Refinable y) +filterDEqnSolution_static :: ( WithField ℝ Manifold x, SimpleSpace (Needle x) + , Refinable y ) => DifferentialEqn x y -> PointsWeb x (Shade' y) -> Option (PointsWeb x (Shade' y)) filterDEqnSolution_static f = localFocusWeb >>> Hask.traverse `id` \((x,shy), ngbs) -> if null ngbs @@ -483,7 +491,8 @@ =<< intersectShade's ( propagateDEqnSolution_loc f ((x,shy), NE.fromList ngbs) ) -filterDEqnSolutions_static :: (WithField ℝ Manifold x, Refinable y) +filterDEqnSolutions_static :: ( WithField ℝ Manifold x, SimpleSpace (Needle x) + , Refinable y ) => DifferentialEqn x y -> PointsWeb x (ConvexSet y) -> Option (PointsWeb x (ConvexSet y)) filterDEqnSolutions_static f = localFocusWeb >>> Hask.traverse `id` \((x, shy@(ConvexSet hull _)), ngbs) -> if null ngbs @@ -516,7 +525,7 @@ filterDEqnSolutions_adaptive :: ∀ x y badness - . (WithField ℝ Manifold x, Refinable y, badness ~ ℝ) + . (WithField ℝ Manifold x, SimpleSpace (Needle x), Refinable y, badness ~ ℝ) => MetricChoice x -- ^ Scalar product on the domain, for regularising the web. -> DifferentialEqn x y -> (x -> Shade' y -> badness) @@ -540,7 +549,9 @@ smallBadnessGradient, largeBadnessGradient :: ℝ (smallBadnessGradient, largeBadnessGradient) = ( badnessGradRated!!(n`div`4), badnessGradRated!!(n*3`div`4) ) - where n = length badnessGradRated + where n = case length badnessGradRated of + 0 -> error "No neighbours available for badness-grading." + l -> l badnessGradRated = sort [ ngBad / bad | ( LocalWebInfo { _thisNodeData @@ -606,7 +617,7 @@ totalAge = maximum $ _solverNodeAge . _thisNodeData . fst <$> preproc'd errTgtModulation = (1-) . (`mod'`1) . negate . sqrt $ fromIntegral totalAge - badness x = badness' x . (shadeNarrowness %~ (^* errTgtModulation)) + badness x = badness' x . (shadeNarrowness %~ (scaleNorm errTgtModulation)) retraceBonds :: WebLocally x (WebLocally x (OldAndNew (x, SolverNodeState y))) -> [((x, [Needle x]), SolverNodeState y)] @@ -626,7 +637,7 @@ , (xN,_) <- oldAndNew nnWeb ] l -> [(xN.-~.x, ngb^.thisNodeId) | (xN,_) <- l] ] - possibleConflicts = [ metricSq locMetr v + possibleConflicts = [ normSq locMetr v | (v,nnId)<-neighbourCandidates , nnId > myId ] , isOld || null possibleConflicts @@ -634,14 +645,14 @@ ] where focused = oldAndNew' $ locWeb^.thisNodeData^.thisNodeData knownNgbs = snd <$> locWeb^.nodeNeighbours - oldMinDistSq = minimum [ metricSq locMetr vOld + oldMinDistSq = minimum [ normSq locMetr vOld | (_,ngb) <- knownNgbs , let Option (Just vOld) = ngb^.thisNodeCoord .-~. xOld ] -iterateFilterDEqn_adaptive :: (WithField ℝ Manifold x, Refinable y) +iterateFilterDEqn_adaptive :: (WithField ℝ Manifold x, SimpleSpace (Needle x), Refinable y) => MetricChoice x -- ^ Scalar product on the domain, for regularising the web. -> DifferentialEqn x y -> (x -> Shade' y -> ℝ) -- ^ Badness function for local results. diff --git a/Data/SimplicialComplex.hs b/Data/SimplicialComplex.hs --- a/Data/SimplicialComplex.hs +++ b/Data/SimplicialComplex.hs @@ -39,7 +39,6 @@ , simplexVertices, simplexVertices' -- * Simplicial complexes , Triangulation - , singleSimplex -- * Triangulation-builder monad , TriangT , evalTriangT, runTriangT, doTriangT, getTriang @@ -51,10 +50,7 @@ , distinctSimplices, NeighbouringSimplices -- ** Building triangulations , disjointTriangulation - , disjointSimplex , mixinTriangulation - , introVertToTriang - , webinateTriang -- * Misc util , HaskMonad, liftInTriangT, unliftInTriangT , Nat, Zero, One, Two, Three, Succ @@ -70,7 +66,7 @@ import Data.Semigroup import Data.Ord (comparing) -import Data.LinearMap.Category +import Math.LinearMap.Category import Data.Tagged import Data.Manifold.Types.Primitive ((^), empty) @@ -150,15 +146,6 @@ fmap f (TriangSkeleton sk vs) = TriangSkeleton (f<$>sk) vs deriving instance (Show x) => Show (Triangulation n x) --- | Consider a single simplex as a simplicial complex, consisting only of --- this simplex and its faces. -singleSimplex :: ∀ n x . KnownNat n => Simplex n x -> Triangulation n x -singleSimplex (ZS x) = TriangVertices $ pure (x, []) -singleSimplex (x :<| s) - = runIdentity . execTriangT insX $ TriangSkeleton (singleSimplex s) mempty - where insX :: ∀ t . TriangT t n x Identity () - insX = introVertToTriang x [SimplexIT 0] >> return() - nTopSplxs :: Triangulation n' x -> Int nTopSplxs (TriangVertices vs) = Arr.length vs nTopSplxs (TriangSkeleton _ vs) = Arr.length vs @@ -404,12 +391,7 @@ | k <- take (nTopSplxs t) [nTopSplxs tr ..] ] , tr <> t ) -disjointSimplex :: ∀ t m n x . (KnownNat n, HaskMonad m) - => Simplex n x -> TriangT t n x m (SimplexIT t n x) -disjointSimplex s = TriangT $ \tr -> return ( SimplexIT $ nTopSplxs tr - , tr <> singleSimplex s ) - -- | Import a triangulation like with 'disjointTriangulation', -- together with references to some of its subsimplices. mixinTriangulation :: ∀ t m f k n x . ( KnownNat n, KnownNat k @@ -425,59 +407,8 @@ t' = fmap (fmap tgetSimplexIT) t -webinateTriang :: ∀ t m n x . (HaskMonad m, KnownNat n) - => SimplexIT t Z x -> SimplexIT t n x -> TriangT t (S n) x m (SimplexIT t (S n) x) -webinateTriang ptt@(SimplexIT pt) bst@(SimplexIT bs) = do - existsReady <- lookupSimplexCone ptt bst - case existsReady of - Option (Just ext) -> return ext - _ -> TriangT $ \(TriangSkeleton sk cnn) - -> let resi = Arr.length cnn - res = SimplexIT $ Arr.length cnn :: SimplexIT t (S n) x - in case sk of - TriangVertices vs -> return - $ ( res - , TriangSkeleton (TriangVertices - $ vs Arr.// [ (pt, second (resi:) $ vs Arr.! pt) - , (bs, second (resi:) $ vs Arr.! bs) ] - ) $ Arr.snoc cnn (freeTuple$->$(pt, bs), []) ) - TriangSkeleton _ cnn' - -> let (cnbs,_) = cnn' Arr.! bs - in do (cnws,sk') <- unsafeRunTriangT ( do - cnws <- forM cnbs $ \j -> do - kt@(SimplexIT k) <- webinateTriang ptt (SimplexIT j) - addUplink' res kt - return k - addUplink' res bst - return cnws - ) sk - let snocer = (freeSnoc cnws bs, []) - return $ (res, TriangSkeleton sk' $ Arr.snoc cnn snocer) - where addUplink' :: SimplexIT t (S n) x -> SimplexIT t n x -> TriangT t n x m () - addUplink' (SimplexIT i) (SimplexIT j) = TriangT - $ \sk -> pure ((), case sk of - TriangVertices vs - -> let (v,ul) = vs Arr.! j - in TriangVertices $ vs Arr.// [(j, (v, i:ul))] - TriangSkeleton skd us - -> let (b,tl) = us Arr.! j - in TriangSkeleton skd $ us Arr.// [(j, (b, i:tl))] - ) - - -introVertToTriang :: ∀ t m n x . (HaskMonad m, KnownNat n) - => x -> [SimplexIT t n x] -> TriangT t (S n) x m (SimplexIT t Z x) -introVertToTriang v glues = do - j <- fmap (\(Option(Just k)) -> SimplexIT k) . onSkeleton . TriangT - $ return . tVertSnoc - mapM_ (webinateTriang j) glues - return j - where tVertSnoc :: Triangulation Z x -> (Int, Triangulation Z x) - tVertSnoc (TriangVertices vs) - = (Arr.length vs, TriangVertices $ vs `Arr.snoc` (v,[])) - diff --git a/Data/VectorSpace/FiniteDimensional.hs b/Data/VectorSpace/FiniteDimensional.hs deleted file mode 100644 --- a/Data/VectorSpace/FiniteDimensional.hs +++ /dev/null @@ -1,360 +0,0 @@ --- | --- Module : Data.VectorSpace.FiniteDimensional --- Copyright : (c) Justus Sagemüller 2015 --- License : GPL v3 --- --- Maintainer : (@) sagemueller $ geo.uni-koeln.de --- Stability : experimental --- Portability : portable --- -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE TupleSections #-} -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE PolyKinds #-} -{-# LANGUAGE UndecidableInstances #-} -{-# LANGUAGE StandaloneDeriving #-} -{-# LANGUAGE ConstraintKinds #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE UnicodeSyntax #-} - - - - -module Data.VectorSpace.FiniteDimensional ( - FiniteDimensional(..) - , SmoothScalar - , FinVecArrRep(..), concreteArrRep, (⊕), splitArrRep - ) where - - - - -import Prelude hiding ((^)) - -import Data.AffineSpace -import Data.VectorSpace -import Data.LinearMap -import Data.Basis -import Data.MemoTrie -import Data.Tagged -import Data.Void - -import Control.Applicative - -import Data.Manifold.Types.Primitive -import Data.CoNat -import Data.Embedding - -import Control.Arrow - -import qualified Data.Vector as Arr -import qualified Numeric.LinearAlgebra.HMatrix as HMat - - - - --- | Constraint that a space's scalars need to fulfill so it can be used for efficient linear algebra. --- Fulfilled pretty much only by the basic real and complex floating-point types. -type SmoothScalar s = ( VectorSpace s, HMat.Numeric s, HMat.Field s - , Num(HMat.Vector s), HMat.Indexable(HMat.Vector s)s - , HMat.Normed(HMat.Vector s) ) - --- | Many linear algebra operations are best implemented via packed, dense 'HMat.Matrix'es. --- For one thing, that makes common general vector operations quite efficient, --- in particular on high-dimensional spaces. --- More importantly, @hmatrix@ offers linear facilities --- such as inverse and eigenbasis transformations, which aren't available in the --- @vector-space@ library yet. But the classes from that library are strongly preferrable --- to plain matrices and arrays, conceptually. --- --- The 'FiniteDimensional' class is used to convert between both representations. --- It would be nice not to have the requirement of finite dimension on 'HerMetric', --- but it's probably not feasible to get rid of it in forseeable time. --- --- Instead of the run-time 'dimension' information, we would rather have a compile-time --- @type Dimension v :: Nat@, but type-level naturals are not mature enough yet. This --- will almost certainly change in the future. -class (HasBasis v, HasTrie (Basis v), SmoothScalar (Scalar v)) => FiniteDimensional v where - dimension :: Tagged v Int - basisIndex :: Tagged v (Basis v -> Int) - -- | Index must be in @[0 .. dimension-1]@, otherwise this is undefined. - indexBasis :: Tagged v (Int -> Basis v) - completeBasis :: Tagged v [Basis v] - completeBasis = liftA2 (\dim f -> f <$> [0 .. dim - 1]) dimension indexBasis - - completeBasisValues :: [v] - completeBasisValues = defCBVs - where defCBVs :: ∀ v . FiniteDimensional v => [v] - defCBVs = basisValue <$> cb - where Tagged cb = completeBasis :: Tagged v [Basis v] - - asPackedVector :: v -> HMat.Vector (Scalar v) - asPackedVector v = HMat.fromList $ snd <$> decompose v - - asPackedMatrix :: (FiniteDimensional w, Scalar w ~ Scalar v) - => (v :-* w) -> HMat.Matrix (Scalar v) - asPackedMatrix = defaultAsPackedMatrix - where defaultAsPackedMatrix :: forall v w s . - (FiniteDimensional v, FiniteDimensional w, s~Scalar v, s~Scalar w) - => (v :-* w) -> HMat.Matrix s - defaultAsPackedMatrix m = HMat.fromColumns $ asPackedVector . atBasis m <$> cb - where (Tagged cb) = completeBasis :: Tagged v [Basis v] - - fromPackedVector :: HMat.Vector (Scalar v) -> v - fromPackedVector v = result - where result = recompose $ zip cb (HMat.toList v) - cb = witness completeBasis result - - fromPackedMatrix :: (FiniteDimensional w, Scalar w ~ Scalar v) - => HMat.Matrix (Scalar v) -> (v :-* w) - fromPackedMatrix = defaultFromPackedMatrix - where defaultFromPackedMatrix :: forall v w s . - (FiniteDimensional v, FiniteDimensional w, s~Scalar v, s~Scalar w) - => HMat.Matrix s -> (v :-* w) - defaultFromPackedMatrix m = linear $ fromPackedVector . HMat.app m . asPackedVector - -instance (SmoothScalar k) => FiniteDimensional (ZeroDim k) where - dimension = Tagged 0 - basisIndex = Tagged absurd - indexBasis = Tagged $ const undefined - completeBasis = Tagged [] - asPackedVector Origin = HMat.fromList [] - fromPackedVector _ = Origin -instance FiniteDimensional ℝ where - dimension = Tagged 1 - basisIndex = Tagged $ \() -> 0 - indexBasis = Tagged $ \0 -> () - completeBasis = Tagged [()] - asPackedVector x = HMat.fromList [x] - asPackedMatrix f = HMat.asColumn . asPackedVector $ atBasis f () - fromPackedVector v = v HMat.! 0 -instance (FiniteDimensional a, FiniteDimensional b, Scalar a~Scalar b) - => FiniteDimensional (a,b) where - dimension = tupDim - where tupDim :: ∀ a b.(FiniteDimensional a,FiniteDimensional b)=>Tagged(a,b)Int - tupDim = Tagged $ da+db - where (Tagged da)=dimension::Tagged a Int; (Tagged db)=dimension::Tagged b Int - basisIndex = basId - where basId :: ∀ a b . (FiniteDimensional a, FiniteDimensional b) - => Tagged (a,b) (Either (Basis a) (Basis b) -> Int) - basId = Tagged basId' - where basId' (Left ba) = basIda ba - basId' (Right bb) = da + basIdb bb - (Tagged da) = dimension :: Tagged a Int - (Tagged basIda) = basisIndex :: Tagged a (Basis a->Int) - (Tagged basIdb) = basisIndex :: Tagged b (Basis b->Int) - indexBasis = basId - where basId :: ∀ a b . (FiniteDimensional a, FiniteDimensional b) - => Tagged (a,b) (Int -> Either (Basis a) (Basis b)) - basId = Tagged basId' - where basId' i | i < da = Left $ basIda i - | otherwise = Right . basIdb $ i - da - (Tagged da) = dimension :: Tagged a Int - (Tagged basIda) = indexBasis :: Tagged a (Int->Basis a) - (Tagged basIdb) = indexBasis :: Tagged b (Int->Basis b) - completeBasis = cb - where cb :: ∀ a b . (FiniteDimensional a, FiniteDimensional b) - => Tagged (a,b) [Either (Basis a) (Basis b)] - cb = Tagged $ map Left cba ++ map Right cbb - where (Tagged cba) = completeBasis :: Tagged a [Basis a] - (Tagged cbb) = completeBasis :: Tagged b [Basis b] - asPackedVector (a,b) = HMat.vjoin [asPackedVector a, asPackedVector b] - fromPackedVector = fPV - where fPV :: ∀ a b . (FiniteDimensional a, FiniteDimensional b, Scalar a~Scalar b) - => HMat.Vector (Scalar a) -> (a,b) - fPV v = (fromPackedVector l, fromPackedVector r) - where (Tagged da) = dimension :: Tagged a Int - (Tagged db) = dimension :: Tagged b Int - l = HMat.subVector 0 da v - r = HMat.subVector da db v - -instance (FiniteDimensional y, FiniteDimensional x) => AdditiveGroup (x⊗y) where - zeroV = DensTensProd $ (0 HMat.>< 0) [] - negateV (DensTensProd v) = DensTensProd $ negate v - DensTensProd v ^+^ DensTensProd w - | HMat.size v == (0,0) = DensTensProd w - | HMat.size w == (0,0) = DensTensProd v - | otherwise = DensTensProd $ v + w - -instance (FiniteDimensional y, FiniteDimensional x) => VectorSpace (x⊗y) where - type Scalar (x⊗y) = Scalar y - μ *^ DensTensProd v = DensTensProd $ HMat.scale μ v - -instance (FiniteDimensional y, FiniteDimensional x) => InnerSpace (x⊗y) where - DensTensProd v <.> DensTensProd w - | HMat.size v == (0,0) = 0 - | HMat.size w == (0,0) = 0 - | otherwise = HMat.flatten v `HMat.dot` HMat.flatten w - -instance (FiniteDimensional y, FiniteDimensional x) => HasBasis (x⊗y) where - type Basis (x⊗y) = (Basis x, Basis y) - basisValue = bvt - where bvt :: ∀ x y . (FiniteDimensional x, FiniteDimensional y) - => (Basis x, Basis y) -> x ⊗ y - bvt (bx,by) = DensTensProd $ HMat.assoc (nx,ny) 0 [((i,j),1)] - where Tagged nx = dimension :: Tagged x Int - Tagged ny = dimension :: Tagged y Int - Tagged i = ($bx) <$> basisIndex :: Tagged x Int - Tagged j = ($by) <$> basisIndex :: Tagged y Int - decompose = dct - where dct :: ∀ x y . (FiniteDimensional x, FiniteDimensional y) - => x ⊗ y -> [((Basis x, Basis y), Scalar y)] - dct (DensTensProd m) = zip [(i,j) | i <- cbx, j <- cby] - (HMat.toList $ HMat.flatten m) - where Tagged cbx = completeBasis :: Tagged x [Basis x] - Tagged cby = completeBasis :: Tagged y [Basis y] - decompose' = dct - where dct :: ∀ x y . (FiniteDimensional x, FiniteDimensional y) - => x ⊗ y -> (Basis x, Basis y) -> Scalar y - dct (DensTensProd m) (bi, bj) = m `HMat.atIndex` (bxi bi, byj bj) - where Tagged bxi = basisIndex :: Tagged x (Basis x -> Int) - Tagged byj = basisIndex :: Tagged y (Basis y -> Int) - -instance (FiniteDimensional a, FiniteDimensional b, Scalar a ~ Scalar b) - => FiniteDimensional (a⊗b) where - dimension = tensDim - where tensDim :: ∀ a b.(FiniteDimensional a,FiniteDimensional b)=>Tagged(a⊗b)Int - tensDim = Tagged $ da*db - where (Tagged da)=dimension::Tagged a Int; (Tagged db)=dimension::Tagged b Int - basisIndex = basId - where basId :: ∀ a b . (FiniteDimensional a, FiniteDimensional b) - => Tagged (a⊗b) ((Basis a, Basis b) -> Int) - basId = Tagged basId' - where basId' (ba,bb) = db*basIda ba + basIdb bb - (Tagged db) = dimension :: Tagged b Int - (Tagged basIda) = basisIndex :: Tagged a (Basis a->Int) - (Tagged basIdb) = basisIndex :: Tagged b (Basis b->Int) - indexBasis = basId - where basId :: ∀ a b . (FiniteDimensional a, FiniteDimensional b) - => Tagged (a⊗b) (Int -> (Basis a, Basis b)) - basId = Tagged basId' - where basId' i = let (ia,ib) = i`divMod`db - in (basIda ia, basIdb ib) - (Tagged db) = dimension :: Tagged b Int - (Tagged basIda) = indexBasis :: Tagged a (Int->Basis a) - (Tagged basIdb) = indexBasis :: Tagged b (Int->Basis b) - completeBasis = cb - where cb :: ∀ a b . (FiniteDimensional a, FiniteDimensional b) - => Tagged (a⊗b) [(Basis a, Basis b)] - cb = Tagged $ [(ba,bb) | ba<-cba, bb<-cbb] - where (Tagged cba) = completeBasis :: Tagged a [Basis a] - (Tagged cbb) = completeBasis :: Tagged b [Basis b] - asPackedVector (DensTensProd m) = HMat.flatten m - fromPackedVector = fPV - where fPV :: ∀ a b . (FiniteDimensional a, FiniteDimensional b, Scalar a~Scalar b) - => HMat.Vector (Scalar a) -> (a⊗b) - fPV v = DensTensProd $ HMat.reshape db v - where (Tagged db) = dimension :: Tagged b Int - - - - -instance (SmoothScalar x, KnownNat n) => FiniteDimensional (FreeVect n x) where - dimension = natTagPænultimate - basisIndex = Tagged getInRange - indexBasis = Tagged InRange - asPackedVector (FreeVect arr) = Arr.convert arr - fromPackedVector arr = FreeVect (Arr.convert arr) - -- asPackedMatrix = _ -- could be done quite efficiently here! - - - --- | Semantically the same as @'Tagged' tag refvs@, but directly uses the --- packed-vector array representation. --- --- The tag should really be kind-polymorphic, but at least GHC-7.8 doesn't quite --- handle the associated types of the manifold classes then. -newtype FinVecArrRep (tag :: * -> *) refvs scalar - = FinVecArrRep { getFinVecArrRep :: HMat.Vector scalar } - -instance (SmoothScalar s) => AffineSpace (FinVecArrRep t b s) where - type Diff (FinVecArrRep t b s) = FinVecArrRep t b s - (.-.) = (^-^) - (.+^) = (^+^) - -instance (SmoothScalar s) => AdditiveGroup (FinVecArrRep t b s) where - zeroV = FinVecArrRep $ 0 HMat.|> [] - negateV (FinVecArrRep v) = FinVecArrRep $ negate v - FinVecArrRep v ^+^ FinVecArrRep w - | HMat.size v == 0 = FinVecArrRep w - | HMat.size w == 0 = FinVecArrRep v - | otherwise = FinVecArrRep $ v + w - -instance (SmoothScalar s) => VectorSpace (FinVecArrRep t b s) where - type Scalar (FinVecArrRep t b s) = s - μ *^ FinVecArrRep v = FinVecArrRep $ HMat.scale μ v - -instance (SmoothScalar s) => InnerSpace (FinVecArrRep t b s) where - FinVecArrRep v <.> FinVecArrRep w - | HMat.size v == 0 = 0 - | HMat.size w == 0 = 0 - | otherwise = v`HMat.dot`w - -concreteArrRep :: (SmoothScalar s, FiniteDimensional r, Scalar r ~ s) - => Isomorphism (->) r (FinVecArrRep t r s) -concreteArrRep = Isomorphism (FinVecArrRep . asPackedVector) - (fromPackedVector . getFinVecArrRep) - -(⊕) :: ∀ t s v w . ( SmoothScalar s, FiniteDimensional v, FiniteDimensional w - , Scalar v ~ s, Scalar w ~ s ) - => FinVecArrRep t v s -> FinVecArrRep t w s -> FinVecArrRep t (v,w) s -FinVecArrRep v ⊕ FinVecArrRep w - | HMat.size v + HMat.size w == 0 = FinVecArrRep v - | HMat.size v == 0 = FinVecArrRep $ HMat.vjoin [HMat.konst 0 nv, w] - | HMat.size w == 0 = FinVecArrRep $ HMat.vjoin [v, HMat.konst 0 nw] - | otherwise = FinVecArrRep $ HMat.vjoin [v,w] - where Tagged nv = dimension :: Tagged v Int - Tagged nw = dimension :: Tagged w Int - -splitArrRep :: ∀ t s v w . ( SmoothScalar s, FiniteDimensional v, FiniteDimensional w - , Scalar v ~ s, Scalar w ~ s ) - => FinVecArrRep t (v,w) s -> (FinVecArrRep t v s, FinVecArrRep t w s) -splitArrRep (FinVecArrRep vw) - | HMat.size vw == 0 = (FinVecArrRep vw, FinVecArrRep vw) - | otherwise = ( FinVecArrRep $ HMat.subVector 0 nv vw - , FinVecArrRep $ HMat.subVector nv nw vw ) - where Tagged nv = dimension :: Tagged v Int - Tagged nw = dimension :: Tagged w Int - - -instance (SmoothScalar s, FiniteDimensional r, Scalar r ~ s) - => HasBasis (FinVecArrRep t r s) where - type Basis (FinVecArrRep t r s) = Basis r - basisValue = (concreteArrRep$->$) . basisValue - decompose = decompose . (concreteArrRep$<-$) - decompose' = decompose' . (concreteArrRep$<-$) - -instance (SmoothScalar s, FiniteDimensional r, Scalar r ~ s) - => FiniteDimensional (FinVecArrRep t r s) where - dimension = d - where d :: ∀ t r s . FiniteDimensional r => Tagged (FinVecArrRep t r s) Int - d = Tagged n - where Tagged n = dimension :: Tagged r Int - indexBasis = d - where d :: ∀ t r s . FiniteDimensional r => Tagged (FinVecArrRep t r s) (Int -> Basis r) - d = Tagged n - where Tagged n = indexBasis :: Tagged r (Int -> Basis r) - basisIndex = d - where d :: ∀ t r s . FiniteDimensional r => Tagged (FinVecArrRep t r s) (Basis r -> Int) - d = Tagged n - where Tagged n = basisIndex :: Tagged r (Basis r -> Int) - asPackedVector = apv - where apv :: ∀ t r s . (FiniteDimensional r, SmoothScalar s) - => FinVecArrRep t r s -> HMat.Vector s - apv (FinVecArrRep v) - | HMat.size v == 0 = HMat.konst 0 n - | otherwise = v - where Tagged n = dimension :: Tagged r Int - fromPackedVector = FinVecArrRep - - -instance (NaturallyEmbedded m r, FiniteDimensional r, s ~ Scalar r) - => NaturallyEmbedded m (FinVecArrRep t r s) where - embed = (concreteArrRep$<-$) . embed - coEmbed = coEmbed . (concreteArrRep$->$) - - diff --git a/images/examples/DiffableFunction-plots/Hann-window.png b/images/examples/DiffableFunction-plots/Hann-window.png new file mode 100644 Binary files /dev/null and b/images/examples/DiffableFunction-plots/Hann-window.png differ diff --git a/images/examples/DiffableFunction-plots/safe-sqrt.png b/images/examples/DiffableFunction-plots/safe-sqrt.png new file mode 100644 Binary files /dev/null and b/images/examples/DiffableFunction-plots/safe-sqrt.png differ diff --git a/images/examples/TreesAndWebs/2D-cartesiandisk.png b/images/examples/TreesAndWebs/2D-cartesiandisk.png new file mode 100644 Binary files /dev/null and b/images/examples/TreesAndWebs/2D-cartesiandisk.png differ diff --git a/images/examples/TreesAndWebs/2D-normaldistrib.png b/images/examples/TreesAndWebs/2D-normaldistrib.png new file mode 100644 Binary files /dev/null and b/images/examples/TreesAndWebs/2D-normaldistrib.png differ diff --git a/images/examples/TreesAndWebs/2D-scatter.png b/images/examples/TreesAndWebs/2D-scatter.png new file mode 100644 Binary files /dev/null and b/images/examples/TreesAndWebs/2D-scatter.png differ diff --git a/images/examples/TreesAndWebs/2D-scatter_twig-environs.png b/images/examples/TreesAndWebs/2D-scatter_twig-environs.png new file mode 100644 Binary files /dev/null and b/images/examples/TreesAndWebs/2D-scatter_twig-environs.png differ diff --git a/manifolds.cabal b/manifolds.cabal --- a/manifolds.cabal +++ b/manifolds.cabal @@ -1,5 +1,5 @@ Name: manifolds -Version: 0.2.3.0 +Version: 0.3.0.0 Category: Math Synopsis: Coordinate-free hypersurfaces Description: Manifolds, a generalisation of the notion of “smooth curves” or surfaces, @@ -31,6 +31,8 @@ Cabal-Version: >=1.10 Extra-Doc-Files: images/examples/*.png, images/examples/ShadeCombinations/2Dconvolution-skewed.png + images/examples/TreesAndWebs/*.png + images/examples/DiffableFunction-plots/*.png Source-Repository head type: git @@ -40,9 +42,11 @@ Build-Depends: base>=4.5 && < 6 , transformers , vector-space>=0.8 + , free-vector-spaces>=0.1.1 + , linear , MemoTrie , vector - , hmatrix >= 0.16 && < 0.18 + , linearmap-category > 0.1 && < 0.2 , containers , comonad , semigroups @@ -51,7 +55,7 @@ , deepseq , microlens >= 0.4 && <= 0.5, microlens-th , trivial-constraint >= 0.4 - , constrained-categories >= 0.2.3 && < 0.3 + , constrained-categories >= 0.2.3 && < 0.3.1 other-extensions: FlexibleInstances , TypeFamilies , FlexibleContexts @@ -69,9 +73,9 @@ Data.Manifold.Web Data.Manifold.DifferentialEquation Data.SimplicialComplex - Data.LinearMap.HerMetric Data.Function.Differentiable Data.Manifold.Types + Data.Manifold.Types.Stiefel Data.Manifold.Griddable Data.Manifold.Riemannian Other-modules: Data.List.FastNub @@ -80,10 +84,8 @@ Data.Manifold.Cone Data.CoNat Data.Embedding - Data.LinearMap.Category Data.Function.Differentiable.Data Data.Function.Affine - Data.VectorSpace.FiniteDimensional Control.Monad.Trans.OuterMaybe Util.Associate Util.LtdShow