diff --git a/Data/LinearMap/HerMetric.hs b/Data/LinearMap/HerMetric.hs new file mode 100644 --- /dev/null +++ b/Data/LinearMap/HerMetric.hs @@ -0,0 +1,269 @@ +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE TupleSections #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE ScopedTypeVariables #-} + + + + +module Data.LinearMap.HerMetric ( + -- * Metric operator types + HerMetric, HerMetric' + -- * Evaluating metrics + , metricSq, metricSq', metric, metric', metrics, metrics' + -- * Defining metrics by projectors + , projector, projector' + -- * Utility + , adjoint + , transformMetric, transformMetric' + , dualiseMetric, dualiseMetric' + , HasMetric(..) + , (^<.>) + , metriScale, metriScale' + ) where + + + + +import Prelude hiding ((^)) + +import Data.VectorSpace +import Data.LinearMap +import Data.Basis +import Data.MemoTrie + +import Control.Applicative + +import Data.Manifold.Types + + +infixr 7 <.>^, ^<.> + + +-- | 'HerMetric' is a portmanteau of /Hermitian/ and /metric/ (in the sense as +-- used in e.g. general relativity – though those particular ones aren't positive +-- definite and thus not really metrics). +-- +-- Mathematically, there are two directly equivalent ways to describe such a metric: +-- as a bilinear mapping of two vectors to a scalar, or as a linear mapping from +-- a vector space to its dual space. We choose the latter, though you can always +-- as well think of metrics as “quadratic dual vectors”. +-- +-- Yet other possible interpretations of this type include /density matrix/ (as in +-- quantum mechanics), /standard range of statistical fluctuations/, and /volume element/. +newtype HerMetric v = HerMetric { getHerMetric :: v :-* DualSpace v } + + +instance HasMetric v => AdditiveGroup (HerMetric v) where + zeroV = HerMetric zeroV + negateV (HerMetric m) = HerMetric $ negateV m + HerMetric m ^+^ HerMetric n = HerMetric $ m ^+^ n +instance HasMetric v => VectorSpace (HerMetric v) where + type Scalar (HerMetric v) = Scalar v + s *^ (HerMetric m) = HerMetric $ s *^ m + +-- | A metric on the dual space; equivalent to a linear mapping from the dual space +-- to the original vector space. +-- +-- Prime-versions of the functions in this module target those dual-space metrics, so +-- we can avoid some explicit handling of double-dual spaces. +newtype HerMetric' v = HerMetric' { dualMetric :: DualSpace v :-* v } +instance (HasMetric v) => AdditiveGroup (HerMetric' v) where + zeroV = HerMetric' zeroV + negateV (HerMetric' m) = HerMetric' $ negateV m + HerMetric' m ^+^ HerMetric' n = HerMetric' $ m ^+^ n +instance (HasMetric v) => VectorSpace (HerMetric' v) where + type Scalar (HerMetric' v) = Scalar v + s *^ (HerMetric' m) = HerMetric' $ s *^ m + + +-- | A metric on @v@ that simply yields the squared overlap of a vector with the +-- given dual-space reference. +-- +-- It will perhaps be the most common way of defining 'HerMetric' values to start +-- with such dual-space vectors and superimpose the projectors using the 'VectorSpace' +-- instance; e.g. @'projector' (1,0) '^+^' 'projector' (0,2)@ yields a hermitian operator +-- describing the ellipsoid span of the vectors /e/₀ and 2⋅/e/₁. +-- Metrics generated this way are positive definite if no negative coefficients have +-- been introduced with the '*^' scaling operator or with '^-^'. +projector :: HasMetric v => DualSpace v -> HerMetric v +projector u = HerMetric (linear $ \v -> u ^* (u<.>^v)) + +projector' :: HasMetric v => v -> HerMetric' v +projector' v = HerMetric' . linear $ \u -> v ^* (v^<.>u) + + + +-- | Evaluate a vector through a metric. For the canonical metric on a Hilbert space, +-- this will be simply 'magnitudeSq'. +metricSq :: HasMetric v => HerMetric v -> v -> Scalar v +metricSq (HerMetric m) v = lapply m v <.>^ v + +metricSq' :: HasMetric v => HerMetric' v -> DualSpace v -> Scalar v +metricSq' (HerMetric' m) u = lapply m u ^<.> u + +-- | Evaluate a vector's “magnitude” through a metric. This assumes an actual +-- mathematical metric, i.e. positive definite – otherwise the internally used +-- square root may get negative arguments (though it can still produce results if the +-- scalars are complex; however, complex spaces aren't supported yet). +metric :: (HasMetric v, Floating (Scalar v)) => HerMetric v -> v -> Scalar v +metric (HerMetric m) v = sqrt $ lapply m v <.>^ v + +metric' :: (HasMetric v, Floating (Scalar v)) => HerMetric' v -> DualSpace v -> Scalar v +metric' (HerMetric' m) u = sqrt $ lapply m u ^<.> u + +metriScale :: (HasMetric v, Floating (Scalar v)) => HerMetric v -> v -> v +metriScale m v = metric m v *^ v + +metriScale' :: (HasMetric v, Floating (Scalar v)) + => HerMetric' v -> DualSpace v -> DualSpace v +metriScale' m v = metric' m v *^ v + + +-- | Square-sum over the metrics for each dual-space vector. +-- +-- @ +-- metrics m vs ≡ sqrt . sum $ metricSq m '<$>' vs +-- @ +metrics :: (HasMetric v, Floating (Scalar v)) => HerMetric v -> [v] -> Scalar v +metrics m vs = sqrt . sum $ metricSq m <$> vs + +metrics' :: (HasMetric v, Floating (Scalar v)) => HerMetric' v -> [DualSpace v] -> Scalar v +metrics' m vs = sqrt . sum $ metricSq' m <$> vs + + +transformMetric :: (HasMetric v, HasMetric w, Scalar v ~ Scalar w) + => (w :-* v) -> HerMetric v -> HerMetric w +transformMetric t (HerMetric m) = HerMetric $ adjoint t *.* m *.* t + +transformMetric' :: ( HasMetric v, HasMetric w, Scalar v ~ Scalar w ) + => (v :-* w) -> HerMetric' v -> HerMetric' w +transformMetric' t (HerMetric' m) + = HerMetric' $ t *.* m *.* adjoint t + +dualiseMetric :: (HasMetric v, HasMetric (DualSpace v)) + => HerMetric (DualSpace v) -> HerMetric' v +dualiseMetric (HerMetric m) = HerMetric' $ linear doubleDual' *.* m + +dualiseMetric' :: (HasMetric v, HasMetric (DualSpace v)) + => HerMetric' v -> HerMetric (DualSpace v) +dualiseMetric' (HerMetric' m) = HerMetric $ linear doubleDual *.* m + + +-- | While the main purpose of this class is to express 'HerMetric', it's actually +-- all about dual spaces. +class ( HasBasis v, VectorSpace (Scalar v), HasTrie (Basis v) + , VectorSpace (DualSpace v), HasBasis (DualSpace v) + , Scalar v ~ Scalar (DualSpace v), Basis v ~ Basis (DualSpace v) ) + => HasMetric v where + + -- | @'DualSpace' v@ is isomorphic to the space of linear functionals on @v@, i.e. + -- @v ':-*' 'Scalar' v@. + -- Typically (for all Hilbert- / 'InnerSpace's) this is in turn isomorphic to @v@ + -- itself, which will be rather more efficient (hence the distinction between a + -- vector space and its dual is often neglected or reduced to “column vs row + -- vectors”). + -- Mathematically though, it makes sense to keep the concepts apart, even if ultimately + -- @'DualSpace' v ~ v@ (which needs not /always/ be the case, though!). + type DualSpace v :: * + type DualSpace v = v + + -- | Apply a dual space vector (aka linear functional) to a vector. + (<.>^) :: DualSpace v -> v -> Scalar v + + -- | Interpret a functional as a dual-space vector. Like 'linear', this /assumes/ + -- (completely unchecked) that the supplied function is linear. + functional :: (v -> Scalar v) -> DualSpace v + + -- | While isomorphism between a space and its dual isn't generally canonical, + -- the /double-dual/ space should be canonically isomorphic in pretty much + -- all relevant cases. Indeed, it is recommended that they are the very same type; + -- the tuple instance actually assumes this to be able to offer an efficient + -- implementation (namely, 'id') of the isomorphisms. + doubleDual :: HasMetric (DualSpace v) => v -> DualSpace (DualSpace v) + doubleDual' :: HasMetric (DualSpace v) => DualSpace (DualSpace v) -> v + + + +-- | Simple flipped version of '<.>^'. +(^<.>) :: HasMetric v => v -> DualSpace v -> Scalar v +ket ^<.> bra = bra <.>^ ket + +instance (VectorSpace k) => HasMetric (ZeroDim k) where + Origin<.>^Origin = zeroV + functional _ = Origin + doubleDual = id; doubleDual'= id +instance HasMetric Double where + (<.>^) = (<.>) + functional f = f 1 + doubleDual = id; doubleDual'= id +instance ( HasMetric v, HasMetric w, Scalar v ~ Scalar w + , HasMetric (DualSpace v), DualSpace (DualSpace v) ~ v + , HasMetric (DualSpace w), DualSpace (DualSpace w) ~ w + ) => HasMetric (v,w) where + type DualSpace (v,w) = (DualSpace v, DualSpace w) + (v,w)<.>^(v',w') = v<.>^v' ^+^ w<.>^w' + functional f = (functional $ f . (,zeroV), functional $ f . (zeroV,)) + doubleDual = id; doubleDual'= id + + + + + +-- | Transpose a linear operator. Contrary to popular belief, this does not +-- just inverse the direction of mapping between the spaces, but also switch to +-- their duals. +adjoint :: (HasMetric v, HasMetric w, Scalar w ~ Scalar v) + => (v :-* w) -> DualSpace w :-* DualSpace v +adjoint m = linear $ \w -> functional $ \v + -> w <.>^lapply m v + + + +metrConst :: (HasMetric v, v ~ DualSpace v, Num (Scalar v)) => Scalar v -> HerMetric v +metrConst = HerMetric . linear . (*^) + +instance (HasMetric v, v ~ DualSpace v, Num (Scalar v)) => Num (HerMetric v) where + fromInteger = metrConst . fromInteger + (+) = (^+^) + negate = negateV + + -- | This does /not/ work correctly if the metrics don't share an eigenbasis! + HerMetric m * HerMetric n = HerMetric $ m *.* n + + -- | Undefined, though it could actually be done. + abs = error "abs undefined for HerMetric" + signum = error "signum undefined for HerMetric" + + +metrNumFun :: (HasMetric v, v ~ Scalar v, v ~ DualSpace v, Num v) + => (v -> v) -> HerMetric v -> HerMetric v +metrNumFun f (HerMetric m) = HerMetric . linear . (*^) . f $ lapply m 1 + +instance (HasMetric v, v ~ Scalar v, v ~ DualSpace v, Fractional v) + => Fractional (HerMetric v) where + fromRational = metrConst . fromRational + recip = metrNumFun recip + +instance (HasMetric v, v ~ Scalar v, v ~ DualSpace v, Floating v) + => Floating (HerMetric v) where + pi = metrConst pi + sqrt = metrNumFun sqrt + exp = metrNumFun exp + log = metrNumFun log + sin = metrNumFun sin + cos = metrNumFun cos + tan = metrNumFun tan + asin = metrNumFun asin + acos = metrNumFun acos + atan = metrNumFun atan + sinh = metrNumFun sinh + cosh = metrNumFun cosh + asinh = metrNumFun asinh + atanh = metrNumFun atanh + acosh = metrNumFun acosh diff --git a/Data/Manifold.hs b/Data/Manifold.hs --- a/Data/Manifold.hs +++ b/Data/Manifold.hs @@ -29,7 +29,7 @@ {-# LANGUAGE RecordWildCards #-} -module Data.Manifold where +module Data.Manifold (module Data.Manifold, module Data.Manifold.Types) where import Data.List import Data.Maybe @@ -41,6 +41,7 @@ import Data.Basis import Data.Complex hiding (magnitude) import Data.Void +import Data.Manifold.Types import qualified Prelude @@ -180,8 +181,6 @@ -type EuclidSpace v = (HasBasis v, EqFloating(Scalar v), Eq v) -type EqFloating f = (Eq f, Ord f, Floating f) -- | A chart is a homeomorphism from a connected, open subset /Q/ ⊂ /M/ of @@ -368,18 +367,18 @@ data CntnFuncValue d c = CntnFuncValue { runCntnFuncValue :: d :--> c } | CntnFuncConst c -instance HasProxy (:-->) where - type ProxyVal (:-->) d c = CntnFuncValue d c +instance HasAgent (:-->) where + type AgentVal (:-->) d c = CntnFuncValue d c alg f = case f $ CntnFuncValue id of CntnFuncValue q -> q CntnFuncConst c -> const__ c f $~ CntnFuncValue g = CntnFuncValue $ f . g f $~ CntnFuncConst c = CntnFuncConst $ f $ c -instance PointProxy CntnFuncValue (:-->) d c where +instance PointAgent CntnFuncValue (:-->) d c where point = CntnFuncConst -instance CartesianProxy (:-->) where +instance CartesianAgent (:-->) where alg1to2 f = case f $ CntnFuncValue id of (CntnFuncConst c₁, CntnFuncConst c₂) -> const__ (c₁, c₂) (CntnFuncConst c₁, CntnFuncValue f₂) @@ -523,10 +522,6 @@ -data GraphWindowSpec = GraphWindowSpec { - lBound, rBound, bBound, tBound :: Double - , xResolution, yResolution :: Int - } finiteGraphContinℝtoℝ :: GraphWindowSpec -> (Double:-->Double) -> [(Double, Double)] finiteGraphContinℝtoℝ (GraphWindowSpec{..}) fc @@ -606,10 +601,6 @@ -data S2 = S2 { ϑParamS2 :: Double -- [0, π[ - , φParamS2 :: Double -- [0, 2π[ - } - -- instance Manifold S2 where -- type TangentSpace S2 = (Double, Double) @@ -633,7 +624,6 @@ -type Endomorphism a = a->a (.:) :: (c->d) -> (a->b->c) -> a->b->d @@ -645,20 +635,8 @@ -type ℝ = Double -instance VectorSpace () where - type Scalar () = ℝ - _ *^ () = () -instance HasBasis () where - type Basis () = Void - basisValue = absurd - decompose () = [] - decompose' () = absurd -instance InnerSpace () where - () <.> () = 0 - class (RealFloat (Metric v), InnerSpace v) => MetricSpace v where type Metric v :: * type Metric v = ℝ @@ -697,6 +675,4 @@ -(^) :: Num a => a -> Int -> a -(^) = (Prelude.^) diff --git a/Data/Manifold/PseudoAffine.hs b/Data/Manifold/PseudoAffine.hs new file mode 100644 --- /dev/null +++ b/Data/Manifold/PseudoAffine.hs @@ -0,0 +1,1015 @@ +-- | +-- Module : Data.Manifold.PseudoAffine +-- Copyright : (c) Justus Sagemüller 2015 +-- License : GPL v3 +-- +-- Maintainer : (@) sagemueller $ geo.uni-koeln.de +-- Stability : experimental +-- Portability : portable +-- +-- This is the second prototype of a manifold class. It appears to give considerable +-- advantages over 'Data.Manifold.Manifold', so that class will probably soon be replaced +-- with the one we define here (though 'PseudoAffine' does not follow the standard notion +-- of a manifold very closely, it should work quite equivalently for pretty much all +-- Haskell types that qualify as manifolds). +-- +-- Manifolds are interesting as objects of various categories, from continuous to +-- diffeomorphic. At the moment, we mainly focus on /region-wise differentiable functions/, +-- which are a promising compromise between flexibility of definition and provability of +-- analytic properties. In particular, they are well-suited for visualisation purposes. + +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE TupleSections #-} +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE PatternGuards #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE RecordWildCards #-} +{-# LANGUAGE CPP #-} + + +module Data.Manifold.PseudoAffine ( + -- * Manifold class + PseudoAffine(..) + -- * Regions within a manifold + , Region + -- * Hierarchy of manifold-categories + , Differentiable + , PWDiffable, RWDiffable + ) where + + + +import Data.List +import Data.Maybe +import Data.Semigroup +import Data.Function (on) +import Data.Fixed + +import Data.VectorSpace +import Data.LinearMap +import Data.LinearMap.HerMetric +import Data.MemoTrie (HasTrie) +import Data.AffineSpace +import Data.Basis +import Data.Complex hiding (magnitude) +import Data.Void +import Data.Tagged +import Data.Manifold.Types + +import qualified Prelude + +import Control.Category.Constrained.Prelude hiding ((^)) +import Control.Arrow.Constrained +import Control.Monad.Constrained +import Data.Foldable.Constrained + + + + +infix 6 .-~. +infixl 6 .+~^ + +-- | 'PseudoAffine' is intended as an alternative class for 'Data.Manifold.Manifold's. +-- The interface is almost identical to the better-known 'AffineSpace' class, but unlike +-- in the mathematical definition of affine spaces we don't require associativity +-- of '.+~^' with '^+^' – except in an asymptotic sense for small vectors. +-- +-- That innocent-looking change makes the class applicable to vastly more general types: +-- while an affine space is basically nothing but a vector space without particularly +-- designated origin, a pseudo-affine space can have nontrivial topology on the global +-- scale, and yet be used in practically the same way as an affine space. At least the +-- usual spheres and tori make good instances, perhaps the class is in fact equivalent to +-- /parallelisable manifolds/. +class PseudoAffine x where + type PseudoDiff x :: * + (.-~.) :: x -> x -> Option (PseudoDiff x) + (.+~^) :: x -> PseudoDiff x -> x + + +type LocallyScalable s x = ( PseudoAffine x, (PseudoDiff x) ~ PseudoDiff x + , HasMetric (PseudoDiff x) + , DualSpace (PseudoDiff x) ~ DualSpace (PseudoDiff x) + , HasMetric (DualSpace (PseudoDiff x)) + , PseudoDiff x ~ DualSpace (DualSpace (PseudoDiff x)) + , s ~ Scalar (PseudoDiff x) + , s ~ Scalar (DualSpace (PseudoDiff x)) ) +type LinearManifold s x = ( PseudoAffine x, PseudoDiff x ~ x + , HasMetric x, HasMetric (DualSpace x) + , DualSpace (DualSpace x) ~ x + , s ~ Scalar x, s ~ Scalar (DualSpace x) ) +type RealDimension r = ( PseudoAffine r, PseudoDiff r ~ r + , HasMetric r, DualSpace r ~ r, Scalar r ~ r + , RealFloat r ) + + + +palerp :: (PseudoAffine x, VectorSpace (PseudoDiff x)) + => x -> x -> Option (Scalar (PseudoDiff x) -> x) +palerp p1 p2 = fmap (\v t -> p1 .+~^ t *^ v) $ p2 .-~. p1 + + + +#define deriveAffine(t) \ +instance PseudoAffine t where { \ + type PseudoDiff t = Diff t; \ + a.-~.b = pure (a.-.b); \ + (.+~^) = (.+^) } + +deriveAffine(Double) +deriveAffine(Rational) + +instance PseudoAffine (ZeroDim k) where + type PseudoDiff (ZeroDim k) = ZeroDim k + Origin .-~. Origin = pure Origin + Origin .+~^ Origin = Origin +instance (PseudoAffine a, PseudoAffine b) => PseudoAffine (a,b) where + type PseudoDiff (a,b) = (PseudoDiff a, PseudoDiff b) + (a,b).-~.(c,d) = liftA2 (,) (a.-~.c) (b.-~.d) + (a,b).+~^(v,w) = (a.+~^v, b.+~^w) +instance (PseudoAffine a, PseudoAffine b, PseudoAffine c) => PseudoAffine (a,b,c) where + type PseudoDiff (a,b,c) = (PseudoDiff a, PseudoDiff b, PseudoDiff c) + (a,b,c).-~.(d,e,f) = liftA3 (,,) (a.-~.d) (b.-~.e) (c.-~.f) + (a,b,c).+~^(v,w,x) = (a.+~^v, b.+~^w, c.+~^x) + + +instance PseudoAffine S¹ where + type PseudoDiff S¹ = ℝ + S¹ φ₁ .-~. S¹ φ₀ + | δφ > pi = pure (δφ - 2*pi) + | δφ < (-pi) = pure (δφ + 2*pi) + | otherwise = pure δφ + where δφ = φ₁ - φ₀ + S¹ φ₀ .+~^ δφ + | φ' < 0 = S¹ $ φ' + tau + | otherwise = S¹ $ φ' + where φ' = (φ₀ + δφ)`mod'`tau + +instance PseudoAffine S² where + type PseudoDiff S² = ℝ² + S² ϑ₁ φ₁ .-~. S² ϑ₀ φ₀ + | ϑ₀ < pi/2 = pure ( ϑ₁*^embed(S¹ φ₁) ^-^ ϑ₀*^embed(S¹ φ₀) ) + | otherwise = pure ( (pi-ϑ₁)*^embed(S¹ φ₁) ^-^ (pi-ϑ₀)*^embed(S¹ φ₀) ) + S² ϑ₀ φ₀ .+~^ δv + | ϑ₀ < pi/2 = sphereFold PositiveHalfSphere $ ϑ₀*^embed(S¹ φ₀) ^+^ δv + | otherwise = sphereFold NegativeHalfSphere $ (pi-ϑ₀)*^embed(S¹ φ₀) ^+^ δv + +sphereFold :: S⁰ -> ℝ² -> S² +sphereFold hfSphere v + | ϑ₀ > pi = S² (inv $ tau - ϑ₀) ((φ₀+pi)`mod'`tau) + | otherwise = S² (inv ϑ₀) φ₀ + where S¹ φ₀ = coEmbed v + ϑ₀ = magnitude v `mod'` tau + inv ϑ = case hfSphere of PositiveHalfSphere -> ϑ + NegativeHalfSphere -> pi - ϑ + + + +tau :: Double +tau = 2 * pi + + + + + +type LinDevPropag d c = HerMetric (PseudoDiff c) -> HerMetric (PseudoDiff d) + +dev_ε_δ :: RealDimension a + => (a -> a) -> LinDevPropag a a +dev_ε_δ f d = let ε = 1 / metric d 1 in projector $ 1 / sqrt (f ε) + +-- | The category of differentiable functions between manifolds over scalar @s@. +-- +-- As you might guess, these offer /automatic differentiation/ of sorts (basically, +-- simple forward AD), but that's in itself is not really the killer feature here. +-- More interestingly, we actually have the (à la Curry-Howard) /proof/ +-- built in: the function /f/ has at /x/₀ derivative /f'ₓ/₀, +-- if, for¹ /ε/>0, there exists /δ/ such that +-- |/f/ /x/ − (/f/ /x/₀ + /x/⋅/f'ₓ/₀)| < /ε/ +-- for all |/x/ − /x/₀| < /δ/. +-- +-- Observe that, though this looks quite similar to the standard definition +-- of differentiability, it is not equivalent thereto – in fact it does +-- not prove any analytic properties at all. To make it equivalent, we need +-- a lower bound on /δ/: simply /δ/ gives us continuity, and for +-- continuous differentiability, /δ/ must grow at least like √/ε/ +-- for small /ε/. Neither of these conditions are enforced by the type system, +-- but we do require them for any allowed values because these proofs are obviously +-- tremendously useful – for instance, you can have a root-finding algorithm +-- and actually be sure you get /all/ solutions correctly, not just /some/ that are +-- (hopefully) the closest to some reference point you'd need to laborously define! +-- +-- Unfortunately however, this also prevents doing any serious algebra etc. with the +-- category, because even something as simple as division necessary introduces singularities +-- where the derivatives must diverge. +-- Not to speak of many trigonometric e.g. trigonometric functions that +-- are undefined on whole regions. The 'PWDiffable' and 'RWDiffable' categories have explicit +-- handling for those issues built in; you may simply use these categories even when +-- you know the result will be smooth in your relevant domain (or must be, for e.g. +-- physics reasons). +-- +-- ¹(The implementation does not deal with /ε/ and /δ/ as difference-bounding +-- reals, but rather as metric tensors that define a boundary by prohibiting the +-- overlap from exceeding one; this makes the concept actually work on general manifolds.) +newtype Differentiable s d c + = Differentiable { runDifferentiable :: + d -> ( c, PseudoDiff d :-* PseudoDiff c, LinDevPropag d c ) } +type (-->) = Differentiable ℝ + + +instance (VectorSpace s) => Category (Differentiable s) where + type Object (Differentiable s) o = LocallyScalable s o + id = Differentiable $ \x -> (x, idL, const zeroV) + Differentiable f . Differentiable g = Differentiable $ + \x -> let (y, g', devg) = g x + (z, f', devf) = f y + devfg δz = let δy = transformMetric f' δz + εy = devf δz + in transformMetric g' εy ^+^ devg δy ^+^ devg εy + in (z, f'*.*g', devfg) + + +instance (VectorSpace s) => Cartesian (Differentiable s) where + type UnitObject (Differentiable s) = ZeroDim s + swap = Differentiable $ \(x,y) -> ((y,x), lSwap, const zeroV) + where lSwap = linear swap + attachUnit = Differentiable $ \x -> ((x, Origin), lAttachUnit, const zeroV) + where lAttachUnit = linear $ \x -> (x, Origin) + detachUnit = Differentiable $ \(x, Origin) -> (x, lDetachUnit, const zeroV) + where lDetachUnit = linear $ \(x, Origin) -> x + regroup = Differentiable $ \(x,(y,z)) -> (((x,y),z), lRegroup, const zeroV) + where lRegroup = linear regroup + regroup' = Differentiable $ \((x,y),z) -> ((x,(y,z)), lRegroup, const zeroV) + where lRegroup = linear regroup' + + +instance (VectorSpace s) => Morphism (Differentiable s) where + Differentiable f *** Differentiable g = Differentiable h + where h (x,y) = ((fx, gy), lPar, devfg) + where (fx, f', devf) = f x + (gy, g', devg) = g y + devfg δs = transformMetric lfst δx + ^+^ transformMetric lsnd δy + where δx = devf $ transformMetric lcofst δs + δy = devg $ transformMetric lcosnd δs + lPar = linear $ lapply f'***lapply g' + lfst = linear fst; lsnd = linear snd + lcofst = linear (,zeroV); lcosnd = linear (zeroV,) + + +instance (VectorSpace s) => PreArrow (Differentiable s) where + terminal = Differentiable $ \_ -> (Origin, zeroV, const zeroV) + fst = Differentiable $ \(x,_) -> (x, lfst, const zeroV) + where lfst = linear fst + snd = Differentiable $ \(_,y) -> (y, lsnd, const zeroV) + where lsnd = linear snd + Differentiable f &&& Differentiable g = Differentiable h + where h x = ((fx, gx), lFanout, devfg) + where (fx, f', devf) = f x + (gx, g', devg) = g x + devfg δs = (devf $ transformMetric lcofst δs) + ^+^ (devg $ transformMetric lcosnd δs) + lFanout = linear $ lapply f'&&&lapply g' + lcofst = linear (,zeroV); lcosnd = linear (zeroV,) + + +instance (VectorSpace s) => WellPointed (Differentiable s) where + unit = Tagged Origin + globalElement x = Differentiable $ \Origin -> (x, zeroV, const zeroV) + const x = Differentiable $ \_ -> (x, zeroV, const zeroV) + + + +type DfblFuncValue s = GenericAgent (Differentiable s) + +instance (VectorSpace s) => HasAgent (Differentiable s) where + alg = genericAlg + ($~) = genericAgentMap +instance (VectorSpace s) => CartesianAgent (Differentiable s) where + alg1to2 = genericAlg1to2 + alg2to1 = genericAlg2to1 + alg2to2 = genericAlg2to2 +instance (VectorSpace s) + => PointAgent (DfblFuncValue s) (Differentiable s) a x where + point = genericPoint + + + +actuallyLinear :: ( LinearManifold s x, LinearManifold s y ) + => (x:-*y) -> Differentiable s x y +actuallyLinear f = Differentiable $ \x -> (lapply f x, f, const zeroV) + +actuallyAffine :: ( LinearManifold s x, LinearManifold s y ) + => y -> (x:-*y) -> Differentiable s x y +actuallyAffine y₀ f = Differentiable $ \x -> (y₀ ^+^ lapply f x, f, const zeroV) + + +dfblFnValsFunc :: ( LocallyScalable s c, LocallyScalable s c', LocallyScalable s d + , v ~ PseudoDiff c, v' ~ PseudoDiff c' + , ε ~ HerMetric v, ε ~ HerMetric v' ) + => (c' -> (c, v':-*v, ε->ε)) -> DfblFuncValue s d c' -> DfblFuncValue s d c +dfblFnValsFunc f = (Differentiable f $~) + +dfblFnValsCombine :: forall d c c' c'' v v' v'' ε ε' ε'' s. + ( LocallyScalable s c, LocallyScalable s c', LocallyScalable s c'' + , LocallyScalable s d + , v ~ PseudoDiff c, v' ~ PseudoDiff c', v'' ~ PseudoDiff c'' + , ε ~ HerMetric v , ε' ~ HerMetric v' , ε'' ~ HerMetric v'', ε~ε', ε~ε'' ) + => ( c' -> c'' -> (c, (v',v''):-*v, ε -> (ε',ε'')) ) + -> DfblFuncValue s d c' -> DfblFuncValue s d c'' -> DfblFuncValue s d c +dfblFnValsCombine cmb (GenericAgent (Differentiable f)) + (GenericAgent (Differentiable g)) + = GenericAgent . Differentiable $ + \d -> let (c', f', devf) = f d + (c'', g', devg) = g d + (c, h', devh) = cmb c' c'' + h'l = h' *.* lcofst; h'r = h' *.* lcosnd + in ( c + , h' *.* linear (lapply f' &&& lapply g') + , \εc -> let εc' = transformMetric h'l εc + εc'' = transformMetric h'r εc + (δc',δc'') = devh εc + in devf εc' ^+^ devg εc'' + ^+^ transformMetric f' δc' + ^+^ transformMetric g' δc'' + ) + where lcofst = linear(,zeroV) + lcosnd = linear(zeroV,) + + + + + +instance (LinearManifold s v, LocallyScalable s a, Floating s) + => AdditiveGroup (DfblFuncValue s a v) where + zeroV = point zeroV + (^+^) = dfblFnValsCombine $ \a b -> (a^+^b, lPlus, const zeroV) + where lPlus = linear $ uncurry (^+^) + negateV = dfblFnValsFunc $ \a -> (negateV a, lNegate, const zeroV) + where lNegate = linear negateV + +instance (RealDimension n, LocallyScalable n a) + => Num (DfblFuncValue n a n) where + fromInteger i = point $ fromInteger i + (+) = dfblFnValsCombine $ \a b -> (a+b, lPlus, const zeroV) + where lPlus = linear $ uncurry (+) + (*) = dfblFnValsCombine $ + \a b -> ( a*b + , linear $ \(da,db) -> a*db + b*da + , \d -> let d¹₂ = sqrt d in (d¹₂,d¹₂) + -- ε δa δb = (a+δa)·(b+δb) - (a·b + (a·δa + b·δb)) + -- = δa·δb + -- so choose δa = δb = √ε + ) + negate = dfblFnValsFunc $ \a -> (negate a, lNegate, const zeroV) + where lNegate = linear negate + abs = dfblFnValsFunc dfblAbs + where dfblAbs a + | a>0 = (a, idL, dev_ε_δ $ \ε -> a + ε/2) + | a<0 = (-a, negateV idL, dev_ε_δ $ \ε -> ε/2 - a) + | otherwise = (0, zeroV, (^/ sqrt 2)) + signum = dfblFnValsFunc dfblSgn + where dfblSgn a + | a>0 = (1, zeroV, dev_ε_δ $ const a) + | a<0 = (-1, zeroV, dev_ε_δ $ \_ -> -a) + | otherwise = (0, zeroV, const $ projector 1) + + + +-- VectorSpace instance is more problematic than you'd think: multiplication +-- requires the allowed-deviation backpropagators to be split as square +-- roots, but the square root of a nontrivial-vector-space metric requires +-- an eigenbasis transform, which we have not implemented yet. +-- +-- instance (LinearManifold s v, LocallyScalable s a, Floating s) +-- => VectorSpace (DfblFuncValue s a v) where +-- type Scalar (DfblFuncValue s a v) = DfblFuncValue s a (Scalar v) +-- (*^) = dfblFnValsCombine $ \μ v -> (μ*^v, lScl, \ε -> (ε ^* sqrt 2, ε ^* sqrt 2)) +-- where lScl = linear $ uncurry (*^) + + +-- | Important special operator needed to compute intersection of 'Region's. +minDblfuncs :: (LocallyScalable s m, RealDimension s) + => Differentiable s m s -> Differentiable s m s -> Differentiable s m s +minDblfuncs (Differentiable f) (Differentiable g) = Differentiable h + where h x + | fx==gx = ( fx, (f'^+^g')^/2 + , \d -> devf d ^+^ devg d + ^+^ transformMetric (f'^-^g') + (projector $ metric d 1) ) + | fx < gx = ( fx, f' + , \d -> devf d + ^+^ transformMetric (f'^-^g') + (projector $ metric d 1 + gx - fx) ) + where (fx, f', devf) = f x + (gx, g', devg) = g x + + + +-- | A pathwise connected subset of a manifold @m@, whose tangent space has scalar @s@. +data Region s m = Region { regionRefPoint :: m + , regionRDef :: PreRegion s m } + +-- | A 'PreRegion' needs to be associated with a certain reference point ('Region' +-- includes that point) to define a connected subset of a manifold. +data PreRegion s m where + GlobalRegion :: PreRegion s m + PreRegion :: (Differentiable s m s) -- A function that is positive at reference point /p/, + -- decreases and crosses zero at the region's + -- boundaries. (If it goes positive again somewhere + -- else, these areas shall /not/ be considered + -- belonging to the (by definition connected) region.) + -> PreRegion s m + +-- | Set-intersection of regions would not be guaranteed to yield a connected result +-- or even have the reference point of one region contained in the other. This +-- combinator assumes (unchecked) that the references are in a connected +-- sub-intersection, which is used as the result. +unsafePreRegionIntersect :: (RealDimension s, LocallyScalable s a) + => PreRegion s a -> PreRegion s a -> PreRegion s a +unsafePreRegionIntersect GlobalRegion r = r +unsafePreRegionIntersect r GlobalRegion = r +unsafePreRegionIntersect (PreRegion ra) (PreRegion rb) = PreRegion $ minDblfuncs ra rb + +-- | Cartesian product of two regions. +regionProd :: (RealDimension s, LocallyScalable s a, LocallyScalable s b) + => Region s a -> Region s b -> Region s (a,b) +regionProd (Region a₀ ra) (Region b₀ rb) = Region (a₀,b₀) (preRegionProd ra rb) + +-- | Cartesian product of two pre-regions. +preRegionProd :: (RealDimension s, LocallyScalable s a, LocallyScalable s b) + => PreRegion s a -> PreRegion s b -> PreRegion s (a,b) +preRegionProd GlobalRegion GlobalRegion = GlobalRegion +preRegionProd GlobalRegion (PreRegion rb) = PreRegion $ rb . snd +preRegionProd (PreRegion ra) GlobalRegion = PreRegion $ ra . fst +preRegionProd (PreRegion ra) (PreRegion rb) = PreRegion $ minDblfuncs (ra.fst) (rb.snd) + + +positivePreRegion, negativePreRegion :: (RealDimension s) => PreRegion s s +positivePreRegion = PreRegion $ Differentiable prr + where prr x = (1 - 1/xp1, (1/xp1²) *^ idL, dev_ε_δ δ ) + -- ε = (1 − 1/(1+x)) + (-δ · 1/(x+1)²) − (1 − 1/(1+x−δ)) + -- = 1/(1+x−δ) − 1/(1+x) − δ · 1/(x+1)² + -- ε·(1+x−δ) = 1 − (1+x−δ)/(1+x) − δ·(1+x-δ)/(x+1)² + -- ε + ε·x − ε·δ = 1 − 1/(1+x) − x/(1+x) + δ/(1+x) − δ/(x+1) + δ²/(x+1)² + -- = 1 − 1/(1+x) − x/(1+x) + δ²/(x+1)² + -- = (1+x − 1 − x)/(1+x) + δ²/(x+1)² + -- 0 = δ² + ε·(x+1)²·δ + ε·(x+1)³ + -- δ = let mph = -ε·(x+1)²/2 + -- in mph + sqrt(mph² - ε·(x+1)³) + where δ ε = let mph = -ε*xp1²/2 + in mph + sqrt(mph^2 - ε * xp1² * xp1) + xp1 = (x+1) + xp1² = xp1 ^ 2 +negativePreRegion = PreRegion $ ppr . ngt + where PreRegion ppr = positivePreRegion + ngt = actuallyLinear $ linear negate + +preRegionToInfFrom, preRegionFromMinInfTo :: RealDimension s => s -> PreRegion s s +preRegionToInfFrom xs = PreRegion $ ppr . trl + where PreRegion ppr = positivePreRegion + trl = actuallyAffine (-xs) idL +preRegionFromMinInfTo xe = PreRegion $ ppr . flp + where PreRegion ppr = positivePreRegion + flp = actuallyAffine (-xe) (linear negate) + +intervalPreRegion :: RealDimension s => (s,s) -> PreRegion s s +intervalPreRegion (lb,rb) = PreRegion $ Differentiable prr + where m = lb + radius; radius = (rb - lb)/2 + prr x = ( 1 - ((x-m)/radius)^2 + , (2*(m-x)/radius^2) *^ idL + , dev_ε_δ $ (*radius) . sqrt ) + + + + +-- | Category of functions that almost everywhere have an open region in +-- which they are continuously differentiable, i.e. /PieceWiseDiff'able/. +newtype PWDiffable s d c + = PWDiffable { + getDfblDomain :: d -> (PreRegion s d, Differentiable s d c) } + + + +instance (RealDimension s) => Category (PWDiffable s) where + type Object (PWDiffable s) o = LocallyScalable s o + id = PWDiffable $ \x -> (GlobalRegion, id) + PWDiffable f . PWDiffable g = PWDiffable h + where h x₀ = case g x₀ of + (GlobalRegion, gr) + -> let (y₀,_,_) = runDifferentiable gr x₀ + in case f y₀ of + (GlobalRegion, fr) -> (GlobalRegion, fr . gr) + (PreRegion ry, fr) + -> ( PreRegion $ ry . gr, fr . gr ) + (PreRegion rx, gr) + -> let (y₀,_,_) = runDifferentiable gr x₀ + in case f y₀ of + (GlobalRegion, fr) -> (PreRegion rx, fr . gr) + (PreRegion ry, fr) + -> ( PreRegion $ minDblfuncs (ry . gr) rx + , fr . gr ) + where (rx, gr) = g x₀ + +globalDiffable :: Differentiable s a b -> PWDiffable s a b +globalDiffable f = PWDiffable $ const (GlobalRegion, f) + +instance (RealDimension s) => EnhancedCat (PWDiffable s) (Differentiable s) where + arr = globalDiffable + +instance (RealDimension s) => 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 ~ PseudoDiff c, v' ~ PseudoDiff 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 ~ PseudoDiff c, v' ~ PseudoDiff c', v'' ~ PseudoDiff 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 (LinearManifold s 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 = (PWDiffable rcipPW $~) + where rcipPW a₀ + | a₀<0 = (negativePreRegion, Differentiable negp) + | otherwise = (positivePreRegion, Differentiable posp) + negp x = (x'¹, (- x'¹^2) *^ idL, dev_ε_δ δ) + -- ε = 1/x − δ/x² − 1/(x+δ) + -- ε·x + ε·δ = 1 + δ/x − δ/x − δ²/x² − 1 + -- = -δ²/x² + -- 0 = δ² + ε·x²·δ + ε·x³ + -- δ = let mph = -ε·x²/2 in mph + sqrt (mph² − ε·x³) + where δ ε = let mph = -ε*x^2/2 in mph + sqrt (mph^2 - ε*x^3) + x'¹ = recip x + posp x = (x'¹, (- x'¹^2) *^ idL, dev_ε_δ δ) + where δ ε = let mph = -ε*x^2/2 in mph + sqrt (mph^2 + ε*x^3) + x'¹ = recip x + + + + + + +-- | Category of functions that, where defined, have an open region in +-- which they are continuously differentiable. Hence /RegionWiseDiff'able/. +-- Basically these are the partial version of `PWDiffable`. +-- +-- Though the possibility of undefined regions is of course not too nice +-- (we don't need Java to demonstrate this with its everywhere-looming @null@ values...), +-- this category will propably be the “workhorse” for most serious +-- calculus applications, because it contains all the usual trig etc. functions +-- and of course everything algebraic you can do in the reals. +-- +-- The easiest way to define ordinary functions in this category is hence +-- with its 'AgentVal'ues, which have instances of the standard classes 'Num' +-- through 'Floating'. For instance, the following defines the /binary entropy/ +-- as a differentiable function on the interval @]0,1[@: (it will +-- actually /know/ where it's defined and where not! – and I don't mean you +-- need to exhaustively 'isNaN'-check all results...) +-- +-- @ +-- hb :: RWDiffable R R R +-- hb = alg (\\p -> - p * logBase 2 p - (1-p) * logBase 2 (1-p) ) +-- @ +newtype RWDiffable s d c + = RWDiffable { + tryDfblDomain :: d -> (PreRegion s d, Option (Differentiable s d c)) } + +notDefinedHere :: Option (Differentiable s d c) +notDefinedHere = Option Nothing + + + +instance (RealDimension s) => Category (RWDiffable s) where + type Object (RWDiffable s) o = LocallyScalable s o + id = RWDiffable $ \x -> (GlobalRegion, pure id) + RWDiffable f . RWDiffable g = RWDiffable h + where h x₀ = case g x₀ of + (GlobalRegion, Option Nothing) + -> (GlobalRegion, notDefinedHere) + (GlobalRegion, Option (Just gr)) + -> let (y₀,_,_) = runDifferentiable gr x₀ + in case f y₀ of + (GlobalRegion, Option Nothing) + -> (GlobalRegion, notDefinedHere) + (GlobalRegion, Option (Just fr)) + -> (GlobalRegion, pure (fr . gr)) + (PreRegion ry, Option Nothing) + -> ( PreRegion $ ry . gr, Option Nothing ) + (PreRegion ry, Option (Just fr)) + -> ( PreRegion $ ry . gr, pure (fr . gr) ) + (PreRegion rx, Option Nothing) + -> (PreRegion rx, notDefinedHere) + (PreRegion rx, Option (Just gr)) + -> let (y₀,_,_) = runDifferentiable gr x₀ + in case f y₀ of + (GlobalRegion, Option Nothing) + -> (PreRegion rx, notDefinedHere) + (GlobalRegion, Option (Just fr)) + -> (PreRegion rx, pure (fr . gr)) + (PreRegion ry, Option Nothing) + -> ( PreRegion $ minDblfuncs (ry . gr) rx + , notDefinedHere ) + (PreRegion ry, Option (Just fr)) + -> ( PreRegion $ minDblfuncs (ry . gr) rx + , pure (fr . gr) ) + where (rx, gr) = g x₀ + + +globalDiffable' :: Differentiable s a b -> RWDiffable s a b +globalDiffable' f = RWDiffable $ const (GlobalRegion, pure f) + +pwDiffable :: PWDiffable s a b -> RWDiffable s a b +pwDiffable (PWDiffable q) = RWDiffable $ \x₀ -> let (r₀,f₀) = q x₀ in (r₀, pure f₀) + + + +instance (RealDimension s) => EnhancedCat (RWDiffable s) (Differentiable s) where + arr = globalDiffable' +instance (RealDimension s) => EnhancedCat (RWDiffable s) (PWDiffable s) where + arr = pwDiffable + +instance (RealDimension s) => Cartesian (RWDiffable s) where + type UnitObject (RWDiffable s) = ZeroDim s + swap = globalDiffable' swap + attachUnit = globalDiffable' attachUnit + detachUnit = globalDiffable' detachUnit + regroup = globalDiffable' regroup + regroup' = globalDiffable' regroup' + +instance (RealDimension s) => Morphism (RWDiffable s) where + RWDiffable f *** RWDiffable g = RWDiffable h + where h (x,y) = (preRegionProd rfx rgy, liftA2 (***) dff dfg) + where (rfx, dff) = f x + (rgy, dfg) = g y + +instance (RealDimension s) => PreArrow (RWDiffable s) where + RWDiffable f &&& RWDiffable g = RWDiffable h + where h x = (unsafePreRegionIntersect rfx rgx, liftA2 (&&&) dff dfg) + where (rfx, dff) = f x + (rgx, dfg) = g x + terminal = globalDiffable' terminal + fst = globalDiffable' fst + snd = globalDiffable' snd + + +instance (RealDimension s) => WellPointed (RWDiffable s) where + unit = Tagged Origin + globalElement x = RWDiffable $ \Origin -> (GlobalRegion, pure (globalElement x)) + const x = RWDiffable $ \_ -> (GlobalRegion, pure (const x)) + + +type RWDfblFuncValue s = GenericAgent (RWDiffable s) + +instance RealDimension s => HasAgent (RWDiffable s) where + alg = genericAlg + ($~) = genericAgentMap +instance RealDimension s => CartesianAgent (RWDiffable s) where + alg1to2 = genericAlg1to2 + alg2to1 = genericAlg2to1 + alg2to2 = genericAlg2to2 +instance (RealDimension s) + => PointAgent (RWDfblFuncValue s) (RWDiffable s) a x where + point = genericPoint + +grwDfblFnValsFunc + :: ( RealDimension s + , LocallyScalable s c, LocallyScalable s c', LocallyScalable s d + , v ~ PseudoDiff c, v' ~ PseudoDiff c' + , ε ~ HerMetric v, ε ~ HerMetric v' ) + => (c' -> (c, v':-*v, ε->ε)) -> RWDfblFuncValue s d c' -> RWDfblFuncValue s d c +grwDfblFnValsFunc f = (RWDiffable (\_ -> (GlobalRegion, pure (Differentiable f))) $~) + +grwDfblFnValsCombine :: forall d c c' c'' v v' v'' ε ε' ε'' s. + ( LocallyScalable s c, LocallyScalable s c', LocallyScalable s c'' + , LocallyScalable s d, RealDimension s + , v ~ PseudoDiff c, v' ~ PseudoDiff c', v'' ~ PseudoDiff c'' + , ε ~ HerMetric v , ε' ~ HerMetric v' , ε'' ~ HerMetric v'', ε~ε', ε~ε'' ) + => ( c' -> c'' -> (c, (v',v''):-*v, ε -> (ε',ε'')) ) + -> RWDfblFuncValue s d c' -> RWDfblFuncValue s d c'' -> RWDfblFuncValue s d c +grwDfblFnValsCombine cmb (GenericAgent (RWDiffable fpcs)) + (GenericAgent (RWDiffable gpcs)) + = GenericAgent . RWDiffable $ + \d₀ -> let (rc', fmay) = fpcs d₀ + (rc'',gmay) = gpcs d₀ + in (unsafePreRegionIntersect rc' rc'',) $ + case (fmay,gmay) of + (Option(Just(Differentiable f)), Option(Just(Differentiable g))) -> + pure . Differentiable $ \d + -> let (c', f', devf) = f d + (c'',g', devg) = g d + (c, h', devh) = cmb c' c'' + h'l = h' *.* lcofst; h'r = h' *.* lcosnd + in ( c + , h' *.* linear (lapply f' &&& lapply g') + , \εc -> let εc' = transformMetric h'l εc + εc'' = transformMetric h'r εc + (δc',δc'') = devh εc + in devf εc' ^+^ devg εc'' + ^+^ transformMetric f' δc' + ^+^ transformMetric g' δc'' + ) + _ -> notDefinedHere + where lcofst = linear(,zeroV) + lcosnd = linear(zeroV,) + + + +instance (LinearManifold s v, LocallyScalable s a, RealDimension s) + => AdditiveGroup (RWDfblFuncValue s a v) where + zeroV = point zeroV + (^+^) = grwDfblFnValsCombine $ \a b -> (a^+^b, lPlus, const zeroV) + where lPlus = linear $ uncurry (^+^) + negateV = grwDfblFnValsFunc $ \a -> (negateV a, lNegate, const zeroV) + where lNegate = linear negateV + +instance (RealDimension n, LocallyScalable n a) + => Num (RWDfblFuncValue n a n) where + fromInteger i = point $ fromInteger i + (+) = grwDfblFnValsCombine $ \a b -> (a+b, lPlus, const zeroV) + where lPlus = linear $ uncurry (+) + (*) = grwDfblFnValsCombine $ + \a b -> ( a*b + , linear $ \(da,db) -> a*db + b*da + , \d -> let d¹₂ = sqrt d in (d¹₂,d¹₂) + ) + negate = grwDfblFnValsFunc $ \a -> (negate a, lNegate, const zeroV) + where lNegate = linear negate + abs = (RWDiffable absPW $~) + where absPW a₀ + | a₀<0 = (negativePreRegion, pure desc) + | otherwise = (positivePreRegion, pure asc) + desc = actuallyLinear $ linear negate + asc = actuallyLinear idL + signum = (RWDiffable sgnPW $~) + where sgnPW a₀ + | a₀<0 = (negativePreRegion, pure (const 1)) + | otherwise = (positivePreRegion, pure (const $ -1)) + +instance (RealDimension n, LocallyScalable n a) + => Fractional (RWDfblFuncValue n a n) where + fromRational i = point $ fromRational i + recip = (RWDiffable rcipPW $~) + where rcipPW a₀ + | a₀<0 = (negativePreRegion, pure (Differentiable negp)) + | otherwise = (positivePreRegion, pure (Differentiable posp)) + negp x = (x'¹, (- x'¹^2) *^ idL, dev_ε_δ δ) + -- ε = 1/x − δ/x² − 1/(x+δ) + -- ε·x + ε·δ = 1 + δ/x − δ/x − δ²/x² − 1 + -- = -δ²/x² + -- 0 = δ² + ε·x²·δ + ε·x³ + -- δ = let mph = -ε·x²/2 in mph + sqrt (mph² − ε·x³) + where δ ε = let mph = -ε*x^2/2 in mph + sqrt (mph^2 - ε*x^3) + x'¹ = recip x + posp x = (x'¹, (- x'¹^2) *^ idL, dev_ε_δ δ) + where δ ε = let mph = -ε*x^2/2 in mph + sqrt (mph^2 + ε*x^3) + x'¹ = recip x + + + + + +-- Helper for checking ε-estimations in GHCi with dynamic-plot: +-- epsEst (f,f') εsgn δf (ViewXCenter xc) (ViewHeight h) +-- = let δfxc = δf xc +-- in tracePlot $ reverse [ (xc - δ, f xc - δ * f' xc + εsgn*ε) | +-- ε <- [0, h/500 .. h], let δ = δfxc ε] +-- ++ [ (xc + δ, f xc + δ * f' xc + εsgn*ε) | +-- ε <- [0, h/500 .. h], let δ = δfxc ε] +-- Golfed version: +-- epsEst(f,d)s φ(ViewXCenter ξ)(ViewHeight h)=let ζ=φ ξ in tracePlot$[(ξ-δ,f ξ-δ*d ξ+s*abs ε)|ε<-[-h,-0.998*h..h],let δ=ζ(abs ε)*signum ε] + +instance (RealDimension n, LocallyScalable n a) + => Floating (RWDfblFuncValue n a n) where + pi = point pi + + exp = grwDfblFnValsFunc + $ \x -> let ex = exp x + in ( ex, ex *^ idL, dev_ε_δ $ \ε -> acosh(ε/(2*ex) + 1) ) + -- ε = e^(x+δ) − eˣ − eˣ·δ + -- = eˣ·(e^δ − 1 − δ) + -- ≤ eˣ · (e^δ − 1 + e^(-δ) − 1) + -- = eˣ · 2·(cosh(δ) − 1) + -- cosh(δ) ≥ ε/(2·eˣ) + 1 + -- δ ≥ acosh(ε/(2·eˣ) + 1) + log = (RWDiffable lnRW $~) + where lnRW x | x > 0 = (positivePreRegion, pure (Differentiable lnPosR)) + | otherwise = (negativePreRegion, notDefinedHere) + lnPosR x = ( log x, recip x *^ idL, dev_ε_δ $ \ε -> x * sqrt(1 - exp(-ε)) ) + -- ε = ln x + (-δ)/x − ln(x−δ) + -- = ln (x / ((x−δ) · exp(δ/x))) + -- x/e^ε = (x−δ) · exp(δ/x) + -- let γ = δ/x ∈ [0,1[ + -- exp(-ε) = (1−γ) · e^γ + -- ≥ (1−γ) · (1+γ) + -- = 1 − γ² + -- γ ≥ sqrt(1 − exp(-ε)) + -- δ ≥ x · sqrt(1 − exp(-ε)) + + sqrt = (RWDiffable sqrtRW $~) + where sqrtRW x | x > 0 = (positivePreRegion, pure (Differentiable sqrtPosR)) + | otherwise = (negativePreRegion, notDefinedHere) + sqrtPosR x = ( sx, idL ^/ (2*sx), dev_ε_δ $ + \ε -> 2 * (s2 * sqrt sx^3 * sqrt ε + signum (ε*2-sx) * sx * ε) ) + where sx = sqrt x; s2 = sqrt 2 + -- Exact inverse of O(δ²) remainder. + + sin = grwDfblFnValsFunc sinDfb + where sinDfb x = ( sx, cx *^ idL, dev_ε_δ δ ) + where sx = sin x; cx = cos x + δ ε = let δ₀ = sqrt $ 2 * ε / (abs sx + abs cx/3) + in if δ₀ < 1 -- TODO: confirm selection of δ-definition range. + then δ₀ + else max 1 $ (ε - abs sx - 1) / cos x + -- When sin x ≥ 0, cos x ≥ 0, δ ∈ [0,1[ + -- ε = sin x + δ · cos x − sin(x+δ) + -- = sin x + δ · cos x − sin x · cos δ − cos x · sin δ + -- ≤ sin x + δ · cos x − sin x · (1−δ²/2) − cos x · (δ − δ³/6) + -- = sin x · δ²/2 + cos x · δ³/6 + -- ≤ δ² · (sin x / 2 + cos x / 6) + -- δ ≥ sqrt(2 · ε / (sin x + cos x / 3)) + -- For general δ≥0, + -- ε ≤ δ · cos x + sin x + 1 + -- δ ≥ (ε − sin x − 1) / cos x + cos = sin . (globalDiffable' (actuallyAffine (pi/2) idL) $~) + + sinh x = (exp x - exp (-x))/2 + {- = grwDfblFnValsFunc sinhDfb + where sinhDfb x = ( sx, cx *^ idL, dev_ε_δ δ ) + where sx = sinh x; cx = cosh x + δ ε = undefined -} + -- ε = sinh x + δ · cosh x − sinh(x+δ) + -- = ½ · ( eˣ − e⁻ˣ + δ · (eˣ + e⁻ˣ) − exp(x+δ) + exp(-x−δ) ) + -- = ½·e⁻ˣ · ( e²ˣ − 1 + δ · (e²ˣ + 1) − e²ˣ·e^δ + e^-δ ) + -- = ½ · ( eˣ − e⁻ˣ + δ · (eˣ + e⁻ˣ) − exp(x+δ) + exp(-x−δ) ) + cosh x = (exp x + exp (-x))/2 + tanh x = (exp x - exp (-x)) / (exp x + exp (-x)) + + atan = grwDfblFnValsFunc atanDfb + where atanDfb x = ( atnx, idL ^/ (1+x^2), dev_ε_δ δ ) + where atnx = atan x + c = (atnx*2/pi)^2 + p = 1 + abs x/(2*pi) + δ ε = p * (sqrt ε + ε * c) + -- Semi-empirically obtained: with the epsEst helper, + -- it is observed that this function is (for xc≥0) a lower bound + -- to the arctangent. The growth of the p coefficient makes sense + -- and holds for arbitrarily large xc, because those move us linearly + -- away from the only place where the function is not virtually constant + -- (around 0). + + asin = (RWDiffable asinRW $~) + where asinRW x | x < (-1) = (preRegionFromMinInfTo (-1), notDefinedHere) + | x > 1 = (preRegionToInfFrom 1, notDefinedHere) + | otherwise = (intervalPreRegion (-1,1), pure (Differentiable asinDefdR)) + asinDefdR x = ( asinx, asin'x *^ idL, dev_ε_δ δ ) + where asinx = asin x; asin'x = recip (sqrt $ 1 - x^2) + c = 1 - x^2 + δ ε = sqrt ε * c + -- Empirical, with epsEst upper bound. + + acos = (RWDiffable acosRW $~) + where acosRW x | x < (-1) = (preRegionFromMinInfTo (-1), notDefinedHere) + | x > 1 = (preRegionToInfFrom 1, notDefinedHere) + | otherwise = (intervalPreRegion (-1,1), pure (Differentiable acosDefdR)) + acosDefdR x = ( acosx, acos'x *^ idL, dev_ε_δ δ ) + where acosx = acos x; acos'x = - recip (sqrt $ 1 - x^2) + c = 1 - x^2 + δ ε = sqrt ε * c -- Like for asin – it's just a translation/reflection. + + asinh = grwDfblFnValsFunc asinhDfb + where asinhDfb x = ( asinhx, idL ^/ sqrt(1+x^2), dev_ε_δ δ ) + where asinhx = asinh x + δ ε = abs x * sqrt((1 - exp(-ε))*0.8 + ε^2/(3*abs x)) + sqrt(ε/(abs x+0.5)) + -- Empirical, modified from log function (the area hyperbolic sine + -- resembles two logarithmic lobes), with epsEst-checked lower bound. + + acosh = (RWDiffable acoshRW $~) + where acoshRW x | x > 0 = (positivePreRegion, pure (Differentiable acoshDfb)) + | otherwise = (negativePreRegion, notDefinedHere) + acoshDfb x = ( acosh x, idL ^/ sqrt(x^2 - 2), dev_ε_δ δ ) + where δ ε = (2 - 1/sqrt x) * (s2 * sqrt sx^3 * sqrt(ε/s2) + signum (ε*s2-sx) * sx * ε/s2) + sx = sqrt(x-1) + s2 = sqrt 2 + -- Empirical, modified from sqrt function – the area hyperbolic cosine + -- strongly resembles \x -> sqrt(2 · (x-1)). + + atanh = (RWDiffable atnhRW $~) + where atnhRW x | x < (-1) = (preRegionFromMinInfTo (-1), notDefinedHere) + | x > 1 = (preRegionToInfFrom 1, notDefinedHere) + | otherwise = (intervalPreRegion (-1,1), pure (Differentiable atnhDefdR)) + atnhDefdR x = ( atanh x, recip(1-x^2) *^ idL, dev_ε_δ $ \ε -> sqrt(tanh ε)*(1-abs x) ) + -- Empirical, with epsEst upper bound. + + + + diff --git a/Data/Manifold/Types.hs b/Data/Manifold/Types.hs new file mode 100644 --- /dev/null +++ b/Data/Manifold/Types.hs @@ -0,0 +1,140 @@ +-- | +-- Module : Data.Manifold.Types +-- Copyright : (c) Justus Sagemüller 2015 +-- License : GPL v3 +-- +-- Maintainer : (@) sagemueller $ geo.uni-koeln.de +-- Stability : experimental +-- Portability : portable +-- + + +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE UndecidableInstances #-} +-- {-# LANGUAGE OverlappingInstances #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE TupleSections #-} +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE PatternGuards #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE RecordWildCards #-} + + +module Data.Manifold.Types where + + +import Data.VectorSpace +import Data.AffineSpace +import Data.Basis +import Data.Complex hiding (magnitude) +import Data.Void +import Data.Monoid + +import qualified Prelude + +import Control.Category.Constrained.Prelude hiding ((^)) +import Control.Arrow.Constrained +import Control.Monad.Constrained +import Data.Foldable.Constrained + + + + + + +type EuclidSpace v = (HasBasis v, EqFloating(Scalar v), Eq v) +type EqFloating f = (Eq f, Ord f, Floating f) + + + +data GraphWindowSpec = GraphWindowSpec { + lBound, rBound, bBound, tBound :: Double + , xResolution, yResolution :: Int + } + + + + +data ZeroDim k = Origin deriving(Eq, Show) +instance Monoid (ZeroDim k) where + mempty = Origin + mappend Origin Origin = Origin +instance AdditiveGroup (ZeroDim k) where + zeroV = Origin + Origin ^+^ Origin = Origin + negateV Origin = Origin +instance VectorSpace (ZeroDim k) where + type Scalar (ZeroDim k) = k + _ *^ Origin = Origin +instance HasBasis (ZeroDim k) where + type Basis (ZeroDim k) = Void + basisValue = absurd + decompose Origin = [] + decompose' Origin = absurd + +data S⁰ = PositiveHalfSphere | NegativeHalfSphere deriving(Eq, Show) +newtype S¹ = S¹ { φParamS¹ :: Double -- [-π, π[ + } deriving (Show) +data S² = S² { ϑParamS² :: !Double -- [0, π[ + , φParamS² :: !Double -- [-π, π[ + } deriving (Show) + + +class NaturallyEmbedded m v where + embed :: m -> v + coEmbed :: v -> m + + +instance (VectorSpace y) => NaturallyEmbedded x (x,y) where + embed x = (x, zeroV) + coEmbed (x,_) = x +instance (VectorSpace y, VectorSpace z) => NaturallyEmbedded x ((x,y),z) where + embed x = (embed x, zeroV) + coEmbed (x,_) = coEmbed x + +instance NaturallyEmbedded S⁰ ℝ where + embed PositiveHalfSphere = 1 + embed NegativeHalfSphere = -1 + coEmbed x | x>=0 = PositiveHalfSphere + | otherwise = NegativeHalfSphere +instance NaturallyEmbedded S¹ ℝ² where + embed (S¹ φ) = (cos φ, sin φ) + coEmbed (x,y) = S¹ $ atan2 y x +instance NaturallyEmbedded S² ℝ³ where + embed (S² ϑ φ) = ((cos φ * sin ϑ, sin φ * sin ϑ), cos ϑ) + coEmbed ((x,y),z) = S² (acos $ z/r) (atan2 y x) + where r = sqrt $ x^2 + y^2 + z^2 + + + + + +type Endomorphism a = a->a + + +type ℝ = Double +type ℝ² = (ℝ,ℝ) +type ℝ³ = (ℝ²,ℝ) + +instance VectorSpace () where + type Scalar () = ℝ + _ *^ () = () + +instance HasBasis () where + type Basis () = Void + basisValue = absurd + decompose () = [] + decompose' () = absurd +instance InnerSpace () where + () <.> () = 0 + + + +(^) :: Num a => a -> Int -> a +(^) = (Prelude.^) + diff --git a/manifolds.cabal b/manifolds.cabal --- a/manifolds.cabal +++ b/manifolds.cabal @@ -1,5 +1,5 @@ Name: manifolds -Version: 0.1.0.0 +Version: 0.1.0.2 Category: Math Synopsis: Working with manifolds in a direct, embedding-free way. Description: Manifolds, a generalisation of the notion of \"smooth curves\" or sufaces, @@ -33,6 +33,7 @@ Build-Depends: base>=4.5 && < 6 , transformers , vector-space>=0.8 + , MemoTrie , vector , vector-algorithms , containers @@ -41,7 +42,8 @@ , comonad , semigroups , void - , constrained-categories + , tagged + , constrained-categories >= 0.2 && < 0.3 other-extensions: FlexibleInstances , TypeFamilies , FlexibleContexts @@ -54,8 +56,11 @@ , TupleSections ghc-options: -O2 Exposed-modules: Data.Manifold + Data.Manifold.PseudoAffine + Data.LinearMap.HerMetric -- Data.Manifold.Visualisation.R3.GLUT - Other-modules: Data.List.FastNub + Other-modules: Data.Manifold.Types + Data.List.FastNub Util.Associate Util.LtdShow default-language: Haskell2010