diff --git a/Data/LinearMap/Category.hs b/Data/LinearMap/Category.hs --- a/Data/LinearMap/Category.hs +++ b/Data/LinearMap/Category.hs @@ -70,6 +70,12 @@ 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) diff --git a/Data/LinearMap/HerMetric.hs b/Data/LinearMap/HerMetric.hs --- a/Data/LinearMap/HerMetric.hs +++ b/Data/LinearMap/HerMetric.hs @@ -21,7 +21,7 @@ , toDualWith, fromDualWith , metricSq, metricSq', metric, metric', metrics, metrics' -- * Defining metrics - , projector, projector' + , projector, projector', projectors, projector's , euclideanMetric' -- * Metrics induce inner products , spanHilbertSubspace @@ -35,8 +35,10 @@ , transformMetric, transformMetric', dualCoCoProduct , dualiseMetric, dualiseMetric' , recipMetric, recipMetric' + -- ** Eigenvectors , eigenSpan, eigenSpan' , eigenCoSpan, eigenCoSpan' + , eigenSystem, HasEigenSystem, EigenVector , metriNormalise, metriNormalise' , metriScale', metriScale , adjoint @@ -55,6 +57,7 @@ , Stiefel1(..) , linMapAsTensProd, linMapFromTensProd , covariance + , outerProducts ) where @@ -166,15 +169,25 @@ -- 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 = matrixMetric $ HMat.outer uDecomp uDecomp - where uDecomp = asPackedVector u +projector u = HerMetric . pure $ u ⊗ u projector' :: HasMetric v => v -> HerMetric' v -projector' v = matrixMetric' $ HMat.outer vDecomp vDecomp - where vDecomp = asPackedVector 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 @@ -305,14 +318,17 @@ --- | The eigenbasis of a /positive definite/ metric, with each eigenvector scaled --- to the square root of the eigenvalue. +-- | 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, the original metric is obtained. (This holds even for --- non-Hilbert/Banach spaces, even though the concept of eigenbasis and --- “scaled length” doesn't really makes sense then in the usual way!) +-- 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 @@ -325,18 +341,117 @@ 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 = zipWith (HMat.scale . recip . sqrt) (HMat.toList μs) (HMat.toColumns vsm) + 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 = zipWith (HMat.scale . recip . sqrt) (HMat.toList μs) (HMat.toColumns vsm) + 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 (_, HerMetric Nothing) = (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... @@ -562,13 +677,11 @@ 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 = (sumV *** sumV) . unzip - $ (projector.fst &&& projector.snd) <$> eigenSpan' met +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 = (sumV *** sumV) . unzip - $ (projector'.fst &&& projector'.snd) <$> eigenSpan met +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) @@ -617,15 +730,16 @@ metricAsLength :: HerMetric ℝ -> ℝ metricAsLength m = case metricSq m 1 of - o | o > 0 -> recip o - | o < 0 -> error "Metric fails to be positive definite!" - | o == 0 -> error "Trying to use zero metric as length." + o | o > 0 -> sqrt $ recip o + | o < 0 -> error "Metric fails to be positive definite!" + | o == 0 -> error "Trying to use zero metric as length." + | otherwise -> error "Metric yields NaN." metricFromLength :: ℝ -> HerMetric ℝ metricFromLength = projector . recip metric'AsLength :: HerMetric' ℝ -> ℝ -metric'AsLength = recip . (`metric'`1) -- do we really want `recip` here? +metric'AsLength = sqrt . (`metric'`1) spanHilbertSubspace :: ∀ s v w @@ -678,7 +792,7 @@ | null eigSp = showString "zeroV" | otherwise = showParen (p>5) . foldr1 ((.) . (.(" ^+^ "++))) - $ ((("projector "++).).showsPrec 6)<$>eigSp + $ ((("projector "++).).showsPrec 10)<$>eigSp where eigSp = eigenSpan' m instance (HasMetric v, Scalar v ~ Double, Show v) => Show (HerMetric' v) where @@ -686,7 +800,7 @@ | null eigSp = showString "zeroV" | otherwise = showParen (p>5) . foldr1 ((.) . (.(" ^+^ "++))) - $ ((("projector' "++).).showsPrec 6)<$>eigSp + $ ((("projector' "++).).showsPrec 10)<$>eigSp where eigSp = eigenSpan m @@ -705,3 +819,27 @@ => 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/PseudoAffine.hs b/Data/Manifold/PseudoAffine.hs --- a/Data/Manifold/PseudoAffine.hs +++ b/Data/Manifold/PseudoAffine.hs @@ -32,6 +32,7 @@ {-# LANGUAGE FunctionalDependencies #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE LiberalTypeSynonyms #-} +{-# LANGUAGE DataKinds #-} {-# LANGUAGE GADTs #-} {-# LANGUAGE RankNTypes #-} {-# LANGUAGE TupleSections #-} @@ -64,7 +65,7 @@ -- ** Local functions , LocalLinear, LocalAffine -- * Misc - , palerp + , palerp, LocallyCoercible(..) ) where @@ -78,6 +79,7 @@ import Data.Fixed import Data.VectorSpace +import Data.Embedding import Data.LinearMap import Data.LinearMap.HerMetric import Data.LinearMap.Category @@ -219,6 +221,29 @@ class (PseudoAffine m, LinearManifold (Needle m), Interior m ~ m) => Manifold m instance (PseudoAffine m, LinearManifold (Needle m), Interior m ~ m) => Manifold m + + +-- | Instances of this class must be diffeomorphic manifolds, and even have +-- /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 ξ)) + => LocallyCoercible x ξ where + -- | Must be compatible with the canonical isomorphism on the tangent spaces, + -- i.e. + -- @ + -- locallyTrivialDiffeomorphism (p .+~^ 'fromPackedVector' v) + -- ≡ locallyTrivialDiffeomorphism p .+~^ 'fromPackedVector' v + -- @ + locallyTrivialDiffeomorphism :: x -> ξ + +instance LocallyCoercible ℝ ℝ where locallyTrivialDiffeomorphism = id +instance LocallyCoercible (ℝ,ℝ) (ℝ,ℝ) where locallyTrivialDiffeomorphism = id +instance LocallyCoercible (ℝ,(ℝ,ℝ)) (ℝ,(ℝ,ℝ)) where locallyTrivialDiffeomorphism = id +instance LocallyCoercible ((ℝ,ℝ),ℝ) ((ℝ,ℝ),ℝ) where locallyTrivialDiffeomorphism = id + + + type LocallyScalable s x = ( PseudoAffine x , HasMetric (Needle x) , s ~ Scalar (Needle x) ) @@ -268,8 +293,8 @@ type EuclidSpace x = ( AffineManifold x, InnerSpace (Diff x) , DualSpace (Diff x) ~ Diff x, Floating (Scalar (Diff x)) ) -euclideanMetric :: EuclidSpace x => Tagged x (Metric x) -euclideanMetric = Tagged euclideanMetric' +euclideanMetric :: EuclidSpace x => proxy x -> Metric x +euclideanMetric _ = euclideanMetric' -- | A co-needle can be understood as a “paper stack”, with which you can measure @@ -331,6 +356,14 @@ (.+~^) = (.+^) 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 Semimanifold (ZeroDim k) where @@ -360,6 +393,10 @@ Tagged tb = translateP :: Tagged b (Interior b -> Needle b -> Interior b) 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) => Semimanifold (a,b,c) where type Needle (a,b,c) = (Needle a, Needle b, Needle c) @@ -379,6 +416,18 @@ Tagged tc = translateP :: Tagged c (Interior c -> Needle c -> Interior c) 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 @@ -388,6 +437,11 @@ (.+~^) = (.+^) 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) + instance (HasMetric a, FiniteDimensional b, Scalar a~Scalar b) => Semimanifold (a⊗b) where type Needle (a⊗b) = a ⊗ b diff --git a/Data/Manifold/TreeCover.hs b/Data/Manifold/TreeCover.hs --- a/Data/Manifold/TreeCover.hs +++ b/Data/Manifold/TreeCover.hs @@ -22,9 +22,12 @@ {-# LANGUAGE RankNTypes #-} {-# LANGUAGE TupleSections #-} {-# LANGUAGE ParallelListComp #-} +{-# LANGUAGE MonadComprehensions #-} {-# LANGUAGE UnicodeSyntax #-} {-# LANGUAGE ConstraintKinds #-} {-# LANGUAGE PatternGuards #-} +{-# LANGUAGE PatternSynonyms #-} +{-# LANGUAGE ViewPatterns #-} {-# LANGUAGE LambdaCase #-} {-# LANGUAGE TypeOperators #-} {-# LANGUAGE ScopedTypeVariables #-} @@ -35,26 +38,27 @@ module Data.Manifold.TreeCover ( -- * Shades - Shade(..), Shade'(..), IsShade + Shade(..), pattern(:±), Shade'(..), (|±|), IsShade -- ** Lenses , shadeCtr, shadeExpanse, shadeNarrowness -- ** Construction - , fullShade, fullShade', pointsShades + , fullShade, fullShade', pointsShades, pointsCovers -- ** Evaluation , occlusion -- ** Misc - , factoriseShade, intersectShade's + , factoriseShade, intersectShade's, Refinable, refineShade', convolveShade', coerceShade -- * Shade trees - , ShadeTree(..), fromLeafPoints - -- * Simple view helpers - , onlyNodes, onlyLeaves + , ShadeTree(..), fromLeafPoints, onlyLeaves, indexShadeTree + -- * View helpers + , onlyNodes -- ** Auxiliary types , SimpleTree, Trees, NonEmptyTree, GenericTree(..) -- * Misc , sShSaw, chainsaw, HasFlatView(..), shadesMerge, smoothInterpolate , twigsWithEnvirons, completeTopShading, flexTwigsShading , WithAny(..), Shaded, stiAsIntervalMapping, spanShading - , DifferentialEqn, filterDEqnSolution_static + , constShaded, stripShadedUntopological + , DifferentialEqn, filterDEqnSolution_loc -- ** Triangulation-builders , TriangBuild, doTriangBuild, singleFullSimplex, autoglueTriangulation , AutoTriang, elementaryTriang, breakdownAutoTriang @@ -90,6 +94,7 @@ import Data.Manifold.PseudoAffine import Data.Function.Differentiable import Data.Function.Differentiable.Data +import Data.SetLike.Intersection import Data.Embedding import Data.CoNat @@ -158,6 +163,7 @@ factoriseShade :: ( Manifold x, RealDimension (Scalar (Needle x)) , Manifold y, RealDimension (Scalar (Needle y)) ) => shade (x,y) -> (shade x, shade y) + coerceShade :: (Manifold x, Manifold y, LocallyCoercible x y) => shade x -> shade y instance IsShade Shade where shadeCtr f (Shade c e) = fmap (`Shade`e) $ f c @@ -170,6 +176,8 @@ δinv = recipMetric δ 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) shadeExpanse :: Functor f (->) (->) => (Metric' x -> f (Metric' x)) -> Shade x -> f (Shade x) shadeExpanse f (Shade c e) = fmap (Shade c) $ f e @@ -184,6 +192,8 @@ _ -> 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) shadeNarrowness :: Functor f (->) (->) => (Metric x -> f (Metric x)) -> Shade' x -> f (Shade' x) shadeNarrowness f (Shade' c e) = fmap (Shade' c) $ f e @@ -202,6 +212,24 @@ fullShade' :: WithField ℝ Manifold x => x -> Metric x -> Shade' x fullShade' ctr expa = Shade' ctr expa + +-- | 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 + + +-- | Similar to ':±', but instead of expanding the shade, each vector /restricts/ it. +-- Iff these form a orthogonal basis (in whatever sense applicable), then both +-- methods will be equivalent. +-- +-- 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] + + + subshadeId' :: WithField ℝ Manifold x => x -> NonEmpty (Needle' x) -> x -> (Int, HourglassBulb) subshadeId' c expvs x = case x .-~. c of @@ -215,20 +243,40 @@ --- | Attempt to find a 'Shade' that “covers” the given points. +-- | Attempt to find a 'Shade' that describes the distribution of given points. -- At least in an affine space (and thus locally in any manifold), this can be used to -- estimate the parameters of a normal distribution from which some points were --- sampled. +-- sampled. Note that some points will be “outside” of the shade, +-- as happens for a normal distribution with some statistical likelyhood. +-- (Use 'pointsCovers' if you need to prevent that.) -- -- 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 list result. +-- Hence the result type is a list. pointsShades :: WithField ℝ Manifold x => [x] -> [Shade x] pointsShades = map snd . pointsShades' zeroV +-- | 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 + where guaranteeIn (ps, Shade x₀ ex) + = case ps >>= \p -> let Option (Just v) = p.-~.x₀ + in guard (metric ex' v > 1) >> [(p,projector' v)] + of [] -> Shade x₀ ex + outs -> guaranteeIn ( fst<$>outs + , Shade x₀ + $ ex ^+^ sumV (snd<$>outs) + ^/ fromIntegral (2 * length outs) ) + where ex' = recipMetric ex + pointsShade's :: WithField ℝ Manifold x => [x] -> [Shade' x] pointsShade's = map (\(Shade c e) -> Shade' c $ recipMetric e) . pointsShades +pointsCover's :: WithField ℝ Manifold x => [x] -> [Shade' x] +pointsCover's = map (\(Shade c e) -> Shade' c $ recipMetric e) . pointsCovers + pseudoECM :: WithField ℝ Manifold x => NonEmpty x -> (x, ([x],[x])) pseudoECM (p₀ NE.:| psr) = foldl' ( \(acc, (rb,nr)) (i,p) -> case p.-~.acc of @@ -245,7 +293,7 @@ _ -> pointsShades' minExt inc'd ++ pointsShades' minExt unreachable where (ctr,(inc'd,unreachable)) = pseudoECM $ NE.fromList ps - expa = ( (^+^minExt) . (^/ fromIntegral(length ps)) . sumV . map projector' ) + expa = ( (^+^minExt) . (^/ fromIntegral(length ps)) . projector's ) <$> mapM (.-~.ctr) ps @@ -273,7 +321,7 @@ = Just $ let cc = c₂ .+~^ v ^/ 2 Option (Just cv₁) = c₁.-~.cc Option (Just cv₂) = c₂.-~.cc - in Shade cc . sumV $ [e₁, e₂] ++ (projector'<$>[cv₁, cv₂]) + in Shade cc $ e₁ ^+^ e₂ ^+^ projector's [cv₁, cv₂] | otherwise = Nothing shadesMerge _ shs = shs @@ -365,21 +413,24 @@ uds = directionChoices hs traverseDirectionChoices :: (WithField ℝ Manifold x, Hask.Applicative f) - => ( (Needle' x, ShadeTree x) - -> [(Needle' x, ShadeTree x)] + => ( (Int, (Needle' x, ShadeTree x)) + -> [(Int, (Needle' x, ShadeTree x))] -> f (ShadeTree x) ) -> [DBranch x] -> f [DBranch x] -traverseDirectionChoices f dbs = td [] (dbs >>= - \(DBranch ѧ (Hourglass τ β)) - -> [(ѧ,τ), (negateV ѧ,β)]) - where td pds ((ѧ,t):(v,b):vds) +traverseDirectionChoices f dbs + = td [] . scanLeafNums 0 + $ dbs >>= \(DBranch ѧ (Hourglass τ β)) + -> [(ѧ,τ), (negateV ѧ,β)] + where td pds (ѧt@(_,(ѧ,_)):vb:vds) = liftA3 (\t' b' -> (DBranch ѧ (Hourglass t' b') :)) - (f (ѧ,t) $ pds++(v,b):uds) - (f (v,b) $ pds++(ѧ,t):uds) - $ td ((ѧ,t):(v,b):pds) vds + (f ѧt $ vb:uds) + (f vb $ ѧt:uds) + $ td (ѧt:vb:pds) vds where uds = pds ++ vds td _ _ = pure [] + scanLeafNums _ [] = [] + scanLeafNums i₀ ((v,t):vts) = (i₀, (v,t)) : scanLeafNums (i₀ + nLeaves t) vts instance (NFData x, NFData (Needle' x)) => NFData (ShadeTree x) where @@ -419,14 +470,34 @@ -- | Build a quite nicely balanced tree from a cloud of points, on any real manifold. -- --- Example: https://nbviewer.jupyter.org/github/leftaroundabout/manifolds/blob/master/test/generate-ShadeTrees.ipynb#pseudorandomCloudTree +-- 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 = fromLeafPoints' sShIdPartition +indexShadeTree :: ∀ x . WithField ℝ Manifold x + => ShadeTree x -> Int -> Either Int ([ShadeTree x], x) +indexShadeTree _ i + | i<0 = Left i +indexShadeTree sh@(PlainLeaves lvs) i = case length lvs of + n | i Right ([sh], lvs!!i) + | otherwise -> Left $ i-n +indexShadeTree (DisjointBranches n brs) i + | i (`indexShadeTree`i') + result -> return result + ) (Left i) brs + | otherwise = Left $ i-n +indexShadeTree sh@(OverlappingBranches n _ brs) i + | i foldl (\case + Left i' -> (`indexShadeTree`i') + result -> return result + ) (Left i) (toList brs>>=toList) + | otherwise = Left $ i-n + fromFnGraphPoints :: ∀ x y . (WithField ℝ Manifold x, WithField ℝ Manifold y) => [(x,y)] -> ShadeTree (x,y) fromFnGraphPoints = fromLeafPoints' fg_sShIdPart @@ -521,7 +592,7 @@ trunks :: ∀ x. WithField ℝ Manifold x => ShadeTree x -> [Shade x] -trunks (PlainLeaves lvs) = pointsShades lvs +trunks (PlainLeaves lvs) = pointsCovers lvs trunks (DisjointBranches _ brs) = Hask.foldMap trunks brs trunks (OverlappingBranches _ sh _) = [sh] @@ -558,11 +629,19 @@ in overlappingBranches (fs sh) brs' -intersectShade's :: ∀ y . WithField ℝ Manifold y => [Shade' y] -> Option (Shade' y) -intersectShade's [] = error "Global `Shade'` not implemented, so can't do intersection of zero co-shades." -intersectShade's (sh:shs) = Hask.foldrM inter2 sh shs - where inter2 :: Shade' y -> Shade' y -> Option (Shade' y) - inter2 (Shade' c e) (Shade' ζ η) +-- | 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 + -- | @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 | σ<-[0,1], v<-eigenCoSpan' ae ] + + -- | Specialised intersection operation. If @p@ is in @a@ and @b@, then it is + -- also in @refineShade' a b@. (The converse may not hold.) + refineShade' :: Shade' y -> Shade' y -> Option (Shade' y) + refineShade' (Shade' c e) (Shade' ζ η) | μe < 1 && μη < 1 = return $ Shade' iCtr iExpa | otherwise = empty where [c', ζ'] = [ ctr.+~^linearCombo @@ -588,42 +667,92 @@ μe = rcⰰ<.>^rc μη = rζⰰ<.>^rζ iExpa = (e^+^η)^/2 ^+^ projector rcⰰ^/(1-μe) ^+^ projector rζⰰ^/(1-μη) + + -- | If @p@ is in @a@ (red) and @δ@ is in @b@ (green), + -- then @p.+~^δ@ is in @convolveShade' a b@ (blue). + -- +-- Example: https://nbviewer.jupyter.org/github/leftaroundabout/manifolds/blob/master/test/ShadeCombinations.ipynb#shadeConvolutions +-- +-- <> + 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 + ζ = case filter (>0) corelap of + [] -> const 0 + nzrelap + -> let cre₁ = 1/minimum nzrelap + cre₂ = maximum nzrelap + edgeFactor = sqrt ( (1 + cre₁)^2 + (1 + cre₂)^2 ) + / (sqrt (1 + cre₁^2) + sqrt (1 + cre₂^2)) + in \case + 0 -> 0 + sq -> edgeFactor / (recip sq + 1) + +instance Refinable ℝ where + refineShade' (Shade' cl el) (Shade' cr er) + = case (metricSq el 1, metricSq er 1) of + (0, _) -> return $ Shade' cr er + (_, 0) -> return $ Shade' cl el + (ql,qr) | ql>0, qr>0 + -> let [rl,rr] = sqrt . recip <$> [ql,qr] + b = maximum $ zipWith (-) [cl,cr] [rl,rr] + t = minimum $ zipWith (+) [cl,cr] [rl,rr] + in guard (b> + let cm = (b+t)/2 + rm = (t-b)/2 + in return $ Shade' cm (projector $ recip rm) +-- convolveShade' (Shade' y₀ ey) (Shade' δ₀ eδ) +-- = case (metricSq ey 1, metricSq eδ 1) of +-- (wy,wδ) | wy>0, wδ>0 +-- -> Shade' (y₀.+~^δ₀) +-- ( projector . recip +-- $ recip (sqrt wy) + recip (sqrt wδ) ) +-- (_ , _) -> Shade' y₀ zeroV +instance (Refinable a, Refinable b) => Refinable (a,b) + + +intersectShade's :: ∀ y . Refinable y => NonEmpty (Shade' y) -> Option (Shade' y) +intersectShade's (sh:|shs) = Hask.foldrM refineShade' sh shs -type DifferentialEqn x y = Shade' (x,y) -> Shade' (LocalLinear x y) -filterDEqnSolution_loc :: ∀ x y . (WithField ℝ Manifold x, WithField ℝ Manifold y) - => DifferentialEqn x y -> (Shade' (x,y), [Shade' (x,y)]) - -> Option (Shade' y, LocalLinear x y) -filterDEqnSolution_loc f (shxy@(Shade' (x,y) expa), neighbours) = (,j₀) <$> yc + +type DifferentialEqn x y = Shade (x,y) -> Shade' (LocalLinear x y) + + +filterDEqnSolution_loc :: ∀ x y . (WithField ℝ Manifold x, Refinable y) + => DifferentialEqn x y -> ((x, Shade' y), NonEmpty (x, Shade' y)) + -> Option (Shade' y) +filterDEqnSolution_loc f ((x, shy@(Shade' y expay)), neighbours) = yc where jShade@(Shade' j₀ jExpa) = f shxy - marginδs :: [(Needle x, (Needle y, Metric y))] + [shxy] = pointsCovers [ (xs, ys') + | (xs, Shade' ys yse) <- (x,shy):NE.toList neighbours + , δy <- eigenCoSpan' yse + , ys' <- [ys.+~^δy, ys.-~^δy] ] + [Shade' _ expax] = pointsCover's $ x : (fst<$>NE.toList neighbours) + marginδs :: NonEmpty (Needle x, (Needle y, Metric y)) marginδs = [ (δxm, (δym, expany)) - | Shade' (xn, yn) expan <- neighbours - , let (Option (Just δx)) = x.-~.xn - (expanx, expany) = factoriseMetric expan - (Option (Just yc'n)) - = covariance $ recipMetric' expan - xntoMarg = metriNormalise expanx δx - (Option (Just δxm)) - = (xn .+~^ xntoMarg :: x) .-~. x - (Option (Just δym)) - = (yn .+~^ (yc'n $ xntoMarg) :: y - ) .-~. y + | (xn, Shade' yn expany) <- neighbours + , let (Option (Just δxm)) = xn.-~.x + (Option (Just δym)) = yn.-~.y ] back2Centre :: (Needle x, (Needle y, Metric y)) -> Shade' y back2Centre (δx, (δym, expany)) - = Shade' (y.+~^δyb) . recipMetric - $ recipMetric' expany - ^+^ recipMetric' (applyLinMapMetric jExpa δx') + = convolveShade' + (Shade' y expany) + (Shade' δyb $ applyLinMapMetric jExpa δx') where δyb = δym ^-^ (j₀ $ δx) δx' = toDualWith expax δx yc :: Option (Shade' y) yc = intersectShade's $ back2Centre <$> marginδs - (expax, expay) = factoriseMetric expa xSpan = eigenCoSpan' expax @@ -631,8 +760,8 @@ -- The simple list-yielding version (see rev. b4a427d59ec82889bab2fde39225b14a57b694df -- may well be more efficient than this version via a traversal. twigsWithEnvirons :: ∀ x. WithField ℝ Manifold x - => ShadeTree x -> [(ShadeTree x, [ShadeTree x])] -twigsWithEnvirons = execWriter . traverseTwigsWithEnvirons (writer . (fst&&&pure)) + => ShadeTree x -> [((Int, ShadeTree x), [(Int, ShadeTree x)])] +twigsWithEnvirons = execWriter . traverseTwigsWithEnvirons (writer . (snd.fst&&&pure)) data OuterMaybeT f a = OuterNothing | OuterJust (f a) deriving (Hask.Functor) instance (Hask.Applicative f) => Hask.Applicative (OuterMaybeT f) where @@ -642,91 +771,73 @@ traverseTwigsWithEnvirons :: ∀ x f . (WithField ℝ Manifold x, Hask.Applicative f) - => ((ShadeTree x, [ShadeTree x]) -> f (ShadeTree x)) + => ( ((Int, ShadeTree x), [(Int, ShadeTree x)]) -> f (ShadeTree x)) -> ShadeTree x -> f (ShadeTree x) -traverseTwigsWithEnvirons f = fst . go [] - where go :: [ShadeTree x] -> ShadeTree x -> (f (ShadeTree x), Bool) - go _ (DisjointBranches nlvs djbs) = ( fmap (DisjointBranches nlvs) - $ Hask.traverse (fst . go []) djbs - , False ) - go envi ct@(OverlappingBranches nlvs rob@(Shade robc _) brs) +traverseTwigsWithEnvirons f = fst . go [] . (0,) + where go :: [(Int, ShadeTree x)] -> (Int, ShadeTree x) + -> (f (ShadeTree x), Bool) + go _ (i₀, DisjointBranches nlvs djbs) = ( fmap (DisjointBranches nlvs) + . Hask.traverse (fst . go []) + $ NE.zip ioffs djbs + , False ) + where ioffs = NE.scanl (\i -> (+i) . nLeaves) i₀ djbs + go envi ct@(i₀, (OverlappingBranches nlvs rob@(Shade robc _) brs)) = ( case descentResult of OuterNothing -> f - $ purgeRemotes (ct, Hask.foldMap (twigProximæ robc) envi) + $ purgeRemotes + (ct, Hask.foldMap (\(io,te) + -> first (+io) <$> twigProximæ robc te) envi) OuterJust dR -> fmap (OverlappingBranches nlvs rob . NE.fromList) dR , False ) where descentResult = traverseDirectionChoices tdc $ NE.toList brs - tdc (vy, ty) alts = case go envi'' ty of + tdc (io, (vy, ty)) alts = case go envi'' (i₀+io, ty) of (_, True) -> OuterNothing (down, _) -> OuterJust down - where envi'' = filter (trunks >>> \(Shade ce _:_) + where envi'' = filter (snd >>> trunks >>> \(Shade ce _:_) -> let Option (Just δyenv) = ce.-~.robc qq = vy<.>^δyenv in qq > -1 && qq < 5 ) envi' - ++ map snd alts + ++ map ((+i₀)***snd) alts envi' = approach =<< envi - approach apt@(OverlappingBranches _ (Shade envc _) _) - = twigsaveTrim hither apt + approach (i₀e, apt@(OverlappingBranches _ (Shade envc _) _)) + = first (+i₀e) <$> twigsaveTrim hither apt where Option (Just δxenv) = robc .-~. envc hither (DBranch bdir (Hourglass bdc₁ bdc₂)) - | bdir<.>^δxenv > 0 = [bdc₁] - | otherwise = [bdc₂] + | bdir<.>^δxenv > 0 = [(0 , bdc₁)] + | otherwise = [(nLeaves bdc₁, bdc₂)] approach q = [q] - go envi plvs@(PlainLeaves _) = (f $ purgeRemotes (plvs, envi), True) + go envi plvs@(i₀, (PlainLeaves _)) + = (f $ purgeRemotes (plvs, envi), True) - twigProximæ :: x -> ShadeTree x -> [ShadeTree x] - twigProximæ x₀ (DisjointBranches _ djbs) = Hask.foldMap (twigProximæ x₀) djbs + twigProximæ :: x -> ShadeTree x -> [(Int, ShadeTree x)] + twigProximæ x₀ (DisjointBranches _ djbs) + = Hask.foldMap (\(i₀,st) -> first (+i₀) <$> twigProximæ x₀ st) + $ NE.zip ioffs djbs + where ioffs = NE.scanl (\i -> (+i) . nLeaves) 0 djbs twigProximæ x₀ ct@(OverlappingBranches _ (Shade xb qb) brs) = twigsaveTrim hither ct where Option (Just δxb) = x₀ .-~. xb hither (DBranch bdir (Hourglass bdc₁ bdc₂)) | bdir<.>^δxb > 0 = twigProximæ x₀ bdc₁ - | otherwise = twigProximæ x₀ bdc₂ - twigProximæ _ plainLeaves = [plainLeaves] + | otherwise = first (+nLeaves bdc₁) + <$> twigProximæ x₀ bdc₂ + twigProximæ _ plainLeaves = [(0, plainLeaves)] - twigsaveTrim :: (DBranch x -> [ShadeTree x]) - -> ShadeTree x -> [ShadeTree x] + twigsaveTrim :: (DBranch x -> [(Int,ShadeTree x)]) + -> ShadeTree x -> [(Int,ShadeTree x)] twigsaveTrim f ct@(OverlappingBranches _ _ dbs) - = case Hask.mapM (f >>> noLeaf) dbs of + = case Hask.mapM (\(i₀,dbr) -> noLeaf $ first(+i₀)<$>f dbr) + $ NE.zip ioffs dbs of Just pqe -> Hask.fold pqe - _ -> [ct] - where noLeaf [PlainLeaves _] = empty + _ -> [(0,ct)] + where noLeaf [(_,PlainLeaves _)] = empty noLeaf bqs = pure bqs + ioffs = NE.scanl (\i -> (+i) . sum . fmap nLeaves . toList) 0 dbs - purgeRemotes :: (ShadeTree x, [ShadeTree x]) -> (ShadeTree x, [ShadeTree x]) - purgeRemotes (ctm@(OverlappingBranches _ sm@(Shade xm _) _), candidates) - = (ctm, filter unobscured closeby) - where closeby = filter proximate candidates - proximate (OverlappingBranches _ sh@(Shade xh _) _) - = minusLogOcclusion sh xm * minusLogOcclusion sm xh - < 1024 -- = (2⋅4²)². The four-radius occlusion occurs - -- if two 𝑟-sized shades have just enough space - -- to fit another 𝑟-shade between them; then - -- we don't consider the shades neighbours - -- anymore. A factor √2 for the discrepancy - -- between standard deviation and max distance. - proximate _ = True - unobscured ht@(OverlappingBranches _ (Shade xh _) _) - = all (don'tObscure (xh, onlyLeaves ht)) closeby - don'tObscure (xh,lvsh) (OverlappingBranches _ sb@(Shade xb eb) _) - = vmc⋅vhc >= 0 || vm⋅vh >= 0 - where Option (Just vm) = pbm .-~. xb - Option (Just vh) = pbh .-~. xb - Option (Just vmc) = xm .-~. xb - Option (Just vhc) = xh .-~. xb - [pbm, pbh] = [ maximumBy (comparing $ \l -> - let Option (Just w) = l.-~.xb - in v⋅w ) lvs - | lvs <- [lvsm, lvsh] - | v <- [vhc, vmc] ] - (⋅) :: Needle x -> Needle x -> ℝ - v⋅w = toDualWith mb v <.>^ w - mb = recipMetric eb - don'tObscure _ _ = True - lvsm = onlyLeaves ctm - purgeRemotes xyz = xyz - + purgeRemotes :: ((Int,ShadeTree x), [(Int,ShadeTree x)]) + -> ((Int,ShadeTree x), [(Int,ShadeTree x)]) + purgeRemotes = id -- See 7d1f3a4 for the implementation; this didn't work reliable. completeTopShading :: (WithField ℝ Manifold x, WithField ℝ Manifold y) => x`Shaded`y -> [Shade' (x,y)] @@ -734,7 +845,7 @@ = pointsShade's $ (_topological &&& _untopological) <$> plvs completeTopShading (DisjointBranches _ bqs) = take 1 . completeTopShading =<< NE.toList bqs -completeTopShading t = pointsShade's . map (_topological &&& _untopological) $ onlyLeaves t +completeTopShading t = pointsCover's . map (_topological &&& _untopological) $ onlyLeaves t flexTopShading :: ∀ x y f . ( WithField ℝ Manifold x, WithField ℝ Manifold y , Applicative f (->) (->) ) @@ -769,25 +880,8 @@ => (Shade' (x,y) -> f (x, (Shade' y, LocalLinear x y))) -> x`Shaded`y -> f (x`Shaded`y) flexTwigsShading f = traverseTwigsWithEnvirons locFlex - where locFlex :: ∀ μ . (x`Shaded`y, μ) -> f (x`Shaded`y) - locFlex (lsh, _) = flexTopShading f lsh - -filterDEqnSolution_static :: ∀ x y . (WithField ℝ Manifold x, WithField ℝ Manifold y) - => DifferentialEqn x y - -> x`Shaded`y -> Option (x`Shaded`y) -filterDEqnSolution_static deq tr = traverseTwigsWithEnvirons locSoltn tr - where locSoltn :: (x`Shaded`y, [x`Shaded`y]) -> Option (x`Shaded`y) - locSoltn (local, environs) = do - let enviShades = completeTopShading =<< environs - flexed <- flexTopShading - (\oSh@(Shade' (ox,_) _) -> - (ox,) <$> filterDEqnSolution_loc deq (oSh, enviShades) - ) local - top'@(Shade' (top'x,_) top'exp) - <- intersectShade's $ completeTopShading =<< [local, flexed] - let (_, top'ySh) = factoriseShade top' - j' <- covariance $ recipMetric' top'exp - flexTopShading (const $ pure (top'x, (top'ySh, j'))) flexed + where locFlex :: ∀ μ . ((Int, x`Shaded`y), μ) -> f (x`Shaded`y) + locFlex ((_,lsh), _) = flexTopShading f lsh @@ -1085,7 +1179,7 @@ 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 (Tagged m) = euclideanMetric :: Tagged x (Metric x) + where m = euclideanMetric t breakdownAutoTriang :: ∀ n n' x . (KnownNat n', n ~ S n') => AutoTriang n x -> [Simplex n x] breakdownAutoTriang (AutoTriang t) = doTriangBuild t @@ -1345,6 +1439,12 @@ shadeWithoutAnything :: Shade (x`WithAny`y) -> Shade x shadeWithoutAnything (Shade (WithAny _ b) e) = Shade b e +constShaded :: y -> ShadeTree x -> x`Shaded`y +constShaded y = unsafeFmapTree (WithAny y<$>) id (shadeWithAny y) + +stripShadedUntopological :: x`Shaded`y -> ShadeTree x +stripShadedUntopological = unsafeFmapTree (fmap _topological) id shadeWithoutAnything + -- | This is to 'ShadeTree' as 'Data.Map.Map' is to 'Data.Set.Set'. type x`Shaded`y = ShadeTree (x`WithAny`y) @@ -1380,7 +1480,7 @@ stiAsIntervalMapping :: (x ~ ℝ, y ~ ℝ) => x`Shaded`y -> [(x, ((y, Diff y), Linear ℝ x y))] -stiAsIntervalMapping = twigsWithEnvirons >=> pure.fst >=> completeTopShading >=> pure. +stiAsIntervalMapping = twigsWithEnvirons >=> pure.snd.fst >=> completeTopShading >=> pure. \(Shade' (xloc, yloc) shd) -> ( xloc, ( (yloc, recip $ metric shd (0,1)) , case covariance (recipMetric' shd) of @@ -1404,7 +1504,7 @@ where addYs :: NonEmpty x -> NonEmpty (x`WithAny`y) addYs l = foldr (NE.<|) (fmap ( WithAny ymid) l ) (fmap (`WithAny`xmid) yexamp) - where [xsh@(Shade xmid _)] = pointsShades $ toList l + where [xsh@(Shade xmid _)] = pointsCovers $ toList l Shade ymid yexpa = f xsh yexamp = [ ymid .+~^ σ*^δy | δy <- eigenSpan yexpa, σ <- [-1,1] ] diff --git a/Data/Manifold/Types.hs b/Data/Manifold/Types.hs --- a/Data/Manifold/Types.hs +++ b/Data/Manifold/Types.hs @@ -54,7 +54,7 @@ , Cutplane(..) , fathomCutDistance, sideOfCut -- * Linear mappings - , Linear, denseLinear + , Linear, LocalLinear, denseLinear ) where diff --git a/Data/Manifold/Web.hs b/Data/Manifold/Web.hs new file mode 100644 --- /dev/null +++ b/Data/Manifold/Web.hs @@ -0,0 +1,219 @@ +-- | +-- Module : Data.Manifold.Web +-- Copyright : (c) Justus Sagemüller 2016 +-- 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 DeriveTraversable #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE TupleSections #-} +{-# LANGUAGE ParallelListComp #-} +{-# LANGUAGE UnicodeSyntax #-} +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE PatternGuards #-} +{-# LANGUAGE PatternSynonyms #-} +{-# LANGUAGE ViewPatterns #-} +{-# LANGUAGE LambdaCase #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE LiberalTypeSynonyms #-} +{-# LANGUAGE RecordWildCards #-} +{-# LANGUAGE DataKinds #-} + + +module Data.Manifold.Web where + + +import Data.List hiding (filter, all, elem, sum, foldr1) +import Data.Maybe +import qualified Data.Set as Set +import qualified Data.Vector as Arr +import qualified Data.Vector.Unboxed as UArr +import Data.List.NonEmpty (NonEmpty(..)) +import Data.List.FastNub +import qualified Data.List.NonEmpty as NE +import Data.Semigroup +import Data.Ord (comparing) +import Control.DeepSeq + +import Data.VectorSpace +import Data.AffineSpace +import Data.LinearMap +import Data.LinearMap.HerMetric +import Data.LinearMap.Category +import Data.AffineSpace +import Data.Basis +import Data.Complex hiding (magnitude) +import Data.Void +import Data.Tagged +import Data.Proxy + +import Data.SimplicialComplex +import Data.Manifold.Types +import Data.Manifold.Types.Primitive ((^), empty) +import Data.Manifold.PseudoAffine +import Data.Function.Differentiable +import Data.Function.Differentiable.Data +import Data.Manifold.TreeCover + +import Data.Embedding +import Data.CoNat + +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 Data.Functor.Identity +import Control.Monad.Trans.State +import Control.Monad.Trans.Writer +import Control.Monad.Trans.Maybe +import Control.Monad.Trans.Class +import qualified Data.Foldable as Hask +import Data.Foldable (all, elem, toList, sum, foldr1) +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 +import Control.Monad.Constrained hiding (forM) +import Data.Foldable.Constrained +import Data.Traversable.Constrained (Traversable, traverse) + +import GHC.Generics (Generic) + + +type WebNodeId = Int +type NeighbourRefs = UArr.Vector WebNodeId + +data PointsWeb :: * -> * -> * where + PointsWeb :: { + webNodeRsc :: ShadeTree x + , webNodeAssocData :: Arr.Vector (y, NeighbourRefs) + } -> PointsWeb x y + deriving (Generic, Hask.Functor, Hask.Foldable, Hask.Traversable) + +instance (NFData x, NFData (Needle' x), NFData y) => NFData (PointsWeb x y) + +instance Foldable (PointsWeb x) (->) (->) where + ffoldl = uncurry . Hask.foldl' . curry + foldMap = Hask.foldMap +instance Traversable (PointsWeb x) (PointsWeb x) (->) (->) where + traverse f (PointsWeb rsc asd) + = fmap (PointsWeb rsc . (`Arr.zip`ngss) . Arr.fromList) + . traverse f $ Arr.toList ys + where (ys,ngss) = Arr.unzip asd + + + +fromWebNodes :: ∀ x y . WithField ℝ Manifold x + => (Shade x->Metric x) -> [(x,y)] -> PointsWeb x y +fromWebNodes mf = fromShaded mf . fromLeafPoints . map (uncurry WithAny . swap) + +fromShadeTree_auto :: ∀ x . WithField ℝ Manifold x => ShadeTree x -> PointsWeb x () +fromShadeTree_auto = fromShaded (recipMetric . _shadeExpanse) . constShaded () + +fromShadeTree :: ∀ x . WithField ℝ Manifold x + => (Shade x -> Metric x) -> ShadeTree x -> PointsWeb x () +fromShadeTree mf = fromShaded mf . constShaded () + +fromShaded :: ∀ x y . WithField ℝ Manifold x + => (Shade x -> Metric x) -- ^ Local scalar-product generator. You can always + -- use @'recipMetric' . '_shadeExpanse'@ (but this + -- may give distortions compared to an actual + -- Riemannian metric). + -> (x`Shaded`y) -- ^ Source tree. + -> PointsWeb x y +fromShaded metricf shd = PointsWeb shd' assocData + where shd' = stripShadedUntopological shd + assocData = Hask.foldMap locMesh $ twigsWithEnvirons shd + + locMesh :: ((Int, ShadeTree (x`WithAny`y)), [(Int, ShadeTree (x`WithAny`y))]) + -> Arr.Vector (y, NeighbourRefs) + locMesh ((i₀, locT), neighRegions) = Arr.map findNeighbours locLeaves + where locLeaves = Arr.map (first (+i₀)) . Arr.indexed . Arr.fromList + $ onlyLeaves locT + vicinityLeaves = Hask.foldMap + (\(i₀n, ngbR) -> Arr.map (first (+i₀n)) + . Arr.indexed + . Arr.fromList + $ onlyLeaves ngbR + ) neighRegions + findNeighbours :: (Int, x`WithAny`y) -> (y, NeighbourRefs) + findNeighbours (i, WithAny y x) + = (y, UArr.fromList $ fst<$>execState seek mempty) + where seek = do + Hask.forM_ (locLeaves Arr.++ vicinityLeaves) + $ \(iNgb, WithAny _ xNgb) -> + when (iNgb/=i) `id`do + let (Option (Just v)) = xNgb.-~.x + oldNgbs <- get + when (all (\(_,(_,nw)) -> visibleOverlap nw v) oldNgbs) `id`do + let w = w₀ ^/ (w₀<.>^v) + where w₀ = toDualWith locRieM v + put $ (iNgb, (v,w)) + : [ neighbour + | neighbour@(_,(nv,_))<-oldNgbs + , visibleOverlap w nv + ] + + visibleOverlap :: Needle' x -> Needle x -> Bool + visibleOverlap w v = o < 1 + where o = w<.>^v + + locRieM :: Metric x + locRieM = case pointsCovers . map _topological + $ onlyLeaves locT + ++ Hask.foldMap (onlyLeaves . snd) neighRegions of + [sh₀] -> metricf sh₀ + +indexWeb :: WithField ℝ Manifold x => PointsWeb x y -> WebNodeId -> Option (x,y) +indexWeb (PointsWeb rsc assocD) i + | i>=0, i PointsWeb x y -> [((x,y), (x,y))] +webEdges web@(PointsWeb rsc assoc) = (lookId***lookId) <$> toList allEdges + where allEdges :: Set.Set (WebNodeId,WebNodeId) + allEdges = Hask.foldMap (\(i,(_,ngbs)) + -> Set.fromList [(min i i', max i i') + | i'<-UArr.toList ngbs ] + ) $ Arr.indexed assoc + lookId i | Option (Just xy) <- indexWeb web i = xy + + +localFocusWeb :: WithField ℝ Manifold x => PointsWeb x y -> PointsWeb x ((x,y), [(x,y)]) +localFocusWeb (PointsWeb rsc asd) = PointsWeb rsc asd'' + where asd' = Arr.imap (\i (y,n) -> case indexShadeTree rsc i of + Right (_,x) -> ((x,y),n) ) asd + asd''= Arr.map (\(xy,n) -> + ((xy, [fst (asd' Arr.! j) | j<-UArr.toList n]), n) + ) asd' + + +filterDEqnSolution_static :: (WithField ℝ Manifold 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 + then pure shy + else refineShade' shy + =<< filterDEqnSolution_loc f ((x,shy), NE.fromList ngbs) + + diff --git a/Data/SetLike/Intersection.hs b/Data/SetLike/Intersection.hs new file mode 100644 --- /dev/null +++ b/Data/SetLike/Intersection.hs @@ -0,0 +1,37 @@ + +-- | +-- Module : Data.SetLike.Intersection +-- Copyright : (c) Justus Sagemüller 2016 +-- License : GPL v3 +-- +-- Maintainer : (@) sagemueller $ geo.uni-koeln.de +-- Stability : experimental +-- Portability : portable +-- + + +module Data.SetLike.Intersection where + +import Data.Semigroup +import qualified Data.List.NonEmpty as NE +import Data.List.NonEmpty (NonEmpty(..)) + + +newtype IntersectT s x = IntersectT { getIntersectors :: NonEmpty (s x) } + + +singleIntersect :: s x -> IntersectT s x +singleIntersect = IntersectT . pure + +rmTautologyIntersect :: + (s x -> s x -> Option (s x)) -- ^ Subset-finder + -> IntersectT s x -> IntersectT s x +rmTautologyIntersect smaller (IntersectT isoa) = IntersectT $ rti isoa + where rti (s₀:|ss) = reduce [] ss + where reduce [] [] = s₀:|[] + reduce (sp₀:sp) [] = NE.cons s₀ $ rti (sp₀:|sp) + reduce sp (s₁:sr) = case smaller s₀ s₁ of + Option (Just si) -> rti $ si :| (sp ++ sr) + Option Nothing -> reduce (s₁:sp) sr + + diff --git a/Data/VectorSpace/FiniteDimensional.hs b/Data/VectorSpace/FiniteDimensional.hs --- a/Data/VectorSpace/FiniteDimensional.hs +++ b/Data/VectorSpace/FiniteDimensional.hs @@ -86,6 +86,12 @@ 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 diff --git a/images/examples/ShadeCombinations/2Dconvolution-skewed.png b/images/examples/ShadeCombinations/2Dconvolution-skewed.png new file mode 100644 Binary files /dev/null and b/images/examples/ShadeCombinations/2Dconvolution-skewed.png differ diff --git a/images/examples/cartesiandisk-2d-PointsWeb.png b/images/examples/cartesiandisk-2d-PointsWeb.png new file mode 100644 Binary files /dev/null and b/images/examples/cartesiandisk-2d-PointsWeb.png differ diff --git a/images/examples/normaldistrib-2d-PointsWeb.png b/images/examples/normaldistrib-2d-PointsWeb.png new file mode 100644 Binary files /dev/null and b/images/examples/normaldistrib-2d-PointsWeb.png differ diff --git a/images/examples/normaldistrib-2d-ShadeTree.png b/images/examples/normaldistrib-2d-ShadeTree.png new file mode 100644 Binary files /dev/null and b/images/examples/normaldistrib-2d-ShadeTree.png differ diff --git a/images/examples/simple-2d-PointsWeb.png b/images/examples/simple-2d-PointsWeb.png new file mode 100644 Binary files /dev/null and b/images/examples/simple-2d-PointsWeb.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.0.1 +Version: 0.2.2.0 Category: Math Synopsis: Coordinate-free hypersurfaces Description: Manifolds, a generalisation of the notion of “smooth curves” or surfaces, @@ -29,7 +29,8 @@ Maintainer: (@) sagemueller $ geo.uni-koeln.de Build-Type: Simple Cabal-Version: >=1.10 -Extra-Doc-Files: images/examples/*.png +Extra-Doc-Files: images/examples/*.png, + images/examples/ShadeCombinations/2Dconvolution-skewed.png Source-Repository head type: git @@ -64,6 +65,7 @@ Exposed-modules: Data.Manifold Data.Manifold.PseudoAffine Data.Manifold.TreeCover + Data.Manifold.Web Data.SimplicialComplex Data.LinearMap.HerMetric Data.Function.Differentiable @@ -72,6 +74,7 @@ Data.Manifold.Riemannian Other-modules: Data.List.FastNub Data.Manifold.Types.Primitive + Data.SetLike.Intersection Data.Manifold.Cone Data.CoNat Data.Embedding