diff --git a/Math/LinearMap/Category.hs b/Math/LinearMap/Category.hs --- a/Math/LinearMap/Category.hs +++ b/Math/LinearMap/Category.hs @@ -76,6 +76,10 @@ -- ** The classes , module Data.VectorSpace , LSpace + , DimensionAware (..) + , Dimensional (..) + , StaticDimensional (..) + , Dimension , TensorSpace (..) , LinearSpace (..) -- ** Orthonormal systems @@ -88,7 +92,7 @@ -- ** Tensors with basis decomposition , (.⊗) -- ** Hilbert space operations - , (·), DualSpace, riesz, coRiesz, showsPrecAsRiesz, (.<) + , (·), DualSpace, riesz, sRiesz, coRiesz, showsPrecAsRiesz, (.<) -- ** Standard decompositions , TensorDecomposable(..), RieszDecomposable(..) , tensorDecomposeShowsPrec, rieszDecomposeShowsPrec @@ -97,6 +101,8 @@ , Num'(..) , Fractional' , RealFrac', RealFloat', LinearShowable + -- ** Coercions + , VSCCoercion(..) -- ** Double-dual, scalar-scalar etc. identity , ClosedScalarWitness(..), TrivialTensorWitness(..) , ScalarSpaceWitness(..), DualSpaceWitness(..), LinearManifoldWitness(..) @@ -115,6 +121,7 @@ import Math.LinearMap.Category.Instances import Math.LinearMap.Category.Instances.Deriving import Math.LinearMap.Asserted +import Math.VectorSpace.DimensionAware import Math.VectorSpace.Docile import Math.LinearMap.Category.TensorQuot @@ -417,8 +424,9 @@ -- @ -- v '<.>^' (w |&> 'euclideanNorm') ≡ v '<.>' w -- @ -(|&>) :: LSpace v => DualVector v -> Variance v -> v -dv |&> Norm m = symVSC coerceDoubleDual $ m-+$>dv +(|&>) :: ∀ v . LSpace v => DualVector v -> Variance v -> v +dv |&> Norm m = case dualSpaceWitness @v of + DualSpaceWitness -> m-+$>dv -- | 'spanNorm' / 'spanVariance' are inefficient if the number of vectors diff --git a/Math/LinearMap/Category/Class.hs b/Math/LinearMap/Category/Class.hs --- a/Math/LinearMap/Category/Class.hs +++ b/Math/LinearMap/Category/Class.hs @@ -1,6 +1,6 @@ -- | -- Module : Math.LinearMap.Category.Class --- Copyright : (c) Justus Sagemüller 2016 +-- Copyright : (c) Justus Sagemüller 2016-2022 -- License : GPL v3 -- -- Maintainer : (@) jsag $ hvl.no @@ -13,10 +13,13 @@ {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE FunctionalDependencies #-} {-# LANGUAGE TypeOperators #-} +{-# LANGUAGE NoStarIsType #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE TypeApplications #-} {-# LANGUAGE Rank2Types #-} {-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE InstanceSigs #-} {-# LANGUAGE PatternSynonyms #-} {-# LANGUAGE ViewPatterns #-} {-# LANGUAGE UnicodeSyntax #-} @@ -24,39 +27,60 @@ {-# LANGUAGE StandaloneDeriving #-} {-# LANGUAGE DeriveGeneric #-} {-# LANGUAGE GADTs #-} +{-# LANGUAGE DataKinds #-} {-# LANGUAGE DefaultSignatures #-} {-# LANGUAGE CPP #-} module Math.LinearMap.Category.Class where import Data.VectorSpace +import Math.VectorSpace.DimensionAware import Data.AffineSpace import Prelude () import qualified Prelude as Hask -import Control.Category.Constrained.Prelude +import Control.Category.Constrained.Prelude hiding (type (+)) import Control.Arrow.Constrained import Data.Coerce import Data.Type.Coercion import Data.Tagged +import Data.Proxy(Proxy(..)) +import qualified Data.Vector.Generic as GArr + import Math.Manifold.Core.PseudoAffine import Math.LinearMap.Asserted import Math.VectorSpace.ZeroDimensional import Data.VectorSpace.Free +import Control.Monad.ST (ST) + +import Data.Singletons (sing, withSingI) +#if MIN_VERSION_singletons(3,0,0) +import Prelude.Singletons (SNum(..)) +import Data.Maybe.Singletons (SMaybe(..)) +import GHC.TypeLits.Singletons (withKnownNat, SNat(..)) +#else +import Data.Singletons.Prelude.Num (SNum(..)) +import Data.Singletons.Prelude.Maybe (SMaybe(..)) +import Data.Singletons.TypeLits (withKnownNat, SNat(..)) +#endif import Data.Kind (Type) +import GHC.TypeLits (Nat, type (+), type (*), KnownNat, natVal) import qualified GHC.Generics as Gnrx import GHC.Generics (Generic, (:*:)((:*:))) +import qualified Math.VectorSpace.DimensionAware.Theorems.MaybeNat as Maybe + data ClosedScalarWitness s where ClosedScalarWitness :: (Scalar s ~ s, DualVector s ~ s) => ClosedScalarWitness s data TrivialTensorWitness s w where TrivialTensorWitness :: w ~ TensorProduct s w => TrivialTensorWitness s w -class (Num s, LinearSpace s, FreeVectorSpace s) => Num' s where +class (Num s, LinearSpace s, FreeVectorSpace s, 1`Dimensional`s) + => Num' s where closedScalarWitness :: ClosedScalarWitness s default closedScalarWitness :: (Scalar s ~ s, DualVector s ~ s) => ClosedScalarWitness s closedScalarWitness = ClosedScalarWitness @@ -80,43 +104,52 @@ -- vector spaces, under the requirement that they internally use the same basis -- (if any). Note that this does not mean they also need to have the same inner -- product / dual space. -data VSCCoercion a b where - VSCCoercion :: Coercible a b - => VSCCoercion a b +data VSCCoercion s a b where + VSCCoercion :: (Coercible a b, StaticDimension a ~ StaticDimension b) + => VSCCoercion s a b -getVSCCoercion :: VSCCoercion a b -> Coercion a b +getVSCCoercion :: VSCCoercion s a b -> Coercion a b getVSCCoercion VSCCoercion = Coercion -symVSC :: VSCCoercion a b -> VSCCoercion b a +symVSC :: VSCCoercion s a b -> VSCCoercion s b a symVSC VSCCoercion = VSCCoercion -firstVSC :: VSCCoercion a b -> VSCCoercion (a,c) (b,c) +firstVSC :: VSCCoercion s a b -> VSCCoercion s (a,c) (b,c) firstVSC VSCCoercion = VSCCoercion -secondVSC :: VSCCoercion a b -> VSCCoercion (c,a) (c,b) +secondVSC :: VSCCoercion s a b -> VSCCoercion s (c,a) (c,b) secondVSC VSCCoercion = VSCCoercion -unsafeFollowVSC :: Coercible a b => c a b -> VSCCoercion a b +unsafeFollowVSC :: (Coercible a b, StaticDimension a ~ StaticDimension b) + => c a b -> VSCCoercion s a b unsafeFollowVSC _ = VSCCoercion -unsafeFloutVSC :: Coercible a b => c b a -> VSCCoercion a b +unsafeFloutVSC :: (Coercible a b, StaticDimension a ~ StaticDimension b) + => c b a -> VSCCoercion s a b unsafeFloutVSC _ = VSCCoercion -instance Category VSCCoercion where +instance Category (VSCCoercion s) where + type Object (VSCCoercion s) v = (TensorSpace v, Scalar v ~ s) id = VSCCoercion VSCCoercion . VSCCoercion = VSCCoercion -instance EnhancedCat Coercion VSCCoercion where +instance EnhancedCat Coercion (VSCCoercion s) where arr = getVSCCoercion -instance EnhancedCat (->) VSCCoercion where +instance EnhancedCat (->) (VSCCoercion s) where arr VSCCoercion x = coerce x + +infixr 0 -+$=> +(-+$=>) :: VSCCoercion s a b -> a -> b +VSCCoercion -+$=> x = coerce x + +class (DimensionAware v, PseudoAffine v) => TensorSpace v where -class (VectorSpace v, PseudoAffine v) => TensorSpace v where -- | The internal representation of a 'Tensor' product. -- -- For Euclidean spaces, this is generally constructed by replacing each @s@ -- scalar field in the @v@ vector with an entire @w@ vector. I.e., you have -- then a “nested vector” or, if @v@ is a @DualVector@ / “row vector”, a matrix. type TensorProduct v w :: Type + scalarSpaceWitness :: ScalarSpaceWitness v linearManifoldWitness :: LinearManifoldWitness v zeroTensor :: (TensorSpace w, Scalar w ~ Scalar v) @@ -155,8 +188,21 @@ fzipTensorWith :: ( TensorSpace u, TensorSpace w, TensorSpace x , Scalar u ~ Scalar v, Scalar w ~ Scalar v, Scalar x ~ Scalar v ) => Bilinear ((w,x) -+> u) (v⊗w, v⊗x) (v⊗u) - coerceFmapTensorProduct :: Hask.Functor p - => p v -> VSCCoercion a b -> VSCCoercion (TensorProduct v a) (TensorProduct v b) + tensorUnsafeFromArrayWithOffset :: ∀ w α n m + . ( n`Dimensional`v + , TensorSpace w, m`Dimensional`w, Scalar w ~ Scalar v + , GArr.Vector α (Scalar v) ) + => Int -> α (Scalar v) -> (v⊗w) + tensorUnsafeWriteArrayWithOffset :: ∀ w α σ n m + . ( n`Dimensional`v + , TensorSpace w, m`Dimensional`w, Scalar w ~ Scalar v + , GArr.Vector α (Scalar v) ) + => GArr.Mutable α σ (Scalar v) -> Int -> (v⊗w) -> ST σ () + coerceFmapTensorProduct :: ( Hask.Functor p + , TensorSpace a, Scalar a ~ Scalar v + , TensorSpace b, Scalar b ~ Scalar v ) + => p v -> VSCCoercion (Scalar v) a b + -> Coercion (TensorProduct v a) (TensorProduct v b) -- | “Sanity-check” a vector. This typically amounts to detecting any NaN components, -- which should trigger a @Nothing@ result. Otherwise, the result should be @Just@ -- the input, but may also be optimised / memoised if applicable (i.e. for @@ -176,7 +222,8 @@ data DualSpaceWitness v where DualSpaceWitness :: ( LinearSpace (Scalar v), DualVector (Scalar v) ~ Scalar v , LinearSpace (DualVector v), Scalar (DualVector v) ~ Scalar v - , DualVector (DualVector v) ~ v ) + , DualVector (DualVector v) ~ v + , StaticDimension (DualVector v) ~ StaticDimension v ) => DualSpaceWitness v -- | The class of vector spaces @v@ for which @'LinearMap' s v w@ is well-implemented. @@ -214,7 +261,9 @@ , dualSpaceWitness :: DualSpaceWitness v ) of (ScalarSpaceWitness,DualSpaceWitness) -> arr asTensor >>> fromFlatTensor - coerceDoubleDual :: VSCCoercion v (DualVector (DualVector v)) + -- | This will probably be removed in the future, since infinite-dimensional + -- (e.g. Banach-) spaces may be not isomorphic to their double dual. + coerceDoubleDual :: VSCCoercion (Scalar v) v (DualVector (DualVector v)) coerceDoubleDual = case dualSpaceWitness :: DualSpaceWitness v of DualSpaceWitness -> VSCCoercion @@ -278,6 +327,14 @@ DualSpaceWitness -> bilinearFunction $ \f -> arr asTensor >>> getLinearFunction (fmapTensor-+$>f) >>> arr fromTensor + +instance DimensionAware (ZeroDim s) where + type StaticDimension (ZeroDim s) = 'Just 0 + dimensionalityWitness = IsStaticDimensional +instance 0`Dimensional`ZeroDim s where + unsafeFromArrayWithOffset _ _ = Origin + unsafeWriteArrayWithOffset _ _ _ = return () + instance Num' s => TensorSpace (ZeroDim s) where type TensorProduct (ZeroDim s) v = ZeroDim s scalarSpaceWitness = case closedScalarWitness :: ClosedScalarWitness s of @@ -297,7 +354,9 @@ transposeTensor = const0 fmapTensor = biConst0 fzipTensorWith = biConst0 - coerceFmapTensorProduct _ VSCCoercion = VSCCoercion + tensorUnsafeFromArrayWithOffset _ _ = Tensor Origin + tensorUnsafeWriteArrayWithOffset _ _ (Tensor Origin) = return () + coerceFmapTensorProduct _ VSCCoercion = Coercion wellDefinedVector Origin = Just Origin wellDefinedTensor (Tensor Origin) = Just (Tensor Origin) instance Num' s => LinearSpace (ZeroDim s) where @@ -346,27 +405,34 @@ -- linear mappings, but they also form a useful vector space on their own right. newtype Tensor s v w = Tensor {getTensorProduct :: TensorProduct v w} -asTensor :: VSCCoercion (LinearMap s v w) (Tensor s (DualVector v) w) -asTensor = VSCCoercion -fromTensor :: VSCCoercion (Tensor s (DualVector v) w) (LinearMap s v w) -fromTensor = VSCCoercion +asTensor :: ∀ s v w . LinearSpace v + => VSCCoercion s (LinearMap s v w) (Tensor s (DualVector v) w) +asTensor = case dualSpaceWitness @v of + DualSpaceWitness -> VSCCoercion +fromTensor :: ∀ s v w . LinearSpace v + => VSCCoercion s (Tensor s (DualVector v) w) (LinearMap s v w) +fromTensor = case dualSpaceWitness @v of + DualSpaceWitness -> VSCCoercion asLinearMap :: ∀ s v w . (LinearSpace v, Scalar v ~ s) - => VSCCoercion (Tensor s v w) (LinearMap s (DualVector v) w) + => VSCCoercion s (Tensor s v w) (LinearMap s (DualVector v) w) asLinearMap = case dualSpaceWitness :: DualSpaceWitness v of DualSpaceWitness -> VSCCoercion fromLinearMap :: ∀ s v w . (LinearSpace v, Scalar v ~ s) - => VSCCoercion (LinearMap s (DualVector v) w) (Tensor s v w) + => VSCCoercion s (LinearMap s (DualVector v) w) (Tensor s v w) fromLinearMap = case dualSpaceWitness :: DualSpaceWitness v of DualSpaceWitness -> VSCCoercion -pseudoFmapTensorLHS :: (TensorProduct v w ~ TensorProduct v' w) - => c v v' -> VSCCoercion (Tensor s v w) (Tensor s v' w) +pseudoFmapTensorLHS :: ( TensorProduct v w ~ TensorProduct v' w + , StaticDimension v ~ StaticDimension v' ) + => c v v' -> VSCCoercion s (Tensor s v w) (Tensor s v' w) pseudoFmapTensorLHS _ = VSCCoercion -pseudoPrecomposeLinmap :: (TensorProduct (DualVector v) w ~ TensorProduct (DualVector v') w) - => c v' v -> VSCCoercion (LinearMap s v w) (LinearMap s v' w) +pseudoPrecomposeLinmap + :: ( TensorProduct (DualVector v) w ~ TensorProduct (DualVector v') w + , StaticDimension v ~ StaticDimension v' ) + => c v' v -> VSCCoercion s (LinearMap s v w) (LinearMap s v' w) pseudoPrecomposeLinmap _ = VSCCoercion envTensorLHSCoercion :: ( TensorProduct v w ~ TensorProduct v' w @@ -505,8 +571,8 @@ (DualSpaceWitness, DualSpaceWitness, DualSpaceWitness) -> fromTensor $ (fzipTensorWith$id) $ (asTensor $ f, asTensor $ g) terminal = zeroV - fst = sampleLinearFunction $ fst - snd = sampleLinearFunction $ snd + fst = id ⊕ zeroV + snd = zeroV ⊕ id instance Num' s => EnhancedCat (->) (LinearMap s) where arr m = arr $ applyLinear $ m instance Num' s => EnhancedCat (LinearFunction s) (LinearMap s) where @@ -564,10 +630,42 @@ $ \f (Tensor (uw, vw), Tensor (ux, vx)) -> Tensor ( (fzipTensorWith-+$>f)-+$>(uw,ux) , (fzipTensorWith-+$>f)-+$>(vw,vx) ) + tensorUnsafeFromArrayWithOffset :: ∀ nm w o α + . ( nm`Dimensional`(u,v) + , TensorSpace w, o`Dimensional`w, Scalar w ~ Scalar v + , GArr.Vector α (Scalar u) ) + => Int -> α (Scalar u) -> ((u,v)⊗w) + tensorUnsafeFromArrayWithOffset + = case ( staticDimensionSing @u, dimensionalityWitness @u + , staticDimensionSing @v, dimensionalityWitness @v ) of + ( SJust sn, IsStaticDimensional + ,SJust sm, IsStaticDimensional ) + -> let sno = sn %* dimensionalitySing @w + smo = sm %* dimensionalitySing @w + in withKnownNat sno (withKnownNat smo ( + \i arr -> Tensor ( unsafeFromArrayWithOffset i arr + , unsafeFromArrayWithOffset + (i + fromIntegral (natVal sno)) arr ))) + tensorUnsafeWriteArrayWithOffset :: ∀ nm w o α σ + . ( nm`Dimensional`(u,v) + , TensorSpace w, o`Dimensional`w, Scalar w ~ Scalar v + , GArr.Vector α (Scalar u) ) + => GArr.Mutable α σ (Scalar u) -> Int -> ((u,v)⊗w) -> ST σ () + tensorUnsafeWriteArrayWithOffset + = case ( staticDimensionSing @u, dimensionalityWitness @u + , staticDimensionSing @v, dimensionalityWitness @v ) of + ( SJust sn, IsStaticDimensional + ,SJust sm, IsStaticDimensional ) + -> let sno = sn %* dimensionalitySing @w + smo = sm %* dimensionalitySing @w + in withKnownNat sno (withKnownNat smo ( + \arr i (Tensor (x,y)) -> do + unsafeWriteArrayWithOffset arr i x + unsafeWriteArrayWithOffset arr (i + fromIntegral (natVal sno)) y )) coerceFmapTensorProduct p cab = case ( coerceFmapTensorProduct (fst<$>p) cab , coerceFmapTensorProduct (snd<$>p) cab ) of - (VSCCoercion, VSCCoercion) -> VSCCoercion + (Coercion, Coercion) -> Coercion wellDefinedVector (u,v) = liftA2 (,) (wellDefinedVector u) (wellDefinedVector v) wellDefinedTensor (Tensor (u,v)) = liftA2 ((Tensor.) . (,)) (wellDefinedTensor u) (wellDefinedTensor v) @@ -619,12 +717,13 @@ -> bilinearFunction $ \f (LinearMap (fu, fv)) -> ((composeLinear-+$>f)-+$>asLinearMap $ fu) ⊕ ((composeLinear-+$>f)-+$>asLinearMap $ fv) - applyTensorFunctional = case ( dualSpaceWitness :: DualSpaceWitness u - , dualSpaceWitness :: DualSpaceWitness v ) of + applyTensorFunctional = case ( dualSpaceWitness @u, dualSpaceWitness @v ) of (DualSpaceWitness, DualSpaceWitness) -> bilinearFunction $ \(LinearMap (fu,fv)) (Tensor (tu,tv)) - -> ((applyTensorFunctional-+$>asLinearMap$fu)-+$>tu) - + ((applyTensorFunctional-+$>asLinearMap$fv)-+$>tv) + -> ((applyTensorFunctional + -+$>getVSCCoercion asLinearMap$fu)-+$>tu) + + ((applyTensorFunctional + -+$>getVSCCoercion asLinearMap$fv)-+$>tv) applyTensorLinMap = case ( dualSpaceWitness :: DualSpaceWitness u , dualSpaceWitness :: DualSpaceWitness v ) of (DualSpaceWitness, DualSpaceWitness) -> bilinearFunction`id` @@ -632,6 +731,8 @@ in ( (applyTensorLinMap-+$>uncurryLinearMap.asLinearMap $ fu)-+$>tu ) ^+^ ( (applyTensorLinMap-+$>uncurryLinearMap.asLinearMap $ fv)-+$>tv ) useTupleLinearSpaceComponents r = r + coerceDoubleDual = case ( dualSpaceWitness @u, dualSpaceWitness @v ) of + (DualSpaceWitness, DualSpaceWitness) -> VSCCoercion lfstBlock :: ( LSpace u, LSpace v, LSpace w , Scalar u ~ Scalar v, Scalar v ~ Scalar w ) @@ -644,34 +745,81 @@ -- | @((v'⊗w)+>x) -> ((v+>w)+>x) -argFromTensor :: ∀ s v w x . (LinearSpace v, LinearSpace w, Scalar v ~ s, Scalar w ~ s) - => VSCCoercion (LinearMap s (Tensor s (DualVector v) w) x) +argFromTensor :: ∀ s v w x . ( LinearSpace v, LinearSpace w + , Scalar v ~ s, Scalar w ~ s + , TensorSpace x, Scalar x ~ s + ) + => VSCCoercion s (LinearMap s (Tensor s (DualVector v) w) x) (LinearMap s (LinearMap s v w) x) argFromTensor = case dualSpaceWitness :: DualSpaceWitness v of DualSpaceWitness -> curryLinearMap >>> fromLinearMap >>> coUncurryLinearMap -- | @((v+>w)+>x) -> ((v'⊗w)+>x)@ -argAsTensor :: ∀ s v w x . (LinearSpace v, LinearSpace w, Scalar v ~ s, Scalar w ~ s) - => VSCCoercion (LinearMap s (LinearMap s v w) x) +argAsTensor :: ∀ s v w x . ( LinearSpace v, LinearSpace w + , Scalar v ~ s, Scalar w ~ s + , TensorSpace x, Scalar x ~ s + ) + => VSCCoercion s (LinearMap s (LinearMap s v w) x) (LinearMap s (Tensor s (DualVector v) w) x) argAsTensor = case dualSpaceWitness :: DualSpaceWitness v of DualSpaceWitness -> uncurryLinearMap <<< asLinearMap <<< coCurryLinearMap +tensorDimensionAssoc :: ∀ u v w s r . (TensorSpace u, TensorSpace v, TensorSpace w) + => (( Maybe.ZipWithTimes (Maybe.ZipWithTimes (StaticDimension u) (StaticDimension v)) + (StaticDimension w) + ~ Maybe.ZipWithTimes (StaticDimension u) + (Maybe.ZipWithTimes (StaticDimension v) (StaticDimension w)) + ) => r) -> r +tensorDimensionAssoc φ + = Maybe.zipWithTimesAssoc (staticDimensionSing @u) + (staticDimensionSing @v) + (staticDimensionSing @w) φ + -- | @(u+>(v⊗w)) -> (u+>v)⊗w@ -deferLinearMap :: VSCCoercion (LinearMap s u (Tensor s v w)) (Tensor s (LinearMap s u v) w) -deferLinearMap = VSCCoercion +deferLinearMap :: ∀ s u v w . (TensorSpace u, TensorSpace v, TensorSpace w) + => VSCCoercion s (LinearMap s u (Tensor s v w)) (Tensor s (LinearMap s u v) w) +deferLinearMap + = tensorDimensionAssoc @u @v @w VSCCoercion -- | @(u+>v)⊗w -> u+>(v⊗w)@ -hasteLinearMap :: VSCCoercion (Tensor s (LinearMap s u v) w) (LinearMap s u (Tensor s v w)) -hasteLinearMap = VSCCoercion +hasteLinearMap :: ∀ s u v w . (TensorSpace u, TensorSpace v, TensorSpace w) + => VSCCoercion s (Tensor s (LinearMap s u v) w) (LinearMap s u (Tensor s v w)) +hasteLinearMap = tensorDimensionAssoc @u @v @w VSCCoercion -lassocTensor :: VSCCoercion (Tensor s u (Tensor s v w)) (Tensor s (Tensor s u v) w) -lassocTensor = VSCCoercion -rassocTensor :: VSCCoercion (Tensor s (Tensor s u v) w) (Tensor s u (Tensor s v w)) -rassocTensor = VSCCoercion +lassocTensor :: ∀ s u v w . (TensorSpace u, TensorSpace v, TensorSpace w) + => VSCCoercion s (Tensor s u (Tensor s v w)) (Tensor s (Tensor s u v) w) +lassocTensor = tensorDimensionAssoc @u @v @w VSCCoercion +rassocTensor :: ∀ s u v w . (TensorSpace u, TensorSpace v, TensorSpace w) + => VSCCoercion s (Tensor s (Tensor s u v) w) (Tensor s u (Tensor s v w)) +rassocTensor = tensorDimensionAssoc @u @v @w VSCCoercion + instance ∀ s u v . ( LinearSpace u, TensorSpace v, Scalar u ~ s, Scalar v ~ s ) + => DimensionAware (LinearMap s u v) where + type StaticDimension (LinearMap s u v) + = Maybe.ZipWithTimes (StaticDimension u) (StaticDimension v) + dimensionalityWitness = case (dimensionalityWitness @u, dimensionalityWitness @v) of + (IsStaticDimensional, IsStaticDimensional) + -> withKnownNat (dimensionalitySing @u %* dimensionalitySing @v) + IsStaticDimensional + (IsFlexibleDimensional, _) -> IsFlexibleDimensional + (_, IsFlexibleDimensional) -> IsFlexibleDimensional +instance ∀ s n u m v nm . ( n`Dimensional`u, m`Dimensional`v + , LinearSpace u, TensorSpace v, Scalar u ~ s, Scalar v ~ s + , nm ~ (n*m) ) + => nm`Dimensional`(LinearMap s u v) where + knownDimensionalitySing = dimensionalitySing @u %* dimensionalitySing @v + unsafeFromArrayWithOffset i arr = case dualSpaceWitness @u of + DualSpaceWitness -> case dimensionalityWitness @(DualVector u) of + IsStaticDimensional + -> fromTensor $ unsafeFromArrayWithOffset i arr + unsafeWriteArrayWithOffset arr i lm = case dualSpaceWitness @u of + DualSpaceWitness -> case dimensionalityWitness @(DualVector u) of + IsStaticDimensional + -> unsafeWriteArrayWithOffset arr i (asTensor $ lm) + +instance ∀ s u v . ( LinearSpace u, TensorSpace v, Scalar u ~ s, Scalar v ~ s ) => TensorSpace (LinearMap s u v) where type TensorProduct (LinearMap s u v) w = TensorProduct (DualVector u) (Tensor s v w) scalarSpaceWitness = case ( scalarSpaceWitness :: ScalarSpaceWitness u @@ -724,13 +872,43 @@ ScalarSpaceWitness -> LinearFunction $ \f -> arr deferLinearMap <<< fzipWith (fzipWith f) <<< arr hasteLinearMap *** arr hasteLinearMap + tensorUnsafeFromArrayWithOffset :: ∀ nm w o α + . ( nm`Dimensional`LinearMap s u v + , TensorSpace w, o`Dimensional`w, Scalar w ~ s + , GArr.Vector α s ) + => Int -> α s -> (LinearMap s u v⊗w) + tensorUnsafeFromArrayWithOffset + = case ( dimensionalityWitness @u, staticDimensionSing @u + , dimensionalityWitness @v, staticDimensionSing @v ) of + ( IsStaticDimensional, SJust sn + ,IsStaticDimensional, SJust sm ) + -> withKnownNat (sm%*dimensionalitySing @w) ( + withKnownNat (sn%*(sm%*dimensionalitySing @w)) ( + \i -> arr (deferLinearMap @s @u @v @w) + . unsafeFromArrayWithOffset i)) + tensorUnsafeWriteArrayWithOffset :: ∀ nm w o α σ + . ( nm`Dimensional`LinearMap s u v + , TensorSpace w, o`Dimensional`w, Scalar w ~ s + , GArr.Vector α s ) + => GArr.Mutable α σ s -> Int -> (LinearMap s u v⊗w) -> ST σ () + tensorUnsafeWriteArrayWithOffset + = case ( dimensionalityWitness @u, staticDimensionSing @u + , dimensionalityWitness @v, staticDimensionSing @v ) of + ( IsStaticDimensional, SJust sn + ,IsStaticDimensional, SJust sm ) + -> withKnownNat (sm%*dimensionalitySing @w) ( + withKnownNat (sn%*(sm%*dimensionalitySing @w)) ( + \ar i -> unsafeWriteArrayWithOffset ar i + . arr (hasteLinearMap @s @u @v @w) )) coerceFmapTensorProduct = cftlp dualSpaceWitness - where cftlp :: ∀ a b p . DualSpaceWitness u -> p (LinearMap s u v) -> VSCCoercion a b - -> VSCCoercion (TensorProduct (DualVector u) (Tensor s v a)) + where cftlp :: ∀ a b p . ( TensorSpace a, Scalar a ~ s + , TensorSpace b, Scalar b ~ s ) + => DualSpaceWitness u -> p (LinearMap s u v) -> VSCCoercion s a b + -> Coercion (TensorProduct (DualVector u) (Tensor s v a)) (TensorProduct (DualVector u) (Tensor s v b)) cftlp DualSpaceWitness _ c = coerceFmapTensorProduct ([]::[DualVector u]) - (fmap c :: VSCCoercion (v⊗a) (v⊗b)) + (fmap c :: VSCCoercion s (v⊗a) (v⊗b)) wellDefinedVector = case dualSpaceWitness :: DualSpaceWitness u of DualSpaceWitness -> arr asTensor >>> wellDefinedTensor >>> arr (fmap (getVSCCoercion fromTensor)) wellDefinedTensor @@ -738,8 +916,9 @@ -- | @((u+>v)+>w) -> u⊗(v+>w)@ coCurryLinearMap :: ∀ s u v w . ( LinearSpace u, Scalar u ~ s - , LinearSpace v, Scalar v ~ s ) => - VSCCoercion (LinearMap s (LinearMap s u v) w) (Tensor s u (LinearMap s v w)) + , LinearSpace v, Scalar v ~ s + , TensorSpace w, Scalar w ~ s ) => + VSCCoercion s (LinearMap s (LinearMap s u v) w) (Tensor s u (LinearMap s v w)) coCurryLinearMap = case ( dualSpaceWitness :: DualSpaceWitness u , dualSpaceWitness :: DualSpaceWitness v ) of (DualSpaceWitness, DualSpaceWitness) @@ -747,26 +926,33 @@ -- | @(u⊗(v+>w)) -> (u+>v)+>w@ coUncurryLinearMap :: ∀ s u v w . ( LinearSpace u, Scalar u ~ s - , LinearSpace v, Scalar v ~ s ) => - VSCCoercion (Tensor s u (LinearMap s v w)) (LinearMap s (LinearMap s u v) w) + , LinearSpace v, Scalar v ~ s + , TensorSpace w, Scalar w ~ s ) => + VSCCoercion s (Tensor s u (LinearMap s v w)) (LinearMap s (LinearMap s u v) w) coUncurryLinearMap = case ( dualSpaceWitness :: DualSpaceWitness u , dualSpaceWitness :: DualSpaceWitness v ) of (DualSpaceWitness, DualSpaceWitness) -> fromTensor <<< lassocTensor <<< fmap fromLinearMap -- | @((u⊗v)+>w) -> (u+>(v+>w))@ -curryLinearMap :: ∀ u v w s . ( LinearSpace u, Scalar u ~ s ) - => VSCCoercion (LinearMap s (Tensor s u v) w) (LinearMap s u (LinearMap s v w)) -curryLinearMap = case dualSpaceWitness :: DualSpaceWitness u of - DualSpaceWitness -> (VSCCoercion :: VSCCoercion ((u⊗v)+>w) +curryLinearMap :: ∀ u v w s . ( LinearSpace u, LinearSpace v, TensorSpace w + , Scalar u ~ s , Scalar v ~ s , Scalar w ~ s ) + => VSCCoercion s (LinearMap s (Tensor s u v) w) (LinearMap s u (LinearMap s v w)) +curryLinearMap = case (dualSpaceWitness @u, dualSpaceWitness @v) of + (DualSpaceWitness, DualSpaceWitness) + -> tensorDimensionAssoc @u @(DualVector v) @w + (VSCCoercion :: VSCCoercion s ((u⊗v)+>w) ((DualVector u)⊗(Tensor s (DualVector v) w)) ) >>> fmap fromTensor >>> fromTensor -- | @(u+>(v+>w)) -> ((u⊗v)+>w)@ -uncurryLinearMap :: ∀ u v w s . ( LinearSpace u, Scalar u ~ s ) - => VSCCoercion (LinearMap s u (LinearMap s v w)) (LinearMap s (Tensor s u v) w) -uncurryLinearMap = case dualSpaceWitness :: DualSpaceWitness u of - DualSpaceWitness -> (VSCCoercion :: VSCCoercion +uncurryLinearMap :: ∀ u v w s . ( LinearSpace u, LinearSpace v, TensorSpace w + , Scalar u ~ s , Scalar v ~ s , Scalar w ~ s ) + => VSCCoercion s (LinearMap s u (LinearMap s v w)) (LinearMap s (Tensor s u v) w) +uncurryLinearMap = case (dualSpaceWitness @u, dualSpaceWitness @v) of + (DualSpaceWitness, DualSpaceWitness) + -> tensorDimensionAssoc @u @(DualVector v) @w + (VSCCoercion :: VSCCoercion s ((DualVector u)⊗(Tensor s (DualVector v) w)) ((u⊗v)+>w) ) <<< fmap asTensor <<< asTensor @@ -783,8 +969,9 @@ dualSpaceWitness = case ( dualSpaceWitness :: DualSpaceWitness u , dualSpaceWitness :: DualSpaceWitness v ) of (DualSpaceWitness, DualSpaceWitness) -> DualSpaceWitness - linearId = case dualSpaceWitness :: DualSpaceWitness u of - DualSpaceWitness -> fromTensor . lassocTensor . fromLinearMap . fmap asTensor + linearId = case (dualSpaceWitness @u, dualSpaceWitness @v) of + (DualSpaceWitness, DualSpaceWitness) + -> fromTensor . lassocTensor . fromLinearMap . fmap asTensor . curryLinearMap . fmap fromTensor $ tensorId tensorId = uncurryLinearMap . coUncurryLinearMap . fmap curryLinearMap . coCurryLinearMap . fmap deferLinearMap $ id @@ -813,7 +1000,26 @@ . arr (asTensor . hasteLinearMap) -+$> g useTupleLinearSpaceComponents _ = usingNonTupleTypeAsTupleError + instance ∀ s u v . (TensorSpace u, TensorSpace v, Scalar u ~ s, Scalar v ~ s) + => DimensionAware (Tensor s u v) where + type StaticDimension (Tensor s u v) + = Maybe.ZipWithTimes (StaticDimension u) (StaticDimension v) + dimensionalityWitness = case (dimensionalityWitness @u, dimensionalityWitness @v) of + (IsStaticDimensional, IsStaticDimensional) + -> withKnownNat (dimensionalitySing @u %* dimensionalitySing @v) + IsStaticDimensional + (IsFlexibleDimensional, _) -> IsFlexibleDimensional + (_, IsFlexibleDimensional) -> IsFlexibleDimensional +instance ∀ s n u m v nm . ( n`Dimensional`u, m`Dimensional`v + , TensorSpace u, TensorSpace v, Scalar u ~ s, Scalar v ~ s + , nm ~ (n*m) ) + => nm`Dimensional`(Tensor s u v) where + knownDimensionalitySing = dimensionalitySing @u %* dimensionalitySing @v + unsafeFromArrayWithOffset = tensorUnsafeFromArrayWithOffset + unsafeWriteArrayWithOffset = tensorUnsafeWriteArrayWithOffset + +instance ∀ s u v . (TensorSpace u, TensorSpace v, Scalar u ~ s, Scalar v ~ s) => TensorSpace (Tensor s u v) where type TensorProduct (Tensor s u v) w = TensorProduct u (Tensor s v w) scalarSpaceWitness = case ( scalarSpaceWitness :: ScalarSpaceWitness u @@ -858,12 +1064,41 @@ ScalarSpaceWitness -> LinearFunction $ \f -> arr lassocTensor <<< fzipWith (fzipWith f) <<< arr rassocTensor *** arr rassocTensor - coerceFmapTensorProduct = cftlp - where cftlp :: ∀ a b p . p (Tensor s u v) -> VSCCoercion a b - -> VSCCoercion (TensorProduct u (Tensor s v a)) + tensorUnsafeFromArrayWithOffset :: ∀ nm w o α + . ( nm`Dimensional`Tensor s u v + , TensorSpace w, o`Dimensional`w, Scalar w ~ s + , GArr.Vector α s ) + => Int -> α s -> (Tensor s u v⊗w) + tensorUnsafeFromArrayWithOffset + = case ( dimensionalityWitness @u, staticDimensionSing @u + , dimensionalityWitness @v, staticDimensionSing @v ) of + ( IsStaticDimensional, SJust sn + ,IsStaticDimensional, SJust sm ) + -> withKnownNat (sm%*dimensionalitySing @w) ( + withKnownNat (sn%*(sm%*dimensionalitySing @w)) ( + \i -> arr (lassocTensor @s @u @v @w) + . unsafeFromArrayWithOffset i)) + tensorUnsafeWriteArrayWithOffset :: ∀ nm w o α σ + . ( nm`Dimensional`Tensor s u v + , TensorSpace w, o`Dimensional`w, Scalar w ~ s + , GArr.Vector α s ) + => GArr.Mutable α σ s -> Int -> (Tensor s u v⊗w) -> ST σ () + tensorUnsafeWriteArrayWithOffset + = case ( dimensionalityWitness @u, staticDimensionSing @u + , dimensionalityWitness @v, staticDimensionSing @v ) of + ( IsStaticDimensional, SJust sn + ,IsStaticDimensional, SJust sm ) + -> withKnownNat (sm%*dimensionalitySing @w) ( + withKnownNat (sn%*(sm%*dimensionalitySing @w)) ( + \ar i -> unsafeWriteArrayWithOffset ar i + . arr (rassocTensor @s @u @v @w) )) + coerceFmapTensorProduct :: ∀ a b p . ( TensorSpace a, Scalar a ~ s + , TensorSpace b, Scalar b ~ s ) + => p (Tensor s u v) -> VSCCoercion s a b + -> Coercion (TensorProduct u (Tensor s v a)) (TensorProduct u (Tensor s v b)) - cftlp _ c = coerceFmapTensorProduct ([]::[u]) - (fmap c :: VSCCoercion (v⊗a) (v⊗b)) + coerceFmapTensorProduct _ c = coerceFmapTensorProduct ([]::[u]) + (fmap c :: VSCCoercion s (v⊗a) (v⊗b)) wellDefinedVector = wellDefinedTensor wellDefinedTensor = arr (getVSCCoercion rassocTensor) >>> wellDefinedTensor >>> arr (fmap (getVSCCoercion lassocTensor)) @@ -933,23 +1168,22 @@ fzipWith = case dualSpaceWitness :: DualSpaceWitness v of DualSpaceWitness -> \f -> arr asTensor *** arr asTensor >>> fzipWith f >>> arr fromTensor -instance (TensorSpace v, Scalar v ~ s) - => Functor (Tensor s v) VSCCoercion VSCCoercion where - fmap = crcFmap - where crcFmap :: ∀ s v a b . (TensorSpace v, Scalar v ~ s) - => VSCCoercion a b -> VSCCoercion (Tensor s v a) (Tensor s v b) - crcFmap f = case coerceFmapTensorProduct ([]::[v]) f of - VSCCoercion -> VSCCoercion +instance ∀ v s . (TensorSpace v, Scalar v ~ s) + => Functor (Tensor s v) (VSCCoercion s) (VSCCoercion s) where + fmap :: ∀ a b . ( TensorSpace a, Scalar a ~ s + , TensorSpace b, Scalar b ~ s ) + => VSCCoercion s a b -> VSCCoercion s (Tensor s v a) (Tensor s v b) + fmap f@VSCCoercion = case coerceFmapTensorProduct @v [] f of + Coercion -> VSCCoercion -instance (LinearSpace v, Scalar v ~ s) - => Functor (LinearMap s v) VSCCoercion VSCCoercion where - fmap = crcFmap dualSpaceWitness - where crcFmap :: ∀ s v a b . (LinearSpace v, Scalar v ~ s) - => DualSpaceWitness v -> VSCCoercion a b - -> VSCCoercion (LinearMap s v a) (LinearMap s v b) - crcFmap DualSpaceWitness f - = case coerceFmapTensorProduct ([]::[DualVector v]) f of - VSCCoercion -> VSCCoercion +instance ∀ v s . (LinearSpace v, Scalar v ~ s) + => Functor (LinearMap s v) (VSCCoercion s) (VSCCoercion s) where + fmap :: ∀ a b . ( TensorSpace a, Scalar a ~ s + , TensorSpace b, Scalar b ~ s ) + => VSCCoercion s a b -> VSCCoercion s (LinearMap s v a) (LinearMap s v b) + fmap f@VSCCoercion = case dualSpaceWitness @v of + DualSpaceWitness -> case coerceFmapTensorProduct @(DualVector v) [] f of + Coercion -> VSCCoercion instance Category (LinearFunction s) where type Object (LinearFunction s) v = (TensorSpace v, Scalar v ~ s) @@ -970,7 +1204,7 @@ terminal = const0 instance EnhancedCat (->) (LinearFunction s) where arr = getLinearFunction -instance EnhancedCat (LinearFunction s) VSCCoercion where +instance EnhancedCat (LinearFunction s) (VSCCoercion s) where arr VSCCoercion = LinearFunction coerce instance (LinearSpace w, Num' s, Scalar w ~ s) @@ -984,15 +1218,42 @@ sampleLinearFunctionFn = LinearFunction $ \f -> sampleLinearFunction -+$> f . applyLinear -fromLinearFn :: VSCCoercion (LinearFunction s (LinearFunction s u v) w) - (Tensor s (LinearFunction s v u) w) -fromLinearFn = VSCCoercion +fromLinearFn :: ∀ s u v w . (DimensionAware u, DimensionAware v, DimensionAware w) + => VSCCoercion s (LinearFunction s (LinearFunction s u v) w) + (Tensor s (LinearFunction s v u) w) +fromLinearFn + = Maybe.zipWithTimesCommu (staticDimensionSing @u) (staticDimensionSing @v) VSCCoercion -asLinearFn :: VSCCoercion (Tensor s (LinearFunction s u v) w) +asLinearFn :: ∀ s u v w . (DimensionAware u, DimensionAware v, DimensionAware w) + => VSCCoercion s (Tensor s (LinearFunction s u v) w) (LinearFunction s (LinearFunction s v u) w) -asLinearFn = VSCCoercion +asLinearFn + = Maybe.zipWithTimesCommu (staticDimensionSing @u) (staticDimensionSing @v) VSCCoercion +instance ∀ s u v . ( LinearSpace u, LinearSpace v + , DimensionAware u, DimensionAware v + , Scalar u ~ s, Scalar v ~ s) + => DimensionAware (LinearFunction s u v) where + type StaticDimension (LinearFunction s u v) + = Maybe.ZipWithTimes (StaticDimension u) (StaticDimension v) + dimensionalityWitness = case (dimensionalityWitness @u, dimensionalityWitness @v) of + (IsStaticDimensional, IsStaticDimensional) + -> withKnownNat (dimensionalitySing @u %* dimensionalitySing @v) + IsStaticDimensional + (IsFlexibleDimensional, _) -> IsFlexibleDimensional + (_, IsFlexibleDimensional) -> IsFlexibleDimensional +instance ∀ s n u m v nm . ( n`Dimensional`u, m`Dimensional`v + , LinearSpace u, LinearSpace v, Scalar u ~ s, Scalar v ~ s + , nm ~ (n*m) ) + => nm`Dimensional`(LinearFunction s u v) where + knownDimensionalitySing = dimensionalitySing @u %* dimensionalitySing @v + unsafeFromArrayWithOffset i ar + = applyLinear-+$>(unsafeFromArrayWithOffset i ar :: LinearMap s u v) + unsafeWriteArrayWithOffset ar i + = unsafeWriteArrayWithOffset ar i . (sampleLinearFunction-+$>) + + instance ∀ s u v . (LinearSpace u, LinearSpace v, Scalar u ~ s, Scalar v ~ s) => TensorSpace (LinearFunction s u v) where type TensorProduct (LinearFunction s u v) w = LinearFunction s (LinearFunction s v u) w @@ -1014,7 +1275,7 @@ #if !MIN_VERSION_manifolds_core(0,6,0) BoundarylessWitness #endif - zeroTensor = fromLinearFn $ const0 + zeroTensor = fromLinearFn -+$=> const0 toFlatTensor = case scalarSpaceWitness :: ScalarSpaceWitness u of ScalarSpaceWitness -> fmap (getVSCCoercion fromLinearFn) $ applyDualVector fromFlatTensor = case ( scalarSpaceWitness :: ScalarSpaceWitness u @@ -1025,8 +1286,8 @@ -+$> coCurryLinearMap $ sampleLinearFunction-+$> f . applyLinear in applyLinear $ fromTensor $ t - addTensors t s = fromLinearFn $ (asLinearFn$t)^+^(asLinearFn$s) - subtractTensors t s = fromLinearFn $ (asLinearFn$t)^-^(asLinearFn$s) + addTensors t s = fromLinearFn -+$=> (asLinearFn-+$=>t)^+^(asLinearFn-+$=>s) + subtractTensors t s = fromLinearFn -+$=> (asLinearFn-+$=>t)^-^(asLinearFn-+$=>s) scaleTensor = bilinearFunction $ \μ (Tensor f) -> Tensor $ μ *^ f negateTensor = LinearFunction $ \(Tensor f) -> Tensor $ negateV f tensorProduct = case scalarSpaceWitness :: ScalarSpaceWitness u of @@ -1038,29 +1299,61 @@ -> Tensor s (LinearFunction s u v) w -+> Tensor s w (LinearFunction s u v) tt ScalarSpaceWitness DualSpaceWitness - = LinearFunction $ arr asLinearFn >>> \f + = LinearFunction $ (asLinearFn-+$=>) >>> \f -> (fmapTensor-+$>applyLinear) -+$> fmap fromTensor . rassocTensor $ transposeTensor . fmap transposeTensor -+$> fmap asTensor . coCurryLinearMap $ sampleLinearFunctionFn -+$> f - fmapTensor = bilinearFunction $ \f -> arr asLinearFn - >>> \g -> fromLinearFn $ f . g + fmapTensor = bilinearFunction $ \f -> (asLinearFn-+$=>) + >>> \g -> fromLinearFn -+$=> f . g fzipTensorWith = case scalarSpaceWitness :: ScalarSpaceWitness u of ScalarSpaceWitness -> bilinearFunction $ \f (g,h) - -> fromLinearFn $ f . ((asLinearFn$g)&&&(asLinearFn$h)) - coerceFmapTensorProduct _ VSCCoercion = VSCCoercion + -> fromLinearFn -+$=> + f . ((asLinearFn-+$=>g)&&&(asLinearFn-+$=>h)) + tensorUnsafeFromArrayWithOffset :: ∀ nm w o α + . ( nm`Dimensional`LinearFunction s u v + , TensorSpace w, o`Dimensional`w, Scalar w ~ s + , GArr.Vector α s ) + => Int -> α s -> (LinearFunction s u v⊗w) + tensorUnsafeFromArrayWithOffset + = case ( dimensionalityWitness @u, staticDimensionSing @u + , dimensionalityWitness @v, staticDimensionSing @v ) of + ( IsStaticDimensional, SJust sn + ,IsStaticDimensional, SJust sm ) + -> withKnownNat (sm%*sn) ( + withKnownNat ((sm%*sn)%*dimensionalitySing @w) ( + \i -> (fromLinearFn @s @v @u @w -+$=>) + . (applyLinear-+$>) + . unsafeFromArrayWithOffset i )) + tensorUnsafeWriteArrayWithOffset :: ∀ nm w o α σ + . ( nm`Dimensional`LinearFunction s u v + , TensorSpace w, o`Dimensional`w, Scalar w ~ s + , GArr.Vector α s ) + => GArr.Mutable α σ s -> Int -> (LinearFunction s u v⊗w) -> ST σ () + tensorUnsafeWriteArrayWithOffset + = case ( dimensionalityWitness @u, staticDimensionSing @u + , dimensionalityWitness @v, staticDimensionSing @v ) of + ( IsStaticDimensional, SJust sn + ,IsStaticDimensional, SJust sm ) + -> withKnownNat (sm%*sn) ( + withKnownNat ((sm%*sn)%*dimensionalitySing @w) ( + \ar i -> unsafeWriteArrayWithOffset ar i + . (sampleLinearFunction-+$>) + . (asLinearFn @s @u @v @w -+$=>) + )) + coerceFmapTensorProduct _ VSCCoercion = Coercion wellDefinedVector = arr sampleLinearFunction >>> wellDefinedVector >>> fmap (arr applyLinear) - wellDefinedTensor = arr asLinearFn >>> (. applyLinear) + wellDefinedTensor = (asLinearFn-+$=>) >>> (. applyLinear) >>> getLinearFunction sampleLinearFunction >>> wellDefinedVector - >>> fmap (arr fromLinearFn <<< \m + >>> fmap ((fromLinearFn-+$=>) <<< \m -> sampleLinearFunction >>> getLinearFunction applyLinear m) -exposeLinearFn :: VSCCoercion (LinearMap s (LinearFunction s u v) w) - (LinearFunction s (LinearFunction s u v) w) +exposeLinearFn :: VSCCoercion s (LinearMap s (LinearFunction s u v) w) + (LinearFunction s (LinearFunction s u v) w) exposeLinearFn = VSCCoercion instance (LinearSpace u, LinearSpace v, Scalar u ~ s, Scalar v ~ s) @@ -1068,16 +1361,18 @@ type DualVector (LinearFunction s u v) = LinearFunction s v u dualSpaceWitness = case ( dualSpaceWitness :: DualSpaceWitness u , dualSpaceWitness :: DualSpaceWitness v ) of - (DualSpaceWitness, DualSpaceWitness) -> DualSpaceWitness + (DualSpaceWitness, DualSpaceWitness) + -> Maybe.zipWithTimesCommu (staticDimensionSing @u) (staticDimensionSing @v) + DualSpaceWitness linearId = symVSC exposeLinearFn $ id tensorId = uncurryLinearMap . symVSC exposeLinearFn $ LinearFunction $ \f -> sampleLinearFunction-+$>tensorProduct-+$>f coerceDoubleDual = VSCCoercion - sampleLinearFunction = LinearFunction . arr $ symVSC exposeLinearFn + sampleLinearFunction = LinearFunction . (-+$=>) $ symVSC exposeLinearFn applyDualVector = case scalarSpaceWitness :: ScalarSpaceWitness u of ScalarSpaceWitness -> bilinearFunction $ \f g -> trace . sampleLinearFunction -+$> f . g - applyLinear = bilinearFunction $ \f g -> (exposeLinearFn $ f) -+$> g + applyLinear = bilinearFunction $ \f g -> (exposeLinearFn -+$=> f) -+$> g applyTensorFunctional = atf scalarSpaceWitness dualSpaceWitness where atf :: ∀ w . (LinearSpace w, Scalar w ~ s) => ScalarSpaceWitness u -> DualSpaceWitness w @@ -1085,28 +1380,29 @@ (LinearMap s (LinearFunction s u v) (DualVector w)) (LinearFunction s (Tensor s (LinearFunction s u v) w) s) atf ScalarSpaceWitness DualSpaceWitness = bilinearFunction $ \f g - -> trace -+$> fromTensor $ transposeTensor + -> trace -+$> fromTensor -+$=> transposeTensor -+$> fmap ((exposeLinearFn $ f) . applyLinear) -+$> ( transposeTensor -+$> deferLinearMap - $ fmap transposeTensor + -+$=> fmap transposeTensor -+$> hasteLinearMap - $ transposeTensor + -+$=> transposeTensor -+$> coCurryLinearMap - $ sampleLinearFunctionFn - -+$> asLinearFn $ g ) + -+$=> sampleLinearFunctionFn + -+$> asLinearFn -+$=> g ) applyTensorLinMap = case scalarSpaceWitness :: ScalarSpaceWitness u of ScalarSpaceWitness -> bilinearFunction $ \f g -> contractMapTensor . transposeTensor - -+$> fmap ((asLinearFn $ g) . applyLinear) + -+$> fmap ((asLinearFn-+$=>g) . applyLinear) -+$> ( transposeTensor -+$> deferLinearMap - $ fmap transposeTensor + -+$=> fmap transposeTensor -+$> hasteLinearMap - $ transposeTensor + -+$=> transposeTensor -+$> coCurryLinearMap - $ sampleLinearFunctionFn - -+$> exposeLinearFn . curryLinearMap $ f ) + -+$=> sampleLinearFunctionFn + -+$> exposeLinearFn + -+$=> curryLinearMap -+$=> f ) useTupleLinearSpaceComponents _ = usingNonTupleTypeAsTupleError @@ -1141,6 +1437,18 @@ usingNonTupleTypeAsTupleError :: a usingNonTupleTypeAsTupleError = error "This is not a tuple type, the method should not be callable." +instance ∀ v s . DimensionAware v => DimensionAware (Gnrx.Rec0 v s) where + type StaticDimension (Gnrx.Rec0 v s) = StaticDimension v + dimensionalityWitness = case dimensionalityWitness @v of + IsStaticDimensional -> IsStaticDimensional + IsFlexibleDimensional -> IsFlexibleDimensional +instance ∀ n v s . n`Dimensional`v => n`Dimensional`(Gnrx.Rec0 v s) where + knownDimensionalitySing = dimensionalitySing @v + unsafeFromArrayWithOffset i ar + = coerce (unsafeFromArrayWithOffset @n @v i ar) + unsafeWriteArrayWithOffset i ar + = coerce (unsafeWriteArrayWithOffset @n @v i ar) + instance ∀ v s . TensorSpace v => TensorSpace (Gnrx.Rec0 v s) where type TensorProduct (Gnrx.Rec0 v s) w = TensorProduct v w wellDefinedVector = fmap Gnrx.K1 . wellDefinedVector . Gnrx.unK1 @@ -1150,9 +1458,9 @@ linearManifoldWitness = genericTensorspaceError zeroTensor = pseudoFmapTensorLHS Gnrx.K1 $ zeroTensor toFlatTensor = LinearFunction $ Gnrx.unK1 >>> getLinearFunction toFlatTensor - >>> arr (pseudoFmapTensorLHS Gnrx.K1) + >>> (pseudoFmapTensorLHS Gnrx.K1-+$=>) fromFlatTensor = LinearFunction $ Gnrx.K1 <<< getLinearFunction fromFlatTensor - <<< arr (pseudoFmapTensorLHS Gnrx.unK1) + <<< (pseudoFmapTensorLHS Gnrx.unK1-+$=>) addTensors (Tensor s) (Tensor t) = pseudoFmapTensorLHS Gnrx.K1 $ addTensors (Tensor s) (Tensor t) subtractTensors (Tensor s) (Tensor t) @@ -1167,9 +1475,10 @@ where tT :: ∀ w . (TensorSpace w, Scalar w ~ Scalar v) => (Gnrx.Rec0 v s ⊗ w) -+> (w ⊗ Gnrx.Rec0 v s) tT = LinearFunction - $ arr (VSCCoercion . coerceFmapTensorProduct ([]::[w]) - (VSCCoercion :: VSCCoercion v (Gnrx.Rec0 v s)) . VSCCoercion) - . getLinearFunction transposeTensor . arr (pseudoFmapTensorLHS Gnrx.unK1) + $ arr (Coercion . coerceFmapTensorProduct @w [] + (VSCCoercion :: VSCCoercion (Scalar v) v (Gnrx.Rec0 v s)) + . Coercion) + . getLinearFunction transposeTensor . (pseudoFmapTensorLHS Gnrx.unK1-+$=>) fmapTensor = LinearFunction $ \f -> envTensorLHSCoercion Gnrx.K1 (fmapTensor-+$>f) fzipTensorWith = bilinearFunction $ @@ -1177,14 +1486,43 @@ $ (fzipTensorWith-+$>f) -+$>( pseudoFmapTensorLHS Gnrx.unK1 $ wt , pseudoFmapTensorLHS Gnrx.unK1 $ xt ) + tensorUnsafeFromArrayWithOffset + :: ∀ w m a . ( TensorSpace w, m`Dimensional`w, Scalar w ~ Scalar v + , GArr.Vector a (Scalar v) ) + => Int -> a (Scalar v) -> (Gnrx.Rec0 v s⊗w) + tensorUnsafeFromArrayWithOffset = case dimensionalityWitness @v of + IsFlexibleDimensional -> error "This is impossible, since this can only be evaluated if `v` is static-dimensional." + IsStaticDimensional -> \i ar + -> coerce (tensorUnsafeFromArrayWithOffset @v @w i ar) + tensorUnsafeWriteArrayWithOffset + :: ∀ w m α σ . ( TensorSpace w, m`Dimensional`w, Scalar w ~ Scalar v + , GArr.Vector α (Scalar v) ) + => GArr.Mutable α σ (Scalar v) -> Int -> (Gnrx.Rec0 v s⊗w) -> ST σ () + tensorUnsafeWriteArrayWithOffset = case dimensionalityWitness @v of + IsFlexibleDimensional -> error "This is impossible, since this can only be evaluated if `v` is static-dimensional." + IsStaticDimensional -> \ar -> coerce (tensorUnsafeWriteArrayWithOffset @v @w ar) coerceFmapTensorProduct = cmtp - where cmtp :: ∀ p a b . Hask.Functor p - => p (Gnrx.Rec0 v s) -> VSCCoercion a b - -> VSCCoercion (TensorProduct (Gnrx.Rec0 v s) a) + where cmtp :: ∀ p a b . ( Hask.Functor p + , TensorSpace a, Scalar a ~ Scalar v + , TensorSpace b, Scalar b ~ Scalar v ) + => p (Gnrx.Rec0 v s) -> VSCCoercion (Scalar v) a b + -> Coercion (TensorProduct (Gnrx.Rec0 v s) a) (TensorProduct (Gnrx.Rec0 v s) b) - cmtp p crc = case coerceFmapTensorProduct ([]::[v]) crc of - VSCCoercion -> VSCCoercion + cmtp p crc = case coerceFmapTensorProduct @v [] crc of + Coercion -> Coercion +instance ∀ i c f p . DimensionAware (f p) => DimensionAware (Gnrx.M1 i c f p) where + type StaticDimension (Gnrx.M1 i c f p) = StaticDimension (f p) + dimensionalityWitness = case dimensionalityWitness @(f p) of + IsStaticDimensional -> IsStaticDimensional + IsFlexibleDimensional -> IsFlexibleDimensional +instance ∀ n i c f p . n`Dimensional`f p => n`Dimensional`Gnrx.M1 i c f p where + knownDimensionalitySing = dimensionalitySing @(f p) + unsafeFromArrayWithOffset i ar + = coerce (unsafeFromArrayWithOffset @n @(f p) i ar) + unsafeWriteArrayWithOffset i ar + = coerce (unsafeWriteArrayWithOffset @n @(f p) i ar) + instance ∀ i c f p . TensorSpace (f p) => TensorSpace (Gnrx.M1 i c f p) where type TensorProduct (Gnrx.M1 i c f p) w = TensorProduct (f p) w wellDefinedVector = fmap Gnrx.M1 . wellDefinedVector . Gnrx.unM1 @@ -1194,9 +1532,9 @@ linearManifoldWitness = genericTensorspaceError zeroTensor = pseudoFmapTensorLHS Gnrx.M1 $ zeroTensor toFlatTensor = LinearFunction $ Gnrx.unM1 >>> getLinearFunction toFlatTensor - >>> arr (pseudoFmapTensorLHS Gnrx.M1) + >>> (pseudoFmapTensorLHS Gnrx.M1-+$=>) fromFlatTensor = LinearFunction $ Gnrx.M1 <<< getLinearFunction fromFlatTensor - <<< arr (pseudoFmapTensorLHS Gnrx.unM1) + <<< (pseudoFmapTensorLHS Gnrx.unM1-+$=>) addTensors (Tensor s) (Tensor t) = pseudoFmapTensorLHS Gnrx.M1 $ addTensors (Tensor s) (Tensor t) subtractTensors (Tensor s) (Tensor t) @@ -1211,9 +1549,10 @@ where tT :: ∀ w . (TensorSpace w, Scalar w ~ Scalar (f p)) => (Gnrx.M1 i c f p ⊗ w) -+> (w ⊗ Gnrx.M1 i c f p) tT = LinearFunction - $ arr (VSCCoercion . coerceFmapTensorProduct ([]::[w]) - (VSCCoercion :: VSCCoercion (f p) (Gnrx.M1 i c f p)) . VSCCoercion) - . getLinearFunction transposeTensor . arr (pseudoFmapTensorLHS Gnrx.unM1) + $ arr (Coercion . coerceFmapTensorProduct ([]::[w]) + (VSCCoercion :: VSCCoercion s (f p) (Gnrx.M1 i c f p)) + . Coercion) + . getLinearFunction transposeTensor . (pseudoFmapTensorLHS Gnrx.unM1-+$=>) fmapTensor = LinearFunction $ \f -> envTensorLHSCoercion Gnrx.M1 (fmapTensor-+$>f) fzipTensorWith = bilinearFunction $ @@ -1221,14 +1560,55 @@ $ (fzipTensorWith-+$>f) -+$>( pseudoFmapTensorLHS Gnrx.unM1 $ wt , pseudoFmapTensorLHS Gnrx.unM1 $ xt ) - coerceFmapTensorProduct = cmtp - where cmtp :: ∀ ぴ a b . Hask.Functor ぴ - => ぴ (Gnrx.M1 i c f p) -> VSCCoercion a b - -> VSCCoercion (TensorProduct (Gnrx.M1 i c f p) a) + tensorUnsafeFromArrayWithOffset + :: ∀ w m a . ( TensorSpace w, m`Dimensional`w, Scalar w ~ Scalar (f p) + , GArr.Vector a (Scalar (f p)) ) + => Int -> a (Scalar (f p)) -> (Gnrx.M1 i c f p⊗w) + tensorUnsafeFromArrayWithOffset = case dimensionalityWitness @(f p) of + IsFlexibleDimensional -> error "This is impossible, since this can only be evaluated if `f p` is static-dimensional." + IsStaticDimensional -> \i ar + -> coerce (tensorUnsafeFromArrayWithOffset @(f p) @w i ar) + tensorUnsafeWriteArrayWithOffset + :: ∀ w m α σ . ( TensorSpace w, m`Dimensional`w, Scalar w ~ Scalar (f p) + , GArr.Vector α (Scalar (f p)) ) + => GArr.Mutable α σ (Scalar (f p)) -> Int -> (Gnrx.M1 i c f p⊗w) -> ST σ () + tensorUnsafeWriteArrayWithOffset = case dimensionalityWitness @(f p) of + IsFlexibleDimensional -> error "This is impossible, since this can only be evaluated if `f p` is static-dimensional." + IsStaticDimensional -> \ar -> + coerce (tensorUnsafeWriteArrayWithOffset @(f p) @w ar) + coerceFmapTensorProduct :: ∀ ぴ a b + . (Hask.Functor ぴ, TensorSpace a, Scalar a ~ Scalar (f p) + , TensorSpace b, Scalar b ~ Scalar (f p) ) + => ぴ (Gnrx.M1 i c f p) -> VSCCoercion (Scalar (f p)) a b + -> Coercion (TensorProduct (Gnrx.M1 i c f p) a) (TensorProduct (Gnrx.M1 i c f p) b) - cmtp p crc = case coerceFmapTensorProduct ([]::[f p]) crc of - VSCCoercion -> VSCCoercion + coerceFmapTensorProduct p crc = case coerceFmapTensorProduct ([]::[f p]) crc of + Coercion -> Coercion +instance ∀ f g p . ( DimensionAware (f p), DimensionAware (g p) + , Scalar (f p) ~ Scalar (g p) ) + => DimensionAware ((f:*:g) p) where + type StaticDimension ((f:*:g) p) + = Maybe.ZipWithPlus (StaticDimension (f p)) (StaticDimension (g p)) + dimensionalityWitness = case ( dimensionalityWitness @(f p) + , dimensionalityWitness @(g p) ) of + (IsStaticDimensional, IsStaticDimensional) + -> withKnownNat (dimensionalitySing @(f p) %+ dimensionalitySing @(g p)) + IsStaticDimensional + (IsFlexibleDimensional, _) -> IsFlexibleDimensional + (_, IsFlexibleDimensional) -> IsFlexibleDimensional +instance ∀ n f m g p nm . ( n`Dimensional`(f p), m`Dimensional`(g p) + , Scalar (f p) ~ Scalar (g p) + , nm ~ (n+m) ) + => nm`Dimensional`((f:*:g) p) where + knownDimensionalitySing = dimensionalitySing @(f p) %+ dimensionalitySing @(g p) + unsafeFromArrayWithOffset i ar + = unsafeFromArrayWithOffset i ar + :*: unsafeFromArrayWithOffset (i + dimension @(f p)) ar + unsafeWriteArrayWithOffset ar i (x:*:y) = do + unsafeWriteArrayWithOffset ar i x + unsafeWriteArrayWithOffset ar (i + dimension @(f p)) y + instance ∀ f g p . ( TensorSpace (f p), TensorSpace (g p), Scalar (f p) ~ Scalar (g p) ) => TensorSpace ((f:*:g) p) where type TensorProduct ((f:*:g) p) w = (f p⊗w, g p⊗w) @@ -1260,22 +1640,63 @@ $ \f (Tensor (uw, vw), Tensor (ux, vx)) -> Tensor ( (fzipTensorWith-+$>f)-+$>(uw,ux) , (fzipTensorWith-+$>f)-+$>(vw,vx) ) + tensorUnsafeFromArrayWithOffset + :: ∀ w m α . ( TensorSpace w, m`Dimensional`w, Scalar w ~ Scalar (f p) + , GArr.Vector α (Scalar (f p)) ) + => Int -> α (Scalar (f p)) -> ((f:*:g) p⊗w) + tensorUnsafeFromArrayWithOffset + = case (dimensionalityWitness @(f p), dimensionalityWitness @(g p)) of + (IsFlexibleDimensional, _) -> error "This is impossible, since this can only be evaluated if `f p` is static-dimensional." + (_, IsFlexibleDimensional) -> error "This is impossible, since this can only be evaluated if `g p` is static-dimensional." + (IsStaticDimensional, IsStaticDimensional) + -> withKnownNat (dimensionalitySing @(f p) %+ dimensionalitySing @(g p)) + (\i ar + -> coerce (tensorUnsafeFromArrayWithOffset @(f p, g p) @w i ar) ) + tensorUnsafeWriteArrayWithOffset + :: ∀ w m α σ . ( TensorSpace w, m`Dimensional`w, Scalar w ~ Scalar (f p) + , GArr.Vector α (Scalar (f p)) ) + => GArr.Mutable α σ (Scalar (f p)) -> Int -> ((f:*:g) p⊗w) -> ST σ () + tensorUnsafeWriteArrayWithOffset + = case (dimensionalityWitness @(f p), dimensionalityWitness @(g p)) of + (IsFlexibleDimensional, _) -> error "This is impossible, since this can only be evaluated if `f p` is static-dimensional." + (_, IsFlexibleDimensional) -> error "This is impossible, since this can only be evaluated if `g p` is static-dimensional." + (IsStaticDimensional, IsStaticDimensional) + -> withKnownNat (dimensionalitySing @(f p) %+ dimensionalitySing @(g p)) + (\ar + -> coerce (tensorUnsafeWriteArrayWithOffset @(f p, g p) @w ar) ) coerceFmapTensorProduct p cab = case ( coerceFmapTensorProduct ((\(u:*:_)->u)<$>p) cab , coerceFmapTensorProduct ((\(_:*:v)->v)<$>p) cab ) of - (VSCCoercion, VSCCoercion) -> VSCCoercion + (Coercion, Coercion) -> Coercion wellDefinedVector (u:*:v) = liftA2 (:*:) (wellDefinedVector u) (wellDefinedVector v) wellDefinedTensor (Tensor (u,v)) = liftA2 ((Tensor.) . (,)) (wellDefinedTensor u) (wellDefinedTensor v) +instance ∀ m . ( Semimanifold m, DimensionAware (Needle (VRep m)) + , Scalar (Needle m) ~ Scalar (Needle (VRep m)) ) + => DimensionAware (GenericNeedle m) where + type StaticDimension (GenericNeedle m) = StaticDimension (Needle (VRep m)) + dimensionalityWitness = case dimensionalityWitness @(Needle (VRep m)) of + IsStaticDimensional -> IsStaticDimensional + IsFlexibleDimensional -> IsFlexibleDimensional +instance ∀ n m . ( Semimanifold m, n`Dimensional`Needle (VRep m) + , Scalar (Needle m) ~ Scalar (Needle (VRep m)) ) + => n`Dimensional`GenericNeedle m where + knownDimensionalitySing = dimensionalitySing @(Needle (VRep m)) + unsafeFromArrayWithOffset i ar + = coerce (unsafeFromArrayWithOffset @n @(Needle (VRep m)) i ar) + unsafeWriteArrayWithOffset ar i + = coerce (unsafeWriteArrayWithOffset @n @(Needle (VRep m)) ar i) + instance ∀ m . ( Semimanifold m, TensorSpace (Needle (VRep m)) , Scalar (Needle m) ~ Scalar (Needle (VRep m)) ) => TensorSpace (GenericNeedle m) where type TensorProduct (GenericNeedle m) w = TensorProduct (Needle (VRep m)) w wellDefinedVector = fmap GenericNeedle . wellDefinedVector . getGenericNeedle wellDefinedTensor = arr (fmap . getVSCCoercion $ pseudoFmapTensorLHS GenericNeedle) - . wellDefinedTensor . arr (pseudoFmapTensorLHS getGenericNeedle) + . wellDefinedTensor + . (pseudoFmapTensorLHS getGenericNeedle-+$=>) scalarSpaceWitness = case scalarSpaceWitness :: ScalarSpaceWitness (Needle (VRep m)) of ScalarSpaceWitness -> ScalarSpaceWitness @@ -1290,10 +1711,10 @@ BoundarylessWitness #endif zeroTensor = pseudoFmapTensorLHS GenericNeedle $ zeroTensor - toFlatTensor = LinearFunction $ arr (pseudoFmapTensorLHS GenericNeedle) + toFlatTensor = LinearFunction $ (pseudoFmapTensorLHS GenericNeedle-+$=>) . getLinearFunction toFlatTensor . getGenericNeedle - fromFlatTensor = LinearFunction $ arr (pseudoFmapTensorLHS getGenericNeedle) + fromFlatTensor = LinearFunction $ (pseudoFmapTensorLHS getGenericNeedle-+$=>) >>> getLinearFunction fromFlatTensor >>> GenericNeedle addTensors (Tensor s) (Tensor t) @@ -1310,10 +1731,13 @@ where tT :: ∀ w . (TensorSpace w, Scalar w ~ Scalar (Needle m)) => (GenericNeedle m ⊗ w) -+> (w ⊗ GenericNeedle m) tT = LinearFunction - $ arr (VSCCoercion . coerceFmapTensorProduct ([]::[w]) - (VSCCoercion :: VSCCoercion (Needle (VRep m)) - (GenericNeedle m)) . VSCCoercion) - . getLinearFunction transposeTensor . arr (pseudoFmapTensorLHS getGenericNeedle) + $ arr (Coercion . coerceFmapTensorProduct ([]::[w]) + (VSCCoercion :: VSCCoercion (Scalar (Needle m)) + (Needle (VRep m)) + (GenericNeedle m)) + . Coercion) + . getLinearFunction transposeTensor + . (pseudoFmapTensorLHS getGenericNeedle-+$=>) fmapTensor = LinearFunction $ \f -> envTensorLHSCoercion GenericNeedle (fmapTensor-+$>f) fzipTensorWith = bilinearFunction $ @@ -1321,15 +1745,39 @@ $ (fzipTensorWith-+$>f) -+$>( pseudoFmapTensorLHS getGenericNeedle $ wt , pseudoFmapTensorLHS getGenericNeedle $ xt ) + tensorUnsafeFromArrayWithOffset + :: ∀ w nn α . ( TensorSpace w, nn`Dimensional`w + , Scalar w ~ (Scalar (Needle (VRep m))) + , GArr.Vector α (Scalar (Needle (VRep m))) ) + => Int -> α (Scalar (Needle (VRep m))) + -> (GenericNeedle m⊗w) + tensorUnsafeFromArrayWithOffset + = case dimensionalityWitness @(Needle (VRep m)) of + IsFlexibleDimensional -> error "This is impossible, since this can only be evaluated if `Needle (VRep m)` is static-dimensional." + IsStaticDimensional -> \i ar + -> coerce (tensorUnsafeFromArrayWithOffset @(Needle (VRep m)) @w i ar) + tensorUnsafeWriteArrayWithOffset + :: ∀ w nn α σ . ( TensorSpace w, nn`Dimensional`w + , Scalar w ~ (Scalar (Needle (VRep m))) + , GArr.Vector α (Scalar (Needle (VRep m))) ) + => GArr.Mutable α σ (Scalar (Needle (VRep m))) + -> Int -> (GenericNeedle m⊗w) -> ST σ () + tensorUnsafeWriteArrayWithOffset + = case dimensionalityWitness @(Needle (VRep m)) of + IsFlexibleDimensional -> error "This is impossible, since this can only be evaluated if `Needle (VRep m)` is static-dimensional." + IsStaticDimensional -> \ar + -> coerce (tensorUnsafeWriteArrayWithOffset @(Needle (VRep m)) @w ar) coerceFmapTensorProduct = cmtp - where cmtp :: ∀ p a b . Hask.Functor p - => p (GenericNeedle m) -> VSCCoercion a b - -> VSCCoercion (TensorProduct (GenericNeedle m) a) + where cmtp :: ∀ p a b . ( Hask.Functor p + , TensorSpace a, Scalar a ~ Scalar (Needle (VRep m)) + , TensorSpace b, Scalar b ~ Scalar (Needle (VRep m)) ) + => p (GenericNeedle m) -> VSCCoercion (Scalar a) a b + -> Coercion (TensorProduct (GenericNeedle m) a) (TensorProduct (GenericNeedle m) b) - cmtp p crc = case coerceFmapTensorProduct ([]::[Needle (VRep m)]) crc of - VSCCoercion -> VSCCoercion + cmtp p crc = case coerceFmapTensorProduct @(Needle (VRep m)) [] crc of + Coercion -> Coercion -instance (LinearSpace v, Num (Scalar v)) => LinearSpace (Gnrx.Rec0 v s) where +instance ∀ v s . (LinearSpace v, Num (Scalar v)) => LinearSpace (Gnrx.Rec0 v s) where type DualVector (Gnrx.Rec0 v s) = DualVector v dualSpaceWitness = genericTensorspaceError linearId = pseudoPrecomposeLinmap Gnrx.unK1 @@ -1345,6 +1793,8 @@ applyTensorLinMap = bilinearFunction $ \(LinearMap f) t -> (applyTensorLinMap-+$>LinearMap f)-+$>pseudoFmapTensorLHS Gnrx.unK1 $ t useTupleLinearSpaceComponents _ = usingNonTupleTypeAsTupleError + coerceDoubleDual = case coerceDoubleDual @v of + VSCCoercion -> VSCCoercion instance (LinearSpace (f p), Num (Scalar (f p))) => LinearSpace (Gnrx.M1 i c f p) where type DualVector (Gnrx.M1 i c f p) = DualVector (f p) @@ -1362,6 +1812,8 @@ applyTensorLinMap = bilinearFunction $ \(LinearMap f) t -> (applyTensorLinMap-+$>LinearMap f)-+$>pseudoFmapTensorLHS Gnrx.unM1 $ t useTupleLinearSpaceComponents _ = usingNonTupleTypeAsTupleError + coerceDoubleDual = case coerceDoubleDual @(f p) of + VSCCoercion -> VSCCoercion data GenericTupleDual f g p = GenericTupleDual !(DualVector (f p)) !(DualVector (g p)) deriving (Generic) @@ -1393,6 +1845,36 @@ p.-~.q = Just $ p.-.q (.-~!) = (.-.) + +instance ( DimensionAware (f p), DimensionAware (g p) + , VectorSpace (DualVector (f p)), VectorSpace (DualVector (g p)) + , Scalar (f p) ~ Scalar (g p) + , Scalar (f p) ~ Scalar (DualVector (f p)) + , Scalar (g p) ~ Scalar (DualVector (g p)) ) + => DimensionAware (GenericTupleDual f g p) where + type StaticDimension (GenericTupleDual f g p) + = Maybe.ZipWithPlus (StaticDimension (f p)) (StaticDimension (g p)) + dimensionalityWitness = case ( dimensionalityWitness @(f p) + , dimensionalityWitness @(g p) ) of + (IsStaticDimensional, IsStaticDimensional) + -> withKnownNat (dimensionalitySing @(f p) %+ dimensionalitySing @(g p)) + IsStaticDimensional + (IsFlexibleDimensional, _) -> IsFlexibleDimensional + (_, IsFlexibleDimensional) -> IsFlexibleDimensional +instance ∀ n f m g p nm . + ( n`Dimensional`f p, m`Dimensional`g p + , VectorSpace (DualVector (f p)), VectorSpace (DualVector (g p)) + , Scalar (f p) ~ Scalar (g p) + , Scalar (f p) ~ Scalar (DualVector (f p)) + , Scalar (g p) ~ Scalar (DualVector (g p)) + , nm ~ (n+m) ) + => nm`Dimensional`GenericTupleDual f g p where + knownDimensionalitySing = dimensionalitySing @(f p) %+ dimensionalitySing @(g p) + unsafeFromArrayWithOffset i ar + = coerce (unsafeFromArrayWithOffset @nm @(GenericTupleDual f g p) i ar) + unsafeWriteArrayWithOffset i ar + = coerce (unsafeWriteArrayWithOffset @nm @(GenericTupleDual f g p) i ar) + instance ( LinearSpace (f p), LinearSpace (g p) , VectorSpace (DualVector (f p)), VectorSpace (DualVector (g p)) , Scalar (f p) ~ Scalar (DualVector (f p)) @@ -1474,12 +1956,49 @@ , asTensor $ fx ) , fromTensor $ (fzipTensorWith-+$>f) -+$> ( asTensor $ gw , asTensor $ gx ) ) + tensorUnsafeFromArrayWithOffset + :: ∀ w m α . ( TensorSpace w, m`Dimensional`w, Scalar w ~ Scalar (f p) + , GArr.Vector α (Scalar (f p)) ) + => Int -> α (Scalar (f p)) -> (GenericTupleDual f g p⊗w) + tensorUnsafeFromArrayWithOffset + = case ( dualSpaceWitness @(f p), dualSpaceWitness @(g p) ) of + (DualSpaceWitness, DualSpaceWitness) -> case + ( dimensionalityWitness @(DualVector (f p)) + , dimensionalityWitness @(DualVector (g p)) ) of + (IsFlexibleDimensional, _) + -> error "This is impossible, since this can only be evaluated if `f p` is static-dimensional." + (_, IsFlexibleDimensional) -> error "This is impossible, since this can only be evaluated if `g p` is static-dimensional." + (IsStaticDimensional, IsStaticDimensional) + -> withKnownNat (dimensionalitySing @(DualVector (f p)) + %+ dimensionalitySing @(DualVector (g p))) + (\i ar + -> coerce (tensorUnsafeFromArrayWithOffset + @(DualVector (f p), DualVector (g p)) @w i ar) ) + tensorUnsafeWriteArrayWithOffset + :: ∀ w m α σ . ( TensorSpace w, m`Dimensional`w, Scalar w ~ Scalar (f p) + , GArr.Vector α (Scalar (f p)) ) + => GArr.Mutable α σ (Scalar (f p)) -> Int -> (GenericTupleDual f g p⊗w) + -> ST σ () + tensorUnsafeWriteArrayWithOffset + = case ( dualSpaceWitness @(f p), dualSpaceWitness @(g p) ) of + (DualSpaceWitness, DualSpaceWitness) -> case + ( dimensionalityWitness @(DualVector (f p)) + , dimensionalityWitness @(DualVector (g p)) ) of + (IsFlexibleDimensional, _) + -> error "This is impossible, since this can only be evaluated if `f p` is static-dimensional." + (_, IsFlexibleDimensional) -> error "This is impossible, since this can only be evaluated if `g p` is static-dimensional." + (IsStaticDimensional, IsStaticDimensional) + -> withKnownNat (dimensionalitySing @(DualVector (f p)) + %+ dimensionalitySing @(DualVector (g p))) + (\ar + -> coerce (tensorUnsafeWriteArrayWithOffset + @(DualVector (f p), DualVector (g p)) @w ar) ) coerceFmapTensorProduct p cab = case ( dualSpaceWitness :: DualSpaceWitness (f p) , dualSpaceWitness :: DualSpaceWitness (g p) ) of (DualSpaceWitness, DualSpaceWitness) -> case ( coerceFmapTensorProduct ((\(GenericTupleDual u _)->u)<$>p) cab , coerceFmapTensorProduct ((\(GenericTupleDual _ v)->v)<$>p) cab ) of - (VSCCoercion, VSCCoercion) -> VSCCoercion + (Coercion, Coercion) -> Coercion @@ -1537,14 +2056,33 @@ (DualSpaceWitness, DualSpaceWitness) -> bilinearFunction $ \(LinearMap (fu,fv)) (Tensor (tu,tv)) -> ((applyTensorFunctional-+$>fu)-+$>tu) + ((applyTensorFunctional-+$>fu)-+$>tu) - applyTensorLinMap = case ( dualSpaceWitness :: DualSpaceWitness (f p) - , dualSpaceWitness :: DualSpaceWitness (g p) ) of - (DualSpaceWitness, DualSpaceWitness) -> bilinearFunction`id` + applyTensorLinMap :: ∀ u w . ( LinearSpace u, TensorSpace w + , Scalar u ~ Scalar (g p), Scalar w ~ Scalar (g p) ) + => LinearFunction (Scalar (g p)) + (LinearMap (Scalar (g p)) (Tensor (Scalar (g p)) ((:*:) f g p) u) w) + (LinearFunction (Scalar (g p)) (Tensor (Scalar (g p)) ((:*:) f g p) u) w) + applyTensorLinMap = case ( dualSpaceWitness @(f p) + , dualSpaceWitness @(g p) + , dualSpaceWitness @u ) of + (DualSpaceWitness, DualSpaceWitness, DualSpaceWitness) -> bilinearFunction`id` \(LinearMap (fu,fv)) (Tensor (tu,tv)) - -> ((applyTensorLinMap -+$> uncurryLinearMap . fmap fromTensor $ fu)-+$>tu) - ^+^ ((applyTensorLinMap -+$> uncurryLinearMap . fmap fromTensor $ fv)-+$>tv) + -> ((applyTensorLinMap -+$> uncurryLinearMap -+$=> fmap fromTensor -+$=> fu)-+$>tu) + ^+^ ((applyTensorLinMap -+$> uncurryLinearMap -+$=> fmap fromTensor -+$=> fv)-+$>tv) useTupleLinearSpaceComponents _ = usingNonTupleTypeAsTupleError + coerceDoubleDual = case ( coerceDoubleDual @(f p), dualSpaceWitness @(f p) + , coerceDoubleDual @(g p), dualSpaceWitness @(g p)) of + (VSCCoercion, DualSpaceWitness, VSCCoercion, DualSpaceWitness) -> VSCCoercion +instance ( LinearSpace (f p), LinearSpace (g p) + , VectorSpace (DualVector (f p)), VectorSpace (DualVector (g p)) + , Scalar (f p) ~ Scalar (DualVector (f p)) + , Scalar (g p) ~ Scalar (DualVector (g p)) + , Scalar (DualVector (f p)) ~ Scalar (DualVector (g p)) ) + => LinearSpace (GenericTupleDual f g p) where + type DualVector (GenericTupleDual f g p) = (f:*:g) p + coerceDoubleDual = case ( coerceDoubleDual @(f p), dualSpaceWitness @(f p) + , coerceDoubleDual @(g p), dualSpaceWitness @(g p)) of + (VSCCoercion, DualSpaceWitness, VSCCoercion, DualSpaceWitness) -> VSCCoercion newtype GenericNeedle' m = GenericNeedle' { getGenericNeedle' :: DualVector (Needle (VRep m)) } @@ -1574,6 +2112,26 @@ => PseudoAffine (GenericNeedle' m) where p.-~.q = pure (p^-^q) (.-~!) = (^-^) + + +instance ∀ m . ( Semimanifold m, DimensionAware (DualVector (Needle (VRep m))) + , Scalar (Needle m) ~ Scalar (DualVector (Needle (VRep m))) ) + => DimensionAware (GenericNeedle' m) where + type StaticDimension (GenericNeedle' m) + = StaticDimension (DualVector (Needle (VRep m))) + dimensionalityWitness = case dimensionalityWitness + @(DualVector (Needle (VRep m))) of + IsStaticDimensional -> IsStaticDimensional + IsFlexibleDimensional -> IsFlexibleDimensional +instance ∀ n m . ( Semimanifold m, n`Dimensional`DualVector (Needle (VRep m)) + , Scalar (Needle m) ~ Scalar (DualVector (Needle (VRep m))) ) + => n`Dimensional`GenericNeedle' m where + knownDimensionalitySing = dimensionalitySing @(DualVector (Needle (VRep m))) + unsafeFromArrayWithOffset i ar + = coerce (unsafeFromArrayWithOffset @n @(DualVector (Needle (VRep m))) i ar) + unsafeWriteArrayWithOffset ar + = coerce (unsafeWriteArrayWithOffset @n @(DualVector (Needle (VRep m))) ar) + instance ∀ m . ( Semimanifold m, TensorSpace (DualVector (Needle (VRep m))) , Scalar (Needle m) ~ Scalar (DualVector (Needle (VRep m))) ) => TensorSpace (GenericNeedle' m) where @@ -1581,7 +2139,7 @@ = TensorProduct (DualVector (Needle (VRep m))) w wellDefinedVector = fmap GenericNeedle' . wellDefinedVector . getGenericNeedle' wellDefinedTensor = arr (fmap . getVSCCoercion $ pseudoFmapTensorLHS GenericNeedle') - . wellDefinedTensor . arr (pseudoFmapTensorLHS getGenericNeedle') + . wellDefinedTensor . (pseudoFmapTensorLHS getGenericNeedle'-+$=>) scalarSpaceWitness = case scalarSpaceWitness :: ScalarSpaceWitness (DualVector (Needle (VRep m))) of ScalarSpaceWitness -> ScalarSpaceWitness @@ -1596,10 +2154,10 @@ BoundarylessWitness #endif zeroTensor = pseudoFmapTensorLHS GenericNeedle' $ zeroTensor - toFlatTensor = LinearFunction $ arr (pseudoFmapTensorLHS GenericNeedle') + toFlatTensor = LinearFunction $ (pseudoFmapTensorLHS GenericNeedle'-+$=>) . getLinearFunction toFlatTensor . getGenericNeedle' - fromFlatTensor = LinearFunction $ arr (pseudoFmapTensorLHS getGenericNeedle') + fromFlatTensor = LinearFunction $ (pseudoFmapTensorLHS getGenericNeedle'-+$=>) >>> getLinearFunction fromFlatTensor >>> GenericNeedle' addTensors (Tensor s) (Tensor t) @@ -1616,10 +2174,14 @@ where tT :: ∀ w . (TensorSpace w, Scalar w ~ Scalar (Needle m)) => (GenericNeedle' m ⊗ w) -+> (w ⊗ GenericNeedle' m) tT = LinearFunction - $ arr (VSCCoercion . coerceFmapTensorProduct ([]::[w]) - (VSCCoercion :: VSCCoercion (DualVector (Needle (VRep m))) - (GenericNeedle' m)) . VSCCoercion) - . getLinearFunction transposeTensor . arr (pseudoFmapTensorLHS getGenericNeedle') + $ arr (Coercion . coerceFmapTensorProduct ([]::[w]) + (VSCCoercion :: VSCCoercion + (Scalar (Needle m)) + (DualVector (Needle (VRep m))) + (GenericNeedle' m)) + . Coercion) + . getLinearFunction transposeTensor + . (pseudoFmapTensorLHS getGenericNeedle'-+$=>) fmapTensor = LinearFunction $ \f -> envTensorLHSCoercion GenericNeedle' (fmapTensor-+$>f) fzipTensorWith = bilinearFunction $ @@ -1627,14 +2189,40 @@ $ (fzipTensorWith-+$>f) -+$>( pseudoFmapTensorLHS getGenericNeedle' $ wt , pseudoFmapTensorLHS getGenericNeedle' $ xt ) - coerceFmapTensorProduct = cmtp - where cmtp :: ∀ p a b . Hask.Functor p - => p (GenericNeedle' m) -> VSCCoercion a b - -> VSCCoercion (TensorProduct (GenericNeedle' m) a) + tensorUnsafeFromArrayWithOffset + :: ∀ w nn α . ( TensorSpace w, nn`Dimensional`w + , Scalar w ~ (Scalar (DualVector (Needle (VRep m)))) + , GArr.Vector α (Scalar (DualVector (Needle (VRep m)))) ) + => Int -> α (Scalar (DualVector (Needle (VRep m)))) + -> (GenericNeedle' m⊗w) + tensorUnsafeFromArrayWithOffset + = case dimensionalityWitness @(DualVector (Needle (VRep m))) of + IsFlexibleDimensional -> error "This is impossible, since this can only be evaluated if `Needle (VRep m)` is static-dimensional." + IsStaticDimensional -> \i ar + -> coerce (tensorUnsafeFromArrayWithOffset + @(DualVector (Needle (VRep m))) @w i ar) + tensorUnsafeWriteArrayWithOffset + :: ∀ w nn α σ . ( TensorSpace w, nn`Dimensional`w + , Scalar w ~ (Scalar (DualVector (Needle (VRep m)))) + , GArr.Vector α (Scalar (DualVector (Needle (VRep m)))) ) + => GArr.Mutable α σ (Scalar (DualVector (Needle (VRep m)))) + -> Int -> (GenericNeedle' m⊗w) -> ST σ () + tensorUnsafeWriteArrayWithOffset + = case dimensionalityWitness @(DualVector (Needle (VRep m))) of + IsFlexibleDimensional -> error "This is impossible, since this can only be evaluated if `Needle (VRep m)` is static-dimensional." + IsStaticDimensional -> \ar + -> coerce (tensorUnsafeWriteArrayWithOffset + @(DualVector (Needle (VRep m))) @w ar) + coerceFmapTensorProduct :: ∀ p a b + . ( Hask.Functor p + , TensorSpace a, Scalar a ~ Scalar (DualVector (Needle (VRep m))) + , TensorSpace b, Scalar b ~ Scalar (DualVector (Needle (VRep m))) ) + => p (GenericNeedle' m) -> VSCCoercion (Scalar a) a b + -> Coercion (TensorProduct (GenericNeedle' m) a) (TensorProduct (GenericNeedle' m) b) - cmtp p crc = case coerceFmapTensorProduct + coerceFmapTensorProduct p crc = case coerceFmapTensorProduct ([]::[DualVector (Needle (VRep m))]) crc of - VSCCoercion -> VSCCoercion + Coercion -> Coercion instance ∀ s m . ( Num' s @@ -1662,6 +2250,8 @@ -> (applyTensorLinMap-+$>LinearMap f) -+$>pseudoFmapTensorLHS getGenericNeedle $ t useTupleLinearSpaceComponents _ = usingNonTupleTypeAsTupleError + coerceDoubleDual = case coerceDoubleDual @(Needle (VRep m)) of + VSCCoercion -> VSCCoercion instance ∀ s m . ( Num' s , Semimanifold m @@ -1696,3 +2286,5 @@ -> (applyTensorLinMap-+$>LinearMap f) -+$>pseudoFmapTensorLHS getGenericNeedle' $ t useTupleLinearSpaceComponents _ = usingNonTupleTypeAsTupleError + coerceDoubleDual = case coerceDoubleDual @(Needle (VRep m)) of + VSCCoercion -> VSCCoercion diff --git a/Math/LinearMap/Category/Instances.hs b/Math/LinearMap/Category/Instances.hs --- a/Math/LinearMap/Category/Instances.hs +++ b/Math/LinearMap/Category/Instances.hs @@ -13,7 +13,10 @@ {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE TypeOperators #-} {-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE DataKinds #-} {-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE InstanceSigs #-} {-# LANGUAGE StandaloneDeriving #-} {-# LANGUAGE UnicodeSyntax #-} {-# LANGUAGE CPP #-} @@ -23,6 +26,7 @@ module Math.LinearMap.Category.Instances where +import Math.VectorSpace.DimensionAware import Math.LinearMap.Category.Class import Data.VectorSpace @@ -42,6 +46,7 @@ import Data.Coerce import Data.Type.Coercion import Data.Tagged +import Data.Proxy import Data.Foldable (foldl') @@ -54,19 +59,30 @@ import Linear ( V0(V0), V1(V1), V2(V2), V3(V3), V4(V4) , _x, _y, _z, _w ) import Control.Lens ((^.)) +import Control.Monad.ST (ST) import qualified Data.Vector as Arr import qualified Data.Vector.Unboxed as UArr +import qualified Data.Vector.Generic as GArr import Math.LinearMap.Asserted import Math.VectorSpace.ZeroDimensional +import qualified Math.VectorSpace.DimensionAware.Theorems.MaybeNat as Maybe import qualified Test.QuickCheck as QC +import GHC.TypeNats (natVal) import qualified GHC.Exts as GHC import qualified GHC.Generics as GHC +import Data.Singletons (SingI, sing, Sing) +#if MIN_VERSION_singletons(3,0,0) +import GHC.TypeLits.Singletons (withKnownNat) +#else +import Data.Singletons.TypeLits (withKnownNat) +#endif + #if MIN_VERSION_manifolds_core(0,6,0) instance LinearSpace v => Semimanifold (EmptyMfd v) where type Needle (EmptyMfd v) = v @@ -116,7 +132,11 @@ fmapTensor = bilinearFunction $ \f (Tensor t) -> Tensor (f-+$>t); \ fzipTensorWith = bilinearFunction \ $ \(LinearFunction f) -> follow Tensor <<< f <<< flout Tensor *** flout Tensor; \ - coerceFmapTensorProduct _ VSCCoercion = VSCCoercion; \ + tensorUnsafeFromArrayWithOffset i ar \ + = Tensor (unsafeFromArrayWithOffset i ar); \ + tensorUnsafeWriteArrayWithOffset ar i (Tensor v) \ + = unsafeWriteArrayWithOffset ar i v; \ + coerceFmapTensorProduct _ VSCCoercion = Coercion; \ wellDefinedTensor (Tensor w) = Tensor <$> wellDefinedVector w }; \ instance LinearSpace (S) where { \ type DualVector (S) = (S); \ @@ -142,7 +162,31 @@ LinearScalarSpace(Float) LinearScalarSpace(Rational) +{-# INLINE tensorUnsafeFromArrayWithOffsetViaList #-} +tensorUnsafeFromArrayWithOffsetViaList + :: ∀ v w n m α . ( n`Dimensional`v + , m`Dimensional`w + , Scalar v ~ Scalar w + , GArr.Vector α (Scalar v) ) + => ([w] -> TensorProduct v w) -> Int -> α (Scalar v) -> (v⊗w) +tensorUnsafeFromArrayWithOffsetViaList l2v i ar + = Tensor $ l2v [ unsafeFromArrayWithOffset + (i + j * dimension @w) ar + | j <- [0 .. dimension @v - 1] ] +{-# INLINE tensorUnsafeWriteArrayWithOffsetViaList #-} +tensorUnsafeWriteArrayWithOffsetViaList + :: ∀ v w n m α σ . ( n`Dimensional`v + , m`Dimensional`w + , Scalar v ~ Scalar w + , GArr.Vector α (Scalar v) ) + => (TensorProduct v w -> [w]) -> GArr.Mutable α σ (Scalar v) + -> Int -> (v⊗w) -> ST σ () +tensorUnsafeWriteArrayWithOffsetViaList v2l ar i (Tensor t) + = forM_ (zip [0..] $ v2l t) $ \(j, v) + -> unsafeWriteArrayWithOffset ar + (i + j * dimension @w) v + #if MIN_VERSION_manifolds_core(0,6,0) #define FreeLinSpaceInteriorDecls #else @@ -150,13 +194,25 @@ toInterior = pure; fromInterior = id; translateP = Tagged (^+^); #endif -#define FreeLinearSpace(V, LV, tp, tenspl, tenid, dspan, contraction, contraaction) \ +#define FreeLinearSpace( V, d, LV, tp \ + , tenspl, tenid, dspan \ + , contraction, contraaction \ + , frls, tols ) \ instance Num s => Semimanifold (V s) where { \ type Needle (V s) = V s; \ FreeLinSpaceInteriorDecls \ (.+~^) = (^+^) }; \ instance Num s => PseudoAffine (V s) where { \ v.-~.w = pure (v^-^w); (.-~!) = (^-^) }; \ +instance ∀ s . (Num' s, Eq s) => DimensionAware (V s) where { \ + type StaticDimension (V s) = 'Just (d); \ + dimensionalityWitness = IsStaticDimensional }; \ +instance ∀ s . (Num' s, Eq s) => (d)`Dimensional`V (s) where { \ + unsafeFromArrayWithOffset \ + = unsafeFromArrayWithOffsetViaList (frls); \ + unsafeWriteArrayWithOffset \ + = unsafeWriteArrayWithOffsetViaList (tols) \ + }; \ instance ∀ s . (Num' s, Eq s) => TensorSpace (V s) where { \ type TensorProduct (V s) w = V w; \ scalarSpaceWitness = case closedScalarWitness :: ClosedScalarWitness s of{ \ @@ -179,7 +235,11 @@ fzipTensorWith = bilinearFunction $ \ \(LinearFunction f) (Tensor vw, Tensor vx) \ -> Tensor $ liftA2 (curry f) vw vx; \ - coerceFmapTensorProduct _ VSCCoercion = VSCCoercion; \ + tensorUnsafeFromArrayWithOffset \ + = tensorUnsafeFromArrayWithOffsetViaList (frls); \ + tensorUnsafeWriteArrayWithOffset \ + = tensorUnsafeWriteArrayWithOffsetViaList (tols); \ + coerceFmapTensorProduct _ VSCCoercion = Coercion; \ wellDefinedTensor = getTensorProduct >>> Hask.traverse wellDefinedVector \ >>> fmap Tensor }; \ instance ∀ s . (Num' s, Eq s) => LinearSpace (V s) where { \ @@ -207,28 +267,34 @@ applyTensorFunctional = bilinearFunction $ \(LinearMap f) (Tensor t) \ -> sum $ liftA2 (<.>^) f t; \ applyTensorLinMap = bilinearFunction $ \(LinearMap f) (Tensor t) \ - -> foldl' (^+^) zeroV $ liftA2 (arr fromTensor >>> \ + -> foldl' (^+^) zeroV $ liftA2 ((fromTensor-+$=>) >>> \ getLinearFunction . getLinearFunction applyLinear) f t; \ composeLinear = bilinearFunction $ \ \f (LinearMap g) -> LinearMap $ fmap ((applyLinear-+$>f)-+$>) g; \ useTupleLinearSpaceComponents _ = usingNonTupleTypeAsTupleError } -FreeLinearSpace( V0 +FreeLinearSpace( V0, 0 , LinearMap , \(Tensor V0) -> zeroV , \_ -> LinearMap V0 , V0 , LinearMap V0 , \V0 -> zeroV - , \V0 _ -> 0 ) -FreeLinearSpace( V1 + , \V0 _ -> 0 + , \[] -> V0 + , \V0 -> [] + ) +FreeLinearSpace( V1, 1 , LinearMap , \(Tensor (V1 w₀)) -> w₀⊗V1 1 , \w -> LinearMap $ V1 (Tensor $ V1 w) , V1 V1 , LinearMap . V1 . blockVectSpan $ V1 1 , \(V1 (V1 w)) -> w - , \(V1 x) f -> (f$x)^._x ) -FreeLinearSpace( V2 + , \(V1 x) f -> (f$x)^._x + , \[x] -> V1 x + , \(V1 x) -> [x] + ) +FreeLinearSpace( V2, 2 , LinearMap , \(Tensor (V2 w₀ w₁)) -> w₀⊗V2 1 0 ^+^ w₁⊗V2 0 1 @@ -239,8 +305,11 @@ (blockVectSpan $ V2 0 1) , \(V2 (V2 w₀ _) (V2 _ w₁)) -> w₀^+^w₁ - , \(V2 x y) f -> (f$x)^._x + (f$y)^._y ) -FreeLinearSpace( V3 + , \(V2 x y) f -> (f$x)^._x + (f$y)^._y + , \(x:y:[]) -> V2 x y + , \(V2 x y) -> (x:y:[]) + ) +FreeLinearSpace( V3, 3 , LinearMap , \(Tensor (V3 w₀ w₁ w₂)) -> w₀⊗V3 1 0 0 ^+^ w₁⊗V3 0 1 0 @@ -257,8 +326,11 @@ , \(V3 (V3 w₀ _ _) (V3 _ w₁ _) (V3 _ _ w₂)) -> w₀^+^w₁^+^w₂ - , \(V3 x y z) f -> (f$x)^._x + (f$y)^._y + (f$z)^._z ) -FreeLinearSpace( V4 + , \(V3 x y z) f -> (f$x)^._x + (f$y)^._y + (f$z)^._z + , \(x:y:z:[]) -> V3 x y z + , \(V3 x y z) -> x:y:z:[] + ) +FreeLinearSpace( V4, 4 , LinearMap , \(Tensor (V4 w₀ w₁ w₂ w₃)) -> w₀⊗V4 1 0 0 0 ^+^ w₁⊗V4 0 1 0 0 @@ -280,7 +352,10 @@ (V4 _ w₁ _ _) (V4 _ _ w₂ _) (V4 _ _ _ w₃)) -> w₀^+^w₁^+^w₂^+^w₃ - , \(V4 x y z w) f -> (f$x)^._x + (f$y)^._y + (f$z)^._z + (f$w)^._w ) + , \(V4 x y z w) f -> (f$x)^._x + (f$y)^._y + (f$z)^._z + (f$w)^._w + , \(x:y:z:w:[]) -> V4 x y z w + , \(V4 x y z w) -> x:y:z:w:[] + ) @@ -310,6 +385,9 @@ instance (Num' n, UArr.Unbox n) => PseudoAffine (FinSuppSeq n) where v.-~.w = Just $ v.-.w; (.-~!) = (.-.) +instance (Num' n, UArr.Unbox n) => DimensionAware (FinSuppSeq n) where + type StaticDimension (FinSuppSeq n) = 'Nothing + dimensionalityWitness = IsFlexibleDimensional instance (Num' n, UArr.Unbox n) => TensorSpace (FinSuppSeq n) where type TensorProduct (FinSuppSeq n) v = [v] wellDefinedVector (FinSuppSeq v) = FinSuppSeq <$> UArr.mapM wellDefinedVector v @@ -333,7 +411,11 @@ fmapTensor = bilinearFunction $ \f (Tensor a) -> Tensor $ map (f$) a fzipTensorWith = bilinearFunction $ \f (Tensor a, Tensor b) -> Tensor $ zipWith (curry $ arr f) a b - coerceFmapTensorProduct _ VSCCoercion = VSCCoercion + tensorUnsafeFromArrayWithOffset + = notStaticDimensionalContradiction @(FinSuppSeq n) + tensorUnsafeWriteArrayWithOffset + = notStaticDimensionalContradiction @(FinSuppSeq n) + coerceFmapTensorProduct _ VSCCoercion = Coercion wellDefinedTensor (Tensor a) = Tensor <$> Hask.traverse wellDefinedVector a @@ -347,6 +429,9 @@ instance (Num' n, UArr.Unbox n) => PseudoAffine (Sequence n) where v.-~.w = Just $ v.-.w; (.-~!) = (.-.) +instance (Num' n, UArr.Unbox n) => DimensionAware (Sequence n) where + type StaticDimension (Sequence n) = 'Nothing + dimensionalityWitness = IsFlexibleDimensional instance (Num' n, UArr.Unbox n) => TensorSpace (Sequence n) where type TensorProduct (Sequence n) v = [v] wellDefinedVector (SoloChunk n c) = SoloChunk n <$> UArr.mapM wellDefinedVector c @@ -373,14 +458,18 @@ fmapTensor = bilinearFunction $ \f (Tensor a) -> Tensor $ map (f$) a fzipTensorWith = bilinearFunction $ \f (Tensor a, Tensor b) -> Tensor $ zipWith (curry $ arr f) a b - coerceFmapTensorProduct _ VSCCoercion = VSCCoercion + tensorUnsafeFromArrayWithOffset + = notStaticDimensionalContradiction @(Sequence n) + tensorUnsafeWriteArrayWithOffset + = notStaticDimensionalContradiction @(Sequence n) + coerceFmapTensorProduct _ VSCCoercion = Coercion -instance (Num' n, UArr.Unbox n) => LinearSpace (Sequence n) where +instance ∀ n . (Num' n, UArr.Unbox n) => LinearSpace (Sequence n) where type DualVector (Sequence n) = FinSuppSeq n dualSpaceWitness = case closedScalarWitness :: ClosedScalarWitness n of ClosedScalarWitness -> DualSpaceWitness linearId = LinearMap [basisValue i | i<-[0..]] - tensorId = LinearMap [asTensor $ fmap (LinearFunction $ + tensorId = LinearMap [asTensor -+$=> fmap (LinearFunction $ \w -> Tensor $ replicate (i-1) zeroV ++ [w]) $ id | i<-[0..]] applyDualVector = bilinearFunction $ adv Seq.minimumChunkSize where adv _ (FinSuppSeq v) (Seq.SoloChunk o q) @@ -404,12 +493,14 @@ \(LinearMap m) (Tensor t) -> sumV $ zipWith (getLinearFunction . getLinearFunction applyLinear) m t useTupleLinearSpaceComponents _ = usingNonTupleTypeAsTupleError -instance (Num' n, UArr.Unbox n) => LinearSpace (FinSuppSeq n) where + coerceDoubleDual = case scalarSpaceWitness @n of + ScalarSpaceWitness -> VSCCoercion +instance ∀ n . (Num' n, UArr.Unbox n) => LinearSpace (FinSuppSeq n) where type DualVector (FinSuppSeq n) = Sequence n dualSpaceWitness = case closedScalarWitness :: ClosedScalarWitness n of ClosedScalarWitness -> DualSpaceWitness linearId = LinearMap [basisValue i | i<-[0..]] - tensorId = LinearMap [asTensor $ fmap (LinearFunction $ + tensorId = LinearMap [asTensor -+$=> fmap (LinearFunction $ \w -> Tensor $ replicate (i-1) zeroV ++ [w]) $ id | i<-[0..]] applyDualVector = bilinearFunction $ adv Seq.minimumChunkSize where adv _ (Seq.SoloChunk o q) (FinSuppSeq v) @@ -427,6 +518,8 @@ \(LinearMap m) (Tensor t) -> sumV $ zipWith (getLinearFunction . getLinearFunction applyLinear) m t useTupleLinearSpaceComponents _ = usingNonTupleTypeAsTupleError + coerceDoubleDual = case scalarSpaceWitness @n of + ScalarSpaceWitness -> VSCCoercion @@ -471,6 +564,19 @@ instance (TensorSpace v, Scalar v ~ s) => PseudoAffine (SymmetricTensor s v) where (.-~!) = (^-^) p.-~.q = pure (p^-^q) +instance ∀ s v . (Num' s, TensorSpace v, Scalar v ~ s) + => DimensionAware (SymmetricTensor s v) where + type StaticDimension (SymmetricTensor s v) + = Maybe.FmapTriangularNum (StaticDimension v) + dimensionalityWitness = case dimensionalityWitness @v of + IsFlexibleDimensional -> IsFlexibleDimensional + IsStaticDimensional + -> withKnownNat (Maybe.triangularNumSing (dimensionalitySing @v)) + IsStaticDimensional +instance ∀ s v n m . ( Num' s, n`Dimensional`v, TensorSpace v, Scalar v ~ s + , m ~ Maybe.TriangularNum n ) + => m`Dimensional`(SymmetricTensor s v) where + knownDimensionalitySing = Maybe.triangularNumSing $ dimensionalitySing @v instance (Num' s, TensorSpace v, Scalar v ~ s) => TensorSpace (SymmetricTensor s v) where type TensorProduct (SymmetricTensor s v) x = Tensor s v (Tensor s v x) wellDefinedVector (SymTensor t) = SymTensor <$> wellDefinedVector t @@ -491,45 +597,61 @@ tensorProduct = bilinearFunction $ \(SymTensor t) g -> Tensor $ fmap (LinearFunction (⊗g)) $ t transposeTensor = LinearFunction $ \(Tensor f) -> getLinearFunction ( - arr (fmap VSCCoercion) . transposeTensor . arr lassocTensor) f + undefined -- arr (fmap VSCCoercion) + . transposeTensor . arr lassocTensor) f fmapTensor = bilinearFunction $ \f (Tensor t) -> Tensor $ fmap (fmap f) $ t fzipTensorWith = bilinearFunction $ \f (Tensor s, Tensor t) -> Tensor $ fzipWith (fzipWith f) $ (s,t) - coerceFmapTensorProduct _ crc = fmap (fmap crc) + coerceFmapTensorProduct _ crc = undefined -- case fmap (fmap crc) :: VSCCoercion of + -- VSCCoercion -> Coercion wellDefinedTensor (Tensor t) = Tensor <$> wellDefinedVector t -instance (Num' s, LinearSpace v, Scalar v ~ s) => LinearSpace (SymmetricTensor s v) where +instance ∀ s v . (Num' s, LinearSpace v, Scalar v ~ s) + => LinearSpace (SymmetricTensor s v) where type DualVector (SymmetricTensor s v) = SymmetricTensor s (DualVector v) dualSpaceWitness = case ( closedScalarWitness :: ClosedScalarWitness s , dualSpaceWitness :: DualSpaceWitness v ) of (ClosedScalarWitness, DualSpaceWitness) -> DualSpaceWitness linearId = case dualSpaceWitness :: DualSpaceWitness v of - DualSpaceWitness -> LinearMap $ rassocTensor . asTensor - . fmap (unsafeFollowVSC SymTensor . asTensor) $ id - tensorId = LinearMap $ asTensor . fmap asTensor . curryLinearMap - . fmap asTensor - . curryLinearMap - . fmap (unsafeFollowVSC $ \t -> Tensor $ rassocTensor $ t) - $ id + DualSpaceWitness -> LinearMap undefined -- $ rassocTensor . asTensor + -- . fmap (unsafeFollowVSC SymTensor . asTensor) $ id + tensorId = LinearMap undefined -- $ asTensor . fmap asTensor . curryLinearMap + -- . fmap asTensor + -- . curryLinearMap + -- . fmap (unsafeFollowVSC $ \t -> Tensor $ rassocTensor $ t) + -- $ id applyLinear = case dualSpaceWitness :: DualSpaceWitness v of DualSpaceWitness -> bilinearFunction $ \(LinearMap f) (SymTensor t) -> (getLinearFunction applyLinear $ fromTensor . deferLinearMap . asLinearMap $ f) $ t applyDualVector = bilinearFunction $ \(SymTensor f) (SymTensor v) -> getLinearFunction - (getLinearFunction applyDualVector $ fromTensor $ f) v - applyTensorFunctional = case dualSpaceWitness :: DualSpaceWitness v of - DualSpaceWitness -> bilinearFunction $ \(LinearMap f) (Tensor t) + (getLinearFunction applyDualVector $ fromTensor -+$=> f) v + applyTensorFunctional :: ∀ u . (LinearSpace u, Scalar u ~ s) + => LinearFunction s + (LinearMap s (SymmetricTensor s v) (DualVector u)) + (LinearFunction s (Tensor s (SymmetricTensor s v) u) s) + applyTensorFunctional = case (dualSpaceWitness @v, dualSpaceWitness @u) of + (DualSpaceWitness, DualSpaceWitness) + -> bilinearFunction $ \(LinearMap f) (Tensor t) -> getLinearFunction (getLinearFunction applyTensorFunctional - $ fromTensor . fmap fromTensor $ f) t - applyTensorLinMap = case dualSpaceWitness :: DualSpaceWitness v of - DualSpaceWitness -> bilinearFunction $ \(LinearMap (Tensor f)) (Tensor t) + $ fromTensor . fmap fromTensor -+$=> f) t + applyTensorLinMap :: ∀ u w . ( LinearSpace u, Scalar u ~ s + , TensorSpace w, Scalar w ~ s ) + => LinearFunction s + (LinearMap s (Tensor s (SymmetricTensor s v) u) w) + (LinearFunction s (Tensor s (SymmetricTensor s v) u) w) + applyTensorLinMap = case (dualSpaceWitness @v, dualSpaceWitness @u) of + (DualSpaceWitness, DualSpaceWitness) + -> bilinearFunction $ \(LinearMap (Tensor f)) (Tensor t) -> getLinearFunction (getLinearFunction applyTensorLinMap $ uncurryLinearMap . fmap (uncurryLinearMap . fromTensor . fmap fromTensor) - $ LinearMap f) t + -+$=> LinearMap f) t useTupleLinearSpaceComponents _ = usingNonTupleTypeAsTupleError + coerceDoubleDual = case (dualSpaceWitness @v, scalarSpaceWitness @s) of + (DualSpaceWitness, ScalarSpaceWitness) -> VSCCoercion @@ -546,7 +668,7 @@ type v⊗〃+>w = LinearMap (Scalar v) (SymmetricTensor (Scalar v) v) w currySymBilin :: LinearSpace v => (v⊗〃+>w) -+> (v+>(v+>w)) -currySymBilin = LinearFunction . arr $ fmap fromTensor . fromTensor . VSCCoercion +currySymBilin = undefined -- LinearFunction . arr $ fmap fromTensor . fromTensor . VSCCoercion @@ -619,6 +741,42 @@ => InnerSpace (LinearMap ℝ ℝ v) where LinearMap f <.> LinearMap g = f<.>g +instance ( TensorSpace u, TensorSpace v, TensorSpace w + , Num s, Scalar u ~ s, Scalar v ~ s, Scalar w ~ s + , InnerSpace (Tensor s u w), InnerSpace (Tensor s v w) ) + => InnerSpace (Tensor s (u,v) w) where + Tensor (uw,vw) <.> Tensor (uw',vw') = uw<.>uw' + vw<.>vw' +instance ( LinearSpace u, LinearSpace v, TensorSpace w + , Num s, Scalar u ~ s, Scalar v ~ s, Scalar w ~ s + , InnerSpace (LinearMap s u w), InnerSpace (LinearMap s v w) ) + => InnerSpace (LinearMap s (u,v) w) where + (<.>) = case (dualSpaceWitness @u, dualSpaceWitness @v) of + (DualSpaceWitness, DualSpaceWitness) + -> \(LinearMap (uw,vw)) (LinearMap (uw',vw')) + -> (asLinearMap$uw)<.>(asLinearMap$uw') + + (asLinearMap$vw)<.>(asLinearMap$vw') + +instance ( TensorSpace u, TensorSpace v, TensorSpace w + , Num s, Scalar u ~ s, Scalar v ~ s, Scalar w ~ s + , InnerSpace (Tensor s u (Tensor s v w)) ) + => InnerSpace (Tensor s (Tensor s u v) w) where + s <.> t = (rassocTensor$s)<.>(rassocTensor$t) +instance ( LinearSpace u, TensorSpace v, TensorSpace w + , Num s, Scalar u ~ s, Scalar v ~ s, Scalar w ~ s + , InnerSpace (LinearMap s u (Tensor s v w)) ) + => InnerSpace (Tensor s (LinearMap s u v) w) where + s <.> t = (hasteLinearMap$s)<.>(hasteLinearMap$t) +instance ( LinearSpace u, LinearSpace v, TensorSpace w + , Num s, Scalar u ~ s, Scalar v ~ s, Scalar w ~ s + , InnerSpace (LinearMap s u (LinearMap s v w)) ) + => InnerSpace (LinearMap s (Tensor s u v) w) where + s <.> t = (curryLinearMap$s)<.>(curryLinearMap$t) +instance ( LinearSpace u, LinearSpace v, TensorSpace w + , Num s, Scalar u ~ s, Scalar v ~ s, Scalar w ~ s + , InnerSpace (Tensor s u (LinearMap s v w)) ) + => InnerSpace (LinearMap s (LinearMap s u v) w) where + s <.> t = (coCurryLinearMap$s)<.>(coCurryLinearMap$t) + instance (Show v) => Show (Tensor ℝ ℝ v) where showsPrec p (Tensor t) = showParen (p>9) $ ("Tensor "++) . showsPrec 10 t @@ -641,3 +799,51 @@ FreeArbitrarySpace(V2) FreeArbitrarySpace(V3) FreeArbitrarySpace(V4) + +instance ( QC.Arbitrary (Tensor s u w), QC.Arbitrary (Tensor s v w) + , Scalar u ~ s, Scalar v ~ s, Scalar w ~ s ) + => QC.Arbitrary (Tensor s (u,v) w) where + arbitrary = Tensor <$> QC.arbitrary + shrink (Tensor t) = Tensor <$> QC.shrink t + +instance ( LinearSpace u, LinearSpace v, TensorSpace w + , QC.Arbitrary (LinearMap s u w), QC.Arbitrary (LinearMap s v w) + , Scalar u ~ s, Scalar v ~ s, Scalar w ~ s ) + => QC.Arbitrary (LinearMap s (u,v) w) where + arbitrary = case (dualSpaceWitness @u, dualSpaceWitness @v) of + (DualSpaceWitness, DualSpaceWitness) -> LinearMap <$> do + (,) <$> (arr fromLinearMap <$> QC.arbitrary) + <*> (arr fromLinearMap <$> QC.arbitrary) + shrink = case (dualSpaceWitness @u, dualSpaceWitness @v) of + (DualSpaceWitness, DualSpaceWitness) -> \(LinearMap (x,y)) -> LinearMap <$> do + (x',y') <- QC.shrink (asLinearMap $ x, asLinearMap $ y) + return (fromLinearMap $ x', fromLinearMap $ y') + +instance ( TensorSpace u, TensorSpace v, TensorSpace w + , QC.Arbitrary (u⊗(v⊗w)) + , Scalar u ~ s, Scalar v ~ s, Scalar w ~ s ) + => QC.Arbitrary (Tensor s (Tensor s u v) w) where + arbitrary = arr lassocTensor <$> QC.arbitrary + shrink (Tensor t) = arr lassocTensor <$> QC.shrink (Tensor t) + +instance ( LinearSpace u, LinearSpace v, TensorSpace w + , QC.Arbitrary (u+>(v+>w)) + , Scalar u ~ s, Scalar v ~ s, Scalar w ~ s ) + => QC.Arbitrary (LinearMap s (Tensor s u v) w) where + arbitrary = arr uncurryLinearMap <$> QC.arbitrary + shrink f = arr uncurryLinearMap <$> QC.shrink (curryLinearMap $ f) + +instance ( LinearSpace u, TensorSpace v, TensorSpace w + , QC.Arbitrary (u+>(v⊗w)) + , Scalar u ~ s, Scalar v ~ s, Scalar w ~ s ) + => QC.Arbitrary (Tensor s (LinearMap s u v) w) where + arbitrary = arr deferLinearMap <$> QC.arbitrary + shrink (Tensor t) = arr deferLinearMap <$> QC.shrink (LinearMap t) + +instance ( LinearSpace u, LinearSpace v, TensorSpace w + , QC.Arbitrary (u⊗(v+>w)) + , Scalar u ~ s, Scalar v ~ s, Scalar w ~ s ) + => QC.Arbitrary (LinearMap s (LinearMap s u v) w) where + arbitrary = arr coUncurryLinearMap <$> QC.arbitrary + shrink f = arr coUncurryLinearMap <$> QC.shrink (coCurryLinearMap $ f) + diff --git a/Math/LinearMap/Category/Instances/Deriving.hs b/Math/LinearMap/Category/Instances/Deriving.hs --- a/Math/LinearMap/Category/Instances/Deriving.hs +++ b/Math/LinearMap/Category/Instances/Deriving.hs @@ -11,11 +11,15 @@ {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE TypeOperators #-} +{-# LANGUAGE NoStarIsType #-} {-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FunctionalDependencies #-} {-# LANGUAGE AllowAmbiguousTypes #-} {-# LANGUAGE TypeApplications #-} +{-# LANGUAGE DataKinds #-} {-# LANGUAGE Rank2Types #-} {-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE InstanceSigs #-} {-# LANGUAGE UnicodeSyntax #-} {-# LANGUAGE StandaloneDeriving #-} {-# LANGUAGE DeriveGeneric #-} @@ -40,6 +44,7 @@ import Math.LinearMap.Category.Class import Math.LinearMap.Category.Instances +import Math.VectorSpace.DimensionAware import Math.VectorSpace.Docile import Data.VectorSpace @@ -49,27 +54,45 @@ import Data.Tree (Forest) import Data.MemoTrie import Data.Hashable +import Data.Void import Prelude () import qualified Prelude as Hask -import Control.Category.Constrained.Prelude +import Control.Category.Constrained.Prelude hiding (type (+)) import Control.Arrow.Constrained +import Control.Monad.ST (ST) import Data.Coerce import Data.Type.Coercion import Data.Tagged +import Data.Proxy import qualified Data.Kind as Kind import Data.Traversable (traverse) import Data.Default.Class +import qualified Data.Vector.Generic as GArr +import qualified Data.Vector.Generic.Mutable as GMArr + import Math.Manifold.Core.PseudoAffine import Math.LinearMap.Asserted import Math.VectorSpace.ZeroDimensional import Data.VectorSpace.Free import GHC.Generics (Generic) +import GHC.TypeLits (Nat, KnownNat, natVal, type (+), type (*)) +#if MIN_VERSION_singletons(3,0,0) +import GHC.TypeLits.Singletons (withKnownNat) +import Data.Maybe.Singletons (FromJust) +import Prelude.Singletons +#else +import Data.Singletons.TypeLits (withKnownNat) +import Data.Singletons.Prelude +#endif + (SingI, sing, withSingI, SMaybe(..), SNum(..)) +import qualified Math.VectorSpace.DimensionAware.Theorems.MaybeNat as Maybe + import Language.Haskell.TH import Language.Haskell.TH.Syntax (Name(..), OccName(..) #if MIN_VERSION_template_haskell(2,17,0) @@ -117,10 +140,23 @@ makeLinearSpaceFromBasis v = makeLinearSpaceFromBasis' def $ deQuantifyType v -data LinearSpaceFromBasisDerivationConfig = LinearSpaceFromBasisDerivationConfig +data LinearSpaceFromBasisDerivationConfig + = LinearSpaceFromBasisDerivationConfig + { _treatBasisAsFinite :: Bool + } instance Default LinearSpaceFromBasisDerivationConfig where def = LinearSpaceFromBasisDerivationConfig + { _treatBasisAsFinite = False } + +requireExtensions :: [Extension] -> Q () +requireExtensions reqExts = do + exts <- extsEnabled + forM_ reqExts $ \re -> do + if re`elem`exts + then return () + else reportError $ "This macro requires -X"++show re++"." + -- | More general version of 'makeLinearSpaceFromBasis', that can be used with -- parameterised types. makeLinearSpaceFromBasis' :: LinearSpaceFromBasisDerivationConfig @@ -129,17 +165,15 @@ Specificity #endif ], Cxt, Type) -> DecsQ -makeLinearSpaceFromBasis' _ cxtv = do +makeLinearSpaceFromBasis' config cxtv = do (cxt,v) <- do (_, cxt', v') <- cxtv return (pure cxt', pure v') - exts <- extsEnabled - if not $ all (`elem`exts) [TypeFamilies, ScopedTypeVariables, TypeApplications] - then reportError "This macro requires -XTypeFamilies, -XScopedTypeVariables and -XTypeApplications." - else pure () + requireExtensions [ TypeFamilies, MultiParamTypeClasses + , ScopedTypeVariables, TypeApplications ] - sequence + sequence ( [ InstanceD Nothing <$> cxt <*> [t|Semimanifold $v|] <*> [d| type instance Needle $v = $v #if !MIN_VERSION_manifolds_core(0,6,0) @@ -161,6 +195,15 @@ $(varP '(.+^)) = (^+^) $(varP '(.-.)) = (^-^) |] + , if _treatBasisAsFinite config + then InstanceD Nothing <$> cxt <*> [t|DimensionAware $v|] <*> [d| + type instance StaticDimension $v = Cardinality (Basis $v) + $(varP 'dimensionalityWitness) = IsStaticDimensional + |] + else InstanceD Nothing <$> cxt <*> [t|DimensionAware $v|] <*> [d| + type instance StaticDimension $v = 'Nothing + $(varP 'dimensionalityWitness) = IsFlexibleDimensional + |] , InstanceD Nothing <$> cxt <*> [t|TensorSpace $v|] <*> [d| type instance TensorProduct $v w = Basis $v :->: w $(varP 'wellDefinedVector) = \v @@ -191,6 +234,16 @@ $(varP 'fzipTensorWith) = bilinearFunction $ \(LinearFunction f) (Tensor tv, Tensor tw) -> Tensor $ liftA2 (curry f) tv tw + $(varP 'tensorUnsafeFromArrayWithOffset) + = \i ar -> Tensor . trie + $ \bv -> let w = unsafeFromArrayWithOffset + (i + dimensionOf w * lookupBasisIndex bv) + ar + in w + $(varP 'tensorUnsafeWriteArrayWithOffset) + = \ar i (Tensor t) -> forM_ (zip [0..] enumBasis) $ \(j,bv) -> do + let w = untrie t bv + unsafeWriteArrayWithOffset ar (i + dimensionOf w * j) w $(varP 'coerceFmapTensorProduct) = \_ VSCCoercion -> error "Cannot yet coerce tensors defined from a `HasBasis` instance. This would require `RoleAnnotations` on `:->:`. Cf. https://gitlab.haskell.org/ghc/ghc/-/issues/8177" |] @@ -243,7 +296,20 @@ $(varP 'useTupleLinearSpaceComponents) = \_ -> usingNonTupleTypeAsTupleError |] - ] + ] ++ if _treatBasisAsFinite config then [do + dim <- pure . VarT <$> newName "n" + InstanceD Nothing <$> ((:)<$>[t|($dim)#(Basis $v)|]<*>cxt) + <*> [t|Dimensional $dim $v|] <*> [d| + $(varP 'unsafeFromArrayWithOffset) + = \i ar -> recompose + [ (b, ar GArr.! (i+j)) + | (j, b) <- zip [0..] enumBasis ] + $(varP 'unsafeWriteArrayWithOffset) + = \ar i v -> forM_ (zip [0..] enumBasis) $ \(j,b) -> + GMArr.unsafeWrite ar (i+j) $ decompose' v b + |]] + else [] + ) data FiniteDimensionalFromBasisDerivationConfig = FiniteDimensionalFromBasisDerivationConfig @@ -263,7 +329,7 @@ #endif ], Cxt, Type) -> DecsQ makeFiniteDimensionalFromBasis' _ cxtv = do - generalInsts <- makeLinearSpaceFromBasis' def cxtv + generalInsts <- makeLinearSpaceFromBasis' def{_treatBasisAsFinite=True} cxtv (cxt,v) <- do (_, cxt', v') <- cxtv return (pure cxt', pure v') @@ -406,7 +472,80 @@ (.-~!) = (^-^) p.-~.q = pure (p^-^q) +type family Cardinality b :: Maybe Nat + +type instance Cardinality Void = 'Just 0 +type instance Cardinality () = 'Just 1 +type instance Cardinality (Either a b) + = Maybe.ZipWithPlus (Cardinality a) (Cardinality b) +type instance Cardinality (a,b) + = Maybe.ZipWithTimes (Cardinality a) (Cardinality b) + +class (KnownNat n, KnownCardinality b, Cardinality b ~ 'Just n) => n#b where + enumBasis :: [b] + lookupBasisIndex :: b -> Int + +instance 0#Void where + enumBasis = [] + lookupBasisIndex = absurd +instance 1#() where + enumBasis = [()] + lookupBasisIndex () = 0 +instance (n#a, m#b, KnownNat nm, nm~(n+m)) => nm # Either a b where + enumBasis = (Left<$>enumBasis)++(Right<$>enumBasis) + lookupBasisIndex (Left bl) = lookupBasisIndex bl + lookupBasisIndex (Right br) = fromIntegral (natVal @n Proxy) + lookupBasisIndex br +instance (n#a, m#b, KnownNat nm, nm~(n*m)) => nm # (a,b) where + enumBasis = (,)<$>enumBasis<*>enumBasis + lookupBasisIndex (ba,bb) = fromIntegral (natVal @m Proxy) * lookupBasisIndex ba + + lookupBasisIndex bb + +type FiniteCardinality b = FromJust (Cardinality b) + +data CardinalityWitness b where + FiniteCardinality :: n#b => CardinalityWitness b + NonfiniteCardinality :: Cardinality b ~ 'Nothing => CardinalityWitness b + +class KnownCardinality b where + cardinalityWitness :: CardinalityWitness b + +instance KnownCardinality Void where cardinalityWitness = FiniteCardinality +instance KnownCardinality () where cardinalityWitness = FiniteCardinality +instance ∀ a b . (KnownCardinality a, KnownCardinality b) + => KnownCardinality (Either a b) where + cardinalityWitness = case (cardinalityWitness @a, cardinalityWitness @b) of + (FiniteCardinality, FiniteCardinality) + -> withKnownNat (sing @(FiniteCardinality a)%+sing @(FiniteCardinality b)) + FiniteCardinality + (NonfiniteCardinality, _) -> NonfiniteCardinality + (_, NonfiniteCardinality) -> NonfiniteCardinality +instance ∀ a b . (KnownCardinality a, KnownCardinality b) + => KnownCardinality (a,b) where + cardinalityWitness = case (cardinalityWitness @a, cardinalityWitness @b) of + (FiniteCardinality, FiniteCardinality) + -> withKnownNat (sing @(FiniteCardinality a)%*sing @(FiniteCardinality b)) + FiniteCardinality + (NonfiniteCardinality, _) -> NonfiniteCardinality + (_, NonfiniteCardinality) -> NonfiniteCardinality + +instance (HasBasis v, KnownCardinality (Basis v)) + => DimensionAware (DualVectorFromBasis v) where + type StaticDimension (DualVectorFromBasis v) = Cardinality (Basis v) + dimensionalityWitness = case cardinalityWitness @(Basis v) of + NonfiniteCardinality -> IsFlexibleDimensional + FiniteCardinality -> IsStaticDimensional +instance ( HasBasis v, n#Basis v ) + => n`Dimensional`DualVectorFromBasis v where + unsafeFromArrayWithOffset i ar + = recompose [ (b, ar GArr.! (i+j)) + | (j, b) <- zip [0..] enumBasis ] + unsafeWriteArrayWithOffset ar i v + = forM_ (zip [0..] enumBasis) $ \(j,b) -> do + GMArr.unsafeWrite ar (i+j) $ decompose' v b + + instance ∀ v . ( HasBasis v, Num' (Scalar v) + , KnownCardinality (Basis v) , Scalar (Scalar v) ~ Scalar v , HasTrie (Basis v) , Eq v ) @@ -439,6 +578,24 @@ fzipTensorWith = bilinearFunction $ \(LinearFunction f) (Tensor tv, Tensor tw) -> Tensor $ liftA2 (curry f) tv tw + tensorUnsafeFromArrayWithOffset :: ∀ w m α + . (m`Dimensional`w, Scalar w ~ Scalar v, GArr.Vector α (Scalar v)) + => Int -> α (Scalar v) + -> Tensor (Scalar v) (DualVectorFromBasis v) w + tensorUnsafeFromArrayWithOffset i ar + = case cardinalityWitness @(Basis v) of + FiniteCardinality -> Tensor . trie + $ lookupBasisIndex >>> \j + -> unsafeFromArrayWithOffset (i + j * dimension @w) ar + tensorUnsafeWriteArrayWithOffset :: ∀ w m α σ + . (m`Dimensional`w, Scalar w ~ Scalar v, GArr.Vector α (Scalar v)) + => GArr.Mutable α σ (Scalar v) -> Int + -> Tensor (Scalar v) (DualVectorFromBasis v) w -> ST σ () + tensorUnsafeWriteArrayWithOffset ar i + = case cardinalityWitness @(Basis v) of + FiniteCardinality -> \(Tensor t) + -> forM_ (zip [0..] enumBasis) $ \(j, bv) -> + unsafeWriteArrayWithOffset ar (i + j * dimension @w) $ untrie t bv coerceFmapTensorProduct _ VSCCoercion = error "Cannot yet coerce tensors defined from a `HasBasis` instance. This would require `RoleAnnotations` on `:->:`. Cf. https://gitlab.haskell.org/ghc/ghc/-/issues/8177" @@ -446,6 +603,9 @@ -- | Do not manually instantiate this class. It is used internally -- by 'makeLinearSpaceFromBasis'. class ( HasBasis v, Num' (Scalar v) + , KnownCardinality (Basis v) + , StaticDimension v ~ (Cardinality (Basis v)) + , SingI (StaticDimension v) , LinearSpace v, DualVector v ~ DualVectorFromBasis v) => BasisGeneratedSpace v where proveTensorProductIsTrie @@ -507,6 +667,7 @@ | (i, Tensor fi) <- enumerate f ] ) useTupleLinearSpaceComponents _ = usingNonTupleTypeAsTupleError + coerceDoubleDual = VSCCoercion zipWith' :: (a -> b -> c) -> [a] -> [b] -> ([c], [b]) @@ -621,9 +782,13 @@ #if !MIN_VERSION_manifolds_core(0,6,0) , Semimanifold v, Interior v ~ v #endif + , StaticDimension v ~ StaticDimension (VectorSpaceImplementation v) ) => AbstractTensorSpace v where + sameDimensionalInAbstraction + :: n`Dimensional`VectorSpaceImplementation v + => (n`Dimensional`v => ρ) -> ρ abstractTensorProductsCoercion - :: VSCCoercion (TensorProduct v w) + :: Coercion (TensorProduct v w) (TensorProduct (VectorSpaceImplementation v) w) class ( AbstractTensorSpace v, LinearSpace (VectorSpaceImplementation v) @@ -640,31 +805,39 @@ = case dualSpaceWitness @(VectorSpaceImplementation v) of DualSpaceWitness -> scalarsSameInAbstraction @v φ -abstractDualVectorCoercion :: ∀ a - . VSCCoercion (AbstractDualVector a (VectorSpaceImplementation a)) - (DualVector (VectorSpaceImplementation a)) -abstractDualVectorCoercion = VSCCoercion +abstractDualVectorCoercion :: ∀ a s + . ( LinearSpace (VectorSpaceImplementation a) + , Scalar (VectorSpaceImplementation a) ~ s ) + => VSCCoercion s (AbstractDualVector a (VectorSpaceImplementation a)) + (DualVector (VectorSpaceImplementation a)) +abstractDualVectorCoercion = case dualSpaceWitness @(VectorSpaceImplementation a) of + DualSpaceWitness -> VSCCoercion -abstractTensorsCoercion :: ∀ a c w +abstractDualTensorsCoercion :: ∀ a s c w . ( AbstractVectorSpace a, LinearSpace c , c ~ VectorSpaceImplementation a, TensorSpace w ) - => VSCCoercion (AbstractDualVector a c⊗w) (DualVector c⊗w) -abstractTensorsCoercion = VSCCoercion + => VSCCoercion s (AbstractDualVector a c⊗w) (DualVector c⊗w) +abstractDualTensorsCoercion = case dualSpaceWitness @c of + DualSpaceWitness -> VSCCoercion -abstractLinmapCoercion :: ∀ a c w +abstractLinmapCoercion :: ∀ a s c w . ( AbstractLinearSpace a, LinearSpace c , c ~ VectorSpaceImplementation a, TensorSpace w ) - => VSCCoercion (AbstractDualVector a c+>w) (DualVector c+>w) + => VSCCoercion s (AbstractDualVector a c+>w) (DualVector c+>w) abstractLinmapCoercion = case ( dualSpaceWitness @c , abstractTensorProductsCoercion @a @w ) of - (DualSpaceWitness, VSCCoercion) -> VSCCoercion + (DualSpaceWitness, Coercion) -> VSCCoercion -coerceLinearMapCodomain :: ∀ v w x . ( LinearSpace v, Coercible w x ) +coerceLinearMapCodomain :: ∀ v w x + . ( LinearSpace v + , TensorSpace w, TensorSpace x, Coercible w x + , Scalar w ~ Scalar v, Scalar x ~ Scalar w + , StaticDimension w ~ StaticDimension x ) => (v+>w) -> (v+>x) coerceLinearMapCodomain = case dualSpaceWitness @v of DualSpaceWitness -> \(LinearMap m) -> LinearMap $ (coerceFmapTensorProduct ([]::[DualVector v]) - (VSCCoercion :: VSCCoercion w x) $ m) + (VSCCoercion :: VSCCoercion (Scalar v) w x) $ m) instance (Show (DualVector c)) => Show (AbstractDualVector a c) where showsPrec p (AbstractDualVector_ φ) = showParen (p>10) @@ -711,6 +884,30 @@ instance ∀ a c . ( AbstractLinearSpace a, VectorSpaceImplementation a ~ c , TensorSpace (DualVector c) ) + => DimensionAware (AbstractDualVector a c) where + type StaticDimension (AbstractDualVector a c) = StaticDimension c + dimensionalityWitness = case dimensionalityWitness @c of + IsStaticDimensional -> IsStaticDimensional + IsFlexibleDimensional -> IsFlexibleDimensional +instance ∀ n a c . ( AbstractLinearSpace a + , VectorSpaceImplementation a ~ c + , n`Dimensional`c + , TensorSpace (DualVector c) ) + => n`Dimensional`AbstractDualVector a c where + knownDimensionalitySing = dimensionalitySing @c + unsafeFromArrayWithOffset i + = scalarsSameInAbstraction @a ( + case (dualSpaceWitness @c, dimensionalityWitness @(DualVector c)) of + (DualSpaceWitness, IsStaticDimensional) + -> AbstractDualVector_ . unsafeFromArrayWithOffset i ) + unsafeWriteArrayWithOffset ar i + = scalarsSameInAbstraction @a ( + case (dualSpaceWitness @c, dimensionalityWitness @(DualVector c)) of + (DualSpaceWitness, IsStaticDimensional) + -> \(AbstractDualVector_ v) -> unsafeWriteArrayWithOffset ar i v ) + +instance ∀ a c . ( AbstractLinearSpace a, VectorSpaceImplementation a ~ c + , TensorSpace (DualVector c) ) => TensorSpace (AbstractDualVector a c) where type TensorProduct (AbstractDualVector a c) w = TensorProduct (DualVector c) w @@ -773,9 +970,10 @@ where tt :: ∀ w . ( TensorSpace w, Scalar w ~ Scalar a , Scalar (DualVector c) ~ Scalar a ) => (AbstractDualVector a c ⊗ w) -+> (w ⊗ AbstractDualVector a c) - tt = case coerceFmapTensorProduct @w [] + tt = case dualSpaceWitness @c of + DualSpaceWitness -> case coerceFmapTensorProduct @w [] (VSCCoercion @(DualVector c) @(AbstractDualVector a c)) of - VSCCoercion -> coerce (transposeTensor @(DualVector c) @w) + Coercion -> coerce (transposeTensor @(DualVector c) @w) fmapTensor = ft where ft :: ∀ w x . ( TensorSpace w, Scalar w ~ Scalar a , TensorSpace x, Scalar x ~ Scalar a ) @@ -791,7 +989,25 @@ (AbstractDualVector a c ⊗ u) ft = scalarsSameInAbstractionAndDuals @a (coerce $ fzipTensorWith @(DualVector c) @u @w @x) - coerceFmapTensorProduct _ = coerceFmapTensorProduct ([]::[DualVector c]) + tensorUnsafeFromArrayWithOffset :: ∀ w m α + . ( m`Dimensional`w, TensorSpace w, Scalar w ~ Scalar a + , GArr.Vector α (Scalar a) ) + => Int -> α (Scalar a) -> Tensor (Scalar a) (AbstractDualVector a c) w + tensorUnsafeFromArrayWithOffset i + = case dimensionalityWitness @(DualVector c) of + IsStaticDimensional -> scalarsSameInAbstractionAndDuals @a + (coerce . tensorUnsafeFromArrayWithOffset @(DualVector c) @w i) + tensorUnsafeWriteArrayWithOffset :: ∀ w m α σ + . ( m`Dimensional`w, TensorSpace w, Scalar w ~ Scalar a + , GArr.Vector α (Scalar a) ) + => GArr.Mutable α σ (Scalar a) -> Int + -> Tensor (Scalar a) (AbstractDualVector a c) w -> ST σ () + tensorUnsafeWriteArrayWithOffset ar + = case dimensionalityWitness @(DualVector c) of + IsStaticDimensional -> scalarsSameInAbstractionAndDuals @a + (coerce (tensorUnsafeWriteArrayWithOffset @(DualVector c) @w ar)) + coerceFmapTensorProduct _ = scalarsSameInAbstractionAndDuals @a + (coerceFmapTensorProduct @(DualVector c) []) witnessAbstractDualVectorTensorSpacyness :: ∀ a c r . ( AbstractLinearSpace a, VectorSpaceImplementation a ~ c @@ -824,8 +1040,8 @@ -> witnessAbstractDualVectorTensorSpacyness @a ( let LinearMap ida = linearId :: (DualVector c ⊗ w) +> (DualVector c ⊗ w) in LinearMap $ - symVSC (abstractTensorProductsCoercion @a - @(DualVector w ⊗ (AbstractDualVector a c⊗w)) ) + sym (abstractTensorProductsCoercion @a + @(DualVector w ⊗ (AbstractDualVector a c⊗w)) ) . coerceFmapTensorProduct ([]::[c ⊗ DualVector w]) (VSCCoercion @(DualVector c ⊗ w) @(AbstractDualVector a c ⊗ w)) $ ida ) @@ -844,7 +1060,7 @@ LinearFunction $ \f -> (applyTensorFunctional @(DualVector c) -+$> abstractLinmapCoercion @a $ f) - . LinearFunction (abstractTensorsCoercion @a $) + . LinearFunction (abstractDualTensorsCoercion @a $) ) applyTensorLinMap = atlm where atlm :: ∀ u w . ( LinearSpace u, Scalar u ~ Scalar a @@ -859,9 +1075,10 @@ @a @((Tensor (Scalar a) (DualVector u) w)) $ coerce f) :: (DualVector c⊗u)+>w in (applyTensorLinMap @(DualVector c)-+$>f') - . LinearFunction (abstractTensorsCoercion @a $) + . LinearFunction (abstractDualTensorsCoercion @a $) ) useTupleLinearSpaceComponents = \_ -> usingNonTupleTypeAsTupleError + coerceDoubleDual = VSCCoercion instance ∀ a c . ( AbstractLinearSpace a, VectorSpaceImplementation a ~ c , FiniteDimensional a, FiniteDimensional c @@ -883,14 +1100,14 @@ => (AbstractDualVector a c +> w) -> (SubBasis (AbstractDualVector a c), DList w) dclm = case (dualFinitenessWitness @c, abstractTensorProductsCoercion @a @w) of - (DualFinitenessWitness DualSpaceWitness, VSCCoercion) + (DualFinitenessWitness DualSpaceWitness, Coercion) -> coerce (decomposeLinMap @(DualVector c) @w) decomposeLinMapWithin = scalarsSameInAbstraction @a dclm where dclm :: ∀ w . (LSpace w, Scalar w ~ Scalar c) => SubBasis (AbstractDualVector a c) -> (AbstractDualVector a c +> w) -> Either (SubBasis (AbstractDualVector a c), DList w) (DList w) dclm = case (dualFinitenessWitness @c, abstractTensorProductsCoercion @a @w) of - (DualFinitenessWitness DualSpaceWitness, VSCCoercion) + (DualFinitenessWitness DualSpaceWitness, Coercion) -> coerce (decomposeLinMapWithin @(DualVector c) @w) recomposeSB = case dualFinitenessWitness @c of DualFinitenessWitness DualSpaceWitness -> scalarsSameInAbstraction @a @@ -907,13 +1124,13 @@ => SubBasis (AbstractDualVector a c) -> [w] -> (AbstractDualVector a c +> w, [w]) rlm = case (dualFinitenessWitness @c, abstractTensorProductsCoercion @a @w) of - (DualFinitenessWitness DualSpaceWitness, VSCCoercion) + (DualFinitenessWitness DualSpaceWitness, Coercion) -> coerce (recomposeLinMap @(DualVector c) @w) recomposeContraLinMap = scalarsSameInAbstraction @a rclm where rclm :: ∀ f w . (LinearSpace w, Scalar w ~ Scalar c, Hask.Functor f) => (f (Scalar w) -> w) -> f a -> AbstractDualVector a c +> w rclm = case (dualFinitenessWitness @c, abstractTensorProductsCoercion @a @w) of - (DualFinitenessWitness DualSpaceWitness, VSCCoercion) -> \f -> + (DualFinitenessWitness DualSpaceWitness, Coercion) -> \f -> (coerce $ recomposeContraLinMap @(DualVector c) @w @f) f . fmap (coerce :: a -> c) recomposeContraLinMapTensor = scalarsSameInAbstraction @a rclmt @@ -928,7 +1145,7 @@ , abstractTensorProductsCoercion @a @(DualVector u) , abstractTensorProductsCoercion @a @(Tensor (Scalar a) (DualVector u) w) ) of - (DualFinitenessWitness DualSpaceWitness, VSCCoercion, VSCCoercion) -> \f -> + (DualFinitenessWitness DualSpaceWitness, Coercion, Coercion) -> \f -> (coerce $ recomposeContraLinMapTensor @(DualVector c) @u @w @f) f . fmap (coerce :: (AbstractDualVector a c+>DualVector u) -> (DualVector c+>DualVector u)) @@ -959,14 +1176,14 @@ tdbc = case (dualSpaceWitness @c, dualSpaceWitness @w) of (DualSpaceWitness, DualSpaceWitness) -> case abstractTensorProductsCoercion @a @(DualVector w) of - VSCCoercion -> coerce (tensorDualBasisCandidates @(DualVector c) @w) - symTensorDualBasisCandidates = scalarsSameInAbstraction @a + Coercion -> coerce (tensorDualBasisCandidates @(DualVector c) @w) + symTensorDualBasisCandidates = case dualSpaceWitness @c of + DualSpaceWitness -> scalarsSameInAbstraction @a ( case ( coerceFmapTensorProduct @c [] (VSCCoercion @a @c) . abstractTensorProductsCoercion @a @a , coerceFmapTensorProduct @(DualVector c) [] - (VSCCoercion @(AbstractDualVector a c) @(DualVector c)) - , dualSpaceWitness @c ) of - (VSCCoercion, VSCCoercion, DualSpaceWitness) + (VSCCoercion @(AbstractDualVector a c) @(DualVector c)) ) of + (Coercion, Coercion) -> coerce (symTensorDualBasisCandidates @(DualVector c)) ) @@ -1010,6 +1227,16 @@ abstractVS_innerProd = scalarsSameInAbstraction @v ( coerce ((<.>) @(VectorSpaceImplementation v)) ) +abstractVS_dimensionalityWitness + :: ∀ v . ( AbstractTensorSpace v + , DimensionAware (VectorSpaceImplementation v) ) + => DimensionalityWitness v +abstractVS_dimensionalityWitness + = case dimensionalityWitness @(VectorSpaceImplementation v) of + IsStaticDimensional + -> sameDimensionalInAbstraction @v IsStaticDimensional + IsFlexibleDimensional -> IsFlexibleDimensional + abstractVS_scalarsSameInAbstraction :: ∀ v ρ . ( AbstractVectorSpace v , Scalar (VectorSpaceImplementation v) ~ Scalar v @@ -1112,16 +1339,17 @@ abstractVS_wellDefinedTensor = scalarsSameInAbstraction @v (case abstractTensorProductsCoercion @v @w of - VSCCoercion -> coerce (wellDefinedTensor @(VectorSpaceImplementation v) @w)) + Coercion -> coerce (wellDefinedTensor @(VectorSpaceImplementation v) @w)) abstractVS_tensorProduct :: ∀ v w . ( AbstractTensorSpace v , TensorSpace w, Scalar w ~ Scalar v ) => Bilinear v w (v ⊗ w) abstractVS_tensorProduct = scalarsSameInAbstraction @v ( case ( abstractTensorProductsCoercion @v @w ) of - VSCCoercion -> coerce (tensorProduct @(VectorSpaceImplementation v) @w) ) + Coercion -> coerce (tensorProduct @(VectorSpaceImplementation v) @w) ) abstractVS_transposeTensor :: ∀ v w . ( AbstractTensorSpace v + , TensorSpace v , TensorSpace w, Scalar w ~ Scalar v ) => (v ⊗ w) -+> (w ⊗ v) abstractVS_transposeTensor @@ -1129,7 +1357,7 @@ ( abstractTensorProductsCoercion @v @w , coerceFmapTensorProduct @w [] (VSCCoercion @(VectorSpaceImplementation v) @(v)) ) of - (VSCCoercion, VSCCoercion) -> scalarsSameInAbstraction @v + (Coercion, Coercion) -> scalarsSameInAbstraction @v (coerce (transposeTensor @(VectorSpaceImplementation v) @w)) ) @@ -1141,7 +1369,7 @@ = scalarsSameInAbstraction @v ( case ( abstractTensorProductsCoercion @v @u , abstractTensorProductsCoercion @v @w ) of - (VSCCoercion, VSCCoercion) + (Coercion, Coercion) -> coerce (fmapTensor @(VectorSpaceImplementation v) @u @w) ) abstractVS_fzipTensorsWith :: ∀ v u w x . ( AbstractTensorSpace v @@ -1153,18 +1381,58 @@ ( case ( abstractTensorProductsCoercion @v @u , abstractTensorProductsCoercion @v @w , abstractTensorProductsCoercion @v @x ) of - (VSCCoercion, VSCCoercion, VSCCoercion) + (Coercion, Coercion, Coercion) -> coerce (fzipTensorWith @(VectorSpaceImplementation v) @u @w @x) ) -abstractVS_coerceFmapTensorProduct :: ∀ v u w p . - ( AbstractTensorSpace v - ) => p v -> VSCCoercion u w -> VSCCoercion (TensorProduct v u) (TensorProduct v w) +abstractVS_unsafeFromArrayWithOffset :: ∀ v n α + . ( AbstractTensorSpace v, n`Dimensional`VectorSpaceImplementation v + , GArr.Vector α (Scalar v) ) + => Int -> α (Scalar v) -> v +abstractVS_unsafeFromArrayWithOffset = scalarsSameInAbstraction @v + (\i -> coerce + . unsafeFromArrayWithOffset @n @(VectorSpaceImplementation v) i ) + +abstractVS_unsafeWriteArrayWithOffset :: ∀ v n α σ + . ( AbstractTensorSpace v, n`Dimensional`VectorSpaceImplementation v + , GArr.Vector α (Scalar v) ) + => GArr.Mutable α σ (Scalar v) -> Int -> v -> ST σ () +abstractVS_unsafeWriteArrayWithOffset = scalarsSameInAbstraction @v + (\ar -> coerce (unsafeWriteArrayWithOffset @n @(VectorSpaceImplementation v) ar)) + +abstractVS_tensorUnsafeFromArrayWithOffset :: ∀ v w n m α + . ( AbstractTensorSpace v, TensorSpace v + , n`Dimensional`VectorSpaceImplementation v + , TensorSpace w, m`Dimensional`w, Scalar w ~ Scalar v + , GArr.Vector α (Scalar v) ) + => Int -> α (Scalar v) -> (v⊗w) +abstractVS_tensorUnsafeFromArrayWithOffset = scalarsSameInAbstraction @v + (\i -> arr (symVSC abstractTensorsCoercion) + . tensorUnsafeFromArrayWithOffset @(VectorSpaceImplementation v) i ) + +abstractVS_tensorUnsafeWriteArrayWithOffset :: ∀ v w n m α σ + . ( AbstractTensorSpace v, TensorSpace v + , n`Dimensional`VectorSpaceImplementation v + , TensorSpace w, m`Dimensional`w, Scalar w ~ Scalar v + , GArr.Vector α (Scalar v) ) + => GArr.Mutable α σ (Scalar v) -> Int -> (v⊗w) -> ST σ () +abstractVS_tensorUnsafeWriteArrayWithOffset = scalarsSameInAbstraction @v + (\ar i -> + tensorUnsafeWriteArrayWithOffset @(VectorSpaceImplementation v) ar i + . arr abstractTensorsCoercion ) + +abstractVS_coerceFmapTensorProduct :: ∀ v u w s p . + ( AbstractTensorSpace v, TensorSpace u, TensorSpace w + , Scalar (VectorSpaceImplementation v) ~ s + , Scalar u ~ s + , Scalar w ~ s + ) => p v -> VSCCoercion s u w + -> Coercion (TensorProduct v u) (TensorProduct v w) abstractVS_coerceFmapTensorProduct _ crc = case ( abstractTensorProductsCoercion @v @u , abstractTensorProductsCoercion @v @w , coerceFmapTensorProduct @(VectorSpaceImplementation v) [] crc ) of - (VSCCoercion, VSCCoercion, VSCCoercion) -> VSCCoercion + (Coercion, Coercion, Coercion) -> Coercion abstractVS_dualSpaceWitness :: ∀ v . (AbstractLinearSpace v , LinearSpace v @@ -1176,16 +1444,23 @@ DualSpaceWitness -> DualSpaceWitness ) -abstractVS_linearId :: ∀ v . ( AbstractLinearSpace v +abstractVS_linearId :: ∀ v . ( AbstractLinearSpace v, TensorSpace v , LinearSpace (VectorSpaceImplementation v) ) => v +> v abstractVS_linearId = case dualSpaceWitness @(VectorSpaceImplementation v) of - DualSpaceWitness -> case coerceFmapTensorProduct + DualSpaceWitness -> scalarsSameInAbstraction @v ( case coerceFmapTensorProduct @(DualVector (VectorSpaceImplementation v)) [] (VSCCoercion @v @(VectorSpaceImplementation v)) of - VSCCoercion -> coerce (linearId @(VectorSpaceImplementation v)) + Coercion -> coerce (linearId @(VectorSpaceImplementation v)) + ) +abstractTensorsCoercion :: ∀ v w s . AbstractTensorSpace v + => VSCCoercion s (v⊗w) (VectorSpaceImplementation v⊗w) +abstractTensorsCoercion = case abstractTensorProductsCoercion @v @w of + Coercion -> VSCCoercion + abstractVS_tensorId :: ∀ v w . ( AbstractLinearSpace v + , TensorSpace v , LinearSpace (VectorSpaceImplementation v) , LinearSpace w, Scalar w ~ Scalar v ) => (v ⊗ w) +> (v ⊗ w) @@ -1193,21 +1468,18 @@ (case (dualSpaceWitness @w, dualSpaceWitness @(VectorSpaceImplementation v)) of (DualSpaceWitness, DualSpaceWitness) -> case coerceFmapTensorProduct @(DualVector w) [] - $ VSCCoercion @(TensorProduct (VectorSpaceImplementation v) w) - @(VectorSpaceImplementation v ⊗ w) - . abstractTensorProductsCoercion @v @w - . VSCCoercion @(v ⊗ w) @(TensorProduct v w) of - VSCCoercion + $ abstractTensorsCoercion @v @w of + Coercion -> case ( coerceFmapTensorProduct @(DualVector (VectorSpaceImplementation v)) [] - (VSCCoercion :: VSCCoercion + (VSCCoercion :: VSCCoercion (Scalar v) (Tensor (Scalar v) (DualVector w) (Tensor (Scalar v) v w)) (Tensor (Scalar v) (DualVector w) (Tensor (Scalar v) (VectorSpaceImplementation v) w))) ) of - VSCCoercion + Coercion -> coerce (tensorId @(VectorSpaceImplementation v) @w) ) @@ -1233,7 +1505,7 @@ => Bilinear (DualVector (v⊗u)) (v⊗u) (Scalar v) abstractVS_applyTensorFunctional = scalarsSameInAbstraction @v (case abstractTensorProductsCoercion @v @u of - VSCCoercion -> coerce (applyTensorFunctional @(VectorSpaceImplementation v) @u)) + Coercion -> coerce (applyTensorFunctional @(VectorSpaceImplementation v) @u)) abstractVS_applyTensorLinMap :: ∀ v u w . ( AbstractLinearSpace v @@ -1243,12 +1515,12 @@ => Bilinear ((v⊗u)+>w) (v⊗u) w abstractVS_applyTensorLinMap = scalarsSameInAbstraction @v ( case abstractTensorProductsCoercion @v @u of - VSCCoercion -> coerce (applyTensorLinMap @(VectorSpaceImplementation v) @u @w) ) + Coercion -> coerce (applyTensorLinMap @(VectorSpaceImplementation v) @u @w) ) abstractSubbasisCoercion :: ∀ v . Coercible (SubBasis v) (SubBasis (VectorSpaceImplementation v)) - => VSCCoercion (SubBasis v) (SubBasis (VectorSpaceImplementation v)) -abstractSubbasisCoercion = VSCCoercion + => Coercion (SubBasis v) (SubBasis (VectorSpaceImplementation v)) +abstractSubbasisCoercion = Coercion precomposeCoercion :: Coercion a b -> Coercion (b -> c) (a -> c) precomposeCoercion Coercion = Coercion @@ -1276,7 +1548,7 @@ ( AbstractLinearSpace v, FiniteDimensional (VectorSpaceImplementation v) , Coercible (SubBasis v) (SubBasis (VectorSpaceImplementation v)) ) => SubBasis v -abstractVS_entireBasis = symVSC (abstractSubbasisCoercion @v) +abstractVS_entireBasis = sym (abstractSubbasisCoercion @v) $ entireBasis @(VectorSpaceImplementation v) abstractVS_enumerateSubBasis :: ∀ v . @@ -1284,7 +1556,7 @@ , Coercible (SubBasis v) (SubBasis (VectorSpaceImplementation v)) ) => SubBasis v -> [v] abstractVS_enumerateSubBasis = precomposeCoercion - (getVSCCoercion $ abstractSubbasisCoercion @v) + (abstractSubbasisCoercion @v) $ coerce (enumerateSubBasis @(VectorSpaceImplementation v)) abstractVS_decomposeLinMap :: ∀ v w . @@ -1295,9 +1567,9 @@ => (v +> w) -> (SubBasis v, DList w) abstractVS_decomposeLinMap = scalarsSameInAbstraction @v ( postcomposeCoercion (firstCoercion $ sym - (getVSCCoercion $ abstractSubbasisCoercion @v)) + (abstractSubbasisCoercion @v)) $ case abstractTensorProductsCoercion @v @w of - VSCCoercion -> ( coerce (decomposeLinMap @(VectorSpaceImplementation v) @w) + Coercion -> ( coerce (decomposeLinMap @(VectorSpaceImplementation v) @w) :: (v +> w) -> ( SubBasis (VectorSpaceImplementation v) , DList w ) ) ) @@ -1308,10 +1580,10 @@ , LSpace w, Scalar w ~ Scalar v ) => SubBasis v -> (v +> w) -> Either (SubBasis v, DList w) (DList w) abstractVS_decomposeLinMapWithin = scalarsSameInAbstraction @v - ( precomposeCoercion (getVSCCoercion $ abstractSubbasisCoercion @v) + ( precomposeCoercion (abstractSubbasisCoercion @v) . postcomposeCoercion (postcomposeCoercion . leftCoercion . firstCoercion $ sym - (getVSCCoercion $ abstractSubbasisCoercion @v)) + (abstractSubbasisCoercion @v)) $ coerce (decomposeLinMapWithin @(VectorSpaceImplementation v) @w) ) @@ -1320,7 +1592,7 @@ , Coercible (SubBasis v) (SubBasis (VectorSpaceImplementation v)) ) => SubBasis v -> [Scalar v] -> (v, [Scalar v]) abstractVS_recomposeSB = scalarsSameInAbstraction @v - ( precomposeCoercion (getVSCCoercion $ abstractSubbasisCoercion @v) + ( precomposeCoercion (abstractSubbasisCoercion @v) $ coerce (recomposeSB @(VectorSpaceImplementation v)) ) @@ -1330,9 +1602,9 @@ , FiniteDimensional w, Scalar w ~ Scalar v ) => SubBasis v -> SubBasis w -> [Scalar v] -> (v ⊗ w, [Scalar v]) abstractVS_recomposeSBTensor = scalarsSameInAbstraction @v - ( precomposeCoercion (getVSCCoercion $ abstractSubbasisCoercion @v) + ( precomposeCoercion (abstractSubbasisCoercion @v) $ case abstractTensorProductsCoercion @v @w of - VSCCoercion -> coerce (recomposeSBTensor @(VectorSpaceImplementation v) @w) + Coercion -> coerce (recomposeSBTensor @(VectorSpaceImplementation v) @w) ) abstractVS_recomposeLinMap :: ∀ v w . ( AbstractLinearSpace v @@ -1341,7 +1613,7 @@ , LSpace w, Scalar w ~ Scalar v ) => SubBasis v -> [w] -> (v +> w, [w]) abstractVS_recomposeLinMap = scalarsSameInAbstraction @v - ( precomposeCoercion (getVSCCoercion $ abstractSubbasisCoercion @v) + ( precomposeCoercion (abstractSubbasisCoercion @v) $ coerce (recomposeLinMap @(VectorSpaceImplementation v) @w) ) @@ -1391,7 +1663,7 @@ => (v ⊗ w) -> (v ⊗ w) -> Bool abstractVS_tensorEquality = scalarsSameInAbstraction @v ( case abstractTensorProductsCoercion @v @w of - VSCCoercion -> coerce (tensorEquality @(VectorSpaceImplementation v) @w) + Coercion -> coerce (tensorEquality @(VectorSpaceImplementation v) @w) ) abstractVS_dualBasisCandidates :: ∀ v . ( AbstractLinearSpace v @@ -1415,7 +1687,7 @@ , abstractTensorProductsCoercion @v @(DualVector w) , abstractTensorProductsCoercion @v @w ) of - (VSCCoercion, VSCCoercion, VSCCoercion) + (VSCCoercion, Coercion, Coercion) -> coerce (tensorDualBasisCandidates @(VectorSpaceImplementation v) @w) ) @@ -1437,7 +1709,7 @@ , coerceFmapTensorProduct @(VectorSpaceImplementation v) [] (VSCCoercion @v @(VectorSpaceImplementation v)) ) of - (VSCCoercion, VSCCoercion, VSCCoercion, VSCCoercion) + (Coercion, Coercion, Coercion, Coercion) -> coerce (symTensorDualBasisCandidates @(VectorSpaceImplementation v)) ) @@ -1447,6 +1719,10 @@ copyNewtypeInstances :: Q Type -> [Name] -> DecsQ copyNewtypeInstances cxtv classes = do + requireExtensions [ TypeFamilies, MultiParamTypeClasses + , ScopedTypeVariables, TypeApplications + , DataKinds ] + (tvbs, cxt, (a,c)) <- do (tvbs', cxt', a') <- deQuantifyType cxtv let extractImplementationType (AppT tc (VarT tvb)) atvbs @@ -1561,6 +1837,19 @@ $(varP '(.-~.)) = \p q -> Just (abstractVS_subvs p q) $(varP '(.-~!)) = abstractVS_subvs |] + "DimensionAware" -> InstanceD Nothing <$> cxt <*> + [t|DimensionAware $a|] <*> [d| + type instance StaticDimension $a = StaticDimension $c + $(varP 'dimensionalityWitness) = abstractVS_dimensionalityWitness + |] + "Dimensional" -> do + dim <- pure . VarT <$> newName "n" + InstanceD Nothing <$> ((:)<$>[t|StaticDimension $c ~ 'Just $dim|] + <*>cxt) <*> + [t|Dimensional $dim $a|] <*> [d| + $(varP 'unsafeFromArrayWithOffset) = abstractVS_unsafeFromArrayWithOffset + $(varP 'unsafeWriteArrayWithOffset) = abstractVS_unsafeWriteArrayWithOffset + |] "TensorSpace" -> InstanceD Nothing <$> cxt <*> [t|TensorSpace $a|] <*> [d| type instance TensorProduct $a w = TensorProduct $c w @@ -1579,12 +1868,18 @@ $(varP 'transposeTensor) = abstractVS_transposeTensor $(varP 'fmapTensor) = abstractVS_fmapTensor $(varP 'fzipTensorWith) = abstractVS_fzipTensorsWith + $(varP 'tensorUnsafeFromArrayWithOffset) + = abstractVS_tensorUnsafeFromArrayWithOffset + $(varP 'tensorUnsafeWriteArrayWithOffset) + = abstractVS_tensorUnsafeWriteArrayWithOffset $(varP 'coerceFmapTensorProduct) = abstractVS_coerceFmapTensorProduct |] "AbstractTensorSpace" -> InstanceD Nothing <$> cxt <*> [t|AbstractTensorSpace $a|] <*> [d| + $(varP 'sameDimensionalInAbstraction) + = \φ -> φ $(varP 'abstractTensorProductsCoercion) - = VSCCoercion + = Coercion |] "LinearSpace" -> InstanceD Nothing <$> cxt <*> [t|LinearSpace $a|] <*> [d| diff --git a/Math/LinearMap/Coercion.hs b/Math/LinearMap/Coercion.hs --- a/Math/LinearMap/Coercion.hs +++ b/Math/LinearMap/Coercion.hs @@ -10,6 +10,9 @@ module Math.LinearMap.Coercion ( VSCCoercion(..) + , (-+$=>) + -- * Conversion between the internal types + , fromLinearMap, asLinearMap, fromTensor, asTensor ) where import Math.LinearMap.Category.Class diff --git a/Math/VectorSpace/DimensionAware.hs b/Math/VectorSpace/DimensionAware.hs new file mode 100644 --- /dev/null +++ b/Math/VectorSpace/DimensionAware.hs @@ -0,0 +1,269 @@ +-- | +-- Module : Math.VectorSpace.DimensionAware +-- Copyright : (c) Justus Sagemüller 2022 +-- License : GPL v3 +-- +-- Maintainer : (@) jsag $ hvl.no +-- Stability : experimental +-- Portability : portable +-- + + +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE DefaultSignatures #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE UnicodeSyntax #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE CPP #-} + +module Math.VectorSpace.DimensionAware where + +import Data.VectorSpace + +import Data.Singletons (SingI, sing, Sing) +#if MIN_VERSION_singletons(3,0,0) +import Prelude.Singletons (SNum(..)) +import Data.Maybe.Singletons +import GHC.TypeLits.Singletons (withKnownNat) +#else +import Data.Singletons.Prelude.Num (SNum(..)) +import Data.Singletons.Prelude.Maybe (SMaybe(..)) +import Data.Singletons.TypeLits (withKnownNat) +#endif + +import qualified Data.Vector.Generic as GArr +import qualified Data.Vector.Generic.Mutable as GMArr +import Control.Monad.ST (ST) + +import Control.Monad + +import GHC.TypeLits +import GHC.Exts (Constraint) +import Data.Proxy (Proxy(..)) + +import Data.Ratio + +import qualified Math.VectorSpace.DimensionAware.Theorems.MaybeNat as Maybe + + +-- | Low-level case distinction between spaces with a dimension that is both fixed +-- and low enough that it makes sense to treat it this way, and more general +-- spaces where this is not feasible. +-- +-- Use this type only when defining instances of 'DimensionAware'. When making +-- decisions based on dimensionality, 'DimensionalityCases' is more convenient. +data DimensionalityWitness v where + IsStaticDimensional :: (n`Dimensional`v) => DimensionalityWitness v + IsFlexibleDimensional :: StaticDimension v ~ 'Nothing => DimensionalityWitness v + + +-- | This class does not really pose any restrictions on a vector space type, but +-- allows it to express its dimension. +-- This is for optimisation purposes only, specifically to allow low-dimensional vectors +-- to be represented efficiently in unboxed arrays / matrices. +class VectorSpace v => DimensionAware v where + -- | If this is `Nothing`, + -- it can mean the dimension is infinite, or just big, or simply unknown / not + -- considered in the implementation. + type StaticDimension v :: Maybe Nat + type StaticDimension v = 'Nothing + + dimensionalityWitness :: DimensionalityWitness v + + +instance DimensionAware Float where + type StaticDimension Float = 'Just 1 + dimensionalityWitness = IsStaticDimensional +instance DimensionAware Double where + type StaticDimension Double = 'Just 1 + dimensionalityWitness = IsStaticDimensional +instance DimensionAware Int where + type StaticDimension Int = 'Just 1 + dimensionalityWitness = IsStaticDimensional +instance DimensionAware Integer where + type StaticDimension Integer = 'Just 1 + dimensionalityWitness = IsStaticDimensional +instance Integral n => DimensionAware (Ratio n) where + type StaticDimension (Ratio n) = 'Just 1 + dimensionalityWitness = IsStaticDimensional + +instance ∀ u v . (DimensionAware u, DimensionAware v, Scalar u ~ Scalar v) + => DimensionAware (u,v) where + type StaticDimension (u,v) = Maybe.ZipWithPlus (StaticDimension u) (StaticDimension v) + dimensionalityWitness = case (dimensionalityWitness @u, dimensionalityWitness @v) of + (IsStaticDimensional, IsStaticDimensional) + -> withKnownNat (dimensionalitySing @u %+ dimensionalitySing @v) + IsStaticDimensional + (IsFlexibleDimensional, _) -> IsFlexibleDimensional + (_, IsFlexibleDimensional) -> IsFlexibleDimensional + + +class (DimensionAware v, StaticDimension v ~ 'Just n) + => n`Dimensional`v | v -> n where + knownDimensionalitySing :: Sing n + {-# INLINE knownDimensionalitySing #-} + default knownDimensionalitySing :: KnownNat n => Sing n + knownDimensionalitySing = sing + -- | Read basis expansion from an array, starting at the specified offset. + -- The array must have at least length @n + offset@, else the behaviour is undefined. + unsafeFromArrayWithOffset :: GArr.Vector α (Scalar v) => Int -> α (Scalar v) -> v + unsafeWriteArrayWithOffset :: GArr.Vector α (Scalar v) + => GArr.Mutable α σ (Scalar v) -> Int -> v -> ST σ () + +-- | Batteries-included version of 'DimensionalityWitness'. +data DimensionalityCases v where + StaticDimensionalCase :: (KnownNat n, n`Dimensional`v) => DimensionalityCases v + FlexibleDimensionalCase :: StaticDimension v ~ 'Nothing => DimensionalityCases v + +#if !MIN_VERSION_singletons(3,0,0) +type family FromJust (a :: Maybe k) :: k where + FromJust ('Just v) = v +#endif + +type Dimension v = FromJust (StaticDimension v) + +#if !MIN_VERSION_singletons(3,0,0) +type family IsJust (a :: Maybe k) :: Bool where + IsJust ('Just _) = 'True + IsJust _ = 'False +#endif + +class DimensionAware v => StaticDimensional v where + dimensionIsStatic :: ∀ r . (∀ n . (KnownNat n, n`Dimensional`v) => r) -> r + +{-# INLINE dimensionalitySing #-} +dimensionalitySing :: ∀ v n . n`Dimensional`v => Sing n +dimensionalitySing = knownDimensionalitySing @n @v + +instance ( DimensionAware v, IsJust (StaticDimension v) ~ 'True ) + => StaticDimensional v where + dimensionIsStatic = case dimensionalityWitness @v of + IsStaticDimensional -> \φ -> withKnownNat (dimensionalitySing @v) φ + +dimensionality :: ∀ v . DimensionAware v => DimensionalityCases v +dimensionality = case dimensionalityWitness @v of + IsStaticDimensional -> withKnownNat (dimensionalitySing @v) StaticDimensionalCase + IsFlexibleDimensional -> FlexibleDimensionalCase + +{-# INLINE dimension #-} +dimension :: ∀ v n a . (n`Dimensional`v, Integral a) => a +dimension = withKnownNat (dimensionalitySing @v) (fromIntegral $ natVal @n Proxy) + +-- | Convenience function. The result does never depend on the runtime input, only +-- on its type. +dimensionOf :: ∀ v n a . (n`Dimensional`v, Integral a) => v -> a +dimensionOf _ = dimension @v + +{-# INLINE unsafeFromArray #-} +-- | Read basis expansion from an array. The array must have length @n@, else the +-- behaviour is undefined. +unsafeFromArray :: ∀ v n α . (n`Dimensional`v, GArr.Vector α (Scalar v)) + => α (Scalar v) -> v +unsafeFromArray = unsafeFromArrayWithOffset 0 + +-- | Read basis expansion from an array, if the size equals the dimension. +fromArray :: ∀ v n α . (n`Dimensional`v, GArr.Vector α (Scalar v)) + => α (Scalar v) -> Maybe v +fromArray ar + | GArr.length ar == dimension @v = Just $ unsafeFromArray ar + | otherwise = Nothing + +{-# INLINE toArray #-} +-- | Write out basis expansion to an array, whose length will always be @n@. +toArray :: ∀ v n α . (n`Dimensional`v, GArr.Vector α (Scalar v)) + => v -> α (Scalar v) +toArray v = GArr.create (do + ar <- GMArr.new $ dimension @v + unsafeWriteArrayWithOffset ar 0 v + return ar + ) + +{-# INLINE staticDimensionSing #-} +staticDimensionSing :: ∀ v . DimensionAware v => Sing (StaticDimension v) +staticDimensionSing = case dimensionalityWitness @v of + IsStaticDimensional -> SJust (dimensionalitySing @v) + IsFlexibleDimensional -> sing + +{-# INLINE scalarUnsafeFromArrayWithOffset #-} +scalarUnsafeFromArrayWithOffset :: (v ~ Scalar v, GArr.Vector α v) + => Int -> α v -> v +scalarUnsafeFromArrayWithOffset i = (`GArr.unsafeIndex`i) + +{-# INLINE scalarUnsafeWriteArrayWithOffset #-} +scalarUnsafeWriteArrayWithOffset :: (v ~ Scalar v, GArr.Vector α v) + => GArr.Mutable α σ v -> Int -> v -> ST σ () +scalarUnsafeWriteArrayWithOffset ar i = GMArr.unsafeWrite ar i + +{-# INLINE unsafeFromArrayWithOffsetViaList #-} +unsafeFromArrayWithOffsetViaList + :: ∀ v n α . (n`Dimensional`v, GArr.Vector α (Scalar v)) + => ([Scalar v] -> v) -> Int -> α (Scalar v) -> v +unsafeFromArrayWithOffsetViaList l2v i + = l2v . GArr.toList . GArr.unsafeSlice i (dimension @v) + +{-# INLINE unsafeWriteArrayWithOffsetViaList #-} +unsafeWriteArrayWithOffsetViaList + :: ∀ v n α σ . (n`Dimensional`v, GArr.Vector α (Scalar v)) + => (v -> [Scalar v]) -> GArr.Mutable α σ (Scalar v) + -> Int -> v -> ST σ () +unsafeWriteArrayWithOffsetViaList v2l ar i + = GMArr.unsafeCopy (GMArr.unsafeSlice i (dimension @v) ar) + <=< GArr.unsafeThaw @(ST σ) @α . GArr.fromList . v2l + +instance 1`Dimensional`Float where + {-# INLINE unsafeFromArrayWithOffset #-} + unsafeFromArrayWithOffset = scalarUnsafeFromArrayWithOffset + {-# INLINE unsafeWriteArrayWithOffset #-} + unsafeWriteArrayWithOffset = scalarUnsafeWriteArrayWithOffset +instance 1`Dimensional`Double where + {-# INLINE unsafeFromArrayWithOffset #-} + unsafeFromArrayWithOffset = scalarUnsafeFromArrayWithOffset + {-# INLINE unsafeWriteArrayWithOffset #-} + unsafeWriteArrayWithOffset = scalarUnsafeWriteArrayWithOffset +instance 1`Dimensional`Int where + {-# INLINE unsafeFromArrayWithOffset #-} + unsafeFromArrayWithOffset = scalarUnsafeFromArrayWithOffset + {-# INLINE unsafeWriteArrayWithOffset #-} + unsafeWriteArrayWithOffset = scalarUnsafeWriteArrayWithOffset +instance 1`Dimensional`Integer where + {-# INLINE unsafeFromArrayWithOffset #-} + unsafeFromArrayWithOffset = scalarUnsafeFromArrayWithOffset + {-# INLINE unsafeWriteArrayWithOffset #-} + unsafeWriteArrayWithOffset = scalarUnsafeWriteArrayWithOffset +instance Integral n => 1`Dimensional`Ratio n where + {-# INLINE unsafeFromArrayWithOffset #-} + unsafeFromArrayWithOffset = scalarUnsafeFromArrayWithOffset + {-# INLINE unsafeWriteArrayWithOffset #-} + unsafeWriteArrayWithOffset = scalarUnsafeWriteArrayWithOffset + + +instance ∀ n u m v nm . ( n`Dimensional`u, m`Dimensional`v + , Scalar u ~ Scalar v + , nm ~ (n+m) ) + => nm`Dimensional`(u,v) where + {-# INLINE knownDimensionalitySing #-} + knownDimensionalitySing = dimensionalitySing @u %+ dimensionalitySing @v + {-# INLINE unsafeFromArrayWithOffset #-} + unsafeFromArrayWithOffset i arr + = ( unsafeFromArrayWithOffset i arr + , unsafeFromArrayWithOffset (i + dimension @u) arr ) + {-# INLINE unsafeWriteArrayWithOffset #-} + unsafeWriteArrayWithOffset arr i (x,y) = do + unsafeWriteArrayWithOffset arr i x + unsafeWriteArrayWithOffset arr (i + dimension @u) y + +notStaticDimensionalContradiction :: ∀ v n r + . (n`Dimensional`v, StaticDimension v ~ 'Nothing) => r +notStaticDimensionalContradiction = undefined + diff --git a/Math/VectorSpace/DimensionAware/Theorems/MaybeNat.hs b/Math/VectorSpace/DimensionAware/Theorems/MaybeNat.hs new file mode 100644 --- /dev/null +++ b/Math/VectorSpace/DimensionAware/Theorems/MaybeNat.hs @@ -0,0 +1,120 @@ +-- | +-- Module : Math.VectorSpace.DimensionAware.Theorems.MaybeNat +-- Copyright : (c) Justus Sagemüller 2022 +-- License : GPL v3 +-- +-- Maintainer : (@) jsag $ hvl.no +-- Stability : experimental +-- Portability : portable +-- +{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE Rank2Types #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE PatternSynonyms #-} +{-# LANGUAGE UnicodeSyntax #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE NoStarIsType #-} +{-# LANGUAGE CPP #-} + +module Math.VectorSpace.DimensionAware.Theorems.MaybeNat where + +#if MIN_VERSION_singletons(3,0,0) +import Prelude.Singletons (SNum(..)) +import Data.Maybe.Singletons +import GHC.TypeLits.Singletons (SNat(..), withKnownNat, sDiv) +#else +import Data.Singletons.Prelude.Num (SNum(..)) +import Data.Singletons.Prelude.Maybe (SMaybe(..)) +import Data.Singletons.TypeLits (SNat(..), withKnownNat, sDiv) +#endif +import Data.Singletons +import qualified Data.Type.Natural as DTN +import GHC.TypeLits +import Unsafe.Coerce + +type family ZipWith (f :: k -> l -> m) (a :: Maybe k) (b :: Maybe l) :: Maybe m where + ZipWith f 'Nothing y = 'Nothing + ZipWith f x 'Nothing = 'Nothing + ZipWith f ('Just x) ('Just y) = 'Just (f x y) + +type family ZipWithPlus (a :: Maybe Nat) (b :: Maybe Nat) :: Maybe Nat where + ZipWithPlus 'Nothing y = 'Nothing + ZipWithPlus x 'Nothing = 'Nothing + ZipWithPlus ('Just x) ('Just y) = 'Just (x+y) + +type family ZipWithTimes (a :: Maybe Nat) (b :: Maybe Nat) :: Maybe Nat where + ZipWithTimes 'Nothing y = 'Nothing + ZipWithTimes x 'Nothing = 'Nothing + ZipWithTimes ('Just x) ('Just y) = 'Just (x*y) + +type family MaybePred (a :: Nat) :: Maybe Nat where + MaybePred 0 = 'Nothing + MaybePred n = 'Just (n-1) + +type family BindMaybePred (a :: Maybe Nat) :: Maybe Nat where + BindMaybePred 'Nothing = 'Nothing + BindMaybePred ('Just n) = MaybePred n + +type TriangularNum (a :: Nat) = (a * (a+1))`Div`2 + +type family FmapTriangularNum (a :: Maybe Nat) where + FmapTriangularNum 'Nothing = 'Nothing + FmapTriangularNum ('Just n) = ('Just (TriangularNum n)) + +justNatSing :: ∀ (n :: Nat) . Sing n -> Sing ('Just n) +justNatSing SNat = sing + +succMaybePredSing :: ∀ n . DTN.SNat n -> Sing (MaybePred (n+1)) +succMaybePredSing s = unsafeCoerce (DTN.withKnownNat s (justNatSing (SNat @n))) + +maybePredSing :: ∀ a . Sing a -> Sing (MaybePred a) +maybePredSing α = withKnownNat α + (case DTN.viewNat (DTN.sNat @a) of + DTN.IsZero -> sing + DTN.IsSucc β -> succMaybePredSing β + ) + +binMaybePredSing :: ∀ a . Sing a -> Sing (BindMaybePred a) +binMaybePredSing SNothing = sing +binMaybePredSing (SJust ν) = maybePredSing ν + +triangularNumSing :: ∀ a . Sing a -> Sing (TriangularNum a) +triangularNumSing α = (α %* (α%+(sing @1)))`sDiv`(sing @2) + +fmapTriangularNumSing :: ∀ a . Sing a -> Sing (FmapTriangularNum a) +fmapTriangularNumSing SNothing = SNothing +fmapTriangularNumSing (SJust α) = SJust (triangularNumSing α) + +zipWithPlusSing :: ∀ a b r . Sing a -> Sing b -> Sing (ZipWithPlus a b) +zipWithPlusSing SNothing _ = sing +zipWithPlusSing _ SNothing = sing +zipWithPlusSing (SJust α) (SJust β) = withKnownNat (α%+β) sing + +zipWithTimesSing :: ∀ a b r . Sing a -> Sing b -> Sing (ZipWithTimes a b) +zipWithTimesSing SNothing _ = sing +zipWithTimesSing _ SNothing = sing +zipWithTimesSing (SJust α) (SJust β) = withKnownNat (α%*β) sing + +zipWithTimesAssoc :: ∀ a b c r . Sing a -> Sing b -> Sing c + -> ((ZipWithTimes a (ZipWithTimes b c) ~ ZipWithTimes (ZipWithTimes a b) c) => r) + -> r +zipWithTimesAssoc SNothing _ _ φ = φ +zipWithTimesAssoc _ SNothing _ φ = φ +zipWithTimesAssoc _ _ SNothing φ = φ +zipWithTimesAssoc (SJust _) (SJust _) (SJust _) φ = φ + +zipWithTimesCommu :: ∀ a b r . Sing a -> Sing b + -> ((ZipWithTimes a b ~ ZipWithTimes b a) => r) -> r +zipWithTimesCommu SNothing _ φ = φ +zipWithTimesCommu _ SNothing φ = φ +zipWithTimesCommu (SJust _) (SJust _) φ = φ + diff --git a/Math/VectorSpace/Docile.hs b/Math/VectorSpace/Docile.hs --- a/Math/VectorSpace/Docile.hs +++ b/Math/VectorSpace/Docile.hs @@ -25,6 +25,7 @@ {-# LANGUAGE RankNTypes #-} {-# LANGUAGE EmptyCase #-} {-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE InstanceSigs #-} {-# LANGUAGE TypeApplications #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE DefaultSignatures #-} @@ -456,7 +457,12 @@ instance ∀ s u v . ( SemiInner u, SemiInner v, Scalar u ~ s, Scalar v ~ s ) => SemiInner (Tensor s u v) where dualBasisCandidates = tensorDualBasisCandidates - tensorDualBasisCandidates = map (second $ arr rassocTensor) + tensorDualBasisCandidates :: ∀ w + . (SemiInner w, Scalar w ~ s) + => [(Int, Tensor s (Tensor s u v) w)] + -> Forest (Int, LinearMap s (Tensor s u v) (DualVector w)) + tensorDualBasisCandidates = case dualSpaceWitness @w of + DualSpaceWitness -> map (second $ arr rassocTensor) >>> tensorDualBasisCandidates >>> map (fmap . second $ arr uncurryLinearMap) @@ -474,7 +480,12 @@ -> [(Int, Tensor s (DualVector u) v)]) >>> tensorDualBasisCandidates >>> coerce - tensorDualBasisCandidates = map (second $ arr hasteLinearMap) + tensorDualBasisCandidates :: ∀ w + . (SemiInner w, Scalar w ~ s) + => [(Int, Tensor s (LinearMap s u v) w)] + -> Forest (Int, LinearMap s (LinearMap s u v) (DualVector w)) + tensorDualBasisCandidates = case dualSpaceWitness @w of + DualSpaceWitness -> map (second $ arr hasteLinearMap) >>> dualBasisCandidates >>> map (fmap . second $ arr coUncurryLinearMap) @@ -922,8 +933,9 @@ oscld = map (sqrt 0.5*)<$>o in concat (sd₀ []) ++ d ++ concat oscld ++ mkSym nw (n-1) (zipWith (.) sds $ (:)<$>oscld) rest - recomposeContraLinMap f tenss - = LinearMap . arr (rassocTensor . asTensor) . rcCLM dualFinitenessWitness f + recomposeContraLinMap f tenss = case dualSpaceWitness @v of + DualSpaceWitness -> + LinearMap . arr (rassocTensor . asTensor) . rcCLM dualFinitenessWitness f $ fmap getSymmetricTensor tenss where rcCLM :: (Hask.Functor f, LinearSpace w, s~Scalar w) => DualFinitenessWitness v @@ -1181,8 +1193,8 @@ instance ( RieszDecomposable x, RieszDecomposable y , Scalar x ~ Scalar y, Scalar (DualVector x) ~ Scalar (DualVector y) ) => RieszDecomposable (x,y) where - rieszDecomposition m = map (first Left) (rieszDecomposition $ fst . m) - ++ map (first Right) (rieszDecomposition $ snd . m) + rieszDecomposition m = map (first Left) (rieszDecomposition $ fmap fst -+$> m) + ++ map (first Right) (rieszDecomposition $ fmap snd -+$> m) instance RieszDecomposable (ZeroDim ℝ) where rieszDecomposition _ = [] @@ -1296,11 +1308,26 @@ b .⊗ w = basisValue b ⊗ w class (FiniteDimensional v, HasBasis v) => TensorDecomposable v where - tensorDecomposition :: v⊗w -> [(Basis v, w)] + tensorDecomposition :: (TensorSpace w, Scalar w ~ Scalar v) + => v⊗w -> [(Basis v, w)] + tensorDecompose' :: (TensorSpace w, Scalar w ~ Scalar v) + => v⊗w -> Basis v -> w showsPrecBasis :: Int -> Basis v -> ShowS +instance ( TensorDecomposable u, TensorSpace v + , HasBasis u, HasBasis v + , Num' s, Scalar u ~ s, Scalar v ~ s + ) => HasBasis (Tensor s u v) where + type Basis (Tensor s u v) = (Basis u, Basis v) + basisValue (bu, bv) = basisValue bu ⊗ basisValue bv + decompose t = [ ((bu,bv),s) + | (bu,v) <- tensorDecomposition t + , (bv,s) <- decompose v ] + decompose' t (bu, bv) = decompose' (tensorDecompose' t bu) bv + instance TensorDecomposable ℝ where tensorDecomposition (Tensor r) = [((), r)] + tensorDecompose' (Tensor r) () = r showsPrecBasis _ = shows instance ∀ x y . ( TensorDecomposable x, TensorDecomposable y , Scalar x ~ Scalar y, Scalar (DualVector x) ~ Scalar (DualVector y) ) @@ -1308,6 +1335,10 @@ tensorDecomposition (Tensor (tx,ty)) = map (first Left) (tensorDecomposition tx) ++ map (first Right) (tensorDecomposition ty) + tensorDecompose' (Tensor (tx,ty)) (Left bx) + = tensorDecompose' tx bx + tensorDecompose' (Tensor (tx,ty)) (Right by) + = tensorDecompose' ty by showsPrecBasis p (Left bx) = showParen (p>9) $ ("Left "++) . showsPrecBasis @x 10 bx showsPrecBasis p (Right by) @@ -1315,9 +1346,11 @@ instance TensorDecomposable (ZeroDim ℝ) where tensorDecomposition _ = [] + tensorDecompose' _ = absurd showsPrecBasis _ = absurd instance TensorDecomposable (V0 ℝ) where tensorDecomposition _ = [] + tensorDecompose' _ b = case b of {} #if MIN_VERSION_free_vector_spaces(0,2,0) showsPrecBasis = showsPrec #else @@ -1326,36 +1359,64 @@ instance TensorDecomposable (V1 ℝ) where #if MIN_VERSION_free_vector_spaces(0,2,0) tensorDecomposition (Tensor (V1 w)) = [(e @0, w)] + tensorDecompose' (Tensor (V1 w)) _ = w showsPrecBasis = showsPrec #else tensorDecomposition (Tensor (V1 w)) = [(ex, w)] + tensorDecompose' (Tensor w) (Mat.E q) = w^.q showsPrecBasis _ (Mat.E q) = (V1"ex"^.q ++) #endif instance TensorDecomposable (V2 ℝ) where #if MIN_VERSION_free_vector_spaces(0,2,0) tensorDecomposition (Tensor (V2 x y)) = [ (e @0, x), (e @1, y) ] + tensorDecompose' (Tensor (V2 x y)) b = case getEuclideanBasisIndex b of + { 0 -> x; 1 -> y } showsPrecBasis = showsPrec #else tensorDecomposition (Tensor (V2 x y)) = [ (ex, x), (ey, y) ] + tensorDecompose' (Tensor w) (Mat.E q) = w^.q showsPrecBasis _ (Mat.E q) = (V2"ex""ey"^.q ++) #endif instance TensorDecomposable (V3 ℝ) where #if MIN_VERSION_free_vector_spaces(0,2,0) tensorDecomposition (Tensor (V3 x y z)) = [ (e @0, x), (e @1, y), (e @2, z) ] + tensorDecompose' (Tensor (V3 x y z)) b = case getEuclideanBasisIndex b of + { 0 -> x; 1 -> y; 2 -> z } showsPrecBasis = showsPrec #else tensorDecomposition (Tensor (V3 x y z)) = [ (ex, x), (ey, y), (ez, z) ] + tensorDecompose' (Tensor w) (Mat.E q) = w^.q showsPrecBasis _ (Mat.E q) = (V3"ex""ey""ez"^.q ++) #endif instance TensorDecomposable (V4 ℝ) where #if MIN_VERSION_free_vector_spaces(0,2,0) tensorDecomposition (Tensor (V4 x y z w)) = [(e @0,x), (e @1,y), (e @2,z), (e @3,w)] + tensorDecompose' (Tensor (V4 x y z w)) b = case getEuclideanBasisIndex b of + { 0 -> x; 1 -> y; 2 -> z; 3 -> w } showsPrecBasis = showsPrec #else tensorDecomposition (Tensor (V4 x y z w)) = [ (ex, x), (ey, y), (ez, z), (ew, w) ] + tensorDecompose' (Tensor w) (Mat.E q) = w^.q showsPrecBasis _ (Mat.E q) = (V4"ex""ey""ez""ew"^.q ++) #endif +instance ∀ u v s + . ( TensorDecomposable u, TensorDecomposable v + , Fractional' s, Scalar u ~ s, Scalar v ~ s + , Scalar (DualVector u) ~ s, Scalar (DualVector v) ~ s ) + => TensorDecomposable (Tensor s u v) where + tensorDecomposition :: ∀ w . (TensorSpace w, Scalar w ~ s) + => (Tensor s u v)⊗w -> [((Basis u, Basis v), w)] + tensorDecomposition (Tensor t) = [ ((bu,bv),w) + | (bu,vw) <- tensorDecomposition @u (Tensor t) + , (bv,w) <- tensorDecomposition @v vw ] + tensorDecompose' :: ∀ w . (TensorSpace w, Scalar w ~ s) + => (Tensor s u v)⊗w -> (Basis u, Basis v) -> w + tensorDecompose' (Tensor t) (bu,bv) + = tensorDecompose' @v (tensorDecompose' @u (Tensor t) bu) bv + showsPrecBasis :: Int -> (Basis u, Basis v) -> ShowS + showsPrecBasis = undefined + tensorDecomposeShowsPrec :: ∀ u v s . ( TensorDecomposable u, FiniteDimensional v, Show v, Scalar u ~ s, Scalar v ~ s ) => Int -> Tensor s u v -> ShowS @@ -1390,6 +1451,12 @@ showsPrec = case (dualSpaceWitness::DualSpaceWitness x, dualSpaceWitness::DualSpaceWitness y) of (DualSpaceWitness, DualSpaceWitness) -> tensorDecomposeShowsPrec + +instance ( TensorDecomposable u + , Scalar u ~ s ) + => Show (Tensor s (Tensor s u v) w) where + showsPrec = case (dualSpaceWitness::DualSpaceWitness u) of + DualSpaceWitness -> undefined (^) :: Num a => a -> Int -> a diff --git a/linearmap-category.cabal b/linearmap-category.cabal --- a/linearmap-category.cabal +++ b/linearmap-category.cabal @@ -2,8 +2,8 @@ -- documentation, see http://haskell.org/cabal/users-guide/ name: linearmap-category -version: 0.5.0.1 -synopsis: Native, complete, matrix-free linear algebra. +version: 0.6.0.0 +synopsis: Native, complete-ish, matrix-free linear algebra. description: The term /numerical linear algebra/ is often used almost synonymous with /matrix modifications/. However, what's interesting for most applications are really just /points in some vector space/ @@ -13,19 +13,31 @@ . This library implements the crucial LA operations like solving linear equations and eigenvalue problems, without requiring - that the vectors are represented in some particular basis. Apart - from conceptual elegance (only operations that are actually + that the vectors are represented in some particular basis. + This appoach offers: + 1. conceptual elegance (only operations that are actually geometrically sensible will typecheck – this is far stronger than just confirming that the dimensions match, as some other libraries - do), this also opens up good optimisation possibilities: the - vectors can be unboxed, use dedicated sparse compression, possibly - carry out the computations on accelerated hardware (GPU etc.). - The spaces can even be infinite-dimensional (e.g. function spaces). + do) + 2. opportunity to type tensors more expressively. E.g. instead of + having a tensor with many dimensions that can easily be confused, + one can have e.g. a space of images and take the tensor product + with a linear batch space, etc.. + 3. it opens up optimisation possibilities: the vectors can be + unboxed, use dedicated sparse compression, possibly carry out the + computations on accelerated hardware (GPU etc.). The spaces can in + principle even be infinite-dimensional (e.g. function spaces). . The linear algebra algorithms in this package only require the vectors to support fundamental operations like addition, scalar - products, double-dual-space coercion and tensor products; none of - this requires a basis representation. + products, double-dual-space coercion and tensor products. These + are expressed by a hierarchy of type classes, none of which requires + a basis representation. + Basis representations are optional to allow storage in matrix-based + backends, but this too is done in a way that allows e.g. taking the + tensor product of a lazy function space with a static-dimensional + matrix space with a low-dimensional channels space, and then only + the inner dimensions will be stored in a packed format. homepage: https://github.com/leftaroundabout/linearmap-family license: GPL-3 license-file: LICENSE @@ -37,6 +49,11 @@ -- extra-source-files: cabal-version: >=1.10 +flag singletons3 + description: Whether to use a version of the singletons package after its split in a core- and base part + default: True + manual: False + library exposed-modules: Math.LinearMap.Category Math.LinearMap.Category.Instances.Deriving @@ -45,6 +62,8 @@ Math.VectorSpace.Dual Math.VectorSpace.MiscUtil.MultiConstraints Math.LinearMap.Category.Derivatives + Math.VectorSpace.DimensionAware.Theorems.MaybeNat + Math.VectorSpace.DimensionAware other-modules: Math.LinearMap.Category.Class Math.LinearMap.Asserted Math.LinearMap.Category.TensorQuot @@ -54,7 +73,8 @@ build-depends: base >=4.8 && <5, vector-space >=0.11 && <0.18, MemoTrie, constrained-categories >=0.3 && <0.5, - containers, vector, + containers, + vector >=0.12 && <0.14, tagged, free-vector-spaces >= 0.1.4 && < 0.3, linear, lens, transformers, @@ -65,7 +85,16 @@ call-stack, template-haskell >=2.12 && <2.20, th-abstraction >=0.4 && <0.5, + ghc-typelits-natnormalise >=0.7 && <0.8, + type-natural >=1.0 && <1.2, QuickCheck >=2.11 && <2.15 + if flag(singletons3) + build-depends: + singletons >=3.0 && <3.1, + singletons-base >=3.0 && <3.1 + else + build-depends: + singletons >=2.7 && <3.0 -- hs-source-dirs: default-language: Haskell2010 @@ -76,7 +105,7 @@ build-depends: base, linearmap-category, vector-space , QuickCheck , manifolds-core - , linear + , linear, vector , constrained-categories , tasty, tasty-quickcheck ghc-options: -threaded "-with-rtsopts -N8 -M2G" diff --git a/test/tasty/test.hs b/test/tasty/test.hs --- a/test/tasty/test.hs +++ b/test/tasty/test.hs @@ -13,7 +13,9 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE TypeApplications #-} +{-# LANGUAGE DataKinds #-} {-# LANGUAGE UnicodeSyntax #-} {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE TemplateHaskell #-} @@ -25,9 +27,12 @@ import Control.Arrow.Constrained import Data.AffineSpace +import Linear.V3 import Linear.V4 import Data.Basis +import Data.Coerce import Math.LinearMap.Category +import Math.VectorSpace.DimensionAware (toArray, fromArray, unsafeFromArray) import Math.Manifold.Core.Types import Math.Manifold.Core.PseudoAffine @@ -35,14 +40,16 @@ import Test.Tasty.QuickCheck import qualified Test.QuickCheck as QC +import qualified Data.Vector.Unboxed as UArr + newtype ℝ⁴ = ℝ⁴ { getℝ⁴ :: V4 ℝ } deriving (Eq, Show) copyNewtypeInstances [t| ℝ⁴ |] [ ''AdditiveGroup, ''AffineSpace, ''VectorSpace , ''Semimanifold, ''PseudoAffine - , ''TensorSpace, ''LinearSpace + , ''DimensionAware, ''Dimensional, ''TensorSpace, ''LinearSpace , ''FiniteDimensional, ''SemiInner, ''InnerSpace ] newtype H¹ℝ⁴ a = H¹ℝ⁴ { getH¹ℝ⁴ :: ((a,a),(a,a)) } @@ -52,7 +59,7 @@ . (RealFloat' a, FiniteDimensional a, SemiInner a) => H¹ℝ⁴ a |] [ ''AdditiveGroup, ''AffineSpace, ''VectorSpace , ''Semimanifold, ''PseudoAffine - , ''TensorSpace, ''LinearSpace + , ''DimensionAware, ''Dimensional, ''TensorSpace, ''LinearSpace , ''FiniteDimensional, ''SemiInner ] derivative₄ :: H¹ℝ⁴ ℝ -> ℝ⁴ @@ -61,7 +68,15 @@ instance InnerSpace (H¹ℝ⁴ ℝ) where H¹ℝ⁴ v <.> H¹ℝ⁴ w = v<.>w + derivative₄ (H¹ℝ⁴ v)<.>derivative₄ (H¹ℝ⁴ w) +instance Arbitrary ℝ⁴ where + arbitrary = ℝ⁴ <$> do + V4 <$> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary +instance Arbitrary w => Arbitrary (Tensor ℝ ℝ⁴ w) where + arbitrary = Tensor <$> do + V4 <$> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary + + newtype ℝ⁵ a = ℝ⁵ { getℝ⁵ :: [ℝ] } deriving (Eq, Show) @@ -152,6 +167,68 @@ . QC.expectFailure $ \v -> (riesz-+$>AbstractDualVector v) ≈≈≈ (H¹ℝ⁴ v :: H¹ℝ⁴ Double) ] + , testGroup "Reading from arrays" + [ testProperty "Scalars" + $ \x -> fromArray (uar [x :: ℝ]) === Just x + , testProperty "Pairs" + $ \x y -> fromArray (uar [x,y :: ℝ]) === Just (x,y) + , testProperty "Nested pairs" + $ \x y ξ υ -> fromArray (uar [x,y,ξ,υ :: ℝ]) === Just ((x,y),(ξ,υ)) + , testProperty "ℝ³" + $ \x y z -> fromArray (uar [x,y,z :: ℝ]) === Just (V3 x y z) + , testProperty "Tensors: (ℝ,ℝ)⊗ℝ³" + $ \x y z ξ υ ζ -> fromArray (uar [x,y,z + ,ξ,υ,ζ :: ℝ]) + === Just (coerce ( V3 x y z + , V3 ξ υ ζ ) :: (ℝ,ℝ)⊗V3 ℝ) + , testProperty "Tensors: ℝ³⊗(ℝ,ℝ)" + $ \x y z ξ υ ζ -> fromArray (uar [x,ξ + ,y,υ + ,z,ζ :: ℝ]) + === Just (coerce (V3 (x,ξ) + (y,υ) + (z,ζ)) :: V3 ℝ⊗(ℝ,ℝ)) + , testProperty "Tensors: (ℝ,ℝ)⊗(ℝ,ℝ)⊗(ℝ,ℝ)" + $ \a b c d e f g h -> fromArray (uar [a,b,c,d,e,f,g,h :: ℝ]) + == Just (coerce (((a,b),(c,d)),((e,f),(g,h))) + :: (ℝ,ℝ)⊗(ℝ,ℝ)⊗(ℝ,ℝ)) + , testProperty "Linear functions: (ℝ,ℝ)-+>ℝ³" + $ \xx xy yx yy zx zy x y + -> (unsafeFromArray (uar [xx,yx,zx + ,xy,yy,zy]) + -+$> (unsafeFromArray (uar [x,y]) :: (ℝ,ℝ))) + === (unsafeFromArray + (uar [ xx*x + xy*y + , yx*x + yy*y + , zx*x + zy*y ]) :: V3 ℝ) + , testProperty "Linear functions: ℝ³-+>(ℝ,ℝ)" + $ \xx xy xz yx yy yz x y z + -> (unsafeFromArray (uar [xx,yx + ,xy,yy + ,xz,yz]) + -+$> (unsafeFromArray (uar [x,y,z]) :: V3 ℝ)) + === (unsafeFromArray + (uar [ xx*x + xy*y + xz*z + , yx*x + yy*y + yz*z ]) :: (ℝ,ℝ)) + -- N.B. this test is sensitive to the computation + -- order, e.g. it fails with xy*y + xx*x + xz*z due to + -- floating-point non-associativity and the exact ===. + ] + , testGroup "Array conversion" + $ let arrayRoundTrip :: ∀ v n . (n`Dimensional`v, Scalar v ~ ℝ, Eq v, Show v) + => v -> QC.Property + arrayRoundTrip v = fromArray (toArray v :: UArr.Vector ℝ) === Just v + in [ testProperty "ℝ" $ arrayRoundTrip @ℝ + , testProperty "(ℝ,ℝ)" $ arrayRoundTrip @(ℝ,ℝ) + , testProperty "ℝ³" $ arrayRoundTrip @(V3 ℝ) + , testProperty "ℝ⁴ (newtype-derived)" $ arrayRoundTrip @ℝ⁴ + , testProperty "ℝ⁵ (basis-derived)" $ arrayRoundTrip @(ℝ⁵ Int) + , testProperty "ℝ³⊗(ℝ,ℝ)" $ arrayRoundTrip @(V3 ℝ⊗(ℝ,ℝ)) + , testProperty "(ℝ,ℝ)⊗ℝ³" $ arrayRoundTrip @((ℝ,ℝ)⊗V3 ℝ) + , testProperty "ℝ³⊗ℝ³⊗ℝ³" $ arrayRoundTrip @(V3 ℝ⊗V3 ℝ⊗V3 ℝ) + , testProperty "ℝ³+>ℝ³" $ arrayRoundTrip @(V3 ℝ+>V3 ℝ) + , testProperty "ℝ³⊗ℝ⁴⊗ℝ⁵" $ arrayRoundTrip @(V3 ℝ⊗ℝ⁴⊗ℝ⁵ Int) + ] ] @@ -160,3 +237,9 @@ v≈≈≈w | magnitudeSq (v^-^w) < (magnitudeSq v + magnitudeSq w)*1e-8 = QC.property True | otherwise = v===w + +uar :: UArr.Unbox a => [a] -> UArr.Vector a +uar = UArr.fromList + +instance QC.Arbitrary s => QC.Arbitrary (V3 s) where + arbitrary = V3 <$> QC.arbitrary <*> QC.arbitrary <*> QC.arbitrary