manifolds 0.1.5.2 → 0.1.6.2

raw patch · 10 files changed

+1275/−1031 lines, 10 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• Data/CoNat.hs +21/−9
• Data/Function/Differentiable.hs +1172/−0
• Data/LinearMap/HerMetric.hs +46/−11
• Data/Manifold/Cone.hs +6/−32
• Data/Manifold/PseudoAffine.hs +7/−966
• Data/Manifold/Riemannian.hs +11/−2
• Data/Manifold/TreeCover.hs +4/−4
• Data/Manifold/Types/Primitive.hs +2/−1
• Data/SimplicialComplex.hs +4/−4
• manifolds.cabal +2/−2

@@ -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

@@ -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+                                      = exit (nLim'-1) (dir'/2) 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) -> l<limR && r>limL )+       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/&#x2080; derivative /f'&#x2093;/&#x2080;,+--   if, for&#xb9; /&#x3b5;/>0, there exists /&#x3b4;/ such that+--   |/f/ /x/ &#x2212; (/f/ /x/&#x2080; + /x/&#x22c5;/f'&#x2093;/&#x2080;)| < /&#x3b5;/+--   for all |/x/ &#x2212; /x/&#x2080;| < /&#x3b4;/.+-- +--   Observe that, though this looks quite similar to the standard definition+--   of differentiability, it is not equivalent thereto &#x2013; in fact it does+--   not prove any analytic properties at all. To make it equivalent, we need+--   a lower bound on /&#x3b4;/: simply /&#x3b4;/ gives us continuity, and for+--   continuous differentiability, /&#x3b4;/ must grow at least like &#x221a;/&#x3b5;/+--   for small /&#x3b5;/. 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 &#x2013; 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).+--   +--   &#xb9;(The implementation does not deal with /&#x3b5;/ and /&#x3b4;/ 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 &#x201c;workhorse&#x201d; 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! &#x2013; and I don't mean you+--   need to exhaustively 'isNaN'-check all results...)+-- +-- @+-- hb :: RWDiffable &#x211d; &#x211d; &#x211d;+-- 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]+++

@@ -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

@@ -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   

@@ -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/&#x2080; derivative /f'&#x2093;/&#x2080;,---   if, for&#xb9; /&#x3b5;/>0, there exists /&#x3b4;/ such that---   |/f/ /x/ &#x2212; (/f/ /x/&#x2080; + /x/&#x22c5;/f'&#x2093;/&#x2080;)| < /&#x3b5;/---   for all |/x/ &#x2212; /x/&#x2080;| < /&#x3b4;/.--- ---   Observe that, though this looks quite similar to the standard definition---   of differentiability, it is not equivalent thereto &#x2013; in fact it does---   not prove any analytic properties at all. To make it equivalent, we need---   a lower bound on /&#x3b4;/: simply /&#x3b4;/ gives us continuity, and for---   continuous differentiability, /&#x3b4;/ must grow at least like &#x221a;/&#x3b5;/---   for small /&#x3b5;/. 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 &#x2013; 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).---   ---   &#xb9;(The implementation does not deal with /&#x3b5;/ and /&#x3b4;/ 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 &#x201c;workhorse&#x201d; 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! &#x2013; and I don't mean you---   need to exhaustively 'isNaN'-check all results...)--- --- @--- hb :: RWDiffable &#x211d; &#x211d; &#x211d;--- 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.-  -  --   

@@ -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

@@ -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')

@@ -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 

@@ -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       

@@ -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