manifolds 0.1.6.2 → 0.1.6.3

raw patch · 10 files changed

+1022/−400 lines, 10 filesbinary-addedPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• Data/Function/Affine.hs +137/−0
• Data/Function/Differentiable.hs +451/−373
• Data/Function/Differentiable/Data.hs +124/−0
• Data/LinearMap/HerMetric.hs +64/−3
• Data/Manifold/PseudoAffine.hs +29/−3
• Data/Manifold/TreeCover.hs +93/−8
• Data/Manifold/Types/Primitive.hs +22/−1
• Data/VectorSpace/FiniteDimensional.hs +99/−11
• images/examples/Friedrichs-mollifier.png binary
• manifolds.cabal +3/−1

@@ -0,0 +1,137 @@+-- |+-- Module      : Data.Function.Affine+-- 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.Affine (+              Affine(..)+            ) where+    +++import Data.List+import Data.Maybe+import Data.Semigroup++import Data.VectorSpace+import Data.LinearMap+import Data.LinearMap.HerMetric+import Data.MemoTrie (HasTrie(..))+import Data.AffineSpace+import Data.Basis+import Data.Void+import Data.Tagged+import Data.Manifold.Types.Primitive+import Data.Manifold.PseudoAffine++import Data.CoNat+import Data.VectorSpace.FiniteDimensional++import qualified Prelude+import qualified Control.Applicative as Hask++import Control.Category.Constrained.Prelude hiding ((^))+import Control.Arrow.Constrained+import Control.Monad.Constrained+import Data.Foldable.Constrained+++++data Affine s d c+   = Affine { affineCoOffset :: d+            , affineOffset :: c+            , affineSlope :: Needle d :-* Needle c+            }++instance (RealDimension s) => EnhancedCat (->) (Affine s) where+  arr (Affine co ao sl) x = ao .+~^ lapply sl (x.-.co)+++instance (MetricScalar s) => Category (Affine s) where+  type Object (Affine s) o = WithField s LinearManifold o+  id = Affine zeroV zeroV idL+  Affine cof aof slf . Affine cog aog slg+      = Affine cog (aof .+~^ lapply slf (aog.-.cof)) (slf*.*slg)++linearAffine :: ( AdditiveGroup d, AdditiveGroup c+                , HasBasis (Needle d), HasTrie (Basis (Needle d)) )+       => (Needle d -> Needle c) -> Affine s d c+linearAffine = Affine zeroV zeroV . linear++instance (MetricScalar s) => Cartesian (Affine s) where+  type UnitObject (Affine s) = ZeroDim s+  swap = linearAffine swap+  attachUnit = linearAffine (, Origin)+  detachUnit = linearAffine fst+  regroup = linearAffine regroup+  regroup' = linearAffine regroup'++instance (MetricScalar s) => Morphism (Affine s) where+  Affine cof aof slf *** Affine cog aog slg+      = Affine (cof,cog) (aof,aog) (linear $ lapply slf *** lapply slg)++instance (MetricScalar s) => PreArrow (Affine s) where+  terminal = linearAffine $ const Origin+  fst = linearAffine fst+  snd = linearAffine snd+  Affine cof aof slf &&& Affine cog aog slg+      = Affine zeroV (aof.-^lapply slf cof, aog.-^lapply slg cog)+                 (linear $ lapply slf &&& lapply slg)++instance (MetricScalar s) => WellPointed (Affine s) where+  unit = Tagged Origin+  globalElement x = Affine zeroV x zeroV+  const x = Affine zeroV x zeroV++++type AffinFuncValue s = GenericAgent (Affine s)++instance (MetricScalar s) => HasAgent (Affine s) where+  alg = genericAlg+  ($~) = genericAgentMap+instance (MetricScalar s) => CartesianAgent (Affine s) where+  alg1to2 = genericAlg1to2+  alg2to1 = genericAlg2to1+  alg2to2 = genericAlg2to2+instance (MetricScalar s)+      => PointAgent (AffinFuncValue s) (Affine s) a x where+  point = genericPoint++++instance (WithField s LinearManifold v, WithField s LinearManifold a)+    => AdditiveGroup (AffinFuncValue s a v) where+  zeroV = GenericAgent $ Affine zeroV zeroV zeroV+  GenericAgent (Affine cof aof slf) ^+^ GenericAgent (Affine cog aog slg)+       = GenericAgent $ Affine (cof^+^cog) (aof^+^aog) (slf^+^slg)+  negateV (GenericAgent (Affine co ao sl))+      = GenericAgent $ Affine (negateV co) (negateV ao) (negateV sl)+++

@@ -28,22 +28,23 @@   module Data.Function.Differentiable (+            -- * Everywhere differentiable functions+              Differentiable+            -- * Region-wise defined diff'able functions+            , RWDiffable+            -- ** Operators for piecewise definition+            -- $definitionRegionOps+            , (?->), (?>), (?<), (?|:), backupRegions             -- * Regions within a manifold-              Region+            , Region             , smoothIndicator-            -- * Hierarchy of manifold-categories-            -- ** Everywhere differentiable functions-            , Differentiable-            -- ** Almost everywhere diff'able funcs-            , PWDiffable-            -- ** Region-wise defined diff'able funcs-            , RWDiffable-            -- * Misc+            -- * Evaluation of differentiable functions             , discretisePathIn             , discretisePathSegs             , continuityRanges             , regionOfContinuityAround             , analyseLocalBehaviour+            , intervalImages             ) where      @@ -61,6 +62,8 @@ import Data.LinearMap.HerMetric import Data.MemoTrie (HasTrie(..)) import Data.AffineSpace+import Data.Function.Differentiable.Data+import Data.Function.Affine import Data.Basis import Data.Complex hiding (magnitude) import Data.Void@@ -85,7 +88,6 @@   - 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.@@ -110,23 +112,36 @@ 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)+continuityRanges nLim δbf (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+         | RealSubray PositiveHalfSphere xl' <- rangeHere+                      = let stepl' = dir/metric (δbf xl') 2+                        in if dir>0+                            then if definedHere then [(max (xl'+stepl') x₀, huge)]+                                                else []+                            else if definedHere && x₀ > xl'+stepl'+                                  then (xl'+stepl',x₀) : go (xl'-stepl') dir+                                  else go (xl'-stepl') dir+         | RealSubray NegativeHalfSphere xr' <- rangeHere+                      = let stepr' = dir/metric (δbf xr') 2+                        in if dir<0+                            then if definedHere then [(-huge, min (xr'-stepr') x₀)]+                                                else []+                            else if definedHere && x₀ < xr'-stepr'+                                  then (x₀,xr'-stepr') : go (xr'+stepr') dir+                                  else go (xr'+stepr') dir          | otherwise  = exit nLim dir x₀-        where (PreRegion (Differentiable r₀), fq₀) = f x₀+        where (rangeHere, fq₀) = f x₀+              (PreRegion (Differentiable r₀)) = genericisePreRegion rangeHere               (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₀)]@@ -156,8 +171,7 @@        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+       xc = 0  discretisePathSegs :: WithField ℝ Manifold y       => Int              -- ^ Maximum number of path segments and/or points per segment.@@ -166,15 +180,16 @@                           --   (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.+      -> RWDiffable ℝ ℝ y -- ^ Path specification. It is recommended that this+                          --   function be limited to a compact interval (e.g. with+                          --   '?>', '?<' and '?->'). For many functions the discretisation+                          --   will even work on an infinite interval: the point density+                          --   is exponentially decreased towards the infinities. But+                          --   this is still pretty bad for performance.       -> ([[(ℝ,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 )+discretisePathSegs nLim (mx,my) f@(RWDiffable ff)+                            = ( map discretise ivsL, map discretise ivsR )+ where (ivsL, ivsR) = continuityRanges nLim mx f        discretise rng@(l,r) = discretisePathIn nLim rng (mx,my) fr         where (_, Option (Just fr)) = ff $ (l+r)/2 @@ -196,21 +211,60 @@                     | otherwise  = pure 0               in ((fx, lapply j 1), epsprop)        _ -> empty- where inRegion GlobalRegion _ = True+ where                                    -- This check shouldn't really be necessary,+                                          -- because the initial value lies by definition+       inRegion GlobalRegion _ = True     -- in its domain.        inRegion (PreRegion (Differentiable rf)) x          | (yr,_,_) <- rf x   = yr>0+       inRegion (RealSubray PositiveHalfSphere xl) x = x>xl+       inRegion (RealSubray NegativeHalfSphere xr) x = x<xr  -- | 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+smoothIndicator (Region _ r₀) = let (PreRegion r) = genericisePreRegion r₀+                                in  r  regionOfContinuityAround :: RWDiffable ℝ q x -> q -> Region ℝ q regionOfContinuityAround (RWDiffable f) q = Region q . fst . f $ q                +intervalImages ::+         Int                         -- ^ Max number of exploration steps per region+      -> (RieMetric ℝ, RieMetric ℝ)  -- ^ Needed resolution in (x,y) direction+      -> RWDiffable ℝ ℝ ℝ            -- ^ Function to investigate+      -> ( [(ℝInterval,ℝInterval)]+         , [(ℝInterval,ℝInterval)] ) -- ^ (XInterval, YInterval) rectangles in which+                                     --   the function graph lies.+intervalImages nLim (mx,my) f@(RWDiffable fd)+                  = (map (id&&&ivimg) domsL, map (id&&&ivimg) domsR)+ where (domsL, domsR) = continuityRanges nLim mx f+       ivimg (xl,xr) = go xl 1 i₀ ∪ go xr (-1) i₀+        where (_, Option (Just fdd@(Differentiable fddd))) = fd xc+              xc = (xl+xr)/2+              i₀ = minimum&&&maximum $ [fdd$xl, fdd$xc, fdd$xr]+              go x dir (a,b)+                 | dir>0 && x>xc   = (a,b)+                 | dir<0 && x<xc   = (a,b)+                 | χ == 0          = (y + (x-xl)*y', y + (x-xr)*y')+                 | y < a+resoHere  = go (x + dir/χ) dir (y,b)+                 | y > b-resoHere  = go (x + dir/χ) dir (a,y)+                 | otherwise       = go (x + safeStep stepOut₀) dir (a,b)+               where (y, j, δε) = fddd x+                     y' = lapply j 1+                     εx = my y+                     resoHere = metricAsLength εx+                     χ = metric (δε εx) 1+                     safeStep s₀+                         | as_devεδ δε (safetyMarg s₀) > abs s₀  = s₀+                         | otherwise                             = safeStep (s₀*0.5)+                     stepOut₀ | y'*dir>0   = 0.5 * (b-y)/y'+                              | otherwise  = -0.5 * (y-a)/y'+                     safetyMarg stp = minimum [y-a, y+stp*y'-a, b-y, b-y-stp*y']+       infixl 3 ∪+       (a,b) ∪ (c,d) = (min a c, max b d) + hugeℝVal :: ℝ hugeℝVal = 1e+100 @@ -219,8 +273,6 @@   -type LinDevPropag d c = Metric c -> Metric d- unsafe_dev_ε_δ :: RealDimension a                 => String -> (a -> a) -> LinDevPropag a a unsafe_dev_ε_δ errHint f d@@ -251,51 +303,14 @@                     = 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 ℝ +genericiseDifferentiable :: (LocallyScalable s d, LocallyScalable s c)+                    => Differentiable s d c -> Differentiable s d c+genericiseDifferentiable (AffinDiffable (Affine x₀ y₀ f))+     = Differentiable $ \x -> (y₀ .+^ lapply f (x.-.x₀), f, const zeroV)+genericiseDifferentiable f = f + instance (MetricScalar s) => Category (Differentiable s) where   type Object (Differentiable s) o = LocallyScalable s o   id = Differentiable $ \x -> (x, idL, const zeroV)@@ -306,10 +321,16 @@                               εy = devf δz                           in transformMetric g' εy ^+^ devg δy ^+^ devg εy            in (z, f'*.*g', devfg)+  AffinDiffable f . AffinDiffable g = AffinDiffable $ f . g+  f . g = genericiseDifferentiable f . genericiseDifferentiable g  +-- instance (RealDimension s) => EnhancedCat (Differentiable s) (Affine s) where+--   arr (Affine co ao sl) = actuallyAffine (ao .-^ lapply sl co) sl+   instance (RealDimension s) => EnhancedCat (->) (Differentiable s) where   arr (Differentiable f) x = let (y,_,_) = f x in y+  arr (AffinDiffable f) x = f $ x  instance (MetricScalar s) => Cartesian (Differentiable s) where   type UnitObject (Differentiable s) = ZeroDim s@@ -337,6 +358,8 @@                 lPar = linear $ lapply f'***lapply g'          lfst = linear fst; lsnd = linear snd          lcofst = linear (,zeroV); lcosnd = linear (zeroV,)+  AffinDiffable f *** AffinDiffable g = AffinDiffable $ f *** g+  f *** g = genericiseDifferentiable f *** genericiseDifferentiable g   instance (MetricScalar s) => PreArrow (Differentiable s) where@@ -353,6 +376,7 @@                            ^+^ (devg $ transformMetric lcosnd δs)                 lFanout = linear $ lapply f'&&&lapply g'          lcofst = linear (,zeroV); lcosnd = linear (zeroV,)+  f &&& g = genericiseDifferentiable f &&& genericiseDifferentiable g   instance (MetricScalar s) => WellPointed (Differentiable s) where@@ -377,15 +401,23 @@   -actuallyLinear :: ( WithField s LinearManifold x, WithField s LinearManifold y )+actuallyLinear :: ( WithField s LinearManifold x, WithField s LinearManifold y, x~y )             => (x:-*y) -> Differentiable s x y-actuallyLinear f = Differentiable $ \x -> (lapply f x, f, const zeroV)+actuallyLinear f = actuallyAffine zeroV f -actuallyAffine :: ( WithField s LinearManifold x, WithField s AffineManifold y )+actuallyAffine :: ( WithField s LinearManifold x+                  , WithField s LinearManifold y -- Really, this should only need `AffineManifold`.+                  , x~y+                  )             => y -> (x:-*Diff y) -> Differentiable s x y-actuallyAffine y₀ f = Differentiable $ \x -> (y₀ .+^ lapply f x, f, const zeroV)+actuallyAffine y₀ f = AffinDiffable $ Affine zeroV y₀ f  +-- affinPoint :: (WithField s LinearManifold c, WithField s LinearManifold d)+--                   => c -> DfblFuncValue s d c+-- affinPoint p = GenericAgent (AffinDiffable (const p))++ dfblFnValsFunc :: ( LocallyScalable s c, LocallyScalable s c', LocallyScalable s d                   , v ~ Needle c, v' ~ Needle c'                   , ε ~ HerMetric v, ε ~ HerMetric v' )@@ -417,6 +449,9 @@                  )  where lcofst = linear(,zeroV)        lcosnd = linear(zeroV,) +dfblFnValsCombine cmb (GenericAgent fa) (GenericAgent ga) +         = dfblFnValsCombine cmb (GenericAgent $ genericiseDifferentiable fa)+                                 (GenericAgent $ genericiseDifferentiable ga)   @@ -425,16 +460,21 @@ 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)+  GenericAgent (AffinDiffable f) ^+^ GenericAgent (AffinDiffable g)+       = let (GenericAgent h) = GenericAgent f ^+^ GenericAgent g+         in GenericAgent $ AffinDiffable h+  α^+^β = dfblFnValsCombine (\a b -> (a^+^b, lPlus, const zeroV)) α β       where lPlus = linear $ uncurry (^+^)-  negateV = dfblFnValsFunc $ \a -> (negateV a, lNegate, const zeroV)+  negateV (GenericAgent (AffinDiffable f))+       = let (GenericAgent h) = negateV $ GenericAgent f+         in GenericAgent $ AffinDiffable h+  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@@ -443,8 +483,7 @@                            --         = δa·δb                            --   so choose δa = δb = √ε                   )-  negate = dfblFnValsFunc $ \a -> (negate a, lNegate, const zeroV)-      where lNegate = linear negate+  negate = negateV   abs = dfblFnValsFunc dfblAbs    where dfblAbs a           | a>0        = (a, idL, unsafe_dev_ε_δ("abs "++show a) $ \ε -> a + ε/2) @@ -496,20 +535,13 @@ 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+genericisePreRegion :: (RealDimension s, LocallyScalable s m)+                          => PreRegion s m -> PreRegion s m+genericisePreRegion GlobalRegion = PreRegion $ const 1+genericisePreRegion (RealSubray PositiveHalfSphere xl) = preRegionToInfFrom' xl+genericisePreRegion (RealSubray NegativeHalfSphere xr) = preRegionFromMinInfTo' xr+genericisePreRegion r = r  -- | 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@@ -519,7 +551,13 @@                   => PreRegion s a -> PreRegion s a -> PreRegion s a unsafePreRegionIntersect GlobalRegion r = r unsafePreRegionIntersect r GlobalRegion = r+unsafePreRegionIntersect (RealSubray PositiveHalfSphere xl) (RealSubray PositiveHalfSphere xl')+                 = RealSubray PositiveHalfSphere $ max xl xl'+unsafePreRegionIntersect (RealSubray NegativeHalfSphere xr) (RealSubray NegativeHalfSphere xr')+                 = RealSubray NegativeHalfSphere $ min xr xr' unsafePreRegionIntersect (PreRegion ra) (PreRegion rb) = PreRegion $ minDblfuncs ra rb+unsafePreRegionIntersect ra rb+   = unsafePreRegionIntersect (genericisePreRegion ra) (genericisePreRegion rb)  -- | Cartesian product of two regions. regionProd :: (RealDimension s, LocallyScalable s a, LocallyScalable s b)@@ -533,10 +571,16 @@ 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)+preRegionProd ra rb = preRegionProd (genericisePreRegion ra) (genericisePreRegion rb)   positivePreRegion, negativePreRegion :: (RealDimension s) => PreRegion s s-positivePreRegion = PreRegion $ Differentiable prr+positivePreRegion = RealSubray PositiveHalfSphere 0+negativePreRegion = RealSubray NegativeHalfSphere 0+++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) δ )@@ -574,17 +618,21 @@                   | otherwise  = ε * x / ((1+ε)/x + ε)               xp1 = (x+1)               xp1² = xp1 ^ 2-negativePreRegion = PreRegion $ ppr . ngt- where PreRegion ppr = positivePreRegion+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+preRegionToInfFrom = RealSubray PositiveHalfSphere+preRegionFromMinInfTo = RealSubray NegativeHalfSphere++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)+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@@ -596,257 +644,93 @@   --- | 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) )+  RWDiffable f . RWDiffable g = RWDiffable h where+   h x₀ = case g x₀ of+           ( rg, Option (Just gr'@(AffinDiffable gr@(Affine cog aog slg))) )+            -> let y₀ = gr $ x₀+               in case f y₀ of+                   (GlobalRegion, Option (Just (AffinDiffable fr)))+                         -> (rg, Option (Just (AffinDiffable (fr.gr))))+                   (GlobalRegion, fhr)+                         -> (rg, fmap (. gr') fhr)+                   (RealSubray diry yl, fhr)+                      -> let hhr = fmap (. gr') fhr+                         in case lapply slg 1 of+                              y' | y'>0 -> ( unsafePreRegionIntersect rg+                                                  $ RealSubray diry (cog + (yl-aog)/y')+                                   -- aog + y' * (xl − cog) = yl+                                   -- xl = cog + (yl − aog)/y'+                                           , hhr )+                                 | y'<0 -> ( unsafePreRegionIntersect rg+                                                  $ RealSubray (otherHalfSphere diry)+                                                               (cog + (yl-aog)/y')+                                           , hhr )+                                 | otherwise -> (rg, hhr)+                   (PreRegion ry, fhr)+                         -> ( PreRegion $ ry . gr', fmap (. gr') fhr )+           (GlobalRegion, Option (Just gr@(Differentiable grd)))+            -> let (y₀,_,_) = grd x₀+               in case f y₀ of+                   (GlobalRegion, Option Nothing)+                         -> (GlobalRegion, notDefinedHere)+                   (GlobalRegion, Option (Just fr))+                         -> (GlobalRegion, pure (fr . gr))+                   (r, Option Nothing) | PreRegion ry <- genericisePreRegion r+                         -> ( PreRegion $ ry . gr, notDefinedHere )+                   (r, Option (Just fr)) | PreRegion ry <- genericisePreRegion r+                         -> ( PreRegion $ ry . gr, pure (fr . gr) )+           (rg@(RealSubray _ _), Option (Just gr@(Differentiable grd)))+            -> let (y₀,_,_) = grd x₀+               in case f y₀ of+                   (GlobalRegion, Option Nothing)+                         -> (rg, notDefinedHere)+                   (GlobalRegion, Option (Just fr))+                         -> (rg, pure (fr . gr))+                   (rf, Option Nothing)+                     | PreRegion rx <- genericisePreRegion rg+                     , PreRegion ry <- genericisePreRegion rf+                         -> ( PreRegion $ minDblfuncs (ry . gr) rx+                            , notDefinedHere )+                   (rf, Option (Just fr))+                     | PreRegion rx <- genericisePreRegion rg+                     , PreRegion ry <- genericisePreRegion rf+                         -> ( PreRegion $ minDblfuncs (ry . gr) rx+                            , pure (fr . gr) )+           (PreRegion rx, Option (Just gr@(Differentiable grd)))+            -> let (y₀,_,_) = grd x₀+               in case f y₀ of+                   (GlobalRegion, Option Nothing)+                         -> (PreRegion rx, notDefinedHere)+                   (GlobalRegion, Option (Just fr))+                         -> (PreRegion rx, pure (fr . gr))+                   (r, Option Nothing) | PreRegion ry <- genericisePreRegion r+                         -> ( PreRegion $ minDblfuncs (ry . gr) rx+                            , notDefinedHere )+                   (r, Option (Just fr)) | PreRegion ry <- genericisePreRegion r+                         -> ( PreRegion $ minDblfuncs (ry . gr) rx+                            , pure (fr . gr) )+           (r, Option Nothing)+            -> (r, notDefinedHere)+             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@@ -880,20 +764,24 @@  data RWDfblFuncValue s d c where   ConstRWDFV :: c -> RWDfblFuncValue s d c+  RWDFV_IdVar :: RWDfblFuncValue s c 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 RWDFV_IdVar = GenericRWDFV id 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+  alg fq = case fq RWDFV_IdVar of     GenericRWDFV f -> f+    ConstRWDFV c -> const c+    RWDFV_IdVar -> id   ($~) = postCompRW instance RealDimension s => CartesianAgent (RWDiffable s) where-  alg1to2 fgq = case fgq (GenericRWDFV id) of+  alg1to2 fgq = case fgq RWDFV_IdVar of     (GenericRWDFV f, GenericRWDFV g) -> f &&& g   alg2to1 fq = case fq (GenericRWDFV fst) (GenericRWDFV snd) of     GenericRWDFV f -> f@@ -924,7 +812,7 @@         \d₀ -> let (rc', fmay) = fpcs d₀                    (rc'',gmay) = gpcs d₀                in (unsafePreRegionIntersect rc' rc'',) $-                    case (fmay,gmay) of+                    case (genericiseDifferentiable<$>fmay, genericiseDifferentiable<$>gmay) of                       (Option(Just(Differentiable f)), Option(Just(Differentiable g))) ->                         pure . Differentiable $ \d                          -> let (c', f', devf) = f d@@ -946,12 +834,54 @@ grwDfblFnValsCombine cmb fv gv         = grwDfblFnValsCombine cmb (genericiseRWDFV fv) (genericiseRWDFV gv) +          +rwDfbl_plus :: ∀ s a v .+        ( WithField s EuclidSpace v, AdditiveGroup v, v ~ Needle (Interior (Needle v))+        , LocallyScalable s a, RealDimension s )+      => RWDiffable s a v -> RWDiffable s a v -> RWDiffable s a v+rwDfbl_plus (RWDiffable f) (RWDiffable g) = RWDiffable h+   where h x₀ = (rh, liftA2 fgplus ff gf)+          where (rf, ff) = f x₀+                (rg, gf) = g x₀+                rh = unsafePreRegionIntersect rf rg+                fgplus :: Differentiable s a v -> Differentiable s a v -> Differentiable s a v+                fgplus (Differentiable fd) (Differentiable gd) = Differentiable hd+                 where hd x = (fx^+^gx, jf^+^jg, \ε -> δf(ε^*4) ^+^ δg(ε^*4))+                        where (fx, jf, δf) = fd x+                              (gx, jg, δg) = gd x+                fgplus (Differentiable fd) (AffinDiffable ga@(Affine cog aog slg))+                                 = Differentiable hd+                 where hd x = (fx^+^gx, jf^+^slg, δf)+                        where (fx, jf, δf) = fd x+                              gx = ga $ x+                fgplus (AffinDiffable fa@(Affine cof aof slf)) (Differentiable gd)+                                 = Differentiable hd+                 where hd x = (fx^+^gx, slf^+^jg, δg)+                        where (gx, jg, δg) = gd x+                              fx = fa $ x+                fgplus (AffinDiffable fa) (AffinDiffable ga) = AffinDiffable ha+                 where (GenericAgent ha) = GenericAgent fa ^+^ GenericAgent ga +rwDfbl_negateV :: ∀ s a v .+        ( WithField s EuclidSpace v, AdditiveGroup v, v ~ Needle (Interior (Needle v))+        , LocallyScalable s a, RealDimension s )+      => RWDiffable s a v -> RWDiffable s a v+rwDfbl_negateV (RWDiffable f) = RWDiffable h+   where h x₀ = (rf, fmap fneg ff)+          where (rf, ff) = f x₀+                fneg :: Differentiable s a v -> Differentiable s a v+                fneg (Differentiable fd) = Differentiable hd+                 where hd x = (negateV fx, negateV jf, δf)+                        where (fx, jf, δf) = fd x+                fneg (AffinDiffable (Affine cof aof slf))+                        = AffinDiffable $ Affine (negateV cof) (negateV aof) (negateV slf)+ 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 RWDFV_IdVar = GenericRWDFV f postCompRW f (GenericRWDFV g) = GenericRWDFV $ f . g  @@ -960,31 +890,68 @@     => AdditiveGroup (RWDfblFuncValue s a v) where   zeroV = point zeroV   ConstRWDFV c₁ ^+^ ConstRWDFV c₂ = ConstRWDFV (c₁^+^c₂)+  ConstRWDFV c₁ ^+^ RWDFV_IdVar = GenericRWDFV $+                               globalDiffable' (actuallyAffine c₁ idL)+  RWDFV_IdVar ^+^ ConstRWDFV c₂ = GenericRWDFV $+                               globalDiffable' (actuallyAffine c₂ idL)   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 (^+^)+  GenericRWDFV f ^+^ GenericRWDFV g = GenericRWDFV $ rwDfbl_plus f g   negateV (ConstRWDFV c) = ConstRWDFV (negateV c)-  negateV v = grwDfblFnValsFunc (\a -> (negateV a, lNegate, const zeroV)) v-      where lNegate = linear negateV+  negateV RWDFV_IdVar = GenericRWDFV $ globalDiffable' (actuallyLinear $ linear negateV)+  negateV (GenericRWDFV f) = GenericRWDFV $ rwDfbl_negateV f  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₁ * RWDFV_IdVar = GenericRWDFV $+                               globalDiffable' (actuallyLinear $ linear (c₁*))+  RWDFV_IdVar * ConstRWDFV c₂ = GenericRWDFV $+                               globalDiffable' (actuallyLinear $ linear (*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+  f*g = genericiseRWDFV f ⋅ genericiseRWDFV g+   where (⋅) :: ∀ n a . (RealDimension n, LocallyScalable n a)+           => RWDfblFuncValue n a n -> RWDfblFuncValue n a n -> RWDfblFuncValue n a n +         GenericRWDFV (RWDiffable fpcs) ⋅ GenericRWDFV (RWDiffable gpcs)+           = GenericRWDFV . RWDiffable $+               \d₀ -> let (rc₁, fmay) = fpcs d₀+                          (rc₂,gmay) = gpcs d₀+                      in (unsafePreRegionIntersect rc₁ rc₂, mulDi <$> fmay <*> gmay)+          where mulDi :: Differentiable n a n -> Differentiable n a n -> Differentiable n a n+                mulDi (AffinDiffable f@(Affine _ _ slf)) (AffinDiffable g@(Affine _ _ slg))+                   = let f' = lapply slf 1; g' = lapply slg 1+                     in case f'*g' of+                          0 -> AffinDiffable undefined+                          f'g' -> Differentiable $+                           \d -> let c₁ = f $ d; c₂ = g $ d+                                 in ( c₁*c₂+                                    , linear.(*)$ c₁*g' + c₂*f'+                                    , unsafe_dev_ε_δ "*" $ sqrt . (/f'g') )+                mulDi (Differentiable f) (Differentiable g)+                   = Differentiable $+                       \d -> let (c₁, slf, devf) = f d+                                 (c₂, slg, devg) = g d+                                 c = c₁*c₂; c₁² = c₁^2; c₂² = c₂^2+                                 h' = c₁*^slg ^+^ c₂*^slf+                                 in ( c+                                    , h'+                                    , \εc -> let rε² = metric εc 1+                                                 c₁worst² = c₁² + recip(1 + c₂²*rε²)+                                                 c₂worst² = c₂² + recip(1 + c₁²*rε²)+                                             in (4*rε²) *^ dualCoCoProduct slf slg+                                                ^+^ devf (εc^*(4*c₂worst²))+                                                ^+^ devg (εc^*(4*c₁worst²))+                    -- TODO: add formal proof for this (or, if necessary, the correct form)+                                        )+                mulDi f g = mulDi (genericiseDifferentiable f) (genericiseDifferentiable g)+                   negate = negateV   abs = (RWDiffable absPW $~)    where absPW a₀@@ -994,8 +961,8 @@          asc = actuallyLinear idL   signum = (RWDiffable sgnPW $~)    where sgnPW a₀-          | a₀<0       = (negativePreRegion, pure (const 1))-          | otherwise  = (positivePreRegion, pure (const $ -1))+          | a₀<0       = (negativePreRegion, pure (const $ -1))+          | otherwise  = (positivePreRegion, pure (const 1))  instance (RealDimension n, LocallyScalable n a)             => Fractional (RWDfblFuncValue n a n) where@@ -1009,35 +976,35 @@                  --           = -δ²/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)+          where δ ε = let mph = -ε*x^2/2+                          δ₀ = mph + sqrt (mph^2 - ε*x^3)+                      in if δ₀ > 0+                           then δ₀+                           else - x -- numerical underflow of εx³ vs mph+                                    --  ≡ ε*x^3 / (2*mph) (Taylor-expansion of the root)                 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)+          where δ ε = let mph = ε*x^2/2+                          δ₀ = sqrt (mph^2 + ε*x^3) - mph+                      in if δ₀>0 then δ₀ else x                 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) )+            in if ex*2 == ex  -- numerical trouble...+                then if x<0 then ( 0, zeroV, unsafe_dev_ε_δ("exp "++show x) $ \ε -> log ε - x )+                            else ( ex, ex*^idL, unsafe_dev_ε_δ("exp "++show x) $ \_ -> 1e-300 )+                else ( ex, ex *^ idL, unsafe_dev_ε_δ("exp "++show x)+                          $ \ε -> case acosh(ε/(2*ex) + 1) of+                                    δ | δ==δ      -> δ+                                      | otherwise -> log ε - x )                  -- ε = e^(x+δ) − eˣ − eˣ·δ                   --   = eˣ·(e^δ − 1 − δ)                   --   ≤ eˣ · (e^δ − 1 + e^(-δ) − 1)@@ -1139,14 +1106,14 @@   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))+                δ ε = abs x * sqrt((1 - exp(-ε))*0.8 + ε^2/(3*abs x + 1)) + 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) δ )+  acosh = postCompRW . RWDiffable $ \x -> if x>1+                                   then (preRegionToInfFrom 1, pure (Differentiable acoshDfb))+                                   else (preRegionFromMinInfTo 1, notDefinedHere)+   where acoshDfb x = ( acosh x, idL ^/ sqrt(x^2 - 1), 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@@ -1160,7 +1127,118 @@    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.   ++++-- $definitionRegionOps+-- Because the agents of 'RWDiffable' aren't really values in /Hask/, you can't use+-- the standard comparison operators on them, nor the built-in syntax of guards+-- or if-statements.+-- +-- However, because this category allows functions to be undefined in some region,+-- such decisions can be faked quite well: '?->' restricts a function to+-- some region, by simply marking it undefined outside¹, and '?|:' replaces these+-- regions with values from another function.+-- +-- Example: define a function that is compactly supported on the interval ]-1,1[,+-- i.e. exactly zero everywhere outside.+--+-- @+-- Graphics.Dynamic.Plot.R2> plotWindow [diffableFnPlot (\\x -> -1 '?<' x '?<' 1 '?->' exp(1/(x^2 - 1)) '?|:' 0)]+-- @+-- +-- <<images/examples/Friedrichs-mollifier.png>>+-- +-- ¹ Note that it may not be necessary to restrict explicitly: for instance if a+-- square root appears somewhere in an expression, then the expression is automatically+-- restricted so that the root has a positive argument!   +infixr 4 ?->+-- | Require the LHS to be defined before considering the RHS as result.+--   This works analogously to the standard `Control.Applicative.Applicative` method+-- +--   @+--   ('Control.Applicative.*>') :: Maybe a -> Maybe b -> Maybe b+--   Just _ 'Control.Applicative.*>' a = a+--   _      'Control.Applicative.*>' a = Nothing+--   @+(?->) :: (RealDimension n, LocallyScalable n a, LocallyScalable n b, LocallyScalable n c)+      => RWDfblFuncValue n c a -> RWDfblFuncValue n c b -> RWDfblFuncValue n c b+ConstRWDFV _ ?-> f = f+RWDFV_IdVar ?-> f = f+GenericRWDFV (RWDiffable r) ?-> ConstRWDFV c = GenericRWDFV (RWDiffable s)+ where s x₀ = case r x₀ of+                (rd, Option (Just q)) -> (rd, return $ const c)+                (rd, Option Nothing) -> (rd, empty)+GenericRWDFV (RWDiffable f) ?-> GenericRWDFV (RWDiffable g) = GenericRWDFV (RWDiffable h)+ where h x₀ = case f x₀ of+                (rf, Option (Just _)) | (rg, q) <- g x₀+                        -> (unsafePreRegionIntersect rf rg, q)+                (rf, Option Nothing) -> (rf, empty)+c ?-> f = c ?-> genericiseRWDFV f++positiveRegionalId :: RealDimension n => RWDiffable n n n+positiveRegionalId = RWDiffable $ \x₀ ->+       if x₀ > 0 then (positivePreRegion, pure . AffinDiffable $ id)+                 else (negativePreRegion, notDefinedHere)++infixl 5 ?> , ?<+-- | Return the RHS, if it is less than the LHS.+--   (Really the purpose is just to compare the values, but returning one of them+--   allows chaining of comparison operators like in Python.)+--   Note that less-than comparison is <http://www.paultaylor.eu/ASD/ equivalent>+--   to less-or-equal comparison, because there is no such thing as equality.+(?>) :: (RealDimension n, LocallyScalable n a)+           => RWDfblFuncValue n a n -> RWDfblFuncValue n a n -> RWDfblFuncValue n a n+a ?> b = (positiveRegionalId $~ a-b) ?-> b++-- | Return the RHS, if it is greater than the LHS.+(?<) :: (RealDimension n, LocallyScalable n a)+           => RWDfblFuncValue n a n -> RWDfblFuncValue n a n -> RWDfblFuncValue n a n+ConstRWDFV a ?< RWDFV_IdVar = GenericRWDFV . RWDiffable $+       \x₀ -> if a < x₀ then (preRegionToInfFrom a, pure . AffinDiffable $ id)+                        else (preRegionFromMinInfTo a, notDefinedHere)+RWDFV_IdVar ?< ConstRWDFV a = GenericRWDFV . RWDiffable $+       \x₀ -> if x₀ < a then (preRegionFromMinInfTo a, pure . AffinDiffable $ const a)+                        else (preRegionToInfFrom a, notDefinedHere)+a ?< b = (positiveRegionalId $~ b-a) ?-> b++infixl 3 ?|:+-- | Try the LHS, if it is undefined use the RHS. This works analogously to+--   the standard `Control.Applicative.Alternative` method+-- +--   @+--   ('Control.Applicative.<|>') :: Maybe a -> Maybe a -> Maybe a+--   Just x 'Control.Applicative.<|>' _ = Just x+--   _      'Control.Applicative.<|>' a = a+--   @+-- +--  Basically a weaker and agent-ised version of 'backupRegions'.+(?|:) :: (RealDimension n, LocallyScalable n a, LocallyScalable n b)+      => RWDfblFuncValue n a b -> RWDfblFuncValue n a b -> RWDfblFuncValue n a b+ConstRWDFV c ?|: _ = ConstRWDFV c+RWDFV_IdVar ?|: _ = RWDFV_IdVar+GenericRWDFV (RWDiffable f) ?|: ConstRWDFV c = GenericRWDFV (RWDiffable h)+ where h x₀ = case f x₀ of+                (rd, Option (Just q)) -> (rd, Option (Just q))+                (rd, Option Nothing) -> (rd, Option . Just $ const c)+GenericRWDFV (RWDiffable f) ?|: GenericRWDFV (RWDiffable g) = GenericRWDFV (RWDiffable h)+ where h x₀ = case f x₀ of+                (rf, Option (Just q)) -> (rf, pure q)+                (rf, Option Nothing) | (rg, q) <- g x₀+                        -> (unsafePreRegionIntersect rf rg, q)+c ?|: f = c ?|: genericiseRWDFV f++-- | Replace the regions in which the first function is undefined with values+--   from the second function.+backupRegions :: (RealDimension n, LocallyScalable n a, LocallyScalable n b)+      => RWDiffable n a b -> RWDiffable n a b -> RWDiffable n a b+backupRegions (RWDiffable f) (RWDiffable g) = RWDiffable h+ where h x₀ = case f x₀ of+                (rf, q@(Option (Just _))) -> (rf, q)+                (rf, Option Nothing) | (rg, q) <- g x₀+                        -> (unsafePreRegionIntersect rf rg, q)+   

@@ -0,0 +1,124 @@+{-# LANGUAGE TypeOperators, GADTs, FlexibleContexts #-}++module Data.Function.Differentiable.Data where+++import Data.Semigroup+import Data.Function.Affine+import Data.VectorSpace+import Data.LinearMap+import Data.LinearMap.HerMetric++import Data.Manifold.Types.Primitive+import Data.Manifold.PseudoAffine++++type LinDevPropag d c = Metric c -> Metric d+++-- | 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.)+data Differentiable s d c where+   Differentiable :: ( 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+                              ) )+                  -> Differentiable s d c+   AffinDiffable :: LinearManifold d+               => Affine s d d -> Differentiable s d d+                      -- This should ideally map between two general affine spaces,+                      -- but since the special case of affine functions is mostly relevant+                      -- to optimise the propagation of real intervals, we don't do that.++++++++-- | 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+  RealSubray :: RealDimension s => S⁰ -> s -> PreRegion s s+  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+++++++++-- | 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. 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+

@@ -30,11 +30,12 @@   , productMetric, productMetric'   , metricAsLength, metricFromLength, metric'AsLength   -- * Utility for metrics-  , transformMetric, transformMetric'+  , transformMetric, transformMetric', dualCoCoProduct   , dualiseMetric, dualiseMetric'   , recipMetric, recipMetric'   , eigenSpan, eigenSpan'   , eigenCoSpan, eigenCoSpan'+  , metriNormalise, metriNormalise'   , metriScale', metriScale   , adjoint   , extendMetric@@ -48,6 +49,8 @@   , FiniteDimensional(..)   -- * Misc   , Stiefel1(..)+  , linMapAsTensProd, linMapFromTensProd+  , covariance   ) where      @@ -204,6 +207,15 @@ toDualWith (HerMetric Nothing) = const zeroV toDualWith (HerMetric (Just m)) = fromPackedVector . HMat.app m . asPackedVector +-- | Divide a vector by its own norm, according to metric, i.e. normalise it+--   or &#x201c;project to the metric's boundary&#x201d;.+metriNormalise :: (HasMetric v, Floating (Scalar v)) => HerMetric v -> v -> v+metriNormalise m v = v ^/ metric m v++metriNormalise' :: (HasMetric v, Floating (Scalar v))+                 => HerMetric' v -> DualSpace v -> DualSpace v+metriNormalise' m v = v ^/ metric' m v+ -- | &#x201c;Anti-normalise&#x201d; a vector: /multiply/ with its own norm, according to metric. metriScale :: (HasMetric v, Floating (Scalar v)) => HerMetric v -> v -> v metriScale m v = metric m v *^ v@@ -228,16 +240,32 @@ transformMetric :: (HasMetric v, HasMetric w, Scalar v ~ Scalar w)            => (w :-* v) -> HerMetric v -> HerMetric w transformMetric _ (HerMetric Nothing) = HerMetric Nothing-transformMetric t (HerMetric (Just m)) = matrixMetric $ tmat HMat.<> m HMat.<> HMat.tr tmat+transformMetric t (HerMetric (Just m)) = matrixMetric $ HMat.tr tmat HMat.<> m HMat.<> tmat  where tmat = asPackedMatrix t  transformMetric' :: ( HasMetric v, HasMetric w, Scalar v ~ Scalar w )            => (v :-* w) -> HerMetric' v -> HerMetric' w transformMetric' _ (HerMetric' Nothing) = HerMetric' Nothing transformMetric' t (HerMetric' (Just m))-                      = matrixMetric' $ HMat.tr tmat HMat.<> m HMat.<> tmat+                      = matrixMetric' $ tmat HMat.<> m HMat.<> HMat.tr tmat  where tmat = asPackedMatrix t +-- | This does something vaguely like  @\\s t -> (s⋅t)²@,+--   but without actually requiring an inner product on the covectors.+--   Used for calculating the superaffine term of multiplications in+--   'Differentiable' categories.+dualCoCoProduct :: (HasMetric v, HasMetric w, Scalar v ~ Scalar w)+           => (w :-* v) -> (w :-* v) -> HerMetric w+dualCoCoProduct s t = ( (sArr `HMat.dot` (t²PLUSs² HMat.<\> sArr))+                       * (tArr `HMat.dot` (t²PLUSs² HMat.<\> tArr)) )+                    *^ matrixMetric t²PLUSs²+ where tmat = asPackedMatrix t+       tArr = HMat.flatten tmat+       smat = asPackedMatrix s+       sArr = HMat.flatten smat+       t²PLUSs² = tmat HMat.<> HMat.tr tmat + smat HMat.<> HMat.tr smat++ -- | This doesn't really do anything at all, since @'HerMetric' v@ is essentially a --   synonym for @'HerMetric' ('DualSpace' v)@. dualiseMetric :: HasMetric v => HerMetric (DualSpace v) -> HerMetric' v@@ -266,6 +294,7 @@   isInfinite' :: (Eq a, Num a) => a -> Bool+isInfinite' 0 = False isInfinite' x = x==x*2  @@ -522,6 +551,21 @@         = HerMetric' . Just $ HMat.diagBlock [mv, HMat.konst 0 (dw,dw)]  where (Tagged dw) = dimension :: Tagged w Int ++++covariance :: ∀ v w . (HasMetric v, HasMetric w, Scalar v ~ ℝ, Scalar w ~ ℝ)+          => HerMetric' (v,w) -> Option (v:-*w)+covariance (HerMetric' Nothing) = pure zeroV+covariance (HerMetric' (Just m))+    | isInfinite' detvnm  = empty+    | otherwise           = pure . fromPackedMatrix $+                               wmat HMat.<> m HMat.<> vmat HMat.<> vnorml+ where wmat = asPackedMatrix (linear snd :: (v,w):-*w)+       vmat = asPackedMatrix (linear (id&&&const zeroV) :: v:-*(v,w))+       (vnorml, (detvnm, _)) = HMat.invlndet (HMat.tr vmat HMat.<> m HMat.<> vmat)++ metricAsLength :: HerMetric ℝ -> ℝ metricAsLength m = case metricSq m 1 of    o | o > 0    -> recip o@@ -595,3 +639,20 @@                       . foldr1 ((.) . (.(" ^+^ "++)))                       $ ((("projector' "++).).showsPrec 6)<$>eigSp    where eigSp = eigenSpan m++++++++++linMapAsTensProd :: (FiniteDimensional v, FiniteDimensional w, Scalar v~Scalar w)+                    => v:-*w -> DualSpace v ⊗ w+linMapAsTensProd f = DensTensProd $ asPackedMatrix f++linMapFromTensProd :: (FiniteDimensional v, FiniteDimensional w, Scalar v~Scalar w)+                    => DualSpace v ⊗ w -> v:-*w+linMapFromTensProd (DensTensProd m) = linear $+                         asPackedVector >>> HMat.app m >>> fromPackedVector

@@ -48,7 +48,7 @@ module Data.Manifold.PseudoAffine (             -- * Manifold class               Manifold-            , Semimanifold(..)+            , Semimanifold(..), Needle'             , PseudoAffine(..)             -- * Type definitions             -- ** Metrics@@ -225,8 +225,11 @@ --    --   (Actually, 'LinearManifold' is stronger than 'VectorSpace' at the moment, since --   'HasMetric' requires 'FiniteDimensional'. This might be lifted in the future.)-type LinearManifold x = ( PseudoAffine x, Interior x ~ x, Needle x ~ x, HasMetric x )+type LinearManifold x = ( AffineManifold x, Needle x ~ x, HasMetric x ) +type LinearManifold' x = ( PseudoAffine x, AffineSpace x, Diff x ~ x+                         , Interior x ~ x, Needle x ~ x, HasMetric x )+ -- | Require some constraint on a manifold, and also fix the type of the manifold's --   underlying field. For example, @WithField ℝ 'HilbertSpace' v@ constrains --   @v@ to be a real (i.e., 'Double'-) Hilbert space.@@ -243,7 +246,7 @@  -- | The 'AffineSpace' class plus manifold constraints. type AffineManifold m = ( PseudoAffine m, Interior m ~ m, AffineSpace m-                        , Needle m ~ Diff m, LinearManifold (Diff m) )+                        , Needle m ~ Diff m, LinearManifold' (Diff m) )  -- | A Hilbert space is a /complete/ inner product space. Being a vector space, it is --   also a manifold.@@ -263,6 +266,12 @@ euclideanMetric = Tagged euclideanMetric'  +-- | A co-needle can be understood as a “paper stack”, with which you can measure+--   the length that a needle reaches in a given direction by counting the number+--   of holes punched through them.+type Needle' x = DualSpace (Needle x)++ -- | The word “metric” is used in the sense as in general relativity. Cf. 'HerMetric'. type Metric x = HerMetric (Needle x) type Metric' x = HerMetric' (Needle x)@@ -374,6 +383,23 @@ instance (MetricScalar a, KnownNat n) => PseudoAffine (FreeVect n a) where   a.-~.b = pure (a.-.b) +instance (HasMetric a, FiniteDimensional b, Scalar a~Scalar b) => Semimanifold (a⊗b) where+  type Needle (a⊗b) = a ⊗ b+  fromInterior = id+  toInterior = pure+  translateP = Tagged (.+~^)+  (.+~^) = (^+^)+instance (HasMetric a, FiniteDimensional b, Scalar a~Scalar b) => PseudoAffine (a⊗b) where+  a.-~.b = pure (a^-^b)++instance (HasMetric a, FiniteDimensional b, Scalar a~Scalar b) => Semimanifold (a:-*b) where+  type Needle (a:-*b) = DualSpace a ⊗ b+  fromInterior = id+  toInterior = pure+  translateP = Tagged (.+~^)+  p.+~^n = p ^+^ linMapFromTensProd n+instance (HasMetric a, FiniteDimensional b, Scalar a~Scalar b) => PseudoAffine (a:-*b) where+  a.-~.b = pure . linMapAsTensProd $ a^-^b  instance Semimanifold S⁰ where   type Needle S⁰ = ℝ⁰

@@ -77,6 +77,8 @@ import Data.Manifold.Types import Data.Manifold.Types.Primitive ((^), empty) import Data.Manifold.PseudoAffine+import Data.Function.Differentiable+import Data.Function.Differentiable.Data      import Data.Embedding import Data.CoNat@@ -107,7 +109,7 @@ -- | Possibly / Partially / asymPtotically singular metric. data PSM x = PSM {        psmExpanse :: !(Metric' x)-     , relevantEigenspan :: ![DualSpace (Needle x)]+     , relevantEigenspan :: ![Needle' x]      }         @@ -161,7 +163,7 @@ fullShade' ctr expa = Shade' ctr expa  subshadeId' :: WithField ℝ Manifold x-                   => x -> NonEmpty (DualSpace (Needle x)) -> x -> (Int, HourglassBulb)+                   => x -> NonEmpty (Needle' x) -> x -> (Int, HourglassBulb) subshadeId' c expvs x = case x .-~. c of     Option (Just v) -> let (iu,vl) = maximumBy (comparing $ abs . snd)                                       $ zip [0..] (map (v <.>^) $ NE.toList expvs)@@ -291,7 +293,7 @@                  | OverlappingBranches !Int !(Shade x) (NonEmpty (DBranch x))   deriving (Generic)            -data DBranch' x c = DBranch { boughDirection :: !(DualSpace (Needle x))+data DBranch' x c = DBranch { boughDirection :: !(Needle' x)                             , boughContents :: !(Hourglass c) }   deriving (Generic, Hask.Functor, Hask.Foldable) type DBranch x = DBranch' x (ShadeTree x)@@ -306,11 +308,11 @@     -instance (NFData x, NFData (DualSpace (Needle x))) => NFData (ShadeTree x) where+instance (NFData x, NFData (Needle' x)) => NFData (ShadeTree x) where   rnf (PlainLeaves xs) = rnf xs   rnf (DisjointBranches n bs) = n `seq` rnf (NE.toList bs)   rnf (OverlappingBranches n sh bs) = n `seq` sh `seq` rnf (NE.toList bs)-instance (NFData x, NFData (DualSpace (Needle x))) => NFData (DBranch x)+instance (NFData x, NFData (Needle' x)) => NFData (DBranch x)    -- | Experimental. There might be a more powerful instance possible. instance (AffineManifold x) => Semimanifold (ShadeTree x) where@@ -361,12 +363,28 @@ -- @ --  -- <<images/examples/simple-2d-ShadeTree.png>>-fromLeafPoints :: forall x. WithField ℝ Manifold x => [x] -> ShadeTree x-fromLeafPoints = go zeroV+fromLeafPoints :: ∀ x. WithField ℝ Manifold x => [x] -> ShadeTree x+fromLeafPoints = fromLeafPoints' sShIdPartition++++fromFnGraphPoints :: ∀ x y . (WithField ℝ Manifold x, WithField ℝ Manifold y)+                     => [(x,y)] -> ShadeTree (x,y)+fromFnGraphPoints = fromLeafPoints' fg_sShIdPart+ where fg_sShIdPart :: Shade (x,y) -> [(x,y)] -> NonEmpty (DBranch' (x,y) [(x,y)])+       fg_sShIdPart (Shade c expa) xs+        | b:bs <- [DBranch (v, zeroV) mempty+                    | v <- eigenCoSpan+                           (transformMetric' (linear fst) expa :: Metric' x) ]+                      = sShIdPartition' c xs $ b:|bs++fromLeafPoints' :: ∀ x. WithField ℝ Manifold x =>+    (Shade x -> [x] -> NonEmpty (DBranch' x [x])) -> [x] -> ShadeTree x+fromLeafPoints' sShIdPart = go zeroV  where go :: Metric' x -> [x] -> ShadeTree x        go preShExpa = \xs -> case pointsShades' (preShExpa^/10) xs of                      [] -> mempty-                     [(_,rShade)] -> let trials = sShIdPartition rShade xs+                     [(_,rShade)] -> let trials = sShIdPart rShade xs                                      in case reduce rShade trials of                                          Just redBrchs                                            -> OverlappingBranches@@ -444,6 +462,73 @@   ++intersectShade's :: ∀ y . WithField ℝ Manifold y => [Shade' y] -> Option (Shade' y)+intersectShade's [] = error "Global `Shade'` not implemented, so can't do intersection of zero co-shades."+intersectShade's (sh:shs) = Hask.foldrM inter2 sh shs+ where inter2 :: Shade' y -> Shade' y -> Option (Shade' y)+       inter2 (Shade' c e) (Shade' ζ η)+           | μc > 1 && μζ > 1  = empty+           | otherwise         = return $ Shade' (c.+~^w) (e^+^η)+        where Option (Just c2ζ) = ζ.-~.c+              Option (Just ζ2c) = c.-~.ζ+              ζNearest, cNearest :: y+              ζNearest = c .+~^ metriNormalise e c2ζ+              cNearest = ζ .+~^ metriNormalise η ζ2c+              Option (Just rζ) = ζNearest.-~.ζ+              Option (Just rc) = cNearest.-~.c+              μc = metric e rc+              μζ = metric η rζ+              w = c2ζ ^* (μζ/(μc + μζ))+              -- = (c^*μc + ζ^*μζ)/(μc + μζ) − c+              -- = (c^*μc + ζ^*μζ − c^*(μc+μζ))^/(μc + μζ)+              -- = (ζ^*μζ − c^*μζ)^/(μc + μζ)+              -- = (ζ−c)^*μζ/(μc + μζ)+++++type DifferentialEqn x y = RWDiffable ℝ (x,y) (Needle x :-* Needle y)+++filterDEqnSolution_loc :: ∀ x y . (WithField ℝ Manifold x, WithField ℝ Manifold y)+           => DifferentialEqn x y -> (Shade' (x,y), [Shade' (x,y)]) -> [Shade' (x,y)]+filterDEqnSolution_loc (RWDiffable f) (Shade' (x,y) expa, neighbours) = case f (x,y) of+          (_, Option Nothing) -> []+          (r, Option (Just (Differentiable fl)))+                | (fc, fc', δ) <- fl (x,y)+                   -> let flatMet :: HerMetric (Needle (x,y))+                          flatMet = recipMetric -- this won't work, metric is singular.+                               . transformMetric' (linear $ id &&& lapply fc) +                               $ recipMetric' expax+                          -- fcs = lapply fc' <$> xSpan+                          -- flinRange = δ $ projectors fcs+                          marginδs :: [(Needle x, (Needle y, Metric y))]+                          marginδs = [ (δxm, (δym, expany))+                                     | Shade' (xn, yn) expan <- neighbours+                                     , let (Option (Just δx)) = x.-~.xn+                                           (expanx, expany) = factoriseMetric expan+                                           (Option (Just yc'n))+                                                  = covariance $ recipMetric' expan+                                           xntoMarg = metriNormalise expanx δx+                                           (Option (Just δxm))+                                              = (xn .+~^ xntoMarg :: x) .-~. x+                                           (Option (Just δym))+                                              = (yn .+~^ lapply yc'n xntoMarg :: y+                                                  ) .-~. y+                                     ]+                          ycQuad :: y+                          (Option (Just (Shade' ycQuad _))) = intersectShade's+                                     [ Shade' ycn expany+                                     | (δxm,(δym,expany)) <- marginδs+                                     , let fca :: Needle x:-*Needle y+                                           fca = fc .+~^ lapply fc' ((δxm,δym)^/2)+                                           ycn = y .+~^ (δym ^-^ lapply fca δxm)+                                     ]+                                     :: Option (Shade' y)+                      in [Shade' (x,ycQuad) flatMet]+ where (expax, expay) = factoriseMetric expa+       xSpan = eigenCoSpan' expax       

@@ -40,13 +40,15 @@         , ZeroDim(..), isoAttachZeroDim         , ℝ⁰, ℝ, ℝ², ℝ³         -- * Hyperspheres-        , S⁰(..), S¹(..), S²(..)+        , S⁰(..), otherHalfSphere, S¹(..), S²(..)         -- * Projective spaces         , ℝP¹,  ℝP²(..)         -- * Intervals\/disks\/cones         , D¹(..), D²(..)         , ℝay         , CD¹(..), Cℝay(..)+        -- * Tensor products+        , (⊗)(..)         -- * Utility (deprecated)         , NaturallyEmbedded(..)         , GraphWindowSpec(..), Endomorphism, (^), (^.), EqFloating@@ -61,6 +63,8 @@ import Data.Void import Data.Monoid +import qualified Numeric.LinearAlgebra.HMatrix as HMat+ import Control.Applicative (Const(..), Alternative(..))  import qualified Prelude@@ -93,6 +97,10 @@ instance Monoid (ZeroDim k) where   mempty = Origin   mappend Origin Origin = Origin+instance AffineSpace (ZeroDim k) where+  type Diff (ZeroDim k) = ZeroDim k+  Origin .+^ Origin = Origin+  Origin .-. Origin = Origin instance AdditiveGroup (ZeroDim k) where   zeroV = Origin   Origin ^+^ Origin = Origin@@ -117,6 +125,11 @@ --   therefore change to @ℝ⁰ 'Control.Category.Constrained.+' ℝ⁰@: the disjoint sum of two --   single-point spaces. data S⁰ = PositiveHalfSphere | NegativeHalfSphere deriving(Eq, Show)++otherHalfSphere :: S⁰ -> S⁰+otherHalfSphere PositiveHalfSphere = NegativeHalfSphere+otherHalfSphere NegativeHalfSphere = PositiveHalfSphere+ -- | The unit circle. newtype S¹ = S¹ { φParamS¹ :: Double -- ^ Must be in range @[-π, π[@.                 } deriving (Show)@@ -168,6 +181,14 @@ data Cℝay x = Cℝay { hParamCℝay :: !Double -- ^ Range @[0, ∞[@                    , pParamCℝay :: !x      -- ^ Irrelevant at @h = 0@.                    }+++++-- | Dense tensor product of two vector spaces.+newtype x⊗y = DensTensProd { getDensTensProd :: HMat.Matrix (Scalar y) }++  class NaturallyEmbedded m v where   embed :: m -> v

@@ -27,7 +27,7 @@ module Data.VectorSpace.FiniteDimensional (     FiniteDimensional(..)   , SmoothScalar -  , FinVecArrRep(..), concreteArrRep, (⊗), splitArrRep+  , FinVecArrRep(..), concreteArrRep, (⊕), splitArrRep   ) where      @@ -95,7 +95,7 @@    where defaultAsPackedMatrix :: forall v w s .                (FiniteDimensional v, FiniteDimensional w, s~Scalar v, s~Scalar w)                          => (v :-* w) -> HMat.Matrix s-         defaultAsPackedMatrix m = HMat.fromRows $ asPackedVector . atBasis m <$> cb+         defaultAsPackedMatrix m = HMat.fromColumns $ asPackedVector . atBasis m <$> cb           where (Tagged cb) = completeBasis :: Tagged v [Basis v]      fromPackedVector :: HMat.Vector (Scalar v) -> v@@ -103,6 +103,14 @@    where result = recompose $ zip cb (HMat.toList v)          cb = witness completeBasis result +  fromPackedMatrix :: (FiniteDimensional w, Scalar w ~ Scalar v)+                       => HMat.Matrix (Scalar v) -> (v :-* w)+  fromPackedMatrix = defaultFromPackedMatrix+   where defaultFromPackedMatrix :: forall v w s .+               (FiniteDimensional v, FiniteDimensional w, s~Scalar v, s~Scalar w)+                         => HMat.Matrix s -> (v :-* w)+         defaultFromPackedMatrix m = linear $ fromPackedVector . HMat.app m . asPackedVector+   instance (SmoothScalar k) => FiniteDimensional (ZeroDim k) where   dimension = Tagged 0   basisIndex = Tagged absurd@@ -116,16 +124,16 @@   indexBasis = Tagged $ \0 -> ()   completeBasis = Tagged [()]   asPackedVector x = HMat.fromList [x]-  asPackedMatrix f = HMat.asRow . asPackedVector $ atBasis f ()+  asPackedMatrix f = HMat.asColumn . asPackedVector $ atBasis f ()   fromPackedVector v = v HMat.! 0 instance (FiniteDimensional a, FiniteDimensional b, Scalar a~Scalar b)             => FiniteDimensional (a,b) where   dimension = tupDim-   where tupDim :: forall a b.(FiniteDimensional a,FiniteDimensional b)=>Tagged(a,b)Int+   where tupDim :: ∀ a b.(FiniteDimensional a,FiniteDimensional b)=>Tagged(a,b)Int          tupDim = Tagged $ da+db           where (Tagged da)=dimension::Tagged a Int; (Tagged db)=dimension::Tagged b Int   basisIndex = basId-   where basId :: forall a b . (FiniteDimensional a, FiniteDimensional b)+   where basId :: ∀ a b . (FiniteDimensional a, FiniteDimensional b)                      => Tagged (a,b) (Either (Basis a) (Basis b) -> Int)          basId = Tagged basId'           where basId' (Left ba) = basIda ba@@ -134,7 +142,7 @@                 (Tagged basIda) = basisIndex :: Tagged a (Basis a->Int)                 (Tagged basIdb) = basisIndex :: Tagged b (Basis b->Int)   indexBasis = basId-   where basId :: forall a b . (FiniteDimensional a, FiniteDimensional b)+   where basId :: ∀ a b . (FiniteDimensional a, FiniteDimensional b)                      => Tagged (a,b) (Int -> Either (Basis a) (Basis b))          basId = Tagged basId'           where basId' i | i < da     = Left $ basIda i@@ -143,14 +151,14 @@                 (Tagged basIda) = indexBasis :: Tagged a (Int->Basis a)                 (Tagged basIdb) = indexBasis :: Tagged b (Int->Basis b)   completeBasis = cb-   where cb :: forall a b . (FiniteDimensional a, FiniteDimensional b)+   where cb :: ∀ a b . (FiniteDimensional a, FiniteDimensional b)                      => Tagged (a,b) [Either (Basis a) (Basis b)]          cb = Tagged $ map Left cba ++ map Right cbb           where (Tagged cba) = completeBasis :: Tagged a [Basis a]                 (Tagged cbb) = completeBasis :: Tagged b [Basis b]   asPackedVector (a,b) = HMat.vjoin [asPackedVector a, asPackedVector b]   fromPackedVector = fPV-   where fPV :: forall a b . (FiniteDimensional a, FiniteDimensional b, Scalar a~Scalar b)+   where fPV :: ∀ a b . (FiniteDimensional a, FiniteDimensional b, Scalar a~Scalar b)                      => HMat.Vector (Scalar a) -> (a,b)          fPV v = (fromPackedVector l, fromPackedVector r)           where (Tagged da) = dimension :: Tagged a Int@@ -158,7 +166,87 @@                 l = HMat.subVector 0 da v                 r = HMat.subVector da db v               +instance (FiniteDimensional y, FiniteDimensional x) => AdditiveGroup (x⊗y) where+  zeroV = DensTensProd $ (0 HMat.>< 0) []+  negateV (DensTensProd v) = DensTensProd $ negate v+  DensTensProd v ^+^ DensTensProd w+   | HMat.size v == (0,0)  = DensTensProd w+   | HMat.size w == (0,0)  = DensTensProd v+   | otherwise             = DensTensProd $ v + w++instance (FiniteDimensional y, FiniteDimensional x) => VectorSpace (x⊗y) where+  type Scalar (x⊗y) = Scalar y+  μ *^ DensTensProd v = DensTensProd $ HMat.scale μ v++instance (FiniteDimensional y, FiniteDimensional x) => InnerSpace (x⊗y) where+  DensTensProd v <.> DensTensProd w+   | HMat.size v == (0,0)  = 0+   | HMat.size w == (0,0)  = 0+   | otherwise             = HMat.flatten v `HMat.dot` HMat.flatten w++instance (FiniteDimensional y, FiniteDimensional x) => HasBasis (x⊗y) where+  type Basis (x⊗y) = (Basis x, Basis y)+  basisValue = bvt+   where bvt :: ∀ x y . (FiniteDimensional x, FiniteDimensional y)+                       => (Basis x, Basis y) -> x ⊗ y+         bvt (bx,by) = DensTensProd $ HMat.assoc (nx,ny) 0 [((i,j),1)]+          where Tagged nx = dimension :: Tagged x Int+                Tagged ny = dimension :: Tagged y Int+                Tagged i = ($bx) <$> basisIndex :: Tagged x Int+                Tagged j = ($by) <$> basisIndex :: Tagged y Int+  decompose = dct+   where dct :: ∀ x y . (FiniteDimensional x, FiniteDimensional y)+                       => x ⊗ y -> [((Basis x, Basis y), Scalar y)]+         dct (DensTensProd m) = zip [(i,j) | i <- cbx, j <- cby]+                                (HMat.toList $ HMat.flatten m)+          where Tagged cbx = completeBasis :: Tagged x [Basis x]+                Tagged cby = completeBasis :: Tagged y [Basis y]+  decompose' = dct+   where dct :: ∀ x y . (FiniteDimensional x, FiniteDimensional y)+                       => x ⊗ y -> (Basis x, Basis y) -> Scalar y+         dct (DensTensProd m) (bi, bj) = m `HMat.atIndex` (bxi bi, byj bj)+          where Tagged bxi = basisIndex :: Tagged x (Basis x -> Int)+                Tagged byj = basisIndex :: Tagged y (Basis y -> Int)+               +instance (FiniteDimensional a, FiniteDimensional b, Scalar a ~ Scalar b)+                                     => FiniteDimensional (a⊗b) where+  dimension = tensDim+   where tensDim :: ∀ a b.(FiniteDimensional a,FiniteDimensional b)=>Tagged(a⊗b)Int+         tensDim = Tagged $ da*db+          where (Tagged da)=dimension::Tagged a Int; (Tagged db)=dimension::Tagged b Int+  basisIndex = basId+   where basId :: ∀ a b . (FiniteDimensional a, FiniteDimensional b)+                     => Tagged (a⊗b) ((Basis a, Basis b) -> Int)+         basId = Tagged basId'+          where basId' (ba,bb) = db*basIda ba + basIdb bb+                (Tagged db) = dimension :: Tagged b Int+                (Tagged basIda) = basisIndex :: Tagged a (Basis a->Int)+                (Tagged basIdb) = basisIndex :: Tagged b (Basis b->Int)+  indexBasis = basId+   where basId :: ∀ a b . (FiniteDimensional a, FiniteDimensional b)+                     => Tagged (a⊗b) (Int -> (Basis a, Basis b))+         basId = Tagged basId'+          where basId' i = let (ia,ib) = i`divMod`db+                           in (basIda ia, basIdb ib)+                (Tagged db) = dimension :: Tagged b Int+                (Tagged basIda) = indexBasis :: Tagged a (Int->Basis a)+                (Tagged basIdb) = indexBasis :: Tagged b (Int->Basis b)+  completeBasis = cb+   where cb :: ∀ a b . (FiniteDimensional a, FiniteDimensional b)+                     => Tagged (a⊗b) [(Basis a, Basis b)]+         cb = Tagged $ [(ba,bb) | ba<-cba, bb<-cbb]+          where (Tagged cba) = completeBasis :: Tagged a [Basis a]+                (Tagged cbb) = completeBasis :: Tagged b [Basis b]+  asPackedVector (DensTensProd m) = HMat.flatten m+  fromPackedVector = fPV+   where fPV :: ∀ a b . (FiniteDimensional a, FiniteDimensional b, Scalar a~Scalar b)+                     => HMat.Vector (Scalar a) -> (a⊗b)+         fPV v = DensTensProd $ HMat.reshape db v+          where (Tagged db) = dimension :: Tagged b Int   +++   instance (SmoothScalar x, KnownNat n) => FiniteDimensional (FreeVect n x) where   dimension = natTagPænultimate   basisIndex = Tagged getInRange@@ -187,7 +275,7 @@   negateV (FinVecArrRep v) = FinVecArrRep $ negate v   FinVecArrRep v ^+^ FinVecArrRep w    | HMat.size v == 0  = FinVecArrRep w-   | HMat.size w == 0  = FinVecArrRep w+   | HMat.size w == 0  = FinVecArrRep v    | otherwise         = FinVecArrRep $ v + w  instance (SmoothScalar s) => VectorSpace (FinVecArrRep t b s) where@@ -205,10 +293,10 @@ concreteArrRep = Isomorphism (FinVecArrRep     . asPackedVector)                              (fromPackedVector . getFinVecArrRep) -(⊗) :: ∀ t s v w . ( SmoothScalar s, FiniteDimensional v, FiniteDimensional w+(⊕) :: ∀ t s v w . ( SmoothScalar s, FiniteDimensional v, FiniteDimensional w                    , Scalar v ~ s, Scalar w ~ s )           => FinVecArrRep t v s -> FinVecArrRep t w s -> FinVecArrRep t (v,w) s-FinVecArrRep v ⊗ FinVecArrRep w+FinVecArrRep v ⊕ FinVecArrRep w   | HMat.size v + HMat.size w == 0  = FinVecArrRep v   | HMat.size v == 0                = FinVecArrRep $ HMat.vjoin [HMat.konst 0 nv, w]   | HMat.size w == 0                = FinVecArrRep $ HMat.vjoin [v, HMat.konst 0 nw]

binary file changed (absent → 6282 bytes)

@@ -1,5 +1,5 @@ Name:                manifolds-Version:             0.1.6.2+Version:             0.1.6.3 Category:            Math Synopsis:            Coordinate-free hypersurfaces Description:         Manifolds, a generalisation of the notion of “smooth curves” or surfaces,@@ -75,6 +75,8 @@                    Data.CoNat                    Data.Embedding                    Data.LinearMap.Category+                   Data.Function.Differentiable.Data+                   Data.Function.Affine                    Data.VectorSpace.FiniteDimensional                    Util.Associate                    Util.LtdShow