diff --git a/Data/Function/Affine.hs b/Data/Function/Affine.hs new file mode 100644 --- /dev/null +++ b/Data/Function/Affine.hs @@ -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) + + + diff --git a/Data/Function/Differentiable.hs b/Data/Function/Differentiable.hs --- a/Data/Function/Differentiable.hs +++ b/Data/Function/Differentiable.hs @@ -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) -> llimL ) +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 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 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/₀ derivative /f'ₓ/₀, --- if, for¹ /ε/>0, there exists /δ/ such that --- |/f/ /x/ − (/f/ /x/₀ + /x/⋅/f'ₓ/₀)| < /ε/ --- for all |/x/ − /x/₀| < /δ/. --- --- Observe that, though this looks quite similar to the standard definition --- of differentiability, it is not equivalent thereto – in fact it does --- not prove any analytic properties at all. To make it equivalent, we need --- a lower bound on /δ/: simply /δ/ gives us continuity, and for --- continuous differentiability, /δ/ must grow at least like √/ε/ --- for small /ε/. Neither of these conditions are enforced by the type system, --- but we do require them for any allowed values because these proofs are obviously --- tremendously useful – for instance, you can have a root-finding algorithm --- and actually be sure you get /all/ solutions correctly, not just /some/ that are --- (hopefully) the closest to some reference point you'd need to laborously define! --- --- Unfortunately however, this also prevents doing any serious algebra etc. with the --- category, because even something as simple as division necessary introduces singularities --- where the derivatives must diverge. --- Not to speak of many trigonometric e.g. trigonometric functions that --- are undefined on whole regions. The 'PWDiffable' and 'RWDiffable' categories have explicit --- handling for those issues built in; you may simply use these categories even when --- you know the result will be smooth in your relevant domain (or must be, for e.g. --- physics reasons). --- --- ¹(The implementation does not deal with /ε/ and /δ/ as difference-bounding --- reals, but rather as metric tensors that define a boundary by prohibiting the --- overlap from exceeding one; this makes the concept actually work on general manifolds.) -newtype Differentiable s d c - = Differentiable { runDifferentiable :: - d -> ( c -- function value - , Needle d :-* Needle c -- Jacobian - , LinDevPropag d c -- Metric showing how far you can go - -- from x₀ without deviating from the - -- Taylor-1 approximation by more than - -- some error margin - ) } -type (-->) = Differentiable ℝ +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 “workhorse” for most serious --- calculus applications, because it contains all the usual trig etc. functions --- and of course everything algebraic you can do in the reals. --- --- The easiest way to define ordinary functions in this category is hence --- with its 'AgentVal'ues, which have instances of the standard classes 'Num' --- through 'Floating'. For instance, the following defines the /binary entropy/ --- as a differentiable function on the interval @]0,1[@: (it will --- actually /know/ where it's defined and where not! – and I don't mean you --- need to exhaustively 'isNaN'-check all results...) --- --- @ --- hb :: RWDiffable ℝ ℝ ℝ --- hb = alg (\\p -> - p * logBase 2 p - (1-p) * logBase 2 (1-p) ) --- @ -newtype RWDiffable s d c - = RWDiffable { - tryDfblDomain :: d -> (PreRegion s d, Option (Differentiable s d c)) } - -notDefinedHere :: Option (Differentiable s d c) -notDefinedHere = Option Nothing - - - instance (RealDimension s) => Category (RWDiffable s) where type Object (RWDiffable s) o = LocallyScalable s o id = RWDiffable $ \x -> (GlobalRegion, pure id) - RWDiffable f . RWDiffable g = RWDiffable h - where h x₀ = case g x₀ of - (GlobalRegion, Option Nothing) - -> (GlobalRegion, notDefinedHere) - (GlobalRegion, Option (Just gr)) - -> let (y₀,_,_) = runDifferentiable gr x₀ - in case f y₀ of - (GlobalRegion, Option Nothing) - -> (GlobalRegion, notDefinedHere) - (GlobalRegion, Option (Just fr)) - -> (GlobalRegion, pure (fr . gr)) - (PreRegion ry, Option Nothing) - -> ( PreRegion $ ry . gr, notDefinedHere ) - (PreRegion ry, Option (Just fr)) - -> ( PreRegion $ ry . gr, pure (fr . gr) ) - (PreRegion rx, Option Nothing) - -> (PreRegion rx, notDefinedHere) - (PreRegion rx, Option (Just gr)) - -> let (y₀,_,_) = runDifferentiable gr x₀ - in case f y₀ of - (GlobalRegion, Option Nothing) - -> (PreRegion rx, notDefinedHere) - (GlobalRegion, Option (Just fr)) - -> (PreRegion rx, pure (fr . gr)) - (PreRegion ry, Option Nothing) - -> ( PreRegion $ minDblfuncs (ry . gr) rx - , notDefinedHere ) - (PreRegion ry, Option (Just fr)) - -> ( PreRegion $ minDblfuncs (ry . gr) rx - , pure (fr . gr) ) + 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)] +-- @ +-- +-- <> +-- +-- ¹ 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 +-- 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) + diff --git a/Data/Function/Differentiable/Data.hs b/Data/Function/Differentiable/Data.hs new file mode 100644 --- /dev/null +++ b/Data/Function/Differentiable/Data.hs @@ -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/₀ derivative /f'ₓ/₀, +-- if, for¹ /ε/>0, there exists /δ/ such that +-- |/f/ /x/ − (/f/ /x/₀ + /x/⋅/f'ₓ/₀)| < /ε/ +-- for all |/x/ − /x/₀| < /δ/. +-- +-- Observe that, though this looks quite similar to the standard definition +-- of differentiability, it is not equivalent thereto – in fact it does +-- not prove any analytic properties at all. To make it equivalent, we need +-- a lower bound on /δ/: simply /δ/ gives us continuity, and for +-- continuous differentiability, /δ/ must grow at least like √/ε/ +-- for small /ε/. Neither of these conditions are enforced by the type system, +-- but we do require them for any allowed values because these proofs are obviously +-- tremendously useful – for instance, you can have a root-finding algorithm +-- and actually be sure you get /all/ solutions correctly, not just /some/ that are +-- (hopefully) the closest to some reference point you'd need to laborously define! +-- +-- Unfortunately however, this also prevents doing any serious algebra etc. with the +-- category, because even something as simple as division necessary introduces singularities +-- where the derivatives must diverge. +-- Not to speak of many trigonometric e.g. trigonometric functions that +-- are undefined on whole regions. The 'PWDiffable' and 'RWDiffable' categories have explicit +-- handling for those issues built in; you may simply use these categories even when +-- you know the result will be smooth in your relevant domain (or must be, for e.g. +-- physics reasons). +-- +-- ¹(The implementation does not deal with /ε/ and /δ/ as difference-bounding +-- reals, but rather as metric tensors that define a boundary by prohibiting the +-- overlap from exceeding one; this makes the concept actually work on general manifolds.) +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 “workhorse” for most serious +-- calculus applications, because it contains all the usual trig etc. functions +-- and of course everything algebraic you can do in the reals. +-- +-- The easiest way to define ordinary functions in this category is hence +-- with its 'AgentVal'ues, which have instances of the standard classes 'Num' +-- through 'Floating'. For instance, the following defines the /binary entropy/ +-- as a differentiable function on the interval @]0,1[@: (it will +-- actually /know/ where it's defined and where not. And I don't mean you +-- need to exhaustively 'isNaN'-check all results...) +-- +-- @ +-- hb :: RWDiffable ℝ ℝ ℝ +-- hb = alg (\\p -> - p * logBase 2 p - (1-p) * logBase 2 (1-p) ) +-- @ +newtype RWDiffable s d c + = RWDiffable { + tryDfblDomain :: d -> (PreRegion s d, Option (Differentiable s d c)) } + +notDefinedHere :: Option (Differentiable s d c) +notDefinedHere = Option Nothing + diff --git a/Data/LinearMap/HerMetric.hs b/Data/LinearMap/HerMetric.hs --- a/Data/LinearMap/HerMetric.hs +++ b/Data/LinearMap/HerMetric.hs @@ -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 “project to the metric's boundary”. +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 + -- | “Anti-normalise” 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 diff --git a/Data/Manifold/PseudoAffine.hs b/Data/Manifold/PseudoAffine.hs --- a/Data/Manifold/PseudoAffine.hs +++ b/Data/Manifold/PseudoAffine.hs @@ -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⁰ = ℝ⁰ diff --git a/Data/Manifold/TreeCover.hs b/Data/Manifold/TreeCover.hs --- a/Data/Manifold/TreeCover.hs +++ b/Data/Manifold/TreeCover.hs @@ -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 @@ -- @ -- -- <> -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 diff --git a/Data/Manifold/Types/Primitive.hs b/Data/Manifold/Types/Primitive.hs --- a/Data/Manifold/Types/Primitive.hs +++ b/Data/Manifold/Types/Primitive.hs @@ -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 diff --git a/Data/VectorSpace/FiniteDimensional.hs b/Data/VectorSpace/FiniteDimensional.hs --- a/Data/VectorSpace/FiniteDimensional.hs +++ b/Data/VectorSpace/FiniteDimensional.hs @@ -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] diff --git a/images/examples/Friedrichs-mollifier.png b/images/examples/Friedrichs-mollifier.png new file mode 100644 Binary files /dev/null and b/images/examples/Friedrichs-mollifier.png differ diff --git a/manifolds.cabal b/manifolds.cabal --- a/manifolds.cabal +++ b/manifolds.cabal @@ -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