diff --git a/Control/Monad/Trans/OuterMaybe.hs b/Control/Monad/Trans/OuterMaybe.hs new file mode 100644 --- /dev/null +++ b/Control/Monad/Trans/OuterMaybe.hs @@ -0,0 +1,22 @@ +-- | +-- Module : Control.Monad.Trans.OuterMaybe +-- Copyright : (c) Justus Sagemüller 2016 +-- License : GPL v3 +-- +-- Maintainer : (@) sagemueller $ geo.uni-koeln.de +-- Stability : experimental +-- Portability : portable +-- +{-# LANGUAGE DeriveFunctor #-} + +module Control.Monad.Trans.OuterMaybe where + +data OuterMaybeT f a = OuterNothing | OuterJust (f a) deriving (Functor) +instance (Applicative f) => Applicative (OuterMaybeT f) where + pure = OuterJust . pure + OuterJust fs <*> OuterJust xs = OuterJust $ fs <*> xs + _ <*> _ = OuterNothing + + + + diff --git a/Data/CoNat.hs b/Data/CoNat.hs --- a/Data/CoNat.hs +++ b/Data/CoNat.hs @@ -45,7 +45,6 @@ import Data.VectorSpace import Data.AffineSpace import Data.Basis -import Data.AdditiveGroup import qualified Data.List as List import qualified Prelude as Hask hiding(foldl) @@ -60,7 +59,6 @@ import qualified Data.Vector as Arr -import qualified Numeric.LinearAlgebra.HMatrix as HMat import Unsafe.Coerce diff --git a/Data/Embedding.hs b/Data/Embedding.hs --- a/Data/Embedding.hs +++ b/Data/Embedding.hs @@ -29,9 +29,6 @@ module Data.Embedding where -import Data.Tagged -import Data.Semigroup - import qualified Prelude as Hask hiding(foldl) import qualified Control.Applicative as Hask import qualified Control.Monad as Hask diff --git a/Data/Function/Affine.hs b/Data/Function/Affine.hs --- a/Data/Function/Affine.hs +++ b/Data/Function/Affine.hs @@ -37,28 +37,19 @@ -import Data.List -import Data.Maybe import Data.Semigroup import Data.VectorSpace import Data.LinearMap import Data.LinearMap.HerMetric -import Data.MemoTrie (HasTrie(..)) import Data.AffineSpace -import Data.Basis -import Data.Void import Data.Tagged import Data.Manifold.Types.Primitive import Data.Manifold.PseudoAffine -import Data.CoNat -import Data.VectorSpace.FiniteDimensional - import qualified Prelude import qualified Control.Applicative as Hask -import Data.Constraint.Trivial import Control.Category.Constrained.Prelude hiding ((^)) import Control.Category.Constrained.Reified import Control.Arrow.Constrained diff --git a/Data/Function/Differentiable.hs b/Data/Function/Differentiable.hs --- a/Data/Function/Differentiable.hs +++ b/Data/Function/Differentiable.hs @@ -50,33 +50,21 @@ import Data.List -import qualified Data.Vector.Generic as Arr -import qualified Data.Vector import Data.Maybe import Data.Semigroup -import Data.Function (on) import Data.Embedding -import Data.Fixed import Data.VectorSpace import Data.LinearMap import Data.LinearMap.Category import Data.LinearMap.HerMetric -import Data.MemoTrie (HasTrie(..)) import Data.AffineSpace import Data.Function.Differentiable.Data import Data.Function.Affine import Data.Basis -import Data.Complex hiding (magnitude) -import Data.Void import Data.Tagged import Data.Manifold.Types.Primitive import Data.Manifold.PseudoAffine - -import Data.CoNat -import Data.VectorSpace.FiniteDimensional - -import qualified Numeric.LinearAlgebra.HMatrix as HMat import qualified Prelude import qualified Control.Applicative as Hask diff --git a/Data/LinearMap/Category.hs b/Data/LinearMap/Category.hs --- a/Data/LinearMap/Category.hs +++ b/Data/LinearMap/Category.hs @@ -31,15 +31,12 @@ module Data.LinearMap.Category where import Data.Tagged -import Data.Semigroup -import Data.MemoTrie import Data.VectorSpace import Data.LinearMap import Data.VectorSpace.FiniteDimensional import Data.AffineSpace import Data.Basis -import Data.AdditiveGroup import qualified Prelude as Hask hiding(foldl) import qualified Control.Applicative as Hask @@ -139,6 +136,16 @@ 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) diff --git a/Data/LinearMap/HerMetric.hs b/Data/LinearMap/HerMetric.hs --- a/Data/LinearMap/HerMetric.hs +++ b/Data/LinearMap/HerMetric.hs @@ -30,17 +30,19 @@ -- * One-dimensional axes and product spaces , factoriseMetric, factoriseMetric' , productMetric, productMetric' - , metricAsLength, metricFromLength, metric'AsLength + , tryMetricAsLength, metricAsLength, metricFromLength, metric'AsLength -- * Utility for metrics , transformMetric, transformMetric', dualCoCoProduct , dualiseMetric, dualiseMetric' - , recipMetric, recipMetric' + , recipMetric, recipMetric', safeRecipMetric, safeRecipMetric' -- ** Eigenvectors , eigenSpan, eigenSpan' , eigenCoSpan, eigenCoSpan' , eigenSystem, HasEigenSystem, EigenVector + -- ** Scaling operations , metriNormalise, metriNormalise' , metriScale', metriScale + , volumeRatio, euclideanRelativeMetricVolume , adjoint , extendMetric , applyLinMapMetric, applyLinMapMetric' @@ -58,6 +60,7 @@ , linMapAsTensProd, linMapFromTensProd , covariance , outerProducts + , orthogonalComplementSpan ) where @@ -66,10 +69,8 @@ import Data.VectorSpace import Data.LinearMap import Data.Basis -import Data.MemoTrie import Data.Semigroup import Data.Tagged -import Data.Void import qualified Data.List as List import qualified Prelude as Hask @@ -297,20 +298,27 @@ -- 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' (HerMetric Nothing) = singularMetric' -recipMetric' (HerMetric (Just (DenseLinear m))) - | isInfinite' detm = singularMetric' - | otherwise = matrixMetric' minv - where (minv, (detm, _)) = HMat.invlndet m +recipMetric' m | Option (Just m') <- safeRecipMetric' m = m' +recipMetric' _ = singularMetric' -recipMetric :: HasMetric v => HerMetric' v -> HerMetric v -recipMetric (HerMetric' Nothing) = singularMetric -recipMetric (HerMetric' (Just (DenseLinear m))) - | isInfinite' detm = singularMetric - | otherwise = matrixMetric minv +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 @@ -449,7 +457,7 @@ 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, []) + eigenSystem (_, _) = (fmap Stiefel1 completeBasisValues, []) -- | Constraint that a space's scalars need to fulfill so it can be used for 'HerMetric'. @@ -728,11 +736,29 @@ = 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 | 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 ℝ @@ -765,13 +791,27 @@ -- | 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 :: forall s v w +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 @@ -781,8 +821,17 @@ -- 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 diff --git a/Data/List/FastNub.hs b/Data/List/FastNub.hs --- a/Data/List/FastNub.hs +++ b/Data/List/FastNub.hs @@ -6,6 +6,8 @@ import Data.List import Data.Function +import Data.Ord +import Control.Arrow ((&&&)) type FastNub a = (Eq a, Ord a) -- S̶h̶o̶u̶l̶d̶ ̶r̶e̶a̶l̶l̶y̶ ̶b̶e̶ ̶(̶E̶q̶ ̶a̶,̶̶ ̶H̶a̶s̶h̶a̶b̶l̶e̶ ̶a̶)̶ @@ -63,3 +65,10 @@ fis (x:xs) (y:ys) | xy = fis (x:xs) ys | otherwise = x : fis xs ys + + +-- | This function is also defined in "GHC.Exts", but only in a version that requires +-- 𝓞(𝑛⋅log 𝑛) function applications, as opposed to 𝑛 here. +sortWith :: Ord b => (a -> b) -> [a] -> [a] +sortWith f = map snd . sortBy (comparing fst) . map (f &&& id) + diff --git a/Data/Manifold/Cone.hs b/Data/Manifold/Cone.hs --- a/Data/Manifold/Cone.hs +++ b/Data/Manifold/Cone.hs @@ -31,22 +31,12 @@ -import Data.List import qualified Data.Vector.Generic as Arr -import qualified Data.Vector import Data.Maybe import Data.Semigroup -import Data.Function (on) -import Data.Fixed import Data.VectorSpace -import Data.LinearMap import Data.LinearMap.HerMetric -import Data.MemoTrie (HasTrie(..)) -import Data.AffineSpace -import Data.Basis -import Data.Complex hiding (magnitude) -import Data.Void import Data.Tagged import Data.Manifold.Types.Primitive @@ -64,7 +54,6 @@ import Data.Foldable.Constrained import Data.Manifold.PseudoAffine -import Data.Embedding diff --git a/Data/Manifold/DifferentialEquation.hs b/Data/Manifold/DifferentialEquation.hs new file mode 100644 --- /dev/null +++ b/Data/Manifold/DifferentialEquation.hs @@ -0,0 +1,117 @@ +-- | +-- Module : Data.Manifold.DifferentialEquation +-- 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 MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE TupleSections #-} +{-# LANGUAGE ParallelListComp #-} +{-# LANGUAGE UnicodeSyntax #-} +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE PatternGuards #-} +{-# LANGUAGE LambdaCase #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE LiberalTypeSynonyms #-} + + +module Data.Manifold.DifferentialEquation ( + -- * Formulating simple differential eqns. + DifferentialEqn + , constLinearDEqn + , filterDEqnSolution_static, iterateFilterDEqn_static + -- * Cost functions for error bounds + , maxDeviationsGoal + , uncertaintyGoal + , uncertaintyGoal' + , euclideanVolGoal + ) where + + +import Data.List.NonEmpty (NonEmpty(..)) +import qualified Data.List.NonEmpty as NE +import Data.Semigroup + +import Data.VectorSpace +import Data.LinearMap.HerMetric +import Data.LinearMap.Category +import Data.AffineSpace +import Data.Basis + +import Data.Manifold.Types +import Data.Manifold.PseudoAffine +import Data.Function.Differentiable +import Data.Function.Differentiable.Data +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) +import qualified Control.Applicative as Hask +import qualified Control.Monad as Hask hiding(forM_, sequence) +import qualified Data.Foldable as Hask +import qualified Data.Traversable as Hask + +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) + + +constLinearDEqn :: (WithField ℝ LinearManifold x, WithField ℝ LinearManifold y) + => Linear ℝ (DualSpace y) (Linear ℝ 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 + + +-- | A function that variates, relatively speaking, most strongly +-- for arguments around 1. In the zero-limit it approaches a constant +-- (but with arbitrarily large derivative); for η → ∞ the derivative +-- 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. +-- 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 vTgt _ (Shade' _ shy) = goalSensitive η + where η = euclideanRelativeMetricVolume shy / vTgt + +maxDeviationsGoal :: WithField ℝ EuclidSpace y => [Needle y] -> x -> Shade' y -> ℝ +maxDeviationsGoal = uncertaintyGoal . projector's + +uncertaintyGoal :: WithField ℝ EuclidSpace y => Metric' y -> x -> Shade' y -> ℝ +uncertaintyGoal = uncertaintyGoal' . const + +uncertaintyGoal' :: WithField ℝ EuclidSpace 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 + 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 @@ -37,47 +37,24 @@ import Data.List hiding (filter, all, elem, sum) import Data.Maybe -import qualified Data.Map as Map -import qualified Data.Vector as Arr -import Data.List.NonEmpty (NonEmpty(..)) -import Data.List.FastNub -import qualified Data.List.NonEmpty as NE -import Data.Semigroup -import Data.VectorSpace -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 ((^), (^.)) import Data.Manifold.PseudoAffine import Data.Manifold.TreeCover (Shade(..), fullShade, shadeCtr, shadeExpanse) 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.Class import qualified Data.Foldable as Hask import Data.Foldable (all, elem, toList, sum) 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(..), Traversable) import Control.Arrow.Constrained @@ -85,7 +62,6 @@ import Data.Foldable.Constrained import Text.Printf -import GHC.Generics (Generic) data GridAxis m g = GridAxInterval (Shade m) diff --git a/Data/Manifold/PseudoAffine.hs b/Data/Manifold/PseudoAffine.hs --- a/Data/Manifold/PseudoAffine.hs +++ b/Data/Manifold/PseudoAffine.hs @@ -65,17 +65,14 @@ -- ** Local functions , LocalLinear, LocalAffine -- * Misc - , palerp, LocallyCoercible(..) + , alerpB, palerp, palerpB, LocallyCoercible(..) + , ImpliesMetric(..) ) where -import Data.List -import qualified Data.Vector.Generic as Arr -import qualified Data.Vector import Data.Maybe import Data.Semigroup -import Data.Function (on) import Data.Fixed import Data.VectorSpace @@ -83,19 +80,13 @@ import Data.LinearMap import Data.LinearMap.HerMetric import Data.LinearMap.Category -import Data.MemoTrie (HasTrie(..)) import Data.AffineSpace -import Data.Basis -import Data.Complex hiding (magnitude) -import Data.Void import Data.Tagged import Data.Manifold.Types.Primitive import Data.CoNat import Data.VectorSpace.FiniteDimensional -import qualified Numeric.LinearAlgebra.HMatrix as HMat - import qualified Prelude import qualified Control.Applicative as Hask @@ -104,6 +95,7 @@ import Control.Monad.Constrained import Data.Foldable.Constrained +import GHC.Exts (Constraint) @@ -191,6 +183,7 @@ class ( Semimanifold x, Semimanifold (Interior x) , Needle (Interior x) ~ Needle x, Interior (Interior x) ~ Interior x) => PseudoAffine x where + {-# MINIMAL (.-~.) | (.-~!) #-} -- | The path reaching from one point to another. -- Should only yield 'Nothing' if -- @@ -212,6 +205,14 @@ -- manifold”. To adress this problem, these functions basically consider only the -- /interior/ of the space. (.-~.) :: x -> Interior x -> Option (Needle x) + p.-~.q = return $ p.-~!q + + -- | Unsafe version of '.-~.'. If the two points lie in disjoint regions, + -- the behaviour is undefined. + (.-~!) :: x -> Interior x -> Needle x + p.-~!q = case p.-~.q of + Option (Just v) -> v + @@ -320,17 +321,29 @@ -- its end points. -- -- A proper, really well-defined (on global scales) interpolation --- only makes sense on a Riemannian manifold, as geodesics. --- This is a task to be tackled in the future. +-- only makes sense on a Riemannian manifold, as 'Data.Manifold.Riemannian.Geodesic'. palerp :: ∀ x. Manifold x => Interior x -> Interior x -> Option (Scalar (Needle x) -> x) palerp p1 p2 = case (fromInterior p2 :: x) .-~. p1 of Option (Just v) -> return $ \t -> p1 .+~^ t *^ v _ -> empty +-- | Like 'palerp', but actually restricted to the interval between the points, +-- with a signature like 'Data.Manifold.Riemannian.geodesicBetween' +-- rather than 'Data.AffineSpace.alerp'. +palerpB :: ∀ x. WithField ℝ Manifold x => Interior x -> Interior x -> Option (D¹ -> x) +palerpB p1 p2 = case (fromInterior p2 :: x) .-~. p1 of + Option (Just v) -> return $ \(D¹ t) -> p1 .+~^ ((t+1)/2) *^ v + _ -> empty +-- | Like 'alerp', but actually restricted to the interval between the points. +alerpB :: ∀ x. (AffineSpace x, VectorSpace (Diff x), Scalar (Diff x) ~ ℝ) + => x -> x -> D¹ -> x +alerpB p1 p2 = case p2 .-. p1 of + v -> \(D¹ t) -> p1 .+^ ((t+1)/2) *^ v + hugeℝVal :: ℝ hugeℝVal = 1e+100 @@ -582,5 +595,24 @@ + +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 + +instance ImpliesMetric HerMetric' where + type MetricRequirement HerMetric' x = x ~ Needle x + inferMetric' = pure diff --git a/Data/Manifold/Riemannian.hs b/Data/Manifold/Riemannian.hs --- a/Data/Manifold/Riemannian.hs +++ b/Data/Manifold/Riemannian.hs @@ -46,47 +46,27 @@ module Data.Manifold.Riemannian where -import Data.List hiding (filter, all, elem, sum) import Data.Maybe -import qualified Data.Map as Map import qualified Data.Vector as Arr -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.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.Manifold.Types import Data.Manifold.Types.Primitive ((^), empty, embed, coEmbed) import Data.Manifold.PseudoAffine import Data.VectorSpace.FiniteDimensional -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.Class import qualified Data.Foldable as Hask -import Data.Foldable (all, elem, toList, sum) import qualified Data.Traversable as Hask -import Data.Traversable (forM) import qualified Numeric.LinearAlgebra.HMatrix as HMat @@ -96,10 +76,9 @@ import Control.Monad.Constrained hiding (forM) import Data.Foldable.Constrained -import GHC.Generics (Generic) -class PseudoAffine x => Geodesic x where +class Semimanifold x => Geodesic x where geodesicBetween :: x -- ^ Starting point; the interpolation will yield this at -1. -> x -- ^ End point, for +1. diff --git a/Data/Manifold/TreeCover.hs b/Data/Manifold/TreeCover.hs --- a/Data/Manifold/TreeCover.hs +++ b/Data/Manifold/TreeCover.hs @@ -46,9 +46,10 @@ -- ** Evaluation , occlusion -- ** Misc - , factoriseShade, intersectShade's, Refinable, refineShade', convolveShade', coerceShade + , factoriseShade, intersectShade's + , Refinable, subShade', refineShade', convolveShade', coerceShade -- * Shade trees - , ShadeTree(..), fromLeafPoints, onlyLeaves, indexShadeTree + , ShadeTree(..), fromLeafPoints, onlyLeaves, indexShadeTree, positionIndex -- * View helpers , onlyNodes -- ** Auxiliary types @@ -56,9 +57,9 @@ -- * Misc , sShSaw, chainsaw, HasFlatView(..), shadesMerge, smoothInterpolate , twigsWithEnvirons, completeTopShading, flexTwigsShading - , WithAny(..), Shaded, stiAsIntervalMapping, spanShading + , WithAny(..), Shaded, fmapShaded, stiAsIntervalMapping, spanShading , constShaded, stripShadedUntopological - , DifferentialEqn, filterDEqnSolution_loc + , DifferentialEqn, propagateDEqnSolution_loc -- ** Triangulation-builders , TriangBuild, doTriangBuild, singleFullSimplex, autoglueTriangulation , AutoTriang, elementaryTriang, breakdownAutoTriang @@ -78,34 +79,28 @@ 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.SetLike.Intersection +import Data.Manifold.Riemannian import Data.Embedding import Data.CoNat +import Lens.Micro (Lens') + 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.OuterMaybe import Control.Monad.Trans.Class import qualified Data.Foldable as Hask import Data.Foldable (all, elem, toList, sum, foldr1) @@ -153,7 +148,7 @@ class IsShade shade where -- type (*) shade :: *->* -- | Access the center of a 'Shade' or a 'Shade''. - shadeCtr :: Functor f (->) (->) => (Interior x->f (Interior x)) -> shade x -> f (shade x) + shadeCtr :: Lens' (shade x) (Interior x) -- -- | Convert between 'Shade' and 'Shade' (which must be neither singular nor infinite). -- unsafeDualShade :: WithField ℝ Manifold x => shade x -> shade* x -- | Check the statistical likelihood-density of a point being within a shade. @@ -179,7 +174,15 @@ coerceShade (Shade x (HerMetric' δxym)) = Shade (locallyTrivialDiffeomorphism x) (HerMetric' $ unsafeCoerceLinear<$>δxym) -shadeExpanse :: Functor f (->) (->) => (Metric' x -> f (Metric' x)) -> Shade x -> f (Shade x) +instance ImpliesMetric Shade where + type MetricRequirement Shade x = Manifold x + inferMetric' (Shade _ e) = pure e + +instance ImpliesMetric Shade' where + type MetricRequirement Shade' x = Manifold x + inferMetric (Shade' _ e) = pure e + +shadeExpanse :: Lens' (Shade x) (Metric' x) shadeExpanse f (Shade c e) = fmap (Shade c) $ f e instance IsShade Shade' where @@ -195,7 +198,7 @@ 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 :: Lens' (Shade' x) (Metric x) shadeNarrowness f (Shade' c e) = fmap (Shade' c) $ f e instance (AffineManifold x) => Semimanifold (Shade x) where @@ -206,6 +209,32 @@ 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 + geodesicBetween (Shade c e) (Shade ζ η) = pure interp + where ([], sharedSpan) = eigenSystem (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 ] + Option (Just pinterp) = geodesicBetween c ζ + +instance (AffineManifold x) => Semimanifold (Shade' x) where + type Needle (Shade' x) = Diff x + fromInterior = id + toInterior = pure + translateP = Tagged (.+~^) + 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 + geodesicBetween (Shade' c e) (Shade' ζ η) = pure interp + where ([], sharedSpan) = eigenSystem (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 ] + Option (Just pinterp) = geodesicBetween c ζ + fullShade :: WithField ℝ Manifold x => x -> Metric' x -> Shade x fullShade ctr expa = Shade ctr expa @@ -433,6 +462,15 @@ scanLeafNums i₀ ((v,t):vts) = (i₀, (v,t)) : scanLeafNums (i₀ + nLeaves t) vts +indexDBranches :: NonEmpty (DBranch x) -> NonEmpty (DBranch' x (Int, ShadeTree x)) +indexDBranches (DBranch d (Hourglass t b) :| l) -- this could more concisely be written as a traversal + = DBranch d (Hourglass (0,t) (nt,b)) :| ixDBs (nt + nb) l + where nt = nLeaves t; nb = nLeaves b + ixDBs _ [] = [] + ixDBs i₀ (DBranch δ (Hourglass τ β) : l) + = DBranch δ (Hourglass (i₀,τ) (i₀+nτ,β)) : ixDBs (i₀ + nτ + nβ) l + where nτ = nLeaves τ; nβ = nLeaves β + instance (NFData x, NFData (Needle' x)) => NFData (ShadeTree x) where rnf (PlainLeaves xs) = rnf xs rnf (DisjointBranches n bs) = n `seq` rnf (NE.toList bs) @@ -477,6 +515,8 @@ fromLeafPoints = fromLeafPoints' sShIdPartition +-- | The leaves of a shade tree are numbered. For a given index, this function +-- attempts to find the leaf with that ID, within its immediate environment. indexShadeTree :: ∀ x . WithField ℝ Manifold x => ShadeTree x -> Int -> Either Int ([ShadeTree x], x) indexShadeTree _ i @@ -498,6 +538,44 @@ | otherwise = Left $ i-n +-- | “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 + => 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 + -> x -- ^ Position to look up + -> Option (Int, ([ShadeTree x], x)) + -- ^ Index of the leaf near to the query point, the “path” of + -- 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) + | (i,p) <- zip [0..] lvs] of + [] -> empty + l | ((i,p),_) <- minimumBy (comparing snd) l + -> pure (i, ([sh], p)) +positionIndex m (DisjointBranches _ brs) x + = fst . foldl' (\case + (q@(Option (Just _)), i₀) -> const (q, i₀) + (_, i₀) -> \t' -> ( first (+i₀) <$> positionIndex m t' x + , i₀+nLeaves t' ) ) + (empty, 0) + $ brs +positionIndex _ sh@(OverlappingBranches n (Shade c ce) brs) x + | Option (Just vx) <- x.-~.c + = let (_,(i₀,t')) = maximumBy (comparing fst) + [ (σ*ω, t') + | DBranch d (Hourglass t'u t'd) <- NE.toList $ indexDBranches brs + , let ω = d<.>^vx + , (t',σ) <- [(t'u, 1), (t'd, -1)] ] + in ((+i₀) *** first (sh:)) + <$> positionIndex (return $ recipMetric ce) t' x +positionIndex _ _ _ = empty + + + fromFnGraphPoints :: ∀ x y . (WithField ℝ Manifold x, WithField ℝ Manifold y) => [(x,y)] -> ShadeTree (x,y) fromFnGraphPoints = fromLeafPoints' fg_sShIdPart @@ -602,6 +680,17 @@ nLeaves (DisjointBranches n _) = n nLeaves (OverlappingBranches n _ _) = n + +instance ImpliesMetric ShadeTree where + type MetricRequirement ShadeTree x = WithField ℝ Manifold x + inferMetric' (OverlappingBranches _ (Shade _ e) _) = pure e + inferMetric' (PlainLeaves lvs) = case pointsShades lvs of + (Shade _ sh:_) -> pure sh + _ -> empty + inferMetric' (DisjointBranches _ (br:|_)) = inferMetric' br + + + overlappingBranches :: Shade x -> NonEmpty (DBranch x) -> ShadeTree x overlappingBranches shx brs = OverlappingBranches n shx brs where n = sum $ fmap (sum . fmap nLeaves) brs @@ -624,7 +713,7 @@ in DisjointBranches (sum $ nLeaves<$>brs') brs' unsafeFmapTree f fn fs (OverlappingBranches n sh brs) = let brs' = fmap (\(DBranch dir br) - -> DBranch (fn dir) (unsafeFmapTree f fn fs<$>br) + -> DBranch (fn dir) (unsafeFmapTree f fn fs<$>br) ) brs in overlappingBranches (fs sh) brs' @@ -636,37 +725,125 @@ -- @'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 ] + [ ac.+~^σ*^v | σ<-[-1,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 - [ (v, 1 / (1 + metricSq oExpa w)) - | v <- (*^) <$> [-1,1] <*> span - , let p = ctr .+~^ v :: y - Option (Just w) = p.-~.oCtr - ] - | ctr <- [c, ζ ] - | span <- eigenCoSpan'<$> [e, η ] - | (oCtr,oExpa) <- [(ζ,η), (c,e)] - ] - Option (Just c'2ζ') = ζ'.-~.c' - Option (Just c2ζ') = ζ'.-~.c - Option (Just ζ2c') = c'.-~.ζ - μc = metricSq e c2ζ' - μζ = metricSq η ζ2c' - iCtr = c' .+~^ c'2ζ' ^* (μζ/(μc + μζ)) -- weighted mean between c' and ζ'. - Option (Just rc) = c.-~.iCtr - Option (Just rζ) = ζ.-~.iCtr - rcⰰ = toDualWith e rc - rζⰰ = toDualWith η rζ - μe = rcⰰ<.>^rc - μη = rζⰰ<.>^rζ - iExpa = (e^+^η)^/2 ^+^ projector rcⰰ^/(1-μe) ^+^ projector rζⰰ^/(1-μη) + refineShade' (Shade' c₀ (HerMetric (Just e₁))) + (Shade' c₀₂ (HerMetric (Just e₂))) + | Option (Just c₂) <- c₀₂.-~.c₀ + , e₁c₂ <- e₁ $ c₂ + , e₂c₂ <- 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₂eec₂ <- (c₂e₁c₂ + c₂e₂c₂) / α + , [γ₁,γ₂] <- middle . sort + $ quadraticEqnSol c₂e₁c₂ + (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) + , cc' <- cc ^+^ ((γ₁+γ₂)/2)*^c₂ + , rγ <- abs (γ₁ - γ₂) / 2 + , η <- if rγ * c₂eec₂ /= 0 && 1 - rγ^2 * c₂eec₂ > 0 + then sqrt (1 - rγ^2 * c₂eec₂) / (rγ * c₂eec₂) + else 0 + = return $ + Shade' (c₀.+~^cc') + (HerMetric (Just ee) ^+^ projector (ee $ c₂^*η)) + + | otherwise = empty + where σe = e₁^+^e₂ + quadraticEqnSol a b c + | a /= 0 && disc > 0 = [ (σ * sqrt disc - b) / (2*a) + | σ <- [-1, 1] ] + | otherwise = [0] + 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. + -- Let WLOG c₁ = 0, so + -- ⟨x|e₁|x⟩ < 1. + -- cc should minimise the quadratic form + -- β(cc) = ⟨cc−c₁|e₁|cc−c₁⟩ + ⟨cc−c₂|e₂|cc−c₂⟩ + -- = ⟨cc|e₁|cc⟩ + ⟨cc−c₂|e₂|cc−c₂⟩ + -- = ⟨cc|e₁|cc⟩ + ⟨cc|e₂|cc⟩ − 2⋅⟨c₂|e₂|cc⟩ + ⟨c₂|e₂|c₂⟩ + -- It is thus + -- β(cc + δ⋅v) − β cc + -- = ⟨cc + δ⋅v|e₁|cc + δ⋅v⟩ + ⟨cc + δ⋅v|e₂|cc + δ⋅v⟩ − 2⋅⟨c₂|e₂|cc + δ⋅v⟩ + ⟨c₂|e₂|c₂⟩ + -- − ⟨cc|e₁|cc⟩ − ⟨cc|e₂|cc⟩ + 2⋅⟨c₂|e₂|cc⟩ − ⟨c₂|e₂|c₂⟩ + -- = ⟨cc + δ⋅v|e₁|cc + δ⋅v⟩ + ⟨cc + δ⋅v|e₂|cc + δ⋅v⟩ − 2⋅⟨c₂|e₂|δ⋅v⟩ + -- − ⟨cc|e₁|cc⟩ − ⟨cc|e₂|cc⟩ + -- = 2⋅⟨δ⋅v|e₁|cc⟩ + ⟨δ⋅v|e₁|δ⋅v⟩ + 2⋅⟨δ⋅v|e₂|cc⟩ + ⟨δ⋅v|e₂|δ⋅v⟩ − 2⋅⟨c₂|e₂|δ⋅v⟩ + -- = 2⋅δ⋅⟨v|e₁+e₂|cc⟩ − 2⋅δ⋅⟨v|e₂|c₂⟩ + 𝓞(δ²) + -- This should vanish for all v, which is fulfilled by + -- (e₁+e₂)|cc⟩ = e₂|c₂⟩. + -- + -- If we now choose + -- ee = (e₁+e₂) / α + -- then + -- ⟨x−cc|ee|x−cc⟩ ⋅ α + -- = ⟨x−cc|ee|x⟩ ⋅ α − ⟨x−cc|ee|cc⟩ ⋅ α + -- = ⟨x|ee|x−cc⟩ ⋅ α − ⟨x−cc|e₂|c₂⟩ + -- = ⟨x|ee|x⟩ ⋅ α − ⟨x|ee|cc⟩ ⋅ α − ⟨x−cc|e₂|c₂⟩ + -- = ⟨x|e₁+e₂|x⟩ − ⟨x|e₂|c₂⟩ − ⟨x−cc|e₂|c₂⟩ + -- = ⟨x|e₁|x⟩ + ⟨x|e₂|x⟩ − ⟨x|e₂|c₂⟩ − ⟨x−cc|e₂|c₂⟩ + -- < 1 + ⟨x|e₂|x−c₂⟩ − ⟨x−cc|e₂|c₂⟩ + -- = 1 + ⟨x−c₂|e₂|x−c₂⟩ + ⟨c₂|e₂|x−c₂⟩ − ⟨x−cc|e₂|c₂⟩ + -- < 2 + ⟨x−c₂−x+cc|e₂|c₂⟩ + -- = 2 + ⟨cc−c₂|e₂|c₂⟩ + -- Really we want + -- ⟨x−cc|ee|x−cc⟩ ⋅ α < α + -- So choose α = 2 + ⟨cc−c₂|e₂|c₂⟩. + -- + -- The ellipsoid "cc±√ee" captures perfectly the intersection + -- of the boundary of the shades, but it tends to significantly + -- overshoot the interior intersection in perpendicular direction, + -- i.e. in direction of c₂−c₁. E.g. + -- https://github.com/leftaroundabout/manifolds/blob/bc0460b9/manifolds/images/examples/ShadeCombinations/EllipseIntersections.png + -- 1. Really, the relevant points are those where either of the + -- intersector badnesses becomes 1. The intersection shade should + -- be centered between those points. We perform according corrections, + -- but only in c₂ direction, so this can be handled efficiently + -- as a 1D quadratic equation. + -- Consider + -- dⱼ c := ⟨c−cⱼ|eⱼ|c−cⱼ⟩ =! 1 + -- dⱼ (cc + γ⋅c₂) + -- = ⟨cc+γ⋅c₂−cⱼ|eⱼ|cc+γ⋅c₂−cⱼ⟩ + -- = ⟨cc−cⱼ|eⱼ|cc−cⱼ⟩ + 2⋅γ⋅⟨c₂|eⱼ|cc−cⱼ⟩ + γ²⋅⟨c₂|eⱼ|c₂⟩ + -- =! 1 + -- So + -- γⱼ = (- b ± √(b²−4⋅a⋅c)) / 2⋅a + -- where a = ⟨c₂|eⱼ|c₂⟩ + -- b = 2 ⋅ (⟨c₂|eⱼ|cc⟩ − ⟨c₂|eⱼ|cⱼ⟩) + -- c = ⟨cc|eⱼ|cc⟩ − 2⋅⟨cc|eⱼ|cⱼ⟩ + ⟨cⱼ|eⱼ|cⱼ⟩ − 1 + -- The ± sign should be chosen to get the smaller |γ| (otherwise + -- we end up on the wrong side of the shade), i.e. + -- γⱼ = (sgn bⱼ ⋅ √(bⱼ²−4⋅aⱼ⋅cⱼ) − bⱼ) / 2⋅aⱼ + -- 2. Trim the result in that direction to the actual + -- thickness of the lens-shaped intersection: we want + -- ⟨rγ⋅c₂|ee'|rγ⋅c₂⟩ = 1 + -- for a squeezed version of ee, + -- ee' = ee + ee|η⋅c₂⟩⟨η⋅c₂|ee + -- ee' = ee + η² ⋅ ee|c₂⟩⟨c₂|ee + -- ⟨rγ⋅c₂|ee'|rγ⋅c₂⟩ + -- = rγ² ⋅ (⟨c₂|ee|c₂⟩ + η² ⋅ ⟨c₂|ee|c₂⟩²) + -- = rγ² ⋅ ⟨c₂|ee|c₂⟩ + η² ⋅ rγ² ⋅ ⟨c₂|ee|c₂⟩² + -- η² = (1 − rγ²⋅⟨c₂|ee|c₂⟩) / (rγ² ⋅ ⟨c₂|ee|c₂⟩²) + -- η = √(1 − rγ²⋅⟨c₂|ee|c₂⟩) / (rγ ⋅ ⟨c₂|ee|c₂⟩) + -- With ⟨c₂|ee|c₂⟩ = (⟨c₂|e₁|c₂⟩ + ⟨c₂|e₂|c₂⟩)/α. + -- | If @p@ is in @a@ (red) and @δ@ is in @b@ (green), -- then @p.+~^δ@ is in @convolveShade' a b@ (blue). @@ -728,31 +905,31 @@ 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 +propagateDEqnSolution_loc :: ∀ x y . (WithField ℝ Manifold x, Refinable y) + => DifferentialEqn x y -> ((x, Shade' y), NonEmpty (Needle x, Shade' y)) + -> NonEmpty (Shade' y) +propagateDEqnSolution_loc f ((x, shy@(Shade' y _)), neighbours) = ycs where jShade@(Shade' j₀ jExpa) = f shxy [shxy] = pointsCovers [ (xs, ys') - | (xs, Shade' ys yse) <- (x,shy):NE.toList neighbours + | (xs, Shade' ys yse) + <- (x,shy):(first (x.+~^)<$>NE.toList neighbours) , δy <- eigenCoSpan' yse , ys' <- [ys.+~^δy, ys.-~^δy] ] - [Shade' _ expax] = pointsCover's $ x : (fst<$>NE.toList neighbours) + [Shade' _ expax] = pointsCover's $ x : ((x.+~^).fst<$>NE.toList neighbours) marginδs :: NonEmpty (Needle x, (Needle y, Metric y)) marginδs = [ (δxm, (δym, expany)) - | (xn, Shade' yn expany) <- neighbours - , let (Option (Just δxm)) = xn.-~.x - (Option (Just δym)) = yn.-~.y + | (δxm, Shade' yn expany) <- neighbours + , let (Option (Just δym)) = yn.-~.y ] back2Centre :: (Needle x, (Needle y, Metric y)) -> Shade' y back2Centre (δx, (δym, expany)) = convolveShade' (Shade' y expany) - (Shade' δyb $ applyLinMapMetric jExpa δx') + (Shade' δyb $ applyLinMapMetric jExpa (δx'^/(δx'<.>^δx))) where δyb = δym ^-^ (j₀ $ δx) δx' = toDualWith expax δx - yc :: Option (Shade' y) - yc = intersectShade's $ back2Centre <$> marginδs + ycs :: NonEmpty (Shade' y) + ycs = back2Centre <$> marginδs xSpan = eigenCoSpan' expax @@ -763,12 +940,6 @@ => 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 - pure = OuterJust . pure - OuterJust fs <*> OuterJust xs = OuterJust $ fs <*> xs - _ <*> _ = OuterNothing - traverseTwigsWithEnvirons :: ∀ x f . (WithField ℝ Manifold x, Hask.Applicative f) => ( ((Int, ShadeTree x), [(Int, ShadeTree x)]) -> f (ShadeTree x)) @@ -1444,6 +1615,11 @@ stripShadedUntopological :: x`Shaded`y -> ShadeTree x stripShadedUntopological = unsafeFmapTree (fmap _topological) id shadeWithoutAnything + +fmapShaded :: (y -> υ) -> (x`Shaded`y) -> (x`Shaded`υ) +fmapShaded f = unsafeFmapTree (fmap $ \(WithAny y x) -> WithAny (f y) x) + id + (\(Shade yx shx) -> Shade (fmap f yx) shx) -- | This is to 'ShadeTree' as 'Data.Map.Map' is to 'Data.Set.Set'. type x`Shaded`y = ShadeTree (x`WithAny`y) diff --git a/Data/Manifold/Types.hs b/Data/Manifold/Types.hs --- a/Data/Manifold/Types.hs +++ b/Data/Manifold/Types.hs @@ -50,9 +50,12 @@ , D¹(..), D²(..) , ℝay , CD¹(..), Cℝay(..) - -- * Cut-planes + -- * Affine subspaces + -- ** Lines + , Line(..), lineAsPlaneIntersection + -- ** Hyperplanes , Cutplane(..) - , fathomCutDistance, sideOfCut + , fathomCutDistance, sideOfCut, cutPosBetween -- * Linear mappings , Linear, LocalLinear, denseLinear ) where @@ -63,9 +66,7 @@ import Data.MemoTrie (HasTrie(..)) import Data.Basis import Data.Fixed -import Data.Void import Data.Tagged -import Data.Monoid import Data.Semigroup import qualified Numeric.LinearAlgebra.HMatrix as HMat import qualified Data.Vector.Generic as Arr @@ -221,6 +222,8 @@ +data Line x = Line { lineHandle :: x + , lineDirection :: Stiefel1 (Needle' x) } @@ -258,4 +261,16 @@ where fathom v = (cn <.>^ v) / scaleDist scaleDist = metric' met cn + +cutPosBetween :: WithField ℝ Manifold x => Cutplane x -> (x,x) -> Option D¹ +cutPosBetween (Cutplane h (Stiefel1 cn)) (x₀,x₁) + | Option (Just [d₀,d₁]) <- map (cn<.>^) <$> sequenceA [x₀.-~.h, x₁.-~.h] + , d₀*d₁ < 0 + = pure . D¹ $ d₁ / (d₁ - d₀) + | otherwise = empty + + +lineAsPlaneIntersection :: WithField ℝ Manifold x => Line x -> [Cutplane x] +lineAsPlaneIntersection (Line h dir) + = [Cutplane h nrml | nrml <- orthogonalComplementSpan [dir]] 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 @@ -44,7 +44,7 @@ -- * Projective spaces , ℝP¹, ℝP²(..) -- * Intervals\/disks\/cones - , D¹(..), D²(..) + , D¹(..), fromIntv0to1, D²(..) , ℝay , CD¹(..), Cℝay(..) -- * Tensor products @@ -59,7 +59,6 @@ import Data.VectorSpace import Data.AffineSpace import Data.Basis -import Data.Complex hiding (magnitude) import Data.Void import Data.Monoid @@ -67,12 +66,12 @@ import Control.Applicative (Const(..), Alternative(..)) +import Lens.Micro ((^.)) + import qualified Prelude import Control.Category.Constrained.Prelude hiding ((^)) import Control.Arrow.Constrained -import Control.Monad.Constrained -import Data.Foldable.Constrained import Data.Embedding @@ -155,6 +154,10 @@ -- the two points -1 and 1 of 'S⁰', i.e. this is simply a closed interval. newtype D¹ = D¹ { xParamD¹ :: Double -- ^ Range @[-1, 1]@. } +fromIntv0to1 :: ℝ -> D¹ +fromIntv0to1 x | x<0 = D¹ (-1) + | x>1 = D¹ 1 + | otherwise = D¹ $ (x+1)/2 -- | The standard, closed unit disk. Homeomorphic to the cone over 'S¹', but not in the -- the obvious, “flat” way. (And not at all, despite @@ -286,8 +289,4 @@ (^) = (Prelude.^) -infixl 8 ^. -{-# INLINE (^.) #-} -(^.) :: s -> (forall f . Prelude.Functor f => (a->f a) -> s->f s) -> a -o ^. g = getConst (g Const o) diff --git a/Data/Manifold/Web.hs b/Data/Manifold/Web.hs --- a/Data/Manifold/Web.hs +++ b/Data/Manifold/Web.hs @@ -16,7 +16,7 @@ {-# LANGUAGE DeriveTraversable #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE GADTs #-} {-# LANGUAGE RankNTypes #-} @@ -25,17 +25,32 @@ {-# LANGUAGE UnicodeSyntax #-} {-# LANGUAGE ConstraintKinds #-} {-# LANGUAGE PatternGuards #-} -{-# LANGUAGE PatternSynonyms #-} -{-# LANGUAGE ViewPatterns #-} {-# LANGUAGE LambdaCase #-} {-# LANGUAGE TypeOperators #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE LiberalTypeSynonyms #-} -{-# LANGUAGE RecordWildCards #-} -{-# LANGUAGE DataKinds #-} +{-# LANGUAGE TemplateHaskell #-} -module Data.Manifold.Web where +module Data.Manifold.Web ( + -- * The web data type + PointsWeb + -- ** Construction + , fromWebNodes, fromShadeTree_auto, fromShadeTree, fromShaded + -- ** Lookup + , nearestNeighbour, indexWeb, webEdges, toGraph + -- ** Decomposition + , sliceWeb_lin -- , sampleWebAlongLine_lin + -- ** Local environments + , localFocusWeb + -- * Differential equations + -- ** Fixed resolution + , filterDEqnSolution_static, iterateFilterDEqn_static + -- ** Automatic resolution + , filterDEqnSolutions_adaptive, iterateFilterDEqn_adaptive + -- * Misc + , ConvexSet(..), ellipsoid + ) where import Data.List hiding (filter, all, elem, sum, foldr1) @@ -44,49 +59,35 @@ 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.List.FastNub (fastNubBy) import Data.Ord (comparing) +import Data.Semigroup 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.Function (on) +import Data.Fixed (mod') -import Data.SimplicialComplex import Data.Manifold.Types -import Data.Manifold.Types.Primitive ((^), empty) +import Data.Manifold.Types.Primitive import Data.Manifold.PseudoAffine -import Data.Function.Differentiable -import Data.Function.Differentiable.Data import Data.Manifold.TreeCover +import Data.SetLike.Intersection +import Data.Manifold.Riemannian -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 Control.Monad.Trans.List import qualified Data.Foldable as Hask -import Data.Foldable (all, elem, toList, sum, foldr1) +import Data.Foldable (all, toList) import qualified Data.Traversable as Hask import Data.Traversable (forM) - -import qualified Numeric.LinearAlgebra.HMatrix as HMat +import Data.Graph import Control.Category.Constrained.Prelude hiding ((^), all, elem, sum, forM, Foldable(..), foldr1, Traversable, traverse) @@ -95,20 +96,36 @@ import Data.Foldable.Constrained import Data.Traversable.Constrained (Traversable, traverse) +import Control.Comonad (Comonad(..)) +import Lens.Micro ((&), (%~), (^.), (.~)) +import Lens.Micro.TH + import GHC.Generics (Generic) type WebNodeId = Int -type NeighbourRefs = UArr.Vector WebNodeId +data Neighbourhood x = Neighbourhood { + neighbours :: UArr.Vector WebNodeId + , localScalarProduct :: Metric x + } + deriving (Generic) + +instance (NFData x, NFData (HerMetric (Needle 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. +-- +-- †In general, there can be disconnected components, but every connected +-- component is strongly connected. data PointsWeb :: * -> * -> * where PointsWeb :: { webNodeRsc :: ShadeTree x - , webNodeAssocData :: Arr.Vector (y, NeighbourRefs) + , webNodeAssocData :: Arr.Vector (y, Neighbourhood x) } -> PointsWeb x y deriving (Generic, Hask.Functor, Hask.Foldable, Hask.Traversable) -instance (NFData x, NFData (Needle' x), NFData y) => NFData (PointsWeb x y) +instance (NFData x, NFData (HerMetric (Needle x)), NFData (Needle' x), NFData y) => NFData (PointsWeb x y) instance Foldable (PointsWeb x) (->) (->) where ffoldl = uncurry . Hask.foldl' . curry @@ -121,10 +138,18 @@ +type MetricChoice x = Shade x -> Metric x + + fromWebNodes :: ∀ x y . WithField ℝ Manifold x - => (Shade x->Metric x) -> [(x,y)] -> PointsWeb x y + => (MetricChoice x) -> [(x,y)] -> PointsWeb x y fromWebNodes mf = fromShaded mf . fromLeafPoints . map (uncurry WithAny . swap) +fromTopWebNodes :: ∀ x y . WithField ℝ Manifold 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 () @@ -133,33 +158,46 @@ fromShadeTree mf = fromShaded mf . constShaded () fromShaded :: ∀ x y . WithField ℝ Manifold x - => (Shade x -> Metric x) -- ^ Local scalar-product generator. You can always + => (MetricChoice 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 +fromShaded metricf = fromTopShaded metricf . fmapShaded ([],) + +fromTopShaded :: ∀ x y . WithField ℝ Manifold x + => (MetricChoice x) + -> (x`Shaded`([Needle x], y)) -- ^ Source tree, with a priori topology information + -- (needles pointing to already-known neighbour candidates) + -> PointsWeb x y +fromTopShaded 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 + locMesh :: ( (Int, ShadeTree (x`WithAny`([Needle x], y))) + , [(Int, ShadeTree (x`WithAny`([Needle x], y)))]) + -> Arr.Vector (y, Neighbourhood x) + locMesh ((i₀, locT), neighRegions) = Arr.map findNeighbours $ Arr.fromList locLeaves + where locLeaves :: [ (Int, x`WithAny`([Needle x], y)) ] + locLeaves = map (first (+i₀)) . zip [0..] $ onlyLeaves locT + vicinityLeaves :: [(Int, x)] vicinityLeaves = Hask.foldMap - (\(i₀n, ngbR) -> Arr.map (first (+i₀n)) - . Arr.indexed - . Arr.fromList + (\(i₀n, ngbR) -> map ((+i₀n) *** _topological) + . zip [0..] $ 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) -> + findNeighbours :: (Int, x`WithAny`([Needle x], y)) -> (y, Neighbourhood x) + findNeighbours (i, WithAny (vns,y) x) + = (y, Neighbourhood + (UArr.fromList $ fst<$>execState seek mempty) + locRieM ) + where seek :: State [(Int, (Needle x, Needle' x))] () + seek = do + Hask.forM_ ( fastNubBy (comparing fst) + $ map (second _topological) locLeaves + ++ vicinityLeaves ++ aprioriNgbs ) + $ \(iNgb, xNgb) -> when (iNgb/=i) `id`do let (Option (Just v)) = xNgb.-~.x oldNgbs <- get @@ -171,6 +209,12 @@ | neighbour@(_,(nv,_))<-oldNgbs , visibleOverlap w nv ] + aprioriNgbs :: [(Int, x)] + aprioriNgbs = catMaybes + [ getOption $ (second $ const xN) <$> + positionIndex (pure locRieM) shd' xN + | v <- vns + , let xN = x.+~^v :: x ] visibleOverlap :: Needle' x -> Needle x -> Bool visibleOverlap w v = o < 1 @@ -188,32 +232,429 @@ , Right (_,x) <- indexShadeTree rsc i = pure (x, fst (assocD Arr.! i)) | otherwise = empty +unsafeIndexWebData :: PointsWeb x y -> WebNodeId -> y +unsafeIndexWebData (PointsWeb _ asd) i = fst (asd Arr.! i) + webEdges :: ∀ x y . WithField ℝ Manifold x => 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)) + allEdges = Hask.foldMap (\(i,(_, Neighbourhood 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)]) +data InterpolationIv y = InterpolationIv { + _interpolationSegRange :: (ℝ,ℝ) + , _interpolationFunction :: ℝ -> y + } + +type InterpolationSeq y = [InterpolationIv y] + +mkInterpolationSeq_lin :: (x~ℝ, Geodesic y) + => [(x,y)] -> InterpolationSeq y +mkInterpolationSeq_lin [(xψ,yψ), (xω,yω)] + = return $ InterpolationIv + (xψ,xω) + (\x -> let drel = fromIntv0to1 $ (x-xψ)/(xω-xψ) + in yio drel ) + where Option (Just yio) = geodesicBetween yψ yω +mkInterpolationSeq_lin (p₀:p₁:ps) + = mkInterpolationSeq_lin [p₀,p₁] <> mkInterpolationSeq_lin (p₁:ps) +mkInterpolationSeq_lin _ = [] + + +-- | 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) + => PointsWeb x y -> Cutplane x -> [(x,y)] +sliceWeb_lin web = sliceEdgs + where edgs = webEdges web + sliceEdgs cp = [ (xi d, yi d) -- Brute-force search through all edges + | ((x₀,y₀), (x₁,y₁)) <- edgs + , Option (Just d) <- [cutPosBetween cp (x₀,x₁)] + , Option (Just xi) <- [geodesicBetween x₀ x₁] + , Option (Just yi) <- [geodesicBetween y₀ y₁] + ] + +-- sampleWebAlongLine_lin :: ∀ x y . (WithField ℝ Manifold x, Geodesic x, Geodesic y) +-- => PointsWeb x y -> x -> Needle x -> [(x,y)] +-- sampleWebAlongLine_lin web x₀ dir = sampleWebAlongLines_lin web x₀ [(dir, maxBound)] + + +data GridPlanes x = GridPlanes { + _gridPlaneNormal :: Needle' x + , _gridPlaneSpacing :: Needle x + , _gridPlanesCount :: Int + } +data GridSetup x = GridSetup { + _gridStartCorner :: x + , _gridSplitDirs :: [GridPlanes x] + } + +cartesianGrid2D :: (x~ℝ, y~ℝ) => ((x,x), Int) -> ((y,y), Int) -> GridSetup (x,y) +cartesianGrid2D ((x₀,x₁), nx) ((y₀,y₁), ny) + = 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) + => 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) + => 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)] + finalLine ((x₀, GridPlanes _ dir nSpl), verts) + | 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 ] + 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 + +sampleWeb_2Dcartesian_lin :: (x~ℝ, y~ℝ, Geodesic z) + => PointsWeb (x,y) z -> ((x,x),Int) -> ((y,y),Int) -> [(y,[(x,Option z)])] +sampleWeb_2Dcartesian_lin web (xspec@(_,nx)) yspec + = go . sampleWebAlongGrid_lin web $ cartesianGrid2D xspec yspec + where go [] = [] + go l@(((_,y),_):_) = let (ln,l') = splitAt nx l + in (y, map (\((x,_),z) -> (x,z)) ln) : go l' + +sampleEntireWeb_2Dcartesian_lin :: (x~ℝ, y~ℝ, Geodesic z) + => PointsWeb (x,y) z -> Int -> Int -> [(y,[(x,Option z)])] +sampleEntireWeb_2Dcartesian_lin web nx ny + = sampleWeb_2Dcartesian_lin web ((x₀,x₁),nx) ((y₀,y₁),ny) + where x₀ = minimum (fst<$>pts) + x₁ = maximum (fst<$>pts) + y₀ = minimum (snd<$>pts) + y₁ = maximum (snd<$>pts) + pts = fst . fst <$> toList (localFocusWeb web) + +webLocalInfo :: ∀ x y . WithField ℝ Manifold x + => PointsWeb x y -> PointsWeb x (WebLocally x y) +webLocalInfo origWeb = result + where result = wli $ localFocusWeb origWeb + wli (PointsWeb rsc asd) = PointsWeb rsc asd' + where asd' = Arr.imap localInfo asd + localInfo i (((x,y), ngbCo), ngbH) + = ( LocalWebInfo { + _thisNodeCoord = x + , _thisNodeData = y + , _containingWeb = result + , _thisNodeId = i + , _nodeNeighbours = zip (UArr.toList $ neighbours ngbH) ngbCo + , _nodeLocalScalarProduct = localScalarProduct ngbH + , _nodeIsOnBoundary = anyUnopposed (localScalarProduct ngbH) ngbCo + }, ngbH ) + anyUnopposed rieM ngbCo = (`any`ngbCo) $ \(v,_) + -> not $ (`any`ngbCo) $ \(v',_) + -> toDualWith rieM v <.>^ v' < 0 + +localFocusWeb :: WithField ℝ Manifold x + => PointsWeb x y -> PointsWeb x ((x,y), [(Needle 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''= Arr.map (\((x,y),n) -> + (((x,y), [ ( case x'.-~.x of + Option (Just v) -> v + , y') + | j<-UArr.toList (neighbours n) + , let ((x',y'),_) = asd' Arr.! j + ]), n) ) asd' +nearestNeighbour :: WithField ℝ Manifold 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)) + : neighbours + (yProx, _) = asd Arr.! iProx + (_, Neighbourhood neighbourIds locMetr) = asd Arr.! iEst + neighbours = [ (i, (xNgb, metricSq locMetr v)) + | i <- UArr.toList neighbourIds + , let Right (_, xNgb) = indexShadeTree rsc i + Option (Just v) = xNgb.-~.x + ] + Option (Just vEst) = xEst.-~.x + + + +data WebLocally x y = LocalWebInfo { + _thisNodeCoord :: x + , _thisNodeData :: y + , _containingWeb :: PointsWeb x (WebLocally x y) + , _thisNodeId :: WebNodeId + , _nodeNeighbours :: [(WebNodeId, (Needle x, y))] + , _nodeLocalScalarProduct :: Metric x + , _nodeIsOnBoundary :: Bool + } deriving (Generic) +makeLenses ''WebLocally + +instance Hask.Functor (WebLocally x) where + fmap f (LocalWebInfo co dt wb id ng sp bn) + = LocalWebInfo co (f dt) (fmap (fmap f) wb) id (map (second $ second f) ng) sp bn +instance WithField ℝ Manifold x => Comonad (WebLocally x) where + extract = _thisNodeData + duplicate lweb = unsafeIndexWebData deepened $ _thisNodeId lweb + where deepened = webLocalInfo $ _containingWeb lweb + + + + + +toGraph :: WithField ℝ Manifold 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])] + edgs = Arr.toList + . Arr.imap (\i (_, Neighbourhood ngbs _) -> (i, i, UArr.toList ngbs)) + $ webNodeAssocData wb + + + + +data ConvexSet x + = EmptyConvex + | ConvexSet { + convexSetHull :: Shade' x + -- ^ If @p@ is in all intersectors, it must also be in the hull. + , convexSetIntersectors :: [Shade' x] + } + +ellipsoid :: Shade' x -> ConvexSet x +ellipsoid s = ConvexSet s [s] + +intersectors :: ConvexSet x -> Option (NonEmpty (Shade' x)) +intersectors (ConvexSet h []) = pure (h:|[]) +intersectors (ConvexSet _ (i:sts)) = pure (i:|sts) +intersectors _ = empty + +-- | Under intersection. +instance Refinable x => Semigroup (ConvexSet x) where + a<>b = sconcat (a:|[b]) + sconcat csets + | Option (Just allIntersectors) <- sconcat <$> Hask.traverse intersectors csets + , IntersectT ists <- rmTautologyIntersect perfectRefine $ IntersectT allIntersectors + , Option (Just hull') <- intersectShade's ists + = ConvexSet hull' (NE.toList ists) + | otherwise = EmptyConvex + where perfectRefine sh₁ sh₂ + | sh₁`subShade'`sh₂ = pure sh₁ + | sh₂`subShade'`sh₁ = pure sh₂ + | otherwise = empty + + + +itWhileJust :: (a -> Option a) -> a -> [a] +itWhileJust f x | Option (Just y) <- f x = x : itWhileJust f y +itWhileJust _ x = [x] + +dupHead :: NonEmpty a -> NonEmpty a +dupHead (x:|xs) = x:|x:xs + + +iterateFilterDEqn_static :: (WithField ℝ Manifold 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) => 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) + =<< intersectShade's + ( propagateDEqnSolution_loc f ((x,shy), NE.fromList ngbs) ) + +filterDEqnSolutions_static :: (WithField ℝ Manifold 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 + then pure shy + else ((shy<>) . ellipsoid) + <$> intersectShade's + ( propagateDEqnSolution_loc f + ((x,hull), second convexSetHull<$>NE.fromList ngbs) ) + >>= \case EmptyConvex -> empty + c -> pure c + + +data SolverNodeState y = SolverNodeInfo { + _solverNodeStatus :: ConvexSet y + , _solverNodeBadness :: ℝ + , _solverNodeAge :: Int + } +makeLenses ''SolverNodeState + + +type OldAndNew d = (Option d, [d]) + +oldAndNew :: OldAndNew d -> [d] +oldAndNew (Option (Just x), l) = x : l +oldAndNew (_, l) = l + +oldAndNew' :: OldAndNew d -> [(Bool, d)] +oldAndNew' (Option (Just x), l) = (True, x) : fmap (False,) l +oldAndNew' (_, l) = (False,) <$> l + + +filterDEqnSolutions_adaptive :: ∀ x y badness + . (WithField ℝ Manifold x, Refinable y, badness ~ ℝ) + => MetricChoice x -- ^ Scalar product on the domain, for regularising the web. + -> DifferentialEqn x y + -> (x -> Shade' y -> badness) + -> PointsWeb x (SolverNodeState y) + -> Option (PointsWeb x (SolverNodeState y)) +filterDEqnSolutions_adaptive mf f badness' oldState + = fmap (fromTopWebNodes mf . concat . fmap retraceBonds + . Hask.toList . webLocalInfo . webLocalInfo) + $ Hask.traverse (uncurry localChange) preproc'd + where preproc'd :: PointsWeb x ((WebLocally x (SolverNodeState y), [(Shade' y, badness)])) + preproc'd = fmap addPropagation $ webLocalInfo oldState + where addPropagation wl + | null neighbourHulls = (wl, []) + | otherwise = (wl, map (id&&&badness undefined) propFromNgbs) + where propFromNgbs = NE.toList $ propagateDEqnSolution_loc f + ( (thisPos, thisShy), NE.fromList neighbourHulls ) + thisPos = _thisNodeCoord wl :: x + thisShy = convexSetHull . _solverNodeStatus $ _thisNodeData wl + neighbourHulls = second (convexSetHull . _solverNodeStatus) . snd + <$> _nodeNeighbours wl + smallBadnessGradient, largeBadnessGradient :: ℝ + (smallBadnessGradient, largeBadnessGradient) + = ( badnessGradRated!!(n`div`4), badnessGradRated!!(n*3`div`4) ) + where n = length badnessGradRated + badnessGradRated = sort [ ngBad / bad + | ( LocalWebInfo { + _thisNodeData + = SolverNodeInfo _ bad _ + , _nodeNeighbours=ngbs } + , ngbProps) <- Hask.toList preproc'd + , (_, ngBad) <- ngbProps + , ngBad>bad ] + localChange :: WebLocally x (SolverNodeState y) -> [(Shade' y, badness)] + -> Option (OldAndNew (x, SolverNodeState y)) + localChange localInfo@LocalWebInfo{ + _thisNodeCoord = x + , _thisNodeData = SolverNodeInfo + shy@(ConvexSet hull _) prevBadness age + , _nodeNeighbours = ngbs + } + ngbProps + | null ngbs = return (pure (x, SolverNodeInfo shy prevBadness (age+1)), []) + | otherwise = do + let neighbourHulls = second (convexSetHull . _solverNodeStatus) . snd + <$> NE.fromList ngbs + (environAge, unfreshness) + = maximum&&&minimum $ age : (_solverNodeAge . snd . snd <$> ngbs) + case find (\(_, badnessN) + -> badnessN / prevBadness > smallBadnessGradient) + $ ngbProps of + Nothing | age < environAge -- point is an obsolete step-stone; + -> return (empty,empty) -- do not further use it. + _otherwise -> do + shy' <- ((shy<>) . ellipsoid) + <$> intersectShade's (fst <$> NE.fromList ngbProps) + newBadness <- case shy' of + EmptyConvex -> empty + ConvexSet hull' _ -> return $ badness x hull' + let updatedNode = SolverNodeInfo shy' newBadness (age+1) + stepStones <- + if unfreshness < 3 + then return [] + else fmap concat . forM (zip (snd<$>ngbs) ngbProps) + $ \( (vN, SolverNodeInfo (ConvexSet hullN _) + _ ageN) + , (_, nBadnessProp'd) ) -> do + case ageN of + _ | ageN > 0 + , badnessGrad <- nBadnessProp'd / prevBadness + , badnessGrad > largeBadnessGradient -> do + let stepV = vN^/2 + xStep = x .+~^ stepV + shyStep <- intersectShade's $ + propagateDEqnSolution_loc f + ( (xStep, hull) + , NE.cons (negateV stepV, hull) + $ fmap (\(vN',hullN') + -> (vN'^-^stepV, hullN') ) + neighbourHulls ) + return [( xStep + , SolverNodeInfo (ellipsoid shyStep) + (badness xStep shyStep) 1 + )] + _otherwise -> return [] + let updated = (x, updatedNode) + return $ (pure updated, stepStones) + + totalAge = maximum $ _solverNodeAge . _thisNodeData . fst <$> preproc'd + errTgtModulation = (1-) . (`mod'`1) . negate . sqrt $ fromIntegral totalAge + badness x = badness' x . (shadeNarrowness %~ (^* errTgtModulation)) + + retraceBonds :: WebLocally x (WebLocally x (OldAndNew (x, SolverNodeState y))) + -> [((x, [Needle x]), SolverNodeState y)] + retraceBonds locWeb@LocalWebInfo{ _thisNodeId = myId + , _thisNodeCoord = xOld + , _nodeLocalScalarProduct = locMetr } + = [ ( (x, fst<$>neighbourCandidates), snsy) + | (isOld, (x, snsy)) <- focused + , let neighbourCandidates + = [ (v,nnId) + | (_,ngb) <- knownNgbs + , (Option (Just v), nnId) + <- case oldAndNew $ ngb^.thisNodeData of + [] -> [ (xN.-~.x, nnId) + | (nnId, (_,nnWeb)) <- ngb^.nodeNeighbours + , nnId /= myId + , (xN,_) <- oldAndNew nnWeb ] + l -> [(xN.-~.x, ngb^.thisNodeId) | (xN,_) <- l] + ] + possibleConflicts = [ metricSq locMetr v + | (v,nnId)<-neighbourCandidates + , nnId > myId ] + , isOld || null possibleConflicts + || minimum possibleConflicts > oldMinDistSq / 4 + ] + where focused = oldAndNew' $ locWeb^.thisNodeData^.thisNodeData + knownNgbs = snd <$> locWeb^.nodeNeighbours + oldMinDistSq = minimum [ metricSq locMetr vOld + | (_,ngb) <- knownNgbs + , let Option (Just vOld) = ngb^.thisNodeCoord .-~. xOld + ] + + + +iterateFilterDEqn_adaptive :: (WithField ℝ Manifold 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. + -> PointsWeb x (Shade' y) -> [PointsWeb x (Shade' y)] +iterateFilterDEqn_adaptive mf f badness + = map (fmap (convexSetHull . _solverNodeStatus)) + . itWhileJust (filterDEqnSolutions_adaptive mf f badness) + . fmap (\((x,shy),_) -> SolverNodeInfo (ellipsoid shy) + (badness x shy) + 1 + ) + . localFocusWeb + + diff --git a/Data/SimplicialComplex.hs b/Data/SimplicialComplex.hs --- a/Data/SimplicialComplex.hs +++ b/Data/SimplicialComplex.hs @@ -64,22 +64,15 @@ import Data.List hiding (filter, all, elem) import Data.Maybe -import qualified Data.Map as Map import qualified Data.Vector as Arr -import Data.List.NonEmpty (NonEmpty(..)) import Data.List.FastNub import qualified Data.List.NonEmpty as NE import Data.Semigroup import Data.Ord (comparing) -import Data.VectorSpace -import Data.LinearMap import Data.LinearMap.Category -import Data.Void import Data.Tagged -import Data.Proxy -import Data.Manifold.Types import Data.Manifold.Types.Primitive ((^), empty) import Data.Manifold.PseudoAffine @@ -101,7 +94,6 @@ import Control.Monad.Constrained import Data.Foldable.Constrained -import GHC.Generics (Generic) infixr 5 :<|, .<. diff --git a/Util/LtdShow.hs b/Util/LtdShow.hs --- a/Util/LtdShow.hs +++ b/Util/LtdShow.hs @@ -6,7 +6,6 @@ module Util.LtdShow (LtdShow(..)) where import qualified Data.Vector as V -import Data.Vector(fromList, toList, (!), singleton) type Array = V.Vector diff --git a/manifolds.cabal b/manifolds.cabal --- a/manifolds.cabal +++ b/manifolds.cabal @@ -1,5 +1,5 @@ Name: manifolds -Version: 0.2.2.0 +Version: 0.2.3.0 Category: Math Synopsis: Coordinate-free hypersurfaces Description: Manifolds, a generalisation of the notion of “smooth curves” or surfaces, @@ -49,6 +49,7 @@ , void , tagged , deepseq + , microlens >= 0.4 && <= 0.5, microlens-th , trivial-constraint >= 0.4 , constrained-categories >= 0.2.3 && < 0.3 other-extensions: FlexibleInstances @@ -66,6 +67,7 @@ Data.Manifold.PseudoAffine Data.Manifold.TreeCover Data.Manifold.Web + Data.Manifold.DifferentialEquation Data.SimplicialComplex Data.LinearMap.HerMetric Data.Function.Differentiable @@ -82,6 +84,7 @@ Data.Function.Differentiable.Data Data.Function.Affine Data.VectorSpace.FiniteDimensional + Control.Monad.Trans.OuterMaybe Util.Associate Util.LtdShow default-language: Haskell2010