diff --git a/Data/LinearMap/HerMetric.hs b/Data/LinearMap/HerMetric.hs --- a/Data/LinearMap/HerMetric.hs +++ b/Data/LinearMap/HerMetric.hs @@ -28,7 +28,7 @@ -- * One-dimensional axes and product spaces , factoriseMetric, factoriseMetric' , productMetric, productMetric' - , metricAsLength, metric'AsLength + , metricAsLength, metricFromLength, metric'AsLength -- * Utility for metrics , transformMetric, transformMetric' , dualiseMetric, dualiseMetric' @@ -214,7 +214,7 @@ transformMetric :: (HasMetric v, HasMetric w, Scalar v ~ Scalar w) => (w :-* v) -> HerMetric v -> HerMetric w transformMetric _ (HerMetric Nothing) = HerMetric Nothing -transformMetric t (HerMetric (Just m)) = matrixMetric $ HMat.tr tmat HMat.<> m HMat.<> tmat +transformMetric t (HerMetric (Just m)) = matrixMetric $ tmat HMat.<> m HMat.<> HMat.tr tmat where tmat = asPackedMatrix t transformMetric' :: ( HasMetric v, HasMetric w, Scalar v ~ Scalar w ) @@ -515,8 +515,11 @@ metricAsLength :: HerMetric ℝ -> ℝ metricAsLength = recip . (`metric`1) +metricFromLength :: ℝ -> HerMetric ℝ +metricFromLength = projector . recip + metric'AsLength :: HerMetric' ℝ -> ℝ -metric'AsLength = recip . (`metric'`1) +metric'AsLength = recip . (`metric'`1) -- do we really want `recip` here? spanHilbertSubspace :: ∀ s v w diff --git a/Data/Manifold/PseudoAffine.hs b/Data/Manifold/PseudoAffine.hs --- a/Data/Manifold/PseudoAffine.hs +++ b/Data/Manifold/PseudoAffine.hs @@ -52,6 +52,7 @@ , PseudoAffine(..) -- * Regions within a manifold , Region + , smoothIndicator -- * Hierarchy of manifold-categories -- ** Everywhere differentiable functions , Differentiable @@ -71,7 +72,11 @@ , EuclidSpace -- * Misc , palerp - , discretisePath + , discretisePathIn + , discretisePathSegs + , continuousIntervals + , regionOfContinuityAround + , analyseLocalBehaviour ) where @@ -301,23 +306,103 @@ -discretisePath :: WithField ℝ Manifold x +discretisePathIn :: WithField ℝ Manifold x => Int -- ^ Limit the number of steps taken in either direction. Note this will not cap the resolution but /length/ of the discretised path. + -> Region ℝ ℝ -- ^ Parameter interval of interest -> RieMetric x -- ^ Inaccuracy allowance /ε/. -> (Differentiable ℝ ℝ x) -- ^ Path specification. -> [(ℝ,x)] -- ^ Trail of points along the path, such that a linear interpolation deviates nowhere by more as /ε/. -discretisePath nLim m (Differentiable f) - = reverse (tail . take nLim $ traceFwd 0 (-1)) ++ take nLim (traceFwd 0 1) +discretisePathIn nLim (Region xm rLim) m (Differentiable f) + = reverse (tail . take nLim $ traceFwd xm (-1)) ++ take nLim (traceFwd xm 1) where traceFwd x₀ dir - | abs x₀ > 1e+100 = [(x₀, fx₀)] - | otherwise = (x₀, fx₀) : traceFwd xn dir + | rnfn x₀ < 0 = [] + | abs x₀ > hugeℝVal = [(x₀, fx₀)] + | otherwise = (x₀, fx₀) : traceFwd xn dir where (fx₀, _, δx²) = f x₀ εx = m fx₀ χ = metric (δx² εx) 1 xn = x₀ + dir * min (abs x₀+1) (recip χ) + rnfn = case rLim of + GlobalRegion -> const 1 + PreRegion (Differentiable pmbf) -> pmbf >>> \(q,_,_)->q + +discretisePathSegs :: WithField ℝ Manifold x + => Int -- ^ Maximum number of path segments and/or points per segment. + -> RieMetric x -- ^ Inaccuracy allowance /ε/. + -> RWDiffable ℝ ℝ x -- ^ Path specification. + -> [[(ℝ,x)]] -- ^ Trail of points along the path, such that a linear interpolation deviates nowhere by more as /ε/. +discretisePathSegs nLim m (RWDiffable f) = jumpsFwd nLim 0 (True,True) + where jumpsFwd nLim' x₀ (goL,goR) + | abs x₀ > hugeℝVal = [] + | Option Nothing <- fq₀ = error "`discretisePathSegs` not yet implemented for partial functions outside of a null set." + | xr < -hugeℝVal + || xr < hugeℝVal = [pseg] + | not goL = pseg : jumpR + | not goR = pseg : jumpL + | otherwise = pseg : (zip jumpL jumpR >>= \(l,r)->[l,r]) + where (r₀, fq₀) = f x₀ + Option (Just lf) = fq₀ + pseg = first (subtract x₀) <$> + discretisePathIn nLim' (Region x₀ r₀) m (lf . actuallyAffine x₀ idL) + ((xl,_):(xpl,_):_) = pseg + ((xr,_):(xpr,_):_) = reverse pseg + jumpR = jumpsFwd (nLim'-1) (xr*2-xpr) (False,goR) + jumpL = jumpsFwd (nLim'-1) (xl*2-xpl) (goL,False) + + +continuousIntervals :: RWDiffable ℝ ℝ x -> (ℝ,ℝ) -> [(ℝ,ℝ)] +continuousIntervals (RWDiffable f) (xl,xr) = enter xl + where enter x₀ = case f x₀ of + (GlobalRegion, _) -> [(xl,xr)] + (PreRegion r₀, _) -> exit r₀ x₀ + where exit :: Differentiable ℝ ℝ ℝ -> ℝ -> [(ℝ,ℝ)] + exit (Differentiable r) x + | x > xr = [(x₀,xr)] + | y' > 0 = exit (Differentiable r) + (x + metricAsLength (δ (metricFromLength y))) + | -y/y' < 1e-10 = (x₀,x) : enter (x + min 1e-100 (abs x * 1e-8)) + | otherwise = exit (Differentiable r) xn + where (y, y'm, δ) = r x + xn = bisBack $ x - y/y' + where bisBack xq + | ybm > 0 = xbm + | otherwise = bisBack xbm + where (ybm, _, _) = r xbm + xbm = (xq*9 + x)/10 + y' = lapply y'm 1 + +analyseLocalBehaviour :: + RWDiffable ℝ ℝ ℝ + -> ℝ -- ^ /x/₀ value. + -> Option ( (ℝ,ℝ) + , ℝ->Option ℝ ) -- ^ /f/ /x/₀, derivative (i.e. Taylor-1-coefficient), + -- and reverse propagation of /O/ (/δ/²) bound. +analyseLocalBehaviour (RWDiffable f) x₀ = case f x₀ of + (_, Option Nothing) -> Hask.empty + (_, Option (Just (Differentiable fd))) -> return $ + let (fx, j, δf) = fd x₀ + epsprop ε + | ε>0 = case metric (δf $ metricFromLength ε) 1 of + 0 -> Hask.empty + δ' -> return $ recip δ' + | otherwise = pure 0 + in ((fx, lapply j 1), epsprop) +-- | Represent a 'Region' by a smooth function which is positive within the region, +-- and crosses zero at the boundary. +smoothIndicator :: LocallyScalable ℝ q => Region ℝ q -> Differentiable ℝ q ℝ +smoothIndicator (Region _ GlobalRegion) = const 1 +smoothIndicator (Region _ (PreRegion r)) = r +regionOfContinuityAround :: RWDiffable ℝ q x -> q -> Region ℝ q +regionOfContinuityAround (RWDiffable f) q = Region q . fst . f $ q + + + +hugeℝVal :: ℝ +hugeℝVal = 1e+100 + #define deriveAffine(t) \ instance Semimanifold (t) where { \ type Needle (t) = Diff (t); \ @@ -515,7 +600,7 @@ dev_ε_δ :: RealDimension a => (a -> a) -> LinDevPropag a a -dev_ε_δ f d = let ε = 1 / metric d 1 in projector $ 1 / sqrt (f ε) +dev_ε_δ f d = let ε = 1 / metric d 1 in projector $ 1 / f ε -- | The category of differentiable functions between manifolds over scalar @s@. -- diff --git a/manifolds.cabal b/manifolds.cabal --- a/manifolds.cabal +++ b/manifolds.cabal @@ -1,5 +1,5 @@ Name: manifolds -Version: 0.1.5.1 +Version: 0.1.5.2 Category: Math Synopsis: Coordinate-free hypersurfaces Description: Manifolds, a generalisation of the notion of “smooth curves” or surfaces,