diff --git a/Data/CoNat.hs b/Data/CoNat.hs --- a/Data/CoNat.hs +++ b/Data/CoNat.hs @@ -28,7 +28,15 @@ {-# LANGUAGE DataKinds #-} {-# LANGUAGE PolyKinds #-} -module Data.CoNat where +module Data.CoNat ( Nat(..), natToInt, fromNat + , natTagLast, natTagPænultimate, natTagAntepænultimate + , tryToMatchT, tryToMatchTT, tryToMatchTTT + , ftorTryToMatch, ftorTryToMatchT, ftorTryToMatchTT + , KnownNat(..) + , Range(..) + , FreeVect(..), (^)(), freeVector, freeCons, freeSnoc + , replicVector, indices, perfectZipWith, freeRotate + , ) where import Data.Tagged import Data.Semigroup @@ -109,8 +117,8 @@ instance KnownNat Z where theNat = Tagged Z theNatN = Tagged 0 - cozero = pure; cosucc _ = Hask.empty; fCosucc _ = Hask.empty - cozeroT = pure; cosuccT _ = Hask.empty; fCosuccT _ = Hask.empty + cozero = pure; cosucc _ = empty; fCosucc _ = empty + cozeroT = pure; cosuccT _ = empty; fCosuccT _ = empty coNat f _ = f; coNatT f _ = f coInduce s _ = s coInduceT s _ = s @@ -121,13 +129,13 @@ => (∀ j . KnownNat j => b j -> b (S j)) -> b k -> Option (b Z) ttmZ sc nf = case k of Z -> return $ unsafeCoerce nf - S _ -> Hask.empty + S _ -> empty where (Tagged k) = theNat :: Tagged k Nat instance (KnownNat n) => KnownNat (S n) where theNat = natSelfSucc theNatN = natSelfSuccN - cozero _ = Hask.empty; cosucc v = pure v; fCosucc v = v - cozeroT _ = Hask.empty; cosuccT v = pure v; fCosuccT v = v + cozero _ = empty; cosucc v = pure v; fCosucc v = v + cozeroT _ = empty; cosuccT v = pure v; fCosuccT v = v coNat _ f = f; coNatT _ f = f coInduce s f = f $ coInduce s f coInduceT s f = f $ coInduceT s f @@ -138,7 +146,7 @@ => (∀ j . KnownNat j => b j -> b (S j)) -> b k -> Option (b (S n)) ttmS sc nf | k == sn = return $ unsafeCoerce nf | k <= sn = tryToMatch sc $ sc nf - | otherwise = Hask.empty + | otherwise = empty where (Tagged k) = theNatN :: Tagged k Int (Tagged sn) = theNatN :: Tagged (S n) Int @@ -214,7 +222,7 @@ clipToRange :: forall n . KnownNat n => Int -> Option (Range n) clipToRange = \k -> if k < n then Hask.pure $ InRange n - else Hask.empty + else empty where (Tagged n) = theNatN :: Tagged n Int instance KnownNat n => HasTrie (Range n) where @@ -270,7 +278,7 @@ freeVector :: forall l n x . (KnownNat n, Hask.Foldable l) => l x -> Option (FreeVect n x) freeVector c' | List.length c == n = pure . FreeVect $ Arr.fromList c - | otherwise = Hask.empty + | otherwise = empty where (Tagged n) = theNatN :: Tagged n Int c = Hask.toList c' @@ -312,3 +320,7 @@ instance (Monoidal f (->) (->)) => Hask.Applicative (AsHaskFunctor f) where pure x = fmap (const x) . AsHaskFunctor $ pureUnit () AsHaskFunctor fs <*> AsHaskFunctor xs = AsHaskFunctor . fmap (uncurry ($)) $ fzip (fs, xs) + + +empty :: Hask.Alternative m => m a +empty = Hask.empty diff --git a/Data/Function/Differentiable.hs b/Data/Function/Differentiable.hs new file mode 100644 --- /dev/null +++ b/Data/Function/Differentiable.hs @@ -0,0 +1,1172 @@ +-- | +-- Module : Data.Function.Differentiable +-- Copyright : (c) Justus Sagemüller 2015 +-- License : GPL v3 +-- +-- Maintainer : (@) sagemueller $ geo.uni-koeln.de +-- Stability : experimental +-- Portability : portable +-- + +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE LiberalTypeSynonyms #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE TupleSections #-} +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE PatternGuards #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE UnicodeSyntax #-} +{-# LANGUAGE MultiWayIf #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE RecordWildCards #-} +{-# LANGUAGE CPP #-} + + +module Data.Function.Differentiable ( + -- * Regions within a manifold + Region + , smoothIndicator + -- * Hierarchy of manifold-categories + -- ** Everywhere differentiable functions + , Differentiable + -- ** Almost everywhere diff'able funcs + , PWDiffable + -- ** Region-wise defined diff'able funcs + , RWDiffable + -- * Misc + , discretisePathIn + , discretisePathSegs + , continuityRanges + , regionOfContinuityAround + , analyseLocalBehaviour + ) 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 +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 +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 + +import Control.Category.Constrained.Prelude hiding ((^)) +import Control.Arrow.Constrained +import Control.Monad.Constrained +import Data.Foldable.Constrained + + + + + + +discretisePathIn :: WithField ℝ Manifold y + => Int -- ^ Limit the number of steps taken in either direction. Note this will not cap the resolution but /length/ of the discretised path. + -> ℝInterval -- ^ Parameter interval of interest. + -> (RieMetric ℝ, RieMetric y) -- ^ Inaccuracy allowance /ε/. + -> (Differentiable ℝ ℝ y) -- ^ Path specification. + -> [(ℝ,y)] -- ^ Trail of points along the path, such that a linear interpolation deviates nowhere by more as /ε/. +discretisePathIn nLim (xl, xr) (mx,my) (Differentiable f) + = reverse (tail . take nLim $ traceFwd xl xm (-1)) + ++ take nLim (traceFwd xr xm 1) + where traceFwd xlim x₀ dir + | signum (x₀-xlim) == signum dir = [(xlim, fxlim)] + | otherwise = (x₀, fx₀) : traceFwd xlim (x₀+xstep) dir + where (fx₀, jf, δx²) = f x₀ + εx = my fx₀ `extendMetric` lapply jf (metricAsLength $ mx x₀) + χ = metric (δx² εx) 1 + xstep = dir * min (abs x₀+1) (recip χ) + (fxlim, _, _) = f xlim + xm = (xr + xl) / 2 + +type ℝInterval = (ℝ,ℝ) + +continuityRanges :: WithField ℝ Manifold y + => Int -- ^ Max number of exploration steps per region + -> RieMetric ℝ -- ^ Needed resolution of boundaries + -> ℝInterval -- ^ Interval to explore + -> RWDiffable ℝ ℝ y -- ^ Function to investigate + -> ([ℝInterval], [ℝInterval]) -- ^ Subintervals on which the function is guaranteed continuous. +continuityRanges nLim δbf (limL,limR) (RWDiffable f) + | (GlobalRegion, _) <- f xc + = ([], [(-huge,huge)]) + | otherwise = glueMid (go xc (-1)) (go xc 1) + where go x₀ dir + | yq₀ <= abs (lapply jq₀ 1 * step₀) + = go (x₀ + step₀/2) dir + | otherwise = exit nLim dir x₀ + where (PreRegion (Differentiable r₀), fq₀) = f x₀ + (yq₀, jq₀, δyq₀) = r₀ x₀ + step₀ = dir/metric (δbf x₀) 1 + exit _ d xq + | xq < limL = exit 0 d limL + | xq > limR = exit 0 d limR + exit 0 _ xq + | not definedHere = [] + | xq < xc = [(xq,x₀)] + | otherwise = [(x₀,xq)] + exit nLim' dir' xq + | yq₁<0 || as_devεδ δyq yq₁=abs stepp + , resoHere stepp<1 = (if definedHere + then ((min x₀ xq₁, max x₀ xq₁):) + else id) $ go xq₂ dir + | otherwise = exit (nLim'-1) dir xq₁ + where (yq, jq, δyq) = r₀ xq + xq₁ = xq + stepp + xq₂ = xq₁ + stepp + yq₁ = yq + f'x*stepp + yq₂ = yq₁ + f'x*stepp + f'x = lapply jq 1 + stepp | f'x*dir < 0 = -0.9 * abs dir' * yq/f'x + | otherwise = dir' * as_devεδ δyq yq -- TODO: memoise in `exit` recursion + resoHere = metricSq $ δbf xq + resoStep = dir/sqrt(resoHere 1) + definedHere = case fq₀ of + Option (Just _) -> True + Option Nothing -> False + glueMid ((l,le):ls) ((re,r):rs) | le==re = (ls, (l,r):rs) + glueMid l r = (l,r) + huge = exp $ fromIntegral nLim + xc | limL*2 /= limL, limR*2 /= limR = (limR+limL)/2 + | otherwise = max limL . min limR $ 0 + +discretisePathSegs :: WithField ℝ Manifold y + => Int -- ^ Maximum number of path segments and/or points per segment. + -> ( RieMetric ℝ + , RieMetric y ) -- ^ Inaccuracy allowance /δ/ for arguments + -- (mostly relevant for resolution of discontinuity boundaries – + -- consider it a “safety margin from singularities”), + -- and /ε/ for results in the target space. + -> ℝInterval -- ^ Interval of interest. You can make this “infinitely large”. + -> RWDiffable ℝ ℝ y -- ^ Path specification. + -> ([[(ℝ,y)]], [[(ℝ,y)]]) -- ^ Discretised paths: continuous segments in either direction +discretisePathSegs nLim (mx,my) rng@(limL,limR) f@(RWDiffable ff) + = ( map discretise $ trimToRange ivsL + , map discretise $ trimToRange ivsR ) + where (ivsL, ivsR) = continuityRanges nLim mx rng f + trimToRange = map ( \(l,r) -> (max limL l, min limR r) ) + . Data.List.filter ( \(l,r) -> llimL ) + discretise rng@(l,r) = discretisePathIn nLim rng (mx,my) fr + where (_, Option (Just fr)) = ff $ (l+r)/2 + + +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 + (r, Option (Just (Differentiable fd))) + | inRegion r x₀ -> return $ + let (fx, j, δf) = fd x₀ + epsprop ε + | ε>0 = case metric (δf $ metricFromLength ε) 1 of + 0 -> empty + δ' -> return $ recip δ' + | otherwise = pure 0 + in ((fx, lapply j 1), epsprop) + _ -> empty + where inRegion GlobalRegion _ = True + inRegion (PreRegion (Differentiable rf)) x + | (yr,_,_) <- rf x = yr>0 + +-- | 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 + + + + + + +type LinDevPropag d c = Metric c -> Metric d + +unsafe_dev_ε_δ :: RealDimension a + => String -> (a -> a) -> LinDevPropag a a +unsafe_dev_ε_δ errHint f d + = let ε'² = metricSq d 1 + in if ε'²>0 + then let δ = f . sqrt $ recip ε'² + in if δ > 0 + then projector $ recip δ + else error $ "ε-δ propagator function for " + ++errHint++", with ε=" + ++show(sqrt $ recip ε'²) + ++ " gives non-positive δ="++show δ++"." + else zeroV +dev_ε_δ :: RealDimension a + => (a -> a) -> Metric a -> Option (Metric a) +dev_ε_δ f d = let ε'² = metricSq d 1 + in if ε'²>0 + then let δ = f . sqrt $ recip ε'² + in if δ > 0 + then pure . projector $ recip δ + else empty + else pure zeroV + +as_devεδ :: RealDimension a => LinDevPropag a a -> a -> a +as_devεδ ldp ε | ε>0 + , δ'² <- metricSq (ldp . projector $ recip ε) 1 + , δ'² > 0 + = sqrt $ recip δ'² + | otherwise = 0 + +-- | The category of differentiable functions between manifolds over scalar @s@. +-- +-- As you might guess, these offer /automatic differentiation/ of sorts (basically, +-- simple forward AD), but that's in itself is not really the killer feature here. +-- More interestingly, we actually have the (à la Curry-Howard) /proof/ +-- built in: the function /f/ has at /x/₀ derivative /f'ₓ/₀, +-- if, for¹ /ε/>0, there exists /δ/ such that +-- |/f/ /x/ − (/f/ /x/₀ + /x/⋅/f'ₓ/₀)| < /ε/ +-- for all |/x/ − /x/₀| < /δ/. +-- +-- Observe that, though this looks quite similar to the standard definition +-- of differentiability, it is not equivalent thereto – in fact it does +-- not prove any analytic properties at all. To make it equivalent, we need +-- a lower bound on /δ/: simply /δ/ gives us continuity, and for +-- continuous differentiability, /δ/ must grow at least like √/ε/ +-- for small /ε/. Neither of these conditions are enforced by the type system, +-- but we do require them for any allowed values because these proofs are obviously +-- tremendously useful – for instance, you can have a root-finding algorithm +-- and actually be sure you get /all/ solutions correctly, not just /some/ that are +-- (hopefully) the closest to some reference point you'd need to laborously define! +-- +-- Unfortunately however, this also prevents doing any serious algebra etc. with the +-- category, because even something as simple as division necessary introduces singularities +-- where the derivatives must diverge. +-- Not to speak of many trigonometric e.g. trigonometric functions that +-- are undefined on whole regions. The 'PWDiffable' and 'RWDiffable' categories have explicit +-- handling for those issues built in; you may simply use these categories even when +-- you know the result will be smooth in your relevant domain (or must be, for e.g. +-- physics reasons). +-- +-- ¹(The implementation does not deal with /ε/ and /δ/ as difference-bounding +-- reals, but rather as metric tensors that define a boundary by prohibiting the +-- overlap from exceeding one; this makes the concept actually work on general manifolds.) +newtype Differentiable s d c + = Differentiable { runDifferentiable :: + d -> ( c -- function value + , Needle d :-* Needle c -- Jacobian + , LinDevPropag d c -- Metric showing how far you can go + -- from x₀ without deviating from the + -- Taylor-1 approximation by more than + -- some error margin + ) } +type (-->) = Differentiable ℝ + + +instance (MetricScalar s) => Category (Differentiable s) where + type Object (Differentiable s) o = LocallyScalable s o + id = Differentiable $ \x -> (x, idL, const zeroV) + Differentiable f . Differentiable g = Differentiable $ + \x -> let (y, g', devg) = g x + (z, f', devf) = f y + devfg δz = let δy = transformMetric f' δz + εy = devf δz + in transformMetric g' εy ^+^ devg δy ^+^ devg εy + in (z, f'*.*g', devfg) + + +instance (RealDimension s) => EnhancedCat (->) (Differentiable s) where + arr (Differentiable f) x = let (y,_,_) = f x in y + +instance (MetricScalar s) => Cartesian (Differentiable s) where + type UnitObject (Differentiable s) = ZeroDim s + swap = Differentiable $ \(x,y) -> ((y,x), lSwap, const zeroV) + where lSwap = linear swap + attachUnit = Differentiable $ \x -> ((x, Origin), lAttachUnit, const zeroV) + where lAttachUnit = linear $ \x -> (x, Origin) + detachUnit = Differentiable $ \(x, Origin) -> (x, lDetachUnit, const zeroV) + where lDetachUnit = linear $ \(x, Origin) -> x + regroup = Differentiable $ \(x,(y,z)) -> (((x,y),z), lRegroup, const zeroV) + where lRegroup = linear regroup + regroup' = Differentiable $ \((x,y),z) -> ((x,(y,z)), lRegroup, const zeroV) + where lRegroup = linear regroup' + + +instance (MetricScalar s) => Morphism (Differentiable s) where + Differentiable f *** Differentiable g = Differentiable h + where h (x,y) = ((fx, gy), lPar, devfg) + where (fx, f', devf) = f x + (gy, g', devg) = g y + devfg δs = transformMetric lfst δx + ^+^ transformMetric lsnd δy + where δx = devf $ transformMetric lcofst δs + δy = devg $ transformMetric lcosnd δs + lPar = linear $ lapply f'***lapply g' + lfst = linear fst; lsnd = linear snd + lcofst = linear (,zeroV); lcosnd = linear (zeroV,) + + +instance (MetricScalar s) => PreArrow (Differentiable s) where + terminal = Differentiable $ \_ -> (Origin, zeroV, const zeroV) + fst = Differentiable $ \(x,_) -> (x, lfst, const zeroV) + where lfst = linear fst + snd = Differentiable $ \(_,y) -> (y, lsnd, const zeroV) + where lsnd = linear snd + Differentiable f &&& Differentiable g = Differentiable h + where h x = ((fx, gx), lFanout, devfg) + where (fx, f', devf) = f x + (gx, g', devg) = g x + devfg δs = (devf $ transformMetric lcofst δs) + ^+^ (devg $ transformMetric lcosnd δs) + lFanout = linear $ lapply f'&&&lapply g' + lcofst = linear (,zeroV); lcosnd = linear (zeroV,) + + +instance (MetricScalar s) => WellPointed (Differentiable s) where + unit = Tagged Origin + globalElement x = Differentiable $ \Origin -> (x, zeroV, const zeroV) + const x = Differentiable $ \_ -> (x, zeroV, const zeroV) + + + +type DfblFuncValue s = GenericAgent (Differentiable s) + +instance (MetricScalar s) => HasAgent (Differentiable s) where + alg = genericAlg + ($~) = genericAgentMap +instance (MetricScalar s) => CartesianAgent (Differentiable s) where + alg1to2 = genericAlg1to2 + alg2to1 = genericAlg2to1 + alg2to2 = genericAlg2to2 +instance (MetricScalar s) + => PointAgent (DfblFuncValue s) (Differentiable s) a x where + point = genericPoint + + + +actuallyLinear :: ( WithField s LinearManifold x, WithField s LinearManifold y ) + => (x:-*y) -> Differentiable s x y +actuallyLinear f = Differentiable $ \x -> (lapply f x, f, const zeroV) + +actuallyAffine :: ( WithField s LinearManifold x, WithField s AffineManifold y ) + => y -> (x:-*Diff y) -> Differentiable s x y +actuallyAffine y₀ f = Differentiable $ \x -> (y₀ .+^ lapply f x, f, const zeroV) + + +dfblFnValsFunc :: ( LocallyScalable s c, LocallyScalable s c', LocallyScalable s d + , v ~ Needle c, v' ~ Needle c' + , ε ~ HerMetric v, ε ~ HerMetric v' ) + => (c' -> (c, v':-*v, ε->ε)) -> DfblFuncValue s d c' -> DfblFuncValue s d c +dfblFnValsFunc f = (Differentiable f $~) + +dfblFnValsCombine :: forall d c c' c'' v v' v'' ε ε' ε'' s. + ( LocallyScalable s c, LocallyScalable s c', LocallyScalable s c'' + , LocallyScalable s d + , v ~ Needle c, v' ~ Needle c', v'' ~ Needle c'' + , ε ~ HerMetric v , ε' ~ HerMetric v' , ε'' ~ HerMetric v'', ε~ε', ε~ε'' ) + => ( c' -> c'' -> (c, (v',v''):-*v, ε -> (ε',ε'')) ) + -> DfblFuncValue s d c' -> DfblFuncValue s d c'' -> DfblFuncValue s d c +dfblFnValsCombine cmb (GenericAgent (Differentiable f)) + (GenericAgent (Differentiable g)) + = GenericAgent . Differentiable $ + \d -> let (c', f', devf) = f d + (c'', g', devg) = g d + (c, h', devh) = cmb c' c'' + h'l = h' *.* lcofst; h'r = h' *.* lcosnd + in ( c + , h' *.* linear (lapply f' &&& lapply g') + , \εc -> let εc' = transformMetric h'l εc + εc'' = transformMetric h'r εc + (δc',δc'') = devh εc + in devf εc' ^+^ devg εc'' + ^+^ transformMetric f' δc' + ^+^ transformMetric g' δc'' + ) + where lcofst = linear(,zeroV) + lcosnd = linear(zeroV,) + + + + + +instance (WithField s LinearManifold v, LocallyScalable s a, Floating s) + => AdditiveGroup (DfblFuncValue s a v) where + zeroV = point zeroV + (^+^) = dfblFnValsCombine $ \a b -> (a^+^b, lPlus, const zeroV) + where lPlus = linear $ uncurry (^+^) + negateV = dfblFnValsFunc $ \a -> (negateV a, lNegate, const zeroV) + where lNegate = linear negateV + +instance (RealDimension n, LocallyScalable n a) + => Num (DfblFuncValue n a n) where + fromInteger i = point $ fromInteger i + (+) = dfblFnValsCombine $ \a b -> (a+b, lPlus, const zeroV) + where lPlus = linear $ uncurry (+) + (*) = dfblFnValsCombine $ + \a b -> ( a*b + , linear $ \(da,db) -> a*db + b*da + , \d -> let d¹₂ = sqrt d in (d¹₂,d¹₂) + -- ε δa δb = (a+δa)·(b+δb) - (a·b + (a·δa + b·δb)) + -- = δa·δb + -- so choose δa = δb = √ε + ) + negate = dfblFnValsFunc $ \a -> (negate a, lNegate, const zeroV) + where lNegate = linear negate + abs = dfblFnValsFunc dfblAbs + where dfblAbs a + | a>0 = (a, idL, unsafe_dev_ε_δ("abs "++show a) $ \ε -> a + ε/2) + | a<0 = (-a, negateV idL, unsafe_dev_ε_δ("abs "++show a) $ \ε -> ε/2 - a) + | otherwise = (0, zeroV, (^/ sqrt 2)) + signum = dfblFnValsFunc dfblSgn + where dfblSgn a + | a>0 = (1, zeroV, unsafe_dev_ε_δ("signum "++show a) $ const a) + | a<0 = (-1, zeroV, unsafe_dev_ε_δ("signum "++show a) $ \_ -> -a) + | otherwise = (0, zeroV, const $ projector 1) + + + +-- VectorSpace instance is more problematic than you'd think: multiplication +-- requires the allowed-deviation backpropagators to be split as square +-- roots, but the square root of a nontrivial-vector-space metric requires +-- an eigenbasis transform, which we have not implemented yet. +-- +-- instance (WithField s LinearManifold v, LocallyScalable s a, Floating s) +-- => VectorSpace (DfblFuncValue s a v) where +-- type Scalar (DfblFuncValue s a v) = DfblFuncValue s a (Scalar v) +-- (*^) = dfblFnValsCombine $ \μ v -> (μ*^v, lScl, \ε -> (ε ^* sqrt 2, ε ^* sqrt 2)) +-- where lScl = linear $ uncurry (*^) + + +-- | Important special operator needed to compute intersection of 'Region's. +minDblfuncs :: (LocallyScalable s m, RealDimension s) + => Differentiable s m s -> Differentiable s m s -> Differentiable s m s +minDblfuncs (Differentiable f) (Differentiable g) = Differentiable h + where h x + | fx < gx = ( fx, jf + , \d -> devf d ^+^ devg d + ^+^ transformMetric δj + (projector . recip $ recip(metric d 1) + gx - fx) ) + | fx > gx = ( gx, jg + , \d -> devf d ^+^ devg d + ^+^ transformMetric δj + (projector . recip $ recip(metric d 1) + fx - gx) ) + | otherwise = ( fx, (jf^+^jg)^/2 + , \d -> devf d ^+^ devg d + ^+^ transformMetric δj d ) + where (fx, jf, devf) = f x + (gx, jg, devg) = g x + δj = jf ^-^ jg + + +postEndo :: ∀ c a b . (HasAgent c, Object c a, Object c b) + => c a a -> GenericAgent c b a -> GenericAgent c b a +postEndo = genericAgentMap + + +-- | A pathwise connected subset of a manifold @m@, whose tangent space has scalar @s@. +data Region s m = Region { regionRefPoint :: m + , regionRDef :: PreRegion s m } + +-- | A 'PreRegion' needs to be associated with a certain reference point ('Region' +-- includes that point) to define a connected subset of a manifold. +data PreRegion s m where + GlobalRegion :: PreRegion s m + PreRegion :: (Differentiable s m s) -- A function that is positive at reference point /p/, + -- decreases and crosses zero at the region's + -- boundaries. (If it goes positive again somewhere + -- else, these areas shall /not/ be considered + -- belonging to the (by definition connected) region.) + -> PreRegion s m + +-- | Set-intersection of regions would not be guaranteed to yield a connected result +-- or even have the reference point of one region contained in the other. This +-- combinator assumes (unchecked) that the references are in a connected +-- sub-intersection, which is used as the result. +unsafePreRegionIntersect :: (RealDimension s, LocallyScalable s a) + => PreRegion s a -> PreRegion s a -> PreRegion s a +unsafePreRegionIntersect GlobalRegion r = r +unsafePreRegionIntersect r GlobalRegion = r +unsafePreRegionIntersect (PreRegion ra) (PreRegion rb) = PreRegion $ minDblfuncs ra rb + +-- | Cartesian product of two regions. +regionProd :: (RealDimension s, LocallyScalable s a, LocallyScalable s b) + => Region s a -> Region s b -> Region s (a,b) +regionProd (Region a₀ ra) (Region b₀ rb) = Region (a₀,b₀) (preRegionProd ra rb) + +-- | Cartesian product of two pre-regions. +preRegionProd :: (RealDimension s, LocallyScalable s a, LocallyScalable s b) + => PreRegion s a -> PreRegion s b -> PreRegion s (a,b) +preRegionProd GlobalRegion GlobalRegion = GlobalRegion +preRegionProd GlobalRegion (PreRegion rb) = PreRegion $ rb . snd +preRegionProd (PreRegion ra) GlobalRegion = PreRegion $ ra . fst +preRegionProd (PreRegion ra) (PreRegion rb) = PreRegion $ minDblfuncs (ra.fst) (rb.snd) + + +positivePreRegion, negativePreRegion :: (RealDimension s) => PreRegion s s +positivePreRegion = PreRegion $ Differentiable prr + where prr x = ( 1 - 1/xp1 + , (1/xp1²) *^ idL + , unsafe_dev_ε_δ("positivePreRegion@"++show x) δ ) + -- ε = (1 − 1/(1+x)) + (-δ · 1/(x+1)²) − (1 − 1/(1+x−δ)) + -- = 1/(1+x−δ) − 1/(1+x) − δ · 1/(x+1)² + -- + -- ε·(1+x−δ) = 1 − (1+x−δ)/(1+x) − δ·(1+x-δ)/(x+1)² + -- ε·(1+x) − ε·δ = 1 − 1/(1+x) − x/(1+x) + δ/(1+x) + -- − δ/(x+1)² − δ⋅x/(x+1)² + δ²/(x+1)² + -- = 1 − (1+x)/(1+x) + ((x+1) − 1)⋅δ/(x+1)² + -- − δ⋅x/(x+1)² + δ²/(x+1)² + -- = 1 − 1 + x⋅δ/(x+1)² − δ⋅x/(x+1)² + δ²/(x+1)² + -- = δ²/(x+1)² + -- + -- ε·(x+1)⋅(x+1)² − ε·δ⋅(x+1)² = δ² + -- 0 = δ² + ε·(x+1)²·δ − ε·(x+1)³ + -- + -- δ = let μ = ε·(x+1)²/2 -- Exact form + -- in -μ + √(μ² + ε·(x+1)³) -- (not overflow save) + -- + -- Safe approximation for large x: + -- ε = 1/(1+x−δ) − 1/(1+x) − δ · 1/(x+1)² + -- ≤ 1/(1+x−δ) − 1/(1+x) + -- + -- ε⋅(1+x−δ)⋅(1+x) ≤ 1+x − (1+x−δ) = δ + -- + -- δ ≥ ε + ε⋅x − ε⋅δ + ε⋅x + ε⋅x² − ε⋅δ⋅x + -- + -- δ⋅(1 + ε + ε⋅x) ≥ ε + ε⋅x + ε⋅x + ε⋅x² ≥ ε⋅x² + -- + -- δ ≥ ε⋅x²/(1 + ε + ε⋅x) + -- = ε⋅x/(1/x + ε/x + ε) + where δ ε | x<100 = let μ = ε*xp1²/2 + in sqrt(μ^2 + ε * xp1² * xp1) - μ + | otherwise = ε * x / ((1+ε)/x + ε) + xp1 = (x+1) + xp1² = xp1 ^ 2 +negativePreRegion = PreRegion $ ppr . ngt + where PreRegion ppr = positivePreRegion + ngt = actuallyLinear $ linear negate + +preRegionToInfFrom, preRegionFromMinInfTo :: RealDimension s => s -> PreRegion s s +preRegionToInfFrom xs = PreRegion $ ppr . trl + where PreRegion ppr = positivePreRegion + trl = actuallyAffine (-xs) idL +preRegionFromMinInfTo xe = PreRegion $ ppr . flp + where PreRegion ppr = positivePreRegion + flp = actuallyAffine (-xe) (linear negate) + +intervalPreRegion :: RealDimension s => (s,s) -> PreRegion s s +intervalPreRegion (lb,rb) = PreRegion $ Differentiable prr + where m = lb + radius; radius = (rb - lb)/2 + prr x = ( 1 - ((x-m)/radius)^2 + , (2*(m-x)/radius^2) *^ idL + , unsafe_dev_ε_δ("intervalPreRegion@"++show x) $ (*radius) . sqrt ) + + + + +-- | Category of functions that almost everywhere have an open region in +-- which they are continuously differentiable, i.e. /PieceWiseDiff'able/. +newtype PWDiffable s d c + = PWDiffable { + getDfblDomain :: d -> (PreRegion s d, Differentiable s d c) } + + + +instance (RealDimension s) => Category (PWDiffable s) where + type Object (PWDiffable s) o = LocallyScalable s o + id = PWDiffable $ \x -> (GlobalRegion, id) + PWDiffable f . PWDiffable g = PWDiffable h + where h x₀ = case g x₀ of + (GlobalRegion, gr) + -> let (y₀,_,_) = runDifferentiable gr x₀ + in case f y₀ of + (GlobalRegion, fr) -> (GlobalRegion, fr . gr) + (PreRegion ry, fr) + -> ( PreRegion $ ry . gr, fr . gr ) + (PreRegion rx, gr) + -> let (y₀,_,_) = runDifferentiable gr x₀ + in case f y₀ of + (GlobalRegion, fr) -> (PreRegion rx, fr . gr) + (PreRegion ry, fr) + -> ( PreRegion $ minDblfuncs (ry . gr) rx + , fr . gr ) + where (rx, gr) = g x₀ + +globalDiffable :: Differentiable s a b -> PWDiffable s a b +globalDiffable f = PWDiffable $ const (GlobalRegion, f) + +instance (RealDimension s) => EnhancedCat (PWDiffable s) (Differentiable s) where + arr = globalDiffable +instance (RealDimension s) => EnhancedCat (->) (PWDiffable s) where + arr (PWDiffable g) x = let (_,Differentiable f) = g x + (y,_,_) = f x + in y + + +instance (RealDimension s) => Cartesian (PWDiffable s) where + type UnitObject (PWDiffable s) = ZeroDim s + swap = globalDiffable swap + attachUnit = globalDiffable attachUnit + detachUnit = globalDiffable detachUnit + regroup = globalDiffable regroup + regroup' = globalDiffable regroup' + +instance (RealDimension s) => Morphism (PWDiffable s) where + PWDiffable f *** PWDiffable g = PWDiffable h + where h (x,y) = (preRegionProd rfx rgy, dff *** dfg) + where (rfx, dff) = f x + (rgy, dfg) = g y + +instance (RealDimension s) => PreArrow (PWDiffable s) where + PWDiffable f &&& PWDiffable g = PWDiffable h + where h x = (unsafePreRegionIntersect rfx rgx, dff &&& dfg) + where (rfx, dff) = f x + (rgx, dfg) = g x + terminal = globalDiffable terminal + fst = globalDiffable fst + snd = globalDiffable snd + + +instance (RealDimension s) => WellPointed (PWDiffable s) where + unit = Tagged Origin + globalElement x = PWDiffable $ \Origin -> (GlobalRegion, globalElement x) + const x = PWDiffable $ \_ -> (GlobalRegion, const x) + + +type PWDfblFuncValue s = GenericAgent (PWDiffable s) + +instance RealDimension s => HasAgent (PWDiffable s) where + alg = genericAlg + ($~) = genericAgentMap +instance RealDimension s => CartesianAgent (PWDiffable s) where + alg1to2 = genericAlg1to2 + alg2to1 = genericAlg2to1 + alg2to2 = genericAlg2to2 +instance (RealDimension s) + => PointAgent (PWDfblFuncValue s) (PWDiffable s) a x where + point = genericPoint + +gpwDfblFnValsFunc + :: ( RealDimension s + , LocallyScalable s c, LocallyScalable s c', LocallyScalable s d + , v ~ Needle c, v' ~ Needle c' + , ε ~ HerMetric v, ε ~ HerMetric v' ) + => (c' -> (c, v':-*v, ε->ε)) -> PWDfblFuncValue s d c' -> PWDfblFuncValue s d c +gpwDfblFnValsFunc f = (PWDiffable (\_ -> (GlobalRegion, Differentiable f)) $~) + +gpwDfblFnValsCombine :: forall d c c' c'' v v' v'' ε ε' ε'' s. + ( LocallyScalable s c, LocallyScalable s c', LocallyScalable s c'' + , LocallyScalable s d, RealDimension s + , v ~ Needle c, v' ~ Needle c', v'' ~ Needle c'' + , ε ~ HerMetric v , ε' ~ HerMetric v' , ε'' ~ HerMetric v'', ε~ε', ε~ε'' ) + => ( c' -> c'' -> (c, (v',v''):-*v, ε -> (ε',ε'')) ) + -> PWDfblFuncValue s d c' -> PWDfblFuncValue s d c'' -> PWDfblFuncValue s d c +gpwDfblFnValsCombine cmb (GenericAgent (PWDiffable fpcs)) + (GenericAgent (PWDiffable gpcs)) + = GenericAgent . PWDiffable $ + \d₀ -> let (rc', Differentiable f) = fpcs d₀ + (rc'',Differentiable g) = gpcs d₀ + in (unsafePreRegionIntersect rc' rc'',) . Differentiable $ + \d -> let (c', f', devf) = f d + (c'',g', devg) = g d + (c, h', devh) = cmb c' c'' + h'l = h' *.* lcofst; h'r = h' *.* lcosnd + in ( c + , h' *.* linear (lapply f' &&& lapply g') + , \εc -> let εc' = transformMetric h'l εc + εc'' = transformMetric h'r εc + (δc',δc'') = devh εc + in devf εc' ^+^ devg εc'' + ^+^ transformMetric f' δc' + ^+^ transformMetric g' δc'' + ) + where lcofst = linear(,zeroV) + lcosnd = linear(zeroV,) + + +instance (WithField s LinearManifold v, LocallyScalable s a, RealDimension s) + => AdditiveGroup (PWDfblFuncValue s a v) where + zeroV = point zeroV + (^+^) = gpwDfblFnValsCombine $ \a b -> (a^+^b, lPlus, const zeroV) + where lPlus = linear $ uncurry (^+^) + negateV = gpwDfblFnValsFunc $ \a -> (negateV a, lNegate, const zeroV) + where lNegate = linear negateV + +instance (RealDimension n, LocallyScalable n a) + => Num (PWDfblFuncValue n a n) where + fromInteger i = point $ fromInteger i + (+) = gpwDfblFnValsCombine $ \a b -> (a+b, lPlus, const zeroV) + where lPlus = linear $ uncurry (+) + (*) = gpwDfblFnValsCombine $ + \a b -> ( a*b + , linear $ \(da,db) -> a*db + b*da + , \d -> let d¹₂ = sqrt d in (d¹₂,d¹₂) + ) + negate = gpwDfblFnValsFunc $ \a -> (negate a, lNegate, const zeroV) + where lNegate = linear negate + abs = (PWDiffable absPW $~) + where absPW a₀ + | a₀<0 = (negativePreRegion, desc) + | otherwise = (positivePreRegion, asc) + desc = actuallyLinear $ linear negate + asc = actuallyLinear idL + signum = (PWDiffable sgnPW $~) + where sgnPW a₀ + | a₀<0 = (negativePreRegion, const 1) + | otherwise = (positivePreRegion, const $ -1) + +instance (RealDimension n, LocallyScalable n a) + => Fractional (PWDfblFuncValue n a n) where + fromRational i = point $ fromRational i + recip = postEndo . PWDiffable $ \a₀ -> if a₀<0 + then (negativePreRegion, Differentiable negp) + else (positivePreRegion, Differentiable posp) + where negp x = (x'¹, (- x'¹^2) *^ idL, unsafe_dev_ε_δ("1/"++show x) δ) + -- ε = 1/x − δ/x² − 1/(x+δ) + -- ε·x + ε·δ = 1 + δ/x − δ/x − δ²/x² − 1 + -- = -δ²/x² + -- 0 = δ² + ε·x²·δ + ε·x³ + -- δ = let mph = -ε·x²/2 in mph + sqrt (mph² − ε·x³) + where δ ε = let mph = -ε*x^2/2 in mph + sqrt (mph^2 - ε*x^3) + x'¹ = recip x + posp x = (x'¹, (- x'¹^2) *^ idL, unsafe_dev_ε_δ("1/"++show x) δ) + where δ ε = let mph = -ε*x^2/2 in mph + sqrt (mph^2 + ε*x^3) + x'¹ = recip x + + + + + + +-- | Category of functions that, where defined, have an open region in +-- which they are continuously differentiable. Hence /RegionWiseDiff'able/. +-- Basically these are the partial version of `PWDiffable`. +-- +-- Though the possibility of undefined regions is of course not too nice +-- (we don't need Java to demonstrate this with its everywhere-looming @null@ values...), +-- this category will propably be the “workhorse” for most serious +-- calculus applications, because it contains all the usual trig etc. functions +-- and of course everything algebraic you can do in the reals. +-- +-- The easiest way to define ordinary functions in this category is hence +-- with its 'AgentVal'ues, which have instances of the standard classes 'Num' +-- through 'Floating'. For instance, the following defines the /binary entropy/ +-- as a differentiable function on the interval @]0,1[@: (it will +-- actually /know/ where it's defined and where not! – and I don't mean you +-- need to exhaustively 'isNaN'-check all results...) +-- +-- @ +-- hb :: RWDiffable ℝ ℝ ℝ +-- hb = alg (\\p -> - p * logBase 2 p - (1-p) * logBase 2 (1-p) ) +-- @ +newtype RWDiffable s d c + = RWDiffable { + tryDfblDomain :: d -> (PreRegion s d, Option (Differentiable s d c)) } + +notDefinedHere :: Option (Differentiable s d c) +notDefinedHere = Option Nothing + + + +instance (RealDimension s) => Category (RWDiffable s) where + type Object (RWDiffable s) o = LocallyScalable s o + id = RWDiffable $ \x -> (GlobalRegion, pure id) + RWDiffable f . RWDiffable g = RWDiffable h + where h x₀ = case g x₀ of + (GlobalRegion, Option Nothing) + -> (GlobalRegion, notDefinedHere) + (GlobalRegion, Option (Just gr)) + -> let (y₀,_,_) = runDifferentiable gr x₀ + in case f y₀ of + (GlobalRegion, Option Nothing) + -> (GlobalRegion, notDefinedHere) + (GlobalRegion, Option (Just fr)) + -> (GlobalRegion, pure (fr . gr)) + (PreRegion ry, Option Nothing) + -> ( PreRegion $ ry . gr, notDefinedHere ) + (PreRegion ry, Option (Just fr)) + -> ( PreRegion $ ry . gr, pure (fr . gr) ) + (PreRegion rx, Option Nothing) + -> (PreRegion rx, notDefinedHere) + (PreRegion rx, Option (Just gr)) + -> let (y₀,_,_) = runDifferentiable gr x₀ + in case f y₀ of + (GlobalRegion, Option Nothing) + -> (PreRegion rx, notDefinedHere) + (GlobalRegion, Option (Just fr)) + -> (PreRegion rx, pure (fr . gr)) + (PreRegion ry, Option Nothing) + -> ( PreRegion $ minDblfuncs (ry . gr) rx + , notDefinedHere ) + (PreRegion ry, Option (Just fr)) + -> ( PreRegion $ minDblfuncs (ry . gr) rx + , pure (fr . gr) ) + + +globalDiffable' :: Differentiable s a b -> RWDiffable s a b +globalDiffable' f = RWDiffable $ const (GlobalRegion, pure f) + +pwDiffable :: PWDiffable s a b -> RWDiffable s a b +pwDiffable (PWDiffable q) = RWDiffable $ \x₀ -> let (r₀,f₀) = q x₀ in (r₀, pure f₀) + + + +instance (RealDimension s) => EnhancedCat (RWDiffable s) (Differentiable s) where + arr = globalDiffable' +instance (RealDimension s) => EnhancedCat (RWDiffable s) (PWDiffable s) where + arr = pwDiffable + +instance (RealDimension s) => Cartesian (RWDiffable s) where + type UnitObject (RWDiffable s) = ZeroDim s + swap = globalDiffable' swap + attachUnit = globalDiffable' attachUnit + detachUnit = globalDiffable' detachUnit + regroup = globalDiffable' regroup + regroup' = globalDiffable' regroup' + +instance (RealDimension s) => Morphism (RWDiffable s) where + RWDiffable f *** RWDiffable g = RWDiffable h + where h (x,y) = (preRegionProd rfx rgy, liftA2 (***) dff dfg) + where (rfx, dff) = f x + (rgy, dfg) = g y + +instance (RealDimension s) => PreArrow (RWDiffable s) where + RWDiffable f &&& RWDiffable g = RWDiffable h + where h x = (unsafePreRegionIntersect rfx rgx, liftA2 (&&&) dff dfg) + where (rfx, dff) = f x + (rgx, dfg) = g x + terminal = globalDiffable' terminal + fst = globalDiffable' fst + snd = globalDiffable' snd + + +instance (RealDimension s) => WellPointed (RWDiffable s) where + unit = Tagged Origin + globalElement x = RWDiffable $ \Origin -> (GlobalRegion, pure (globalElement x)) + const x = RWDiffable $ \_ -> (GlobalRegion, pure (const x)) + + +data RWDfblFuncValue s d c where + ConstRWDFV :: c -> RWDfblFuncValue s d c + GenericRWDFV :: RWDiffable s d c -> RWDfblFuncValue s d c + +genericiseRWDFV :: (RealDimension s, LocallyScalable s c, LocallyScalable s d) + => RWDfblFuncValue s d c -> RWDfblFuncValue s d c +genericiseRWDFV (ConstRWDFV c) = GenericRWDFV $ const c +genericiseRWDFV v = v + +instance RealDimension s => HasAgent (RWDiffable s) where + type AgentVal (RWDiffable s) d c = RWDfblFuncValue s d c + alg fq = case fq (GenericRWDFV id) of + GenericRWDFV f -> f + ($~) = postCompRW +instance RealDimension s => CartesianAgent (RWDiffable s) where + alg1to2 fgq = case fgq (GenericRWDFV id) of + (GenericRWDFV f, GenericRWDFV g) -> f &&& g + alg2to1 fq = case fq (GenericRWDFV fst) (GenericRWDFV snd) of + GenericRWDFV f -> f + alg2to2 fgq = case fgq (GenericRWDFV fst) (GenericRWDFV snd) of + (GenericRWDFV f, GenericRWDFV g) -> f &&& g +instance (RealDimension s) + => PointAgent (RWDfblFuncValue s) (RWDiffable s) a x where + point = ConstRWDFV + +grwDfblFnValsFunc + :: ( RealDimension s + , LocallyScalable s c, LocallyScalable s c', LocallyScalable s d + , v ~ Needle c, v' ~ Needle c' + , ε ~ HerMetric v, ε ~ HerMetric v' ) + => (c' -> (c, v':-*v, ε->ε)) -> RWDfblFuncValue s d c' -> RWDfblFuncValue s d c +grwDfblFnValsFunc f = (RWDiffable (\_ -> (GlobalRegion, pure (Differentiable f))) $~) + +grwDfblFnValsCombine :: forall d c c' c'' v v' v'' ε ε' ε'' s. + ( LocallyScalable s c, LocallyScalable s c', LocallyScalable s c'' + , LocallyScalable s d, RealDimension s + , v ~ Needle c, v' ~ Needle c', v'' ~ Needle c'' + , ε ~ HerMetric v , ε' ~ HerMetric v' , ε'' ~ HerMetric v'', ε~ε', ε~ε'' ) + => ( c' -> c'' -> (c, (v',v''):-*v, ε -> (ε',ε'')) ) + -> RWDfblFuncValue s d c' -> RWDfblFuncValue s d c'' -> RWDfblFuncValue s d c +grwDfblFnValsCombine cmb (GenericRWDFV (RWDiffable fpcs)) + (GenericRWDFV (RWDiffable gpcs)) + = GenericRWDFV . RWDiffable $ + \d₀ -> let (rc', fmay) = fpcs d₀ + (rc'',gmay) = gpcs d₀ + in (unsafePreRegionIntersect rc' rc'',) $ + case (fmay,gmay) of + (Option(Just(Differentiable f)), Option(Just(Differentiable g))) -> + pure . Differentiable $ \d + -> let (c', f', devf) = f d + (c'',g', devg) = g d + (c, h', devh) = cmb c' c'' + h'l = h' *.* lcofst; h'r = h' *.* lcosnd + in ( c + , h' *.* linear (lapply f' &&& lapply g') + , \εc -> let εc' = transformMetric h'l εc + εc'' = transformMetric h'r εc + (δc',δc'') = devh εc + in devf εc' ^+^ devg εc'' + ^+^ transformMetric f' δc' + ^+^ transformMetric g' δc'' + ) + _ -> notDefinedHere + where lcofst = linear(,zeroV) + lcosnd = linear(zeroV,) +grwDfblFnValsCombine cmb fv gv + = grwDfblFnValsCombine cmb (genericiseRWDFV fv) (genericiseRWDFV gv) + + +postCompRW :: ( RealDimension s + , LocallyScalable s a, LocallyScalable s b, LocallyScalable s c ) + => RWDiffable s b c -> RWDfblFuncValue s a b -> RWDfblFuncValue s a c +postCompRW (RWDiffable f) (ConstRWDFV x) = case f x of + (_, Option (Just fd)) -> ConstRWDFV $ fd $ x +postCompRW f (GenericRWDFV g) = GenericRWDFV $ f . g + + +instance ( WithField s EuclidSpace v, AdditiveGroup v, v ~ Needle (Interior (Needle v)) + , LocallyScalable s a, RealDimension s) + => AdditiveGroup (RWDfblFuncValue s a v) where + zeroV = point zeroV + ConstRWDFV c₁ ^+^ ConstRWDFV c₂ = ConstRWDFV (c₁^+^c₂) + ConstRWDFV c₁ ^+^ GenericRWDFV g = GenericRWDFV $ + globalDiffable' (actuallyAffine c₁ idL) . g + GenericRWDFV f ^+^ ConstRWDFV c₂ = GenericRWDFV $ + globalDiffable' (actuallyAffine c₂ idL) . f + v^+^w = grwDfblFnValsCombine (\a b -> (a^+^b, lPlus, const zeroV)) v w + where lPlus = linear $ uncurry (^+^) + negateV (ConstRWDFV c) = ConstRWDFV (negateV c) + negateV v = grwDfblFnValsFunc (\a -> (negateV a, lNegate, const zeroV)) v + where lNegate = linear negateV + +instance (RealDimension n, LocallyScalable n a) + => Num (RWDfblFuncValue n a n) where + fromInteger i = point $ fromInteger i + (+) = (^+^) + ConstRWDFV c₁ * ConstRWDFV c₂ = ConstRWDFV (c₁*c₂) + ConstRWDFV c₁ * GenericRWDFV g = GenericRWDFV $ + globalDiffable' (actuallyLinear $ linear (c₁*)) . g + GenericRWDFV f * ConstRWDFV c₂ = GenericRWDFV $ + globalDiffable' (actuallyLinear $ linear (*c₂)) . f + v*w = grwDfblFnValsCombine ( + \a b -> ( a*b + , linear $ \(da,db) -> a*db + b*da + , \d -> let d¹₂ = sqrt d in (d¹₂,d¹₂) + ) + ) v w + negate = negateV + abs = (RWDiffable absPW $~) + where absPW a₀ + | a₀<0 = (negativePreRegion, pure desc) + | otherwise = (positivePreRegion, pure asc) + desc = actuallyLinear $ linear negate + asc = actuallyLinear idL + signum = (RWDiffable sgnPW $~) + where sgnPW a₀ + | a₀<0 = (negativePreRegion, pure (const 1)) + | otherwise = (positivePreRegion, pure (const $ -1)) + +instance (RealDimension n, LocallyScalable n a) + => Fractional (RWDfblFuncValue n a n) where + fromRational i = point $ fromRational i + recip = postCompRW . RWDiffable $ \a₀ -> if a₀<0 + then (negativePreRegion, pure (Differentiable negp)) + else (positivePreRegion, pure (Differentiable posp)) + where negp x = (x'¹, (- x'¹^2) *^ idL, unsafe_dev_ε_δ("1/"++show x) δ) + -- ε = 1/x − δ/x² − 1/(x+δ) + -- ε·x + ε·δ = 1 + δ/x − δ/x − δ²/x² − 1 + -- = -δ²/x² + -- 0 = δ² + ε·x²·δ + ε·x³ + -- δ = let mph = -ε·x²/2 in mph + sqrt (mph² − ε·x³) + where δ ε = let mph = -ε*x^2/2 in mph + sqrt (mph^2 - ε*x^3) + x'¹ = recip x + posp x = (x'¹, (- x'¹^2) *^ idL, unsafe_dev_ε_δ("1/"++show x) δ) + where δ ε = let mph = -ε*x^2/2 in mph + sqrt (mph^2 + ε*x^3) + x'¹ = recip x + + + + + +-- Helper for checking ε-estimations in GHCi with dynamic-plot: +-- epsEst (f,f') εsgn δf (ViewXCenter xc) (ViewHeight h) +-- = let δfxc = δf xc +-- in tracePlot $ reverse [ (xc - δ, f xc - δ * f' xc + εsgn*ε) | +-- ε <- [0, h/500 .. h], let δ = δfxc ε] +-- ++ [ (xc + δ, f xc + δ * f' xc + εsgn*ε) | +-- ε <- [0, h/500 .. h], let δ = δfxc ε] +-- Golfed version: +-- epsEst(f,d)s φ(ViewXCenter ξ)(ViewHeight h)=let ζ=φ ξ in tracePlot$[(ξ-δ,f ξ-δ*d ξ+s*abs ε)|ε<-[-h,-0.998*h..h],let δ=ζ(abs ε)*signum ε] + +instance (RealDimension n, LocallyScalable n a) + => Floating (RWDfblFuncValue n a n) where + pi = point pi + + exp = grwDfblFnValsFunc + $ \x -> let ex = exp x + in if ex==0 -- numeric underflow + then ( 0, zeroV, unsafe_dev_ε_δ("exp "++show x) $ \ε -> log ε - x ) + else ( ex, ex *^ idL, unsafe_dev_ε_δ("exp "++show x) $ \ε -> acosh(ε/(2*ex) + 1) ) + -- ε = e^(x+δ) − eˣ − eˣ·δ + -- = eˣ·(e^δ − 1 − δ) + -- ≤ eˣ · (e^δ − 1 + e^(-δ) − 1) + -- = eˣ · 2·(cosh(δ) − 1) + -- cosh(δ) ≥ ε/(2·eˣ) + 1 + -- δ ≥ acosh(ε/(2·eˣ) + 1) + log = postCompRW . RWDiffable $ \x -> if x>0 + then (positivePreRegion, pure (Differentiable lnPosR)) + else (negativePreRegion, notDefinedHere) + where lnPosR x = ( log x, recip x *^ idL, unsafe_dev_ε_δ("log "++show x) $ \ε -> x * sqrt(1 - exp(-ε)) ) + -- ε = ln x + (-δ)/x − ln(x−δ) + -- = ln (x / ((x−δ) · exp(δ/x))) + -- x/e^ε = (x−δ) · exp(δ/x) + -- let γ = δ/x ∈ [0,1[ + -- exp(-ε) = (1−γ) · e^γ + -- ≥ (1−γ) · (1+γ) + -- = 1 − γ² + -- γ ≥ sqrt(1 − exp(-ε)) + -- δ ≥ x · sqrt(1 − exp(-ε)) + + sqrt = postCompRW . RWDiffable $ \x -> if x>0 + then (positivePreRegion, pure (Differentiable sqrtPosR)) + else (negativePreRegion, notDefinedHere) + where sqrtPosR x = ( sx, idL ^/ (2*sx), unsafe_dev_ε_δ("sqrt "++show x) $ + \ε -> 2 * (s2 * sqrt sx^3 * sqrt ε + signum (ε*2-sx) * sx * ε) ) + where sx = sqrt x; s2 = sqrt 2 + -- Exact inverse of O(δ²) remainder. + + sin = grwDfblFnValsFunc sinDfb + where sinDfb x = ( sx, cx *^ idL, unsafe_dev_ε_δ("sin "++show x) δ ) + where sx = sin x; cx = cos x + sx² = sx^2; cx² = cx^2 + sx' = abs sx; cx' = abs cx + sx'³ = sx'*sx²; cx⁴ = cx²*cx² + δ ε = (ε*(1.8 + ε^2/(cx' + (2+40*cx⁴)/ε)) + σ₃³*sx'³)**(1/3) - σ₃*sx' + + σ₂*sqrt ε/(σ₂+cx²) + -- Carefully fine-tuned to give everywhere a good and safe bound. + -- The third root makes it pretty slow too, but since tight + -- deviation bounds can dramatically reduce the number of evaluations + -- needed in the first place, this is probably worthwhile. + σ₂ = 1.4; σ₃ = 1.75; σ₃³ = σ₃^3 + -- Safety margins for overlap between quadratic and cubic model + -- (these aren't naturally compatible to be used both together) + + cos = sin . (globalDiffable' (actuallyAffine (pi/2) idL) $~) + + sinh x = (exp x - exp (-x))/2 + {- = grwDfblFnValsFunc sinhDfb + where sinhDfb x = ( sx, cx *^ idL, unsafe_dev_ε_δ δ ) + where sx = sinh x; cx = cosh x + δ ε = undefined -} + -- ε = sinh x + δ · cosh x − sinh(x+δ) + -- = ½ · ( eˣ − e⁻ˣ + δ · (eˣ + e⁻ˣ) − exp(x+δ) + exp(-x−δ) ) + -- = ½·e⁻ˣ · ( e²ˣ − 1 + δ · (e²ˣ + 1) − e²ˣ·e^δ + e^-δ ) + -- = ½ · ( eˣ − e⁻ˣ + δ · (eˣ + e⁻ˣ) − exp(x+δ) + exp(-x−δ) ) + cosh x = (exp x + exp (-x))/2 + + tanh = grwDfblFnValsFunc tanhDfb + where tanhDfb x = ( tnhx, idL ^/ (cosh x^2), unsafe_dev_ε_δ("tan "++show x) δ ) + where tnhx = tanh x + c = (tnhx*2/pi)^2 + p = 1 + abs x/(2*pi) + δ ε = p * (sqrt ε + ε * c) + -- copied from 'atan' definition. Empirically works safely, in fact + -- with quite a big margin. TODO: find a tighter definition. + + atan = grwDfblFnValsFunc atanDfb + where atanDfb x = ( atnx, idL ^/ (1+x^2), unsafe_dev_ε_δ("atan "++show x) δ ) + where atnx = atan x + c = (atnx*2/pi)^2 + p = 1 + abs x/(2*pi) + δ ε = p * (sqrt ε + ε * c) + -- Semi-empirically obtained: with the epsEst helper, + -- it is observed that this function is (for xc≥0) a lower bound + -- to the arctangent. The growth of the p coefficient makes sense + -- and holds for arbitrarily large xc, because those move us linearly + -- away from the only place where the function is not virtually constant + -- (around 0). + + asin = postCompRW . RWDiffable $ \x -> if + | x < (-1) -> (preRegionFromMinInfTo (-1), notDefinedHere) + | x > 1 -> (preRegionToInfFrom 1, notDefinedHere) + | otherwise -> (intervalPreRegion (-1,1), pure (Differentiable asinDefdR)) + where asinDefdR x = ( asinx, asin'x *^ idL, unsafe_dev_ε_δ("asin "++show x) δ ) + where asinx = asin x; asin'x = recip (sqrt $ 1 - x^2) + c = 1 - x^2 + δ ε = sqrt ε * c + -- Empirical, with epsEst upper bound. + + acos = postCompRW . RWDiffable $ \x -> if + | x < (-1) -> (preRegionFromMinInfTo (-1), notDefinedHere) + | x > 1 -> (preRegionToInfFrom 1, notDefinedHere) + | otherwise -> (intervalPreRegion (-1,1), pure (Differentiable acosDefdR)) + where acosDefdR x = ( acosx, acos'x *^ idL, unsafe_dev_ε_δ("acos "++show x) δ ) + where acosx = acos x; acos'x = - recip (sqrt $ 1 - x^2) + c = 1 - x^2 + δ ε = sqrt ε * c -- Like for asin – it's just a translation/reflection. + + asinh = grwDfblFnValsFunc asinhDfb + where asinhDfb x = ( asinhx, idL ^/ sqrt(1+x^2), unsafe_dev_ε_δ("asinh "++show x) δ ) + where asinhx = asinh x + δ ε = abs x * sqrt((1 - exp(-ε))*0.8 + ε^2/(3*abs x)) + sqrt(ε/(abs x+0.5)) + -- Empirical, modified from log function (the area hyperbolic sine + -- resembles two logarithmic lobes), with epsEst-checked lower bound. + + acosh = postCompRW . RWDiffable $ \x -> if x>0 + then (positivePreRegion, pure (Differentiable acoshDfb)) + else (negativePreRegion, notDefinedHere) + where acoshDfb x = ( acosh x, idL ^/ sqrt(x^2 - 2), unsafe_dev_ε_δ("acosh "++show x) δ ) + where δ ε = (2 - 1/sqrt x) * (s2 * sqrt sx^3 * sqrt(ε/s2) + signum (ε*s2-sx) * sx * ε/s2) + sx = sqrt(x-1) + s2 = sqrt 2 + -- Empirical, modified from sqrt function – the area hyperbolic cosine + -- strongly resembles \x -> sqrt(2 · (x-1)). + + atanh = postCompRW . RWDiffable $ \x -> if + | x < (-1) -> (preRegionFromMinInfTo (-1), notDefinedHere) + | x > 1 -> (preRegionToInfFrom 1, notDefinedHere) + | otherwise -> (intervalPreRegion (-1,1), pure (Differentiable atnhDefdR)) + where atnhDefdR x = ( atanh x, recip(1-x^2) *^ idL, unsafe_dev_ε_δ("atanh "++show x) $ \ε -> sqrt(tanh ε)*(1-abs x) ) + -- Empirical, with epsEst upper bound. + + + + + +isZeroMap :: ∀ v a . (FiniteDimensional v, AdditiveGroup a, Eq a) => (v:-*a) -> Bool +isZeroMap m = all ((==zeroV) . atBasis m) b + where (Tagged b) = completeBasis :: Tagged v [Basis v] + + + diff --git a/Data/LinearMap/HerMetric.hs b/Data/LinearMap/HerMetric.hs --- a/Data/LinearMap/HerMetric.hs +++ b/Data/LinearMap/HerMetric.hs @@ -37,6 +37,7 @@ , eigenCoSpan, eigenCoSpan' , metriScale', metriScale , adjoint + , extendMetric -- * The dual-space class , HasMetric , HasMetric'(..) @@ -121,6 +122,19 @@ metricMatrix' :: Maybe (HMat.Matrix (Scalar v)) } +extendMetric :: (HasMetric v, Scalar v~ℝ) => HerMetric v -> v -> HerMetric v +extendMetric (HerMetric Nothing) _ = HerMetric Nothing +extendMetric (HerMetric (Just m)) v + | isInfinite' detm = HerMetric $ Just m + | isInfinite' detmninv = singularMetric + | otherwise = HerMetric $ Just mn + where -- this could probably be done much more efficiently, with only + -- multiplications, no inverses. + (minv, (detm, _)) = HMat.invlndet m + (mn, (detmninv, _)) = HMat.invlndet (minv + HMat.outer vv vv) + vv = asPackedVector v + + matrixMetric' :: HasMetric v => HMat.Matrix (Scalar v) -> HerMetric' v matrixMetric' = HerMetric' . Just @@ -267,28 +281,24 @@ eigenSpan :: (HasMetric v, Scalar v ~ ℝ) => HerMetric' v -> [v] eigenSpan (HerMetric' Nothing) = [] eigenSpan (HerMetric' (Just m)) = map fromPackedVector eigSpan - where (μs,vsm) = HMat.eigSH' m -- TODO: replace with `eigSH'`, which is unchecked - -- (`HerMetric` is always Hermitian!) + where (μs,vsm) = HMat.eigSH' m eigSpan = zipWith (HMat.scale . sqrt) (HMat.toList μs) (HMat.toColumns vsm) eigenSpan' :: (HasMetric v, Scalar v ~ ℝ) => HerMetric v -> [DualSpace v] eigenSpan' (HerMetric Nothing) = [] eigenSpan' (HerMetric (Just m)) = map fromPackedVector eigSpan - where (μs,vsm) = HMat.eigSH' m -- TODO: replace with `eigSH'`, which is unchecked - -- (`HerMetric` is always Hermitian!) + where (μs,vsm) = HMat.eigSH' m eigSpan = zipWith (HMat.scale . sqrt) (HMat.toList μs) (HMat.toColumns vsm) eigenCoSpan :: (HasMetric v, Scalar v ~ ℝ) => HerMetric' v -> [DualSpace v] eigenCoSpan (HerMetric' Nothing) = [] eigenCoSpan (HerMetric' (Just m)) = map fromPackedVector eigSpan - where (μs,vsm) = HMat.eigSH' m -- TODO: replace with `eigSH'`, which is unchecked - -- (`HerMetric` is always Hermitian!) + where (μs,vsm) = HMat.eigSH' m eigSpan = zipWith (HMat.scale . recip . sqrt) (HMat.toList μs) (HMat.toColumns vsm) eigenCoSpan' :: (HasMetric v, Scalar v ~ ℝ) => HerMetric v -> [v] eigenCoSpan' (HerMetric Nothing) = [] eigenCoSpan' (HerMetric (Just m)) = map fromPackedVector eigSpan - where (μs,vsm) = HMat.eigSH' m -- TODO: replace with `eigSH'`, which is unchecked - -- (`HerMetric` is always Hermitian!) + where (μs,vsm) = HMat.eigSH' m eigSpan = zipWith (HMat.scale . recip . sqrt) (HMat.toList μs) (HMat.toColumns vsm) @@ -459,7 +469,7 @@ normaliseWith :: HasMetric v => HerMetric v -> v -> Option v normaliseWith m v = case metric m v of - 0 -> Hask.empty + 0 -> empty μ -> pure (v ^/ μ) orthonormalPairsWith :: forall v . HasMetric v => HerMetric v -> [v] -> [(v, DualSpace v)] @@ -513,7 +523,10 @@ where (Tagged dw) = dimension :: Tagged w Int metricAsLength :: HerMetric ℝ -> ℝ -metricAsLength = recip . (`metric`1) +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." metricFromLength :: ℝ -> HerMetric ℝ metricFromLength = projector . recip @@ -535,7 +548,7 @@ spanHilbertSubspace met = emb . orthonormalPairsWith met where emb onb' | n'==n = return $ Embedding emb prj . arr identityMatrix - | otherwise = Hask.empty + | otherwise = empty where emb = DenseLinear . HMat.fromColumns $ (asPackedVector . fst) <$> onb prj = DenseLinear . HMat.fromRows $ (asPackedVector . snd) <$> onb n' = length onb' @@ -560,3 +573,25 @@ -- therefore we define this space not as normalised vectors, but rather as all -- vectors modulo scaling by positive factors. newtype Stiefel1 v = Stiefel1 { getStiefel1N :: DualSpace v } + + + + + + + +instance (HasMetric v, Scalar v ~ Double, Show (DualSpace v)) => Show (HerMetric v) where + showsPrec p m + | null eigSp = showString "zeroV" + | otherwise = showParen (p>5) + . foldr1 ((.) . (.(" ^+^ "++))) + $ ((("projector "++).).showsPrec 6)<$>eigSp + where eigSp = eigenSpan' m + +instance (HasMetric v, Scalar v ~ Double, Show v) => Show (HerMetric' v) where + showsPrec p m + | null eigSp = showString "zeroV" + | otherwise = showParen (p>5) + . foldr1 ((.) . (.(" ^+^ "++))) + $ ((("projector' "++).).showsPrec 6)<$>eigSp + where eigSp = eigenSpan m diff --git a/Data/Manifold/Cone.hs b/Data/Manifold/Cone.hs --- a/Data/Manifold/Cone.hs +++ b/Data/Manifold/Cone.hs @@ -127,13 +127,13 @@ type CℝayInterior (ZeroDim ℝ) = ℝ fromCℝayInterior (FinVecArrRep qb) | HMat.size qb == 0 = Cℝay 1 Origin | x <- qb HMat.! 0 = Cℝay (bijectℝtoℝplus x) Origin - toCℝayInterior (Cℝay 0 Origin) = Hask.empty + toCℝayInterior (Cℝay 0 Origin) = empty toCℝayInterior (Cℝay y Origin) = pure . FinVecArrRep $ 1 HMat.|>[bijectℝplustoℝ y] instance ConeSemimfd ℝ where type CℝayInterior ℝ = ℝ² fromCℝayInterior (FinVecArrRep qb) = Cℝay (q'+b') (q'-b') where [q', b'] = HMat.toList $ HMat.cmap ((/2) . bijectℝtoℝplus) qb - toCℝayInterior (Cℝay 0 _) = Hask.empty + toCℝayInterior (Cℝay 0 _) = empty toCℝayInterior (Cℝay h x) = pure . FinVecArrRep . HMat.cmap bijectℝplustoℝ $ HMat.fromList [h+x, h-x] fromCD¹Interior (FinVecArrRep qb) = CD¹ (bijectℝplustoIntv $ q'+b') (q'-b') @@ -152,7 +152,7 @@ fromCD¹Interior xa | x>0 = CD¹ (bijectℝtoIntv x) PositiveHalfSphere | otherwise = CD¹ (-bijectℝtoIntv x) NegativeHalfSphere where x = getFinVecArrRep xa HMat.! 0 - toCD¹Interior (CD¹ 1 _) = Hask.empty + toCD¹Interior (CD¹ 1 _) = empty toCD¹Interior (CD¹ x PositiveHalfSphere) = return . FinVecArrRep . HMat.scalar $ bijectIntvtoℝ x toCD¹Interior (CD¹ x NegativeHalfSphere) @@ -169,7 +169,7 @@ fromCD¹Interior (FinVecArrRep xy) = CD¹ (bijectℝtoIntv r) (S¹ $ atan2 y x) where r = HMat.norm_2 xy [x,y] = HMat.toList xy - toCD¹Interior (CD¹ 1 _) = Hask.empty + toCD¹Interior (CD¹ 1 _) = empty toCD¹Interior (CD¹ r (S¹ φ)) = return . FinVecArrRep . HMat.scale r' $ HMat.fromList [cos φ, sin φ] where r' = bijectIntvtoℝ r @@ -215,7 +215,7 @@ type CℝayInterior (CD¹ x) = (ℝ, ConeVecArr x) fromCℝayInterior i = Cℝay h (embCℝayToCD¹ o) where (Cℝay h o) = simplyCncted_fromCℝayInterior i - toCℝayInterior (Cℝay _ (CD¹ 1 _)) = Hask.empty + toCℝayInterior (Cℝay _ (CD¹ 1 _)) = empty toCℝayInterior (Cℝay h p) = simplyCncted_toCℝayInterior $ Cℝay h (projCD¹ToCℝay p) @@ -244,7 +244,7 @@ cmps = (h - Arr.sum cmps') `Arr.cons` cmps n = fromIntegral $ Arr.length cmps in return $ FinVecArrRep (bijectℝplustoℝ `Arr.map` cmps) -simplyCncted_toCℝayInterior (Cℝay _ _) = Hask.empty +simplyCncted_toCℝayInterior (Cℝay _ _) = empty -- Some essential homeomorphisms @@ -274,28 +274,6 @@ projCD¹ToCℝay :: CD¹ m -> Cℝay m projCD¹ToCℝay (CD¹ h m) = Cℝay (bijectIntvtoℝplus h) m --- instance (WithScalar ℝ PseudoAffine m) => Semimanifold (Cℝay m) where --- type Needle (Cℝay m) = (Needle m, ℝ) --- type Interior (Cℝay m) = (Interior m, ℝ) --- --- fromInterior (im, d) --- | d>38 = Cℝay m d -- from 38 on, the +1 is numerically --- -- insignificant against the exponential. --- | otherwise = cℝay m (log $ exp d + 1) --- -- note that (for the same reason we can shortcut above 38) --- -- such negative arguments will actually yield the value zero. --- -- This means we're actually reaching the “infinitely far” --- -- rim rather quickly. This might be a problem, but normally --- -- shouldn't really matter much. --- -- It would perhaps be better to have homeomorphism that --- -- approaches -1/x in the negative limit, but such a --- -- function doesn't seem as easy to come by. --- where m = fromInterior im --- toInterior (Cℝay m q) --- | q>38 = fmap (,q) im --- | q>0 = fmap (, log $ exp d - 1) im --- | otherwise = Hask.empty --- where im = toInterior m stiefel1Project :: LinearManifold v => DualSpace v -- ^ Must be nonzero. @@ -322,10 +300,6 @@ instance HasUnitSphere ℝ³ where type UnitSphere ℝ³ = S² instance HasUnitSphere (FinVecArrRep t ℝ³ ℝ) where type UnitSphere (FinVecArrRep t ℝ³ ℝ) = S² - --- instance (HasUnitSphere v, v ~ DualSpace v) => NaturallyEmbedded (Stiefel1 v) v where --- embed = embed . unstiefel --- coEmbed = stiefel . coEmbed diff --git a/Data/Manifold/PseudoAffine.hs b/Data/Manifold/PseudoAffine.hs --- a/Data/Manifold/PseudoAffine.hs +++ b/Data/Manifold/PseudoAffine.hs @@ -50,16 +50,6 @@ Manifold , Semimanifold(..) , PseudoAffine(..) - -- * Regions within a manifold - , Region - , smoothIndicator - -- * Hierarchy of manifold-categories - -- ** Everywhere differentiable functions - , Differentiable - -- ** Almost everywhere diff'able funcs - , PWDiffable - -- ** Region-wise defined diff'able funcs - , RWDiffable -- * Type definitions -- ** Metrics , Metric, Metric', euclideanMetric @@ -70,13 +60,9 @@ , WithField , HilbertSpace , EuclidSpace + , LocallyScalable -- * Misc , palerp - , discretisePathIn - , discretisePathSegs - , continuousIntervals - , regionOfContinuityAround - , analyseLocalBehaviour ) where @@ -253,7 +239,7 @@ -- | The 'RealFloat' class plus manifold constraints. type RealDimension r = ( PseudoAffine r, Interior r ~ r, Needle r ~ r , HasMetric r, DualSpace r ~ r, Scalar r ~ r - , RealFloat r ) + , RealFloat r, r ~ ℝ) -- | The 'AffineSpace' class plus manifold constraints. type AffineManifold m = ( PseudoAffine m, Interior m ~ m, AffineSpace m @@ -300,106 +286,11 @@ => 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 - _ -> Hask.empty - + _ -> empty -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 /ε/. -discretisePathIn nLim (Region xm rLim) m (Differentiable f) - = reverse (tail . take nLim $ traceFwd xm (-1)) ++ take nLim (traceFwd xm 1) - where traceFwd x₀ 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 @@ -517,14 +408,14 @@ type Interior D¹ = ℝ fromInterior = D¹ . tanh toInterior (D¹ x) | abs x < 1 = return $ atanh x - | otherwise = Hask.empty + | otherwise = empty translateP = Tagged (+) instance PseudoAffine D¹ where - D¹ 1 .-~. _ = Hask.empty - D¹ (-1) .-~. _ = Hask.empty + D¹ 1 .-~. _ = empty + D¹ (-1) .-~. _ = empty D¹ x .-~. y | abs x < 1 = return $ atanh x - y - | otherwise = Hask.empty + | otherwise = empty instance Semimanifold S² where type Needle S² = ℝ² @@ -586,861 +477,11 @@ - tau :: ℝ tau = 2 * pi toS¹range :: ℝ -> ℝ toS¹range φ = (φ+pi)`mod'`tau - pi - - - - -type LinDevPropag d c = Metric c -> Metric d - -dev_ε_δ :: RealDimension a - => (a -> a) -> LinDevPropag a a -dev_ε_δ f d = let ε = 1 / metric d 1 in projector $ 1 / f ε - --- | The category of differentiable functions between manifolds over scalar @s@. --- --- As you might guess, these offer /automatic differentiation/ of sorts (basically, --- simple forward AD), but that's in itself is not really the killer feature here. --- More interestingly, we actually have the (à la Curry-Howard) /proof/ --- built in: the function /f/ has at /x/₀ derivative /f'ₓ/₀, --- if, for¹ /ε/>0, there exists /δ/ such that --- |/f/ /x/ − (/f/ /x/₀ + /x/⋅/f'ₓ/₀)| < /ε/ --- for all |/x/ − /x/₀| < /δ/. --- --- Observe that, though this looks quite similar to the standard definition --- of differentiability, it is not equivalent thereto – in fact it does --- not prove any analytic properties at all. To make it equivalent, we need --- a lower bound on /δ/: simply /δ/ gives us continuity, and for --- continuous differentiability, /δ/ must grow at least like √/ε/ --- for small /ε/. Neither of these conditions are enforced by the type system, --- but we do require them for any allowed values because these proofs are obviously --- tremendously useful – for instance, you can have a root-finding algorithm --- and actually be sure you get /all/ solutions correctly, not just /some/ that are --- (hopefully) the closest to some reference point you'd need to laborously define! --- --- Unfortunately however, this also prevents doing any serious algebra etc. with the --- category, because even something as simple as division necessary introduces singularities --- where the derivatives must diverge. --- Not to speak of many trigonometric e.g. trigonometric functions that --- are undefined on whole regions. The 'PWDiffable' and 'RWDiffable' categories have explicit --- handling for those issues built in; you may simply use these categories even when --- you know the result will be smooth in your relevant domain (or must be, for e.g. --- physics reasons). --- --- ¹(The implementation does not deal with /ε/ and /δ/ as difference-bounding --- reals, but rather as metric tensors that define a boundary by prohibiting the --- overlap from exceeding one; this makes the concept actually work on general manifolds.) -newtype Differentiable s d c - = Differentiable { runDifferentiable :: - d -> ( c, Needle d :-* Needle c, LinDevPropag d c ) } -type (-->) = Differentiable ℝ - - -instance (MetricScalar s) => Category (Differentiable s) where - type Object (Differentiable s) o = LocallyScalable s o - id = Differentiable $ \x -> (x, idL, const zeroV) - Differentiable f . Differentiable g = Differentiable $ - \x -> let (y, g', devg) = g x - (z, f', devf) = f y - devfg δz = let δy = transformMetric f' δz - εy = devf δz - in transformMetric g' εy ^+^ devg δy ^+^ devg εy - in (z, f'*.*g', devfg) - - -instance (RealDimension s) => EnhancedCat (->) (Differentiable s) where - arr (Differentiable f) x = let (y,_,_) = f x in y - -instance (MetricScalar s) => Cartesian (Differentiable s) where - type UnitObject (Differentiable s) = ZeroDim s - swap = Differentiable $ \(x,y) -> ((y,x), lSwap, const zeroV) - where lSwap = linear swap - attachUnit = Differentiable $ \x -> ((x, Origin), lAttachUnit, const zeroV) - where lAttachUnit = linear $ \x -> (x, Origin) - detachUnit = Differentiable $ \(x, Origin) -> (x, lDetachUnit, const zeroV) - where lDetachUnit = linear $ \(x, Origin) -> x - regroup = Differentiable $ \(x,(y,z)) -> (((x,y),z), lRegroup, const zeroV) - where lRegroup = linear regroup - regroup' = Differentiable $ \((x,y),z) -> ((x,(y,z)), lRegroup, const zeroV) - where lRegroup = linear regroup' - - -instance (MetricScalar s) => Morphism (Differentiable s) where - Differentiable f *** Differentiable g = Differentiable h - where h (x,y) = ((fx, gy), lPar, devfg) - where (fx, f', devf) = f x - (gy, g', devg) = g y - devfg δs = transformMetric lfst δx - ^+^ transformMetric lsnd δy - where δx = devf $ transformMetric lcofst δs - δy = devg $ transformMetric lcosnd δs - lPar = linear $ lapply f'***lapply g' - lfst = linear fst; lsnd = linear snd - lcofst = linear (,zeroV); lcosnd = linear (zeroV,) - - -instance (MetricScalar s) => PreArrow (Differentiable s) where - terminal = Differentiable $ \_ -> (Origin, zeroV, const zeroV) - fst = Differentiable $ \(x,_) -> (x, lfst, const zeroV) - where lfst = linear fst - snd = Differentiable $ \(_,y) -> (y, lsnd, const zeroV) - where lsnd = linear snd - Differentiable f &&& Differentiable g = Differentiable h - where h x = ((fx, gx), lFanout, devfg) - where (fx, f', devf) = f x - (gx, g', devg) = g x - devfg δs = (devf $ transformMetric lcofst δs) - ^+^ (devg $ transformMetric lcosnd δs) - lFanout = linear $ lapply f'&&&lapply g' - lcofst = linear (,zeroV); lcosnd = linear (zeroV,) - - -instance (MetricScalar s) => WellPointed (Differentiable s) where - unit = Tagged Origin - globalElement x = Differentiable $ \Origin -> (x, zeroV, const zeroV) - const x = Differentiable $ \_ -> (x, zeroV, const zeroV) - - - -type DfblFuncValue s = GenericAgent (Differentiable s) - -instance (MetricScalar s) => HasAgent (Differentiable s) where - alg = genericAlg - ($~) = genericAgentMap -instance (MetricScalar s) => CartesianAgent (Differentiable s) where - alg1to2 = genericAlg1to2 - alg2to1 = genericAlg2to1 - alg2to2 = genericAlg2to2 -instance (MetricScalar s) - => PointAgent (DfblFuncValue s) (Differentiable s) a x where - point = genericPoint - - - -actuallyLinear :: ( WithField s LinearManifold x, WithField s LinearManifold y ) - => (x:-*y) -> Differentiable s x y -actuallyLinear f = Differentiable $ \x -> (lapply f x, f, const zeroV) - -actuallyAffine :: ( WithField s LinearManifold x, WithField s LinearManifold y ) - => y -> (x:-*y) -> Differentiable s x y -actuallyAffine y₀ f = Differentiable $ \x -> (y₀ ^+^ lapply f x, f, const zeroV) - - -dfblFnValsFunc :: ( LocallyScalable s c, LocallyScalable s c', LocallyScalable s d - , v ~ Needle c, v' ~ Needle c' - , ε ~ HerMetric v, ε ~ HerMetric v' ) - => (c' -> (c, v':-*v, ε->ε)) -> DfblFuncValue s d c' -> DfblFuncValue s d c -dfblFnValsFunc f = (Differentiable f $~) - -dfblFnValsCombine :: forall d c c' c'' v v' v'' ε ε' ε'' s. - ( LocallyScalable s c, LocallyScalable s c', LocallyScalable s c'' - , LocallyScalable s d - , v ~ Needle c, v' ~ Needle c', v'' ~ Needle c'' - , ε ~ HerMetric v , ε' ~ HerMetric v' , ε'' ~ HerMetric v'', ε~ε', ε~ε'' ) - => ( c' -> c'' -> (c, (v',v''):-*v, ε -> (ε',ε'')) ) - -> DfblFuncValue s d c' -> DfblFuncValue s d c'' -> DfblFuncValue s d c -dfblFnValsCombine cmb (GenericAgent (Differentiable f)) - (GenericAgent (Differentiable g)) - = GenericAgent . Differentiable $ - \d -> let (c', f', devf) = f d - (c'', g', devg) = g d - (c, h', devh) = cmb c' c'' - h'l = h' *.* lcofst; h'r = h' *.* lcosnd - in ( c - , h' *.* linear (lapply f' &&& lapply g') - , \εc -> let εc' = transformMetric h'l εc - εc'' = transformMetric h'r εc - (δc',δc'') = devh εc - in devf εc' ^+^ devg εc'' - ^+^ transformMetric f' δc' - ^+^ transformMetric g' δc'' - ) - where lcofst = linear(,zeroV) - lcosnd = linear(zeroV,) - - - - - -instance (WithField s LinearManifold v, LocallyScalable s a, Floating s) - => AdditiveGroup (DfblFuncValue s a v) where - zeroV = point zeroV - (^+^) = dfblFnValsCombine $ \a b -> (a^+^b, lPlus, const zeroV) - where lPlus = linear $ uncurry (^+^) - negateV = dfblFnValsFunc $ \a -> (negateV a, lNegate, const zeroV) - where lNegate = linear negateV - -instance (RealDimension n, LocallyScalable n a) - => Num (DfblFuncValue n a n) where - fromInteger i = point $ fromInteger i - (+) = dfblFnValsCombine $ \a b -> (a+b, lPlus, const zeroV) - where lPlus = linear $ uncurry (+) - (*) = dfblFnValsCombine $ - \a b -> ( a*b - , linear $ \(da,db) -> a*db + b*da - , \d -> let d¹₂ = sqrt d in (d¹₂,d¹₂) - -- ε δa δb = (a+δa)·(b+δb) - (a·b + (a·δa + b·δb)) - -- = δa·δb - -- so choose δa = δb = √ε - ) - negate = dfblFnValsFunc $ \a -> (negate a, lNegate, const zeroV) - where lNegate = linear negate - abs = dfblFnValsFunc dfblAbs - where dfblAbs a - | a>0 = (a, idL, dev_ε_δ $ \ε -> a + ε/2) - | a<0 = (-a, negateV idL, dev_ε_δ $ \ε -> ε/2 - a) - | otherwise = (0, zeroV, (^/ sqrt 2)) - signum = dfblFnValsFunc dfblSgn - where dfblSgn a - | a>0 = (1, zeroV, dev_ε_δ $ const a) - | a<0 = (-1, zeroV, dev_ε_δ $ \_ -> -a) - | otherwise = (0, zeroV, const $ projector 1) - - - --- VectorSpace instance is more problematic than you'd think: multiplication --- requires the allowed-deviation backpropagators to be split as square --- roots, but the square root of a nontrivial-vector-space metric requires --- an eigenbasis transform, which we have not implemented yet. --- --- instance (WithField s LinearManifold v, LocallyScalable s a, Floating s) --- => VectorSpace (DfblFuncValue s a v) where --- type Scalar (DfblFuncValue s a v) = DfblFuncValue s a (Scalar v) --- (*^) = dfblFnValsCombine $ \μ v -> (μ*^v, lScl, \ε -> (ε ^* sqrt 2, ε ^* sqrt 2)) --- where lScl = linear $ uncurry (*^) - - --- | Important special operator needed to compute intersection of 'Region's. -minDblfuncs :: (LocallyScalable s m, RealDimension s) - => Differentiable s m s -> Differentiable s m s -> Differentiable s m s -minDblfuncs (Differentiable f) (Differentiable g) = Differentiable h - where h x - | fx==gx = ( fx, (f'^+^g')^/2 - , \d -> devf d ^+^ devg d - ^+^ transformMetric (f'^-^g') - (projector $ metric d 1) ) - | fx < gx = ( fx, f' - , \d -> devf d - ^+^ transformMetric (f'^-^g') - (projector $ metric d 1 + gx - fx) ) - where (fx, f', devf) = f x - (gx, g', devg) = g x - - -postEndo :: ∀ c a b . (HasAgent c, Object c a, Object c b) - => c a a -> GenericAgent c b a -> GenericAgent c b a -postEndo = genericAgentMap - - --- | A pathwise connected subset of a manifold @m@, whose tangent space has scalar @s@. -data Region s m = Region { regionRefPoint :: m - , regionRDef :: PreRegion s m } - --- | A 'PreRegion' needs to be associated with a certain reference point ('Region' --- includes that point) to define a connected subset of a manifold. -data PreRegion s m where - GlobalRegion :: PreRegion s m - PreRegion :: (Differentiable s m s) -- A function that is positive at reference point /p/, - -- decreases and crosses zero at the region's - -- boundaries. (If it goes positive again somewhere - -- else, these areas shall /not/ be considered - -- belonging to the (by definition connected) region.) - -> PreRegion s m - --- | Set-intersection of regions would not be guaranteed to yield a connected result --- or even have the reference point of one region contained in the other. This --- combinator assumes (unchecked) that the references are in a connected --- sub-intersection, which is used as the result. -unsafePreRegionIntersect :: (RealDimension s, LocallyScalable s a) - => PreRegion s a -> PreRegion s a -> PreRegion s a -unsafePreRegionIntersect GlobalRegion r = r -unsafePreRegionIntersect r GlobalRegion = r -unsafePreRegionIntersect (PreRegion ra) (PreRegion rb) = PreRegion $ minDblfuncs ra rb - --- | Cartesian product of two regions. -regionProd :: (RealDimension s, LocallyScalable s a, LocallyScalable s b) - => Region s a -> Region s b -> Region s (a,b) -regionProd (Region a₀ ra) (Region b₀ rb) = Region (a₀,b₀) (preRegionProd ra rb) - --- | Cartesian product of two pre-regions. -preRegionProd :: (RealDimension s, LocallyScalable s a, LocallyScalable s b) - => PreRegion s a -> PreRegion s b -> PreRegion s (a,b) -preRegionProd GlobalRegion GlobalRegion = GlobalRegion -preRegionProd GlobalRegion (PreRegion rb) = PreRegion $ rb . snd -preRegionProd (PreRegion ra) GlobalRegion = PreRegion $ ra . fst -preRegionProd (PreRegion ra) (PreRegion rb) = PreRegion $ minDblfuncs (ra.fst) (rb.snd) - - -positivePreRegion, negativePreRegion :: (RealDimension s) => PreRegion s s -positivePreRegion = PreRegion $ Differentiable prr - where prr x = (1 - 1/xp1, (1/xp1²) *^ idL, dev_ε_δ δ ) - -- ε = (1 − 1/(1+x)) + (-δ · 1/(x+1)²) − (1 − 1/(1+x−δ)) - -- = 1/(1+x−δ) − 1/(1+x) − δ · 1/(x+1)² - -- ε·(1+x−δ) = 1 − (1+x−δ)/(1+x) − δ·(1+x-δ)/(x+1)² - -- ε + ε·x − ε·δ = 1 − 1/(1+x) − x/(1+x) + δ/(1+x) − δ/(x+1) + δ²/(x+1)² - -- = 1 − 1/(1+x) − x/(1+x) + δ²/(x+1)² - -- = (1+x − 1 − x)/(1+x) + δ²/(x+1)² - -- 0 = δ² + ε·(x+1)²·δ + ε·(x+1)³ - -- δ = let mph = -ε·(x+1)²/2 - -- in mph + sqrt(mph² - ε·(x+1)³) - where δ ε = let mph = -ε*xp1²/2 - in mph + sqrt(mph^2 - ε * xp1² * xp1) - xp1 = (x+1) - xp1² = xp1 ^ 2 -negativePreRegion = PreRegion $ ppr . ngt - where PreRegion ppr = positivePreRegion - ngt = actuallyLinear $ linear negate - -preRegionToInfFrom, preRegionFromMinInfTo :: RealDimension s => s -> PreRegion s s -preRegionToInfFrom xs = PreRegion $ ppr . trl - where PreRegion ppr = positivePreRegion - trl = actuallyAffine (-xs) idL -preRegionFromMinInfTo xe = PreRegion $ ppr . flp - where PreRegion ppr = positivePreRegion - flp = actuallyAffine (-xe) (linear negate) - -intervalPreRegion :: RealDimension s => (s,s) -> PreRegion s s -intervalPreRegion (lb,rb) = PreRegion $ Differentiable prr - where m = lb + radius; radius = (rb - lb)/2 - prr x = ( 1 - ((x-m)/radius)^2 - , (2*(m-x)/radius^2) *^ idL - , dev_ε_δ $ (*radius) . sqrt ) - - - - --- | Category of functions that almost everywhere have an open region in --- which they are continuously differentiable, i.e. /PieceWiseDiff'able/. -newtype PWDiffable s d c - = PWDiffable { - getDfblDomain :: d -> (PreRegion s d, Differentiable s d c) } - - - -instance (RealDimension s) => Category (PWDiffable s) where - type Object (PWDiffable s) o = LocallyScalable s o - id = PWDiffable $ \x -> (GlobalRegion, id) - PWDiffable f . PWDiffable g = PWDiffable h - where h x₀ = case g x₀ of - (GlobalRegion, gr) - -> let (y₀,_,_) = runDifferentiable gr x₀ - in case f y₀ of - (GlobalRegion, fr) -> (GlobalRegion, fr . gr) - (PreRegion ry, fr) - -> ( PreRegion $ ry . gr, fr . gr ) - (PreRegion rx, gr) - -> let (y₀,_,_) = runDifferentiable gr x₀ - in case f y₀ of - (GlobalRegion, fr) -> (PreRegion rx, fr . gr) - (PreRegion ry, fr) - -> ( PreRegion $ minDblfuncs (ry . gr) rx - , fr . gr ) - where (rx, gr) = g x₀ - -globalDiffable :: Differentiable s a b -> PWDiffable s a b -globalDiffable f = PWDiffable $ const (GlobalRegion, f) - -instance (RealDimension s) => EnhancedCat (PWDiffable s) (Differentiable s) where - arr = globalDiffable -instance (RealDimension s) => EnhancedCat (->) (PWDiffable s) where - arr (PWDiffable g) x = let (_,Differentiable f) = g x - (y,_,_) = f x - in y - - -instance (RealDimension s) => Cartesian (PWDiffable s) where - type UnitObject (PWDiffable s) = ZeroDim s - swap = globalDiffable swap - attachUnit = globalDiffable attachUnit - detachUnit = globalDiffable detachUnit - regroup = globalDiffable regroup - regroup' = globalDiffable regroup' - -instance (RealDimension s) => Morphism (PWDiffable s) where - PWDiffable f *** PWDiffable g = PWDiffable h - where h (x,y) = (preRegionProd rfx rgy, dff *** dfg) - where (rfx, dff) = f x - (rgy, dfg) = g y - -instance (RealDimension s) => PreArrow (PWDiffable s) where - PWDiffable f &&& PWDiffable g = PWDiffable h - where h x = (unsafePreRegionIntersect rfx rgx, dff &&& dfg) - where (rfx, dff) = f x - (rgx, dfg) = g x - terminal = globalDiffable terminal - fst = globalDiffable fst - snd = globalDiffable snd - - -instance (RealDimension s) => WellPointed (PWDiffable s) where - unit = Tagged Origin - globalElement x = PWDiffable $ \Origin -> (GlobalRegion, globalElement x) - const x = PWDiffable $ \_ -> (GlobalRegion, const x) - - -type PWDfblFuncValue s = GenericAgent (PWDiffable s) - -instance RealDimension s => HasAgent (PWDiffable s) where - alg = genericAlg - ($~) = genericAgentMap -instance RealDimension s => CartesianAgent (PWDiffable s) where - alg1to2 = genericAlg1to2 - alg2to1 = genericAlg2to1 - alg2to2 = genericAlg2to2 -instance (RealDimension s) - => PointAgent (PWDfblFuncValue s) (PWDiffable s) a x where - point = genericPoint - -gpwDfblFnValsFunc - :: ( RealDimension s - , LocallyScalable s c, LocallyScalable s c', LocallyScalable s d - , v ~ Needle c, v' ~ Needle c' - , ε ~ HerMetric v, ε ~ HerMetric v' ) - => (c' -> (c, v':-*v, ε->ε)) -> PWDfblFuncValue s d c' -> PWDfblFuncValue s d c -gpwDfblFnValsFunc f = (PWDiffable (\_ -> (GlobalRegion, Differentiable f)) $~) - -gpwDfblFnValsCombine :: forall d c c' c'' v v' v'' ε ε' ε'' s. - ( LocallyScalable s c, LocallyScalable s c', LocallyScalable s c'' - , LocallyScalable s d, RealDimension s - , v ~ Needle c, v' ~ Needle c', v'' ~ Needle c'' - , ε ~ HerMetric v , ε' ~ HerMetric v' , ε'' ~ HerMetric v'', ε~ε', ε~ε'' ) - => ( c' -> c'' -> (c, (v',v''):-*v, ε -> (ε',ε'')) ) - -> PWDfblFuncValue s d c' -> PWDfblFuncValue s d c'' -> PWDfblFuncValue s d c -gpwDfblFnValsCombine cmb (GenericAgent (PWDiffable fpcs)) - (GenericAgent (PWDiffable gpcs)) - = GenericAgent . PWDiffable $ - \d₀ -> let (rc', Differentiable f) = fpcs d₀ - (rc'',Differentiable g) = gpcs d₀ - in (unsafePreRegionIntersect rc' rc'',) . Differentiable $ - \d -> let (c', f', devf) = f d - (c'',g', devg) = g d - (c, h', devh) = cmb c' c'' - h'l = h' *.* lcofst; h'r = h' *.* lcosnd - in ( c - , h' *.* linear (lapply f' &&& lapply g') - , \εc -> let εc' = transformMetric h'l εc - εc'' = transformMetric h'r εc - (δc',δc'') = devh εc - in devf εc' ^+^ devg εc'' - ^+^ transformMetric f' δc' - ^+^ transformMetric g' δc'' - ) - where lcofst = linear(,zeroV) - lcosnd = linear(zeroV,) - - -instance (WithField s LinearManifold v, LocallyScalable s a, RealDimension s) - => AdditiveGroup (PWDfblFuncValue s a v) where - zeroV = point zeroV - (^+^) = gpwDfblFnValsCombine $ \a b -> (a^+^b, lPlus, const zeroV) - where lPlus = linear $ uncurry (^+^) - negateV = gpwDfblFnValsFunc $ \a -> (negateV a, lNegate, const zeroV) - where lNegate = linear negateV - -instance (RealDimension n, LocallyScalable n a) - => Num (PWDfblFuncValue n a n) where - fromInteger i = point $ fromInteger i - (+) = gpwDfblFnValsCombine $ \a b -> (a+b, lPlus, const zeroV) - where lPlus = linear $ uncurry (+) - (*) = gpwDfblFnValsCombine $ - \a b -> ( a*b - , linear $ \(da,db) -> a*db + b*da - , \d -> let d¹₂ = sqrt d in (d¹₂,d¹₂) - ) - negate = gpwDfblFnValsFunc $ \a -> (negate a, lNegate, const zeroV) - where lNegate = linear negate - abs = (PWDiffable absPW $~) - where absPW a₀ - | a₀<0 = (negativePreRegion, desc) - | otherwise = (positivePreRegion, asc) - desc = actuallyLinear $ linear negate - asc = actuallyLinear idL - signum = (PWDiffable sgnPW $~) - where sgnPW a₀ - | a₀<0 = (negativePreRegion, const 1) - | otherwise = (positivePreRegion, const $ -1) - -instance (RealDimension n, LocallyScalable n a) - => Fractional (PWDfblFuncValue n a n) where - fromRational i = point $ fromRational i - recip = postEndo . PWDiffable $ \a₀ -> if a₀<0 - then (negativePreRegion, Differentiable negp) - else (positivePreRegion, Differentiable posp) - where negp x = (x'¹, (- x'¹^2) *^ idL, dev_ε_δ δ) - -- ε = 1/x − δ/x² − 1/(x+δ) - -- ε·x + ε·δ = 1 + δ/x − δ/x − δ²/x² − 1 - -- = -δ²/x² - -- 0 = δ² + ε·x²·δ + ε·x³ - -- δ = let mph = -ε·x²/2 in mph + sqrt (mph² − ε·x³) - where δ ε = let mph = -ε*x^2/2 in mph + sqrt (mph^2 - ε*x^3) - x'¹ = recip x - posp x = (x'¹, (- x'¹^2) *^ idL, dev_ε_δ δ) - where δ ε = let mph = -ε*x^2/2 in mph + sqrt (mph^2 + ε*x^3) - x'¹ = recip x - - - - - - --- | Category of functions that, where defined, have an open region in --- which they are continuously differentiable. Hence /RegionWiseDiff'able/. --- Basically these are the partial version of `PWDiffable`. --- --- Though the possibility of undefined regions is of course not too nice --- (we don't need Java to demonstrate this with its everywhere-looming @null@ values...), --- this category will propably be the “workhorse” for most serious --- calculus applications, because it contains all the usual trig etc. functions --- and of course everything algebraic you can do in the reals. --- --- The easiest way to define ordinary functions in this category is hence --- with its 'AgentVal'ues, which have instances of the standard classes 'Num' --- through 'Floating'. For instance, the following defines the /binary entropy/ --- as a differentiable function on the interval @]0,1[@: (it will --- actually /know/ where it's defined and where not! – and I don't mean you --- need to exhaustively 'isNaN'-check all results...) --- --- @ --- hb :: RWDiffable ℝ ℝ ℝ --- hb = alg (\\p -> - p * logBase 2 p - (1-p) * logBase 2 (1-p) ) --- @ -newtype RWDiffable s d c - = RWDiffable { - tryDfblDomain :: d -> (PreRegion s d, Option (Differentiable s d c)) } - -notDefinedHere :: Option (Differentiable s d c) -notDefinedHere = Option Nothing - - - -instance (RealDimension s) => Category (RWDiffable s) where - type Object (RWDiffable s) o = LocallyScalable s o - id = RWDiffable $ \x -> (GlobalRegion, pure id) - RWDiffable f . RWDiffable g = RWDiffable h - where h x₀ = case g x₀ of - (GlobalRegion, Option Nothing) - -> (GlobalRegion, notDefinedHere) - (GlobalRegion, Option (Just gr)) - -> let (y₀,_,_) = runDifferentiable gr x₀ - in case f y₀ of - (GlobalRegion, Option Nothing) - -> (GlobalRegion, notDefinedHere) - (GlobalRegion, Option (Just fr)) - -> (GlobalRegion, pure (fr . gr)) - (PreRegion ry, Option Nothing) - -> ( PreRegion $ ry . gr, Option Nothing ) - (PreRegion ry, Option (Just fr)) - -> ( PreRegion $ ry . gr, pure (fr . gr) ) - (PreRegion rx, Option Nothing) - -> (PreRegion rx, notDefinedHere) - (PreRegion rx, Option (Just gr)) - -> let (y₀,_,_) = runDifferentiable gr x₀ - in case f y₀ of - (GlobalRegion, Option Nothing) - -> (PreRegion rx, notDefinedHere) - (GlobalRegion, Option (Just fr)) - -> (PreRegion rx, pure (fr . gr)) - (PreRegion ry, Option Nothing) - -> ( PreRegion $ minDblfuncs (ry . gr) rx - , notDefinedHere ) - (PreRegion ry, Option (Just fr)) - -> ( PreRegion $ minDblfuncs (ry . gr) rx - , pure (fr . gr) ) - where (rx, gr) = g x₀ - - -globalDiffable' :: Differentiable s a b -> RWDiffable s a b -globalDiffable' f = RWDiffable $ const (GlobalRegion, pure f) - -pwDiffable :: PWDiffable s a b -> RWDiffable s a b -pwDiffable (PWDiffable q) = RWDiffable $ \x₀ -> let (r₀,f₀) = q x₀ in (r₀, pure f₀) - - - -instance (RealDimension s) => EnhancedCat (RWDiffable s) (Differentiable s) where - arr = globalDiffable' -instance (RealDimension s) => EnhancedCat (RWDiffable s) (PWDiffable s) where - arr = pwDiffable - -instance (RealDimension s) => Cartesian (RWDiffable s) where - type UnitObject (RWDiffable s) = ZeroDim s - swap = globalDiffable' swap - attachUnit = globalDiffable' attachUnit - detachUnit = globalDiffable' detachUnit - regroup = globalDiffable' regroup - regroup' = globalDiffable' regroup' - -instance (RealDimension s) => Morphism (RWDiffable s) where - RWDiffable f *** RWDiffable g = RWDiffable h - where h (x,y) = (preRegionProd rfx rgy, liftA2 (***) dff dfg) - where (rfx, dff) = f x - (rgy, dfg) = g y - -instance (RealDimension s) => PreArrow (RWDiffable s) where - RWDiffable f &&& RWDiffable g = RWDiffable h - where h x = (unsafePreRegionIntersect rfx rgx, liftA2 (&&&) dff dfg) - where (rfx, dff) = f x - (rgx, dfg) = g x - terminal = globalDiffable' terminal - fst = globalDiffable' fst - snd = globalDiffable' snd - - -instance (RealDimension s) => WellPointed (RWDiffable s) where - unit = Tagged Origin - globalElement x = RWDiffable $ \Origin -> (GlobalRegion, pure (globalElement x)) - const x = RWDiffable $ \_ -> (GlobalRegion, pure (const x)) - - -type RWDfblFuncValue s = GenericAgent (RWDiffable s) - -instance RealDimension s => HasAgent (RWDiffable s) where - alg = genericAlg - ($~) = genericAgentMap -instance RealDimension s => CartesianAgent (RWDiffable s) where - alg1to2 = genericAlg1to2 - alg2to1 = genericAlg2to1 - alg2to2 = genericAlg2to2 -instance (RealDimension s) - => PointAgent (RWDfblFuncValue s) (RWDiffable s) a x where - point = genericPoint - -grwDfblFnValsFunc - :: ( RealDimension s - , LocallyScalable s c, LocallyScalable s c', LocallyScalable s d - , v ~ Needle c, v' ~ Needle c' - , ε ~ HerMetric v, ε ~ HerMetric v' ) - => (c' -> (c, v':-*v, ε->ε)) -> RWDfblFuncValue s d c' -> RWDfblFuncValue s d c -grwDfblFnValsFunc f = (RWDiffable (\_ -> (GlobalRegion, pure (Differentiable f))) $~) - -grwDfblFnValsCombine :: forall d c c' c'' v v' v'' ε ε' ε'' s. - ( LocallyScalable s c, LocallyScalable s c', LocallyScalable s c'' - , LocallyScalable s d, RealDimension s - , v ~ Needle c, v' ~ Needle c', v'' ~ Needle c'' - , ε ~ HerMetric v , ε' ~ HerMetric v' , ε'' ~ HerMetric v'', ε~ε', ε~ε'' ) - => ( c' -> c'' -> (c, (v',v''):-*v, ε -> (ε',ε'')) ) - -> RWDfblFuncValue s d c' -> RWDfblFuncValue s d c'' -> RWDfblFuncValue s d c -grwDfblFnValsCombine cmb (GenericAgent (RWDiffable fpcs)) - (GenericAgent (RWDiffable gpcs)) - = GenericAgent . RWDiffable $ - \d₀ -> let (rc', fmay) = fpcs d₀ - (rc'',gmay) = gpcs d₀ - in (unsafePreRegionIntersect rc' rc'',) $ - case (fmay,gmay) of - (Option(Just(Differentiable f)), Option(Just(Differentiable g))) -> - pure . Differentiable $ \d - -> let (c', f', devf) = f d - (c'',g', devg) = g d - (c, h', devh) = cmb c' c'' - h'l = h' *.* lcofst; h'r = h' *.* lcosnd - in ( c - , h' *.* linear (lapply f' &&& lapply g') - , \εc -> let εc' = transformMetric h'l εc - εc'' = transformMetric h'r εc - (δc',δc'') = devh εc - in devf εc' ^+^ devg εc'' - ^+^ transformMetric f' δc' - ^+^ transformMetric g' δc'' - ) - _ -> notDefinedHere - where lcofst = linear(,zeroV) - lcosnd = linear(zeroV,) - - - -instance (WithField s LinearManifold v, LocallyScalable s a, RealDimension s) - => AdditiveGroup (RWDfblFuncValue s a v) where - zeroV = point zeroV - (^+^) = grwDfblFnValsCombine $ \a b -> (a^+^b, lPlus, const zeroV) - where lPlus = linear $ uncurry (^+^) - negateV = grwDfblFnValsFunc $ \a -> (negateV a, lNegate, const zeroV) - where lNegate = linear negateV - -instance (RealDimension n, LocallyScalable n a) - => Num (RWDfblFuncValue n a n) where - fromInteger i = point $ fromInteger i - (+) = grwDfblFnValsCombine $ \a b -> (a+b, lPlus, const zeroV) - where lPlus = linear $ uncurry (+) - (*) = grwDfblFnValsCombine $ - \a b -> ( a*b - , linear $ \(da,db) -> a*db + b*da - , \d -> let d¹₂ = sqrt d in (d¹₂,d¹₂) - ) - negate = grwDfblFnValsFunc $ \a -> (negate a, lNegate, const zeroV) - where lNegate = linear negate - abs = (RWDiffable absPW $~) - where absPW a₀ - | a₀<0 = (negativePreRegion, pure desc) - | otherwise = (positivePreRegion, pure asc) - desc = actuallyLinear $ linear negate - asc = actuallyLinear idL - signum = (RWDiffable sgnPW $~) - where sgnPW a₀ - | a₀<0 = (negativePreRegion, pure (const 1)) - | otherwise = (positivePreRegion, pure (const $ -1)) - -instance (RealDimension n, LocallyScalable n a) - => Fractional (RWDfblFuncValue n a n) where - fromRational i = point $ fromRational i - recip = postEndo . RWDiffable $ \a₀ -> if a₀<0 - then (negativePreRegion, pure (Differentiable negp)) - else (positivePreRegion, pure (Differentiable posp)) - where negp x = (x'¹, (- x'¹^2) *^ idL, dev_ε_δ δ) - -- ε = 1/x − δ/x² − 1/(x+δ) - -- ε·x + ε·δ = 1 + δ/x − δ/x − δ²/x² − 1 - -- = -δ²/x² - -- 0 = δ² + ε·x²·δ + ε·x³ - -- δ = let mph = -ε·x²/2 in mph + sqrt (mph² − ε·x³) - where δ ε = let mph = -ε*x^2/2 in mph + sqrt (mph^2 - ε*x^3) - x'¹ = recip x - posp x = (x'¹, (- x'¹^2) *^ idL, dev_ε_δ δ) - where δ ε = let mph = -ε*x^2/2 in mph + sqrt (mph^2 + ε*x^3) - x'¹ = recip x - - - - - --- Helper for checking ε-estimations in GHCi with dynamic-plot: --- epsEst (f,f') εsgn δf (ViewXCenter xc) (ViewHeight h) --- = let δfxc = δf xc --- in tracePlot $ reverse [ (xc - δ, f xc - δ * f' xc + εsgn*ε) | --- ε <- [0, h/500 .. h], let δ = δfxc ε] --- ++ [ (xc + δ, f xc + δ * f' xc + εsgn*ε) | --- ε <- [0, h/500 .. h], let δ = δfxc ε] --- Golfed version: --- epsEst(f,d)s φ(ViewXCenter ξ)(ViewHeight h)=let ζ=φ ξ in tracePlot$[(ξ-δ,f ξ-δ*d ξ+s*abs ε)|ε<-[-h,-0.998*h..h],let δ=ζ(abs ε)*signum ε] - -instance (RealDimension n, LocallyScalable n a) - => Floating (RWDfblFuncValue n a n) where - pi = point pi - - exp = grwDfblFnValsFunc - $ \x -> let ex = exp x - in ( ex, ex *^ idL, dev_ε_δ $ \ε -> acosh(ε/(2*ex) + 1) ) - -- ε = e^(x+δ) − eˣ − eˣ·δ - -- = eˣ·(e^δ − 1 − δ) - -- ≤ eˣ · (e^δ − 1 + e^(-δ) − 1) - -- = eˣ · 2·(cosh(δ) − 1) - -- cosh(δ) ≥ ε/(2·eˣ) + 1 - -- δ ≥ acosh(ε/(2·eˣ) + 1) - log = postEndo . RWDiffable $ \x -> if x>0 - then (positivePreRegion, pure (Differentiable lnPosR)) - else (negativePreRegion, notDefinedHere) - where lnPosR x = ( log x, recip x *^ idL, dev_ε_δ $ \ε -> x * sqrt(1 - exp(-ε)) ) - -- ε = ln x + (-δ)/x − ln(x−δ) - -- = ln (x / ((x−δ) · exp(δ/x))) - -- x/e^ε = (x−δ) · exp(δ/x) - -- let γ = δ/x ∈ [0,1[ - -- exp(-ε) = (1−γ) · e^γ - -- ≥ (1−γ) · (1+γ) - -- = 1 − γ² - -- γ ≥ sqrt(1 − exp(-ε)) - -- δ ≥ x · sqrt(1 − exp(-ε)) - - sqrt = postEndo . RWDiffable $ \x -> if x>0 - then (positivePreRegion, pure (Differentiable sqrtPosR)) - else (negativePreRegion, notDefinedHere) - where sqrtPosR x = ( sx, idL ^/ (2*sx), dev_ε_δ $ - \ε -> 2 * (s2 * sqrt sx^3 * sqrt ε + signum (ε*2-sx) * sx * ε) ) - where sx = sqrt x; s2 = sqrt 2 - -- Exact inverse of O(δ²) remainder. - - sin = grwDfblFnValsFunc sinDfb - where sinDfb x = ( sx, cx *^ idL, dev_ε_δ δ ) - where sx = sin x; cx = cos x - δ ε = let δ₀ = sqrt $ 2 * ε / (abs sx + abs cx/3) - in if δ₀ < 1 -- TODO: confirm selection of δ-definition range. - then δ₀ - else max 1 $ (ε - abs sx - 1) / cos x - -- When sin x ≥ 0, cos x ≥ 0, δ ∈ [0,1[ - -- ε = sin x + δ · cos x − sin(x+δ) - -- = sin x + δ · cos x − sin x · cos δ − cos x · sin δ - -- ≤ sin x + δ · cos x − sin x · (1−δ²/2) − cos x · (δ − δ³/6) - -- = sin x · δ²/2 + cos x · δ³/6 - -- ≤ δ² · (sin x / 2 + cos x / 6) - -- δ ≥ sqrt(2 · ε / (sin x + cos x / 3)) - -- For general δ≥0, - -- ε ≤ δ · cos x + sin x + 1 - -- δ ≥ (ε − sin x − 1) / cos x - cos = sin . (globalDiffable' (actuallyAffine (pi/2) idL) $~) - - sinh x = (exp x - exp (-x))/2 - {- = grwDfblFnValsFunc sinhDfb - where sinhDfb x = ( sx, cx *^ idL, dev_ε_δ δ ) - where sx = sinh x; cx = cosh x - δ ε = undefined -} - -- ε = sinh x + δ · cosh x − sinh(x+δ) - -- = ½ · ( eˣ − e⁻ˣ + δ · (eˣ + e⁻ˣ) − exp(x+δ) + exp(-x−δ) ) - -- = ½·e⁻ˣ · ( e²ˣ − 1 + δ · (e²ˣ + 1) − e²ˣ·e^δ + e^-δ ) - -- = ½ · ( eˣ − e⁻ˣ + δ · (eˣ + e⁻ˣ) − exp(x+δ) + exp(-x−δ) ) - cosh x = (exp x + exp (-x))/2 - tanh x = (exp x - exp (-x)) / (exp x + exp (-x)) - - atan = grwDfblFnValsFunc atanDfb - where atanDfb x = ( atnx, idL ^/ (1+x^2), dev_ε_δ δ ) - where atnx = atan x - c = (atnx*2/pi)^2 - p = 1 + abs x/(2*pi) - δ ε = p * (sqrt ε + ε * c) - -- Semi-empirically obtained: with the epsEst helper, - -- it is observed that this function is (for xc≥0) a lower bound - -- to the arctangent. The growth of the p coefficient makes sense - -- and holds for arbitrarily large xc, because those move us linearly - -- away from the only place where the function is not virtually constant - -- (around 0). - - asin = postEndo . RWDiffable $ \x -> if - | x < (-1) -> (preRegionFromMinInfTo (-1), notDefinedHere) - | x > 1 -> (preRegionToInfFrom 1, notDefinedHere) - | otherwise -> (intervalPreRegion (-1,1), pure (Differentiable asinDefdR)) - where asinDefdR x = ( asinx, asin'x *^ idL, dev_ε_δ δ ) - where asinx = asin x; asin'x = recip (sqrt $ 1 - x^2) - c = 1 - x^2 - δ ε = sqrt ε * c - -- Empirical, with epsEst upper bound. - - acos = postEndo . RWDiffable $ \x -> if - | x < (-1) -> (preRegionFromMinInfTo (-1), notDefinedHere) - | x > 1 -> (preRegionToInfFrom 1, notDefinedHere) - | otherwise -> (intervalPreRegion (-1,1), pure (Differentiable acosDefdR)) - where acosDefdR x = ( acosx, acos'x *^ idL, dev_ε_δ δ ) - where acosx = acos x; acos'x = - recip (sqrt $ 1 - x^2) - c = 1 - x^2 - δ ε = sqrt ε * c -- Like for asin – it's just a translation/reflection. - - asinh = grwDfblFnValsFunc asinhDfb - where asinhDfb x = ( asinhx, idL ^/ sqrt(1+x^2), dev_ε_δ δ ) - where asinhx = asinh x - δ ε = abs x * sqrt((1 - exp(-ε))*0.8 + ε^2/(3*abs x)) + sqrt(ε/(abs x+0.5)) - -- Empirical, modified from log function (the area hyperbolic sine - -- resembles two logarithmic lobes), with epsEst-checked lower bound. - - acosh = postEndo . RWDiffable $ \x -> if x>0 - then (positivePreRegion, pure (Differentiable acoshDfb)) - else (negativePreRegion, notDefinedHere) - where acoshDfb x = ( acosh x, idL ^/ sqrt(x^2 - 2), dev_ε_δ δ ) - where δ ε = (2 - 1/sqrt x) * (s2 * sqrt sx^3 * sqrt(ε/s2) + signum (ε*s2-sx) * sx * ε/s2) - sx = sqrt(x-1) - s2 = sqrt 2 - -- Empirical, modified from sqrt function – the area hyperbolic cosine - -- strongly resembles \x -> sqrt(2 · (x-1)). - - atanh = postEndo . RWDiffable $ \x -> if - | x < (-1) -> (preRegionFromMinInfTo (-1), notDefinedHere) - | x > 1 -> (preRegionToInfFrom 1, notDefinedHere) - | otherwise -> (intervalPreRegion (-1,1), pure (Differentiable atnhDefdR)) - where atnhDefdR x = ( atanh x, recip(1-x^2) *^ idL, dev_ε_δ $ \ε -> sqrt(tanh ε)*(1-abs x) ) - -- Empirical, with epsEst upper bound. - - - - diff --git a/Data/Manifold/Riemannian.hs b/Data/Manifold/Riemannian.hs --- a/Data/Manifold/Riemannian.hs +++ b/Data/Manifold/Riemannian.hs @@ -69,7 +69,7 @@ import Data.Proxy import Data.Manifold.Types -import Data.Manifold.Types.Primitive ((^), embed, coEmbed) +import Data.Manifold.Types.Primitive ((^), empty, embed, coEmbed) import Data.Manifold.PseudoAffine import Data.VectorSpace.FiniteDimensional @@ -157,7 +157,7 @@ instance Geodesic S⁰ where geodesicBetween PositiveHalfSphere PositiveHalfSphere = return $ const PositiveHalfSphere geodesicBetween NegativeHalfSphere NegativeHalfSphere = return $ const NegativeHalfSphere - geodesicBetween _ _ = Hask.empty + geodesicBetween _ _ = empty instance Geodesic S¹ where geodesicBetween (S¹ φ) (S¹ ϕ) @@ -242,3 +242,12 @@ instance IntervalLike ℝ where toClosedInterval x = D¹ $ tanh x + + + + +class Geodesic m => Riemannian m where + rieMetric :: RieMetric m + +instance Riemannian ℝ where + rieMetric = const m where m = projector 1 diff --git a/Data/Manifold/TreeCover.hs b/Data/Manifold/TreeCover.hs --- a/Data/Manifold/TreeCover.hs +++ b/Data/Manifold/TreeCover.hs @@ -75,7 +75,7 @@ import Data.SimplicialComplex import Data.Manifold.Types -import Data.Manifold.Types.Primitive ((^)) +import Data.Manifold.Types.Primitive ((^), empty) import Data.Manifold.PseudoAffine import Data.Embedding @@ -517,7 +517,7 @@ -- | s' <- getTriangulation $ simplexFaces s ] -- where expandInDir j xs = case sortBy (comparing snd) $ filter ((> -1) . snd) xs_bc of -- ((x, q) : _) | q<0 -> pure x --- _ -> Hask.empty +-- _ -> empty -- where xs_bc = map (\x -> (x, getBaryCoord (emb >-$ x) j)) xs -- (Tagged n) = theNatN :: Tagged n Int @@ -541,7 +541,7 @@ optimalBottomExtension s xs = case filter ((>0).snd) $ zipWith ((. bottomExtendSuitability s) . (,)) [0..] xs of - [] -> Hask.empty + [] -> empty qs -> pure . fst . maximumBy (comparing snd) $ qs @@ -648,7 +648,7 @@ return $ case q of Just(_,is) | s<-bottomExtendSuitability is x, s>0 -> pure s - _ -> Hask.empty + _ -> empty return . fmap sum $ Hask.sequence scores spanSemiOpenSimplex :: ∀ t n n' x . (KnownNat n', WithField ℝ Manifold x, n~S n') 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 @@ -50,6 +50,7 @@ -- * Utility (deprecated) , NaturallyEmbedded(..) , GraphWindowSpec(..), Endomorphism, (^), (^.), EqFloating + , empty ) where @@ -60,7 +61,7 @@ import Data.Void import Data.Monoid -import Control.Applicative (Const(..)) +import Control.Applicative (Const(..), Alternative(..)) import qualified Prelude diff --git a/Data/SimplicialComplex.hs b/Data/SimplicialComplex.hs --- a/Data/SimplicialComplex.hs +++ b/Data/SimplicialComplex.hs @@ -80,7 +80,7 @@ import Data.Proxy import Data.Manifold.Types -import Data.Manifold.Types.Primitive ((^)) +import Data.Manifold.Types.Primitive ((^), empty) import Data.Manifold.PseudoAffine import Data.Embedding @@ -268,7 +268,7 @@ => TriangT t k x m y -> TriangT t n x m (Option y) onSkeleton q@(TriangT qf) = case tryToMatchTTT forgetVolumes q of Option (Just q') -> pure <$> q' - _ -> return Hask.empty + _ -> return empty newtype SimplexIT (t :: *) (n :: Nat) (x :: *) = SimplexIT { tgetSimplexIT' :: Int } @@ -380,7 +380,7 @@ [[iIVert], [jIVert]] <- forM [i,j] $ fmap (filter (not . (`elem` shVerts)) . Hask.toList) . lookSplxVerticesIT return $ pure ((iIVert, jIVert), shBound) - _ -> return Hask.empty + _ -> return empty triangulationBulk :: ∀ t m n k x . (HaskMonad m, KnownNat k, KnownNat n) => TriangT t n x m [Simplex k x] @@ -400,7 +400,7 @@ baseSups :: [SimplexIT t (S k) x] <- lookSupersimplicesIT base return $ case intersect tipSups baseSups of (res:_) -> pure res - _ -> Hask.empty + _ -> empty diff --git a/manifolds.cabal b/manifolds.cabal --- a/manifolds.cabal +++ b/manifolds.cabal @@ -1,5 +1,5 @@ Name: manifolds -Version: 0.1.5.2 +Version: 0.1.6.2 Category: Math Synopsis: Coordinate-free hypersurfaces Description: Manifolds, a generalisation of the notion of “smooth curves” or surfaces, @@ -65,7 +65,7 @@ Data.Manifold.TreeCover Data.SimplicialComplex Data.LinearMap.HerMetric - -- Data.Manifold.Visualisation.R3.GLUT + Data.Function.Differentiable Data.Manifold.Types Data.Manifold.Griddable Data.Manifold.Riemannian