diff --git a/Data/VectorSpace/Free.hs b/Data/VectorSpace/Free.hs --- a/Data/VectorSpace/Free.hs +++ b/Data/VectorSpace/Free.hs @@ -8,10 +8,13 @@ -- Portability : portable -- {-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE CPP #-} {-# LANGUAGE ConstrainedClassMethods #-} +{-# LANGUAGE DefaultSignatures #-} +{-# LANGUAGE ScopedTypeVariables #-} module Data.VectorSpace.Free ( -- * Supported types @@ -32,16 +35,20 @@ -- ** Small , OneDimensional(..) -- ** Free + , FreeVectorSpace(..) , FiniteFreeSpace(..) ) where import Data.AffineSpace import Data.VectorSpace +import Data.Cross +import Data.VectorSpace.Free.Class import Data.VectorSpace.Free.FiniteSupportedSequence (FinSuppSeq) import Data.VectorSpace.Free.Sequence (Sequence) import Data.Basis import Data.MemoTrie +import Data.Void import qualified Linear as L import Linear.V0 @@ -62,6 +69,8 @@ import Foreign.C.Types (CFloat, CDouble) import GHC.Exts (IsList(..)) +import qualified GHC.Generics as Gnrx +import GHC.Generics (Generic, (:*:)(..)) vDecomp :: FoldableWithIndex (L.E v) v => v s -> [(L.E v, s)] vDecomp = ifoldr (\b s l -> (b,s):l) [] @@ -101,7 +110,15 @@ portFinDP(V0) portFinDP(V1) portFinDP(V2) +instance Num s => HasCross2 (V2 s) where + cross2 (V2 x y) = V2 (-y) x + portFinDP(V3) +instance Num s => HasCross3 (V3 s) where + V3 ax ay az `cross3` V3 bx by bz = V3 (ay * bz - az * by) + (az * bx - ax * bz) + (ax * by - ay * bx) + portFinDP(V4) @@ -137,6 +154,10 @@ class (VectorSpace v, Fractional (Scalar v)) => OneDimensional v where -- | Compare the (directed) length of two vectors. (^/^) :: v -> v -> Maybe (Scalar v) + default (^/^) :: ( Generic v, OneDimensional (VRep v) + , Scalar (VRep v) ~ Scalar v ) + => v -> v -> Maybe (Scalar v) + v ^/^ w = (Gnrx.from v :: VRep v) ^/^ Gnrx.from w -- | Unsafe version of '^/^'. (^/!) :: v -> v -> Scalar v v^/!w = case v^/^w of @@ -171,10 +192,24 @@ class (VectorSpace v, Num (Scalar v)) => FiniteFreeSpace v where - {-# MINIMAL freeDimension, toFullUnboxVect, unsafeFromFullUnboxVect #-} freeDimension :: Functor p => p v -> Int + default freeDimension :: (Generic v, FiniteFreeSpace (VRep v)) + => p v -> Int + freeDimension _ = freeDimension ([]::[VRep v]) toFullUnboxVect :: UArr.Unbox (Scalar v) => v -> UArr.Vector (Scalar v) + default toFullUnboxVect + :: ( Generic v, FiniteFreeSpace (VRep v) + , UArr.Unbox (Scalar v) + , Scalar (VRep v) ~ Scalar v ) + => v -> UArr.Vector (Scalar v) + toFullUnboxVect v = toFullUnboxVect (Gnrx.from v :: VRep v) unsafeFromFullUnboxVect :: UArr.Unbox (Scalar v) => UArr.Vector (Scalar v) -> v + default unsafeFromFullUnboxVect + :: ( Generic v, FiniteFreeSpace (VRep v) + , UArr.Unbox (Scalar v) + , Scalar (VRep v) ~ Scalar v ) + => UArr.Vector (Scalar v) -> v + unsafeFromFullUnboxVect v = Gnrx.to (unsafeFromFullUnboxVect v :: VRep v) fromUnboxVect :: UArr.Unbox (Scalar v) => UArr.Vector (Scalar v) -> v fromUnboxVect v = result where result = case UArr.length v of @@ -235,3 +270,36 @@ + +instance FiniteFreeSpace a => FiniteFreeSpace (Gnrx.Rec0 a s) where + freeDimension = freeDimension . fmap Gnrx.unK1 + toFullUnboxVect = toFullUnboxVect . Gnrx.unK1 + unsafeFromFullUnboxVect = Gnrx.K1 . unsafeFromFullUnboxVect + fromUnboxVect = Gnrx.K1 . fromUnboxVect +instance FiniteFreeSpace (f p) => FiniteFreeSpace (Gnrx.M1 i c f p) where + freeDimension = freeDimension . fmap Gnrx.unM1 + toFullUnboxVect = toFullUnboxVect . Gnrx.unM1 + unsafeFromFullUnboxVect = Gnrx.M1 . unsafeFromFullUnboxVect + fromUnboxVect = Gnrx.M1 . fromUnboxVect +instance (FiniteFreeSpace (f p), FiniteFreeSpace (g p), Scalar (f p) ~ Scalar (g p)) + => FiniteFreeSpace ((f :*: g) p) where + freeDimension p = freeDimension (fmap (\(x:*:_)->x) p) + + freeDimension (fmap (\(_:*:y)->y) p) + toFullUnboxVect (u:*:v) = toFullUnboxVect u UArr.++ toFullUnboxVect v + unsafeFromFullUnboxVect uv = u:*:v + where u = unsafeFromFullUnboxVect uv + v = unsafeFromFullUnboxVect $ UArr.drop du uv + du = freeDimension [u] + + + + +instance OneDimensional a => OneDimensional (Gnrx.Rec0 a s) where + Gnrx.K1 v ^/^ Gnrx.K1 w = v ^/^ w + Gnrx.K1 v ^/! Gnrx.K1 w = v ^/! w +instance OneDimensional (f p) => OneDimensional (Gnrx.M1 i c f p) where + Gnrx.M1 v ^/^ Gnrx.M1 w = v ^/^ w + Gnrx.M1 v ^/! Gnrx.M1 w = v ^/! w + + +type VRep v = Gnrx.Rep v Void diff --git a/Data/VectorSpace/Free/Class.hs b/Data/VectorSpace/Free/Class.hs new file mode 100644 --- /dev/null +++ b/Data/VectorSpace/Free/Class.hs @@ -0,0 +1,58 @@ +-- | +-- Module : Data.VectorSpace.Free.Class +-- Copyright : (c) Justus Sagemüller 2016 +-- License : GPL v3 +-- +-- Maintainer : (@) sagemueller $ geo.uni-koeln.de +-- Stability : experimental +-- Portability : portable +-- +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} + +module Data.VectorSpace.Free.Class ( FreeVectorSpace(..) ) where + +import Data.VectorSpace +import Data.Ratio + +import Foreign.C.Types (CFloat, CDouble, CSChar, CShort, CInt, CLong, CLLong, CIntMax) + + +-- | Vector spaces that are spanned by a specific, canonical set of basis vectors. +class (VectorSpace v, Num (Scalar v)) => FreeVectorSpace v where + -- | Element-wise multiplication, equivalent to Matlab's @.*@ operator or + -- @'L.liftI2' (*)@. + (^*^) :: v -> v -> v + -- | Like a monomorphic 'fmap'. Only guaranteed to act on the nonzero entries; + -- whether the function is also applied on zeroes is instance-specific. + vmap :: (Scalar v -> Scalar v) -> v -> v + +instance FreeVectorSpace Float where {(^*^) = (*); vmap = id} +instance FreeVectorSpace Double where {(^*^) = (*); vmap = id} +instance FreeVectorSpace Int where {(^*^) = (*); vmap = id} +instance FreeVectorSpace Integer where {(^*^) = (*); vmap = id} +instance FreeVectorSpace CSChar where {(^*^) = (*); vmap = id} +instance FreeVectorSpace CShort where {(^*^) = (*); vmap = id} +instance FreeVectorSpace CInt where {(^*^) = (*); vmap = id} +instance FreeVectorSpace CLong where {(^*^) = (*); vmap = id} +instance FreeVectorSpace CLLong where {(^*^) = (*); vmap = id} +instance FreeVectorSpace CFloat where {(^*^) = (*); vmap = id} +instance FreeVectorSpace CDouble where {(^*^) = (*); vmap = id} +instance FreeVectorSpace CIntMax where {(^*^) = (*); vmap = id} +instance Integral a => FreeVectorSpace (Ratio a) where {(^*^) = (*); vmap = id} + +instance (FreeVectorSpace v, FreeVectorSpace w, Scalar v ~ Scalar w) + => FreeVectorSpace (v,w) where + (v,w)^*^(v',w') = (v^*^v', w^*^w') + vmap f (v,w) = (vmap f v, vmap f w) +instance ( FreeVectorSpace u, FreeVectorSpace v, FreeVectorSpace w + , Scalar v ~ Scalar u, Scalar v ~ Scalar w ) + => FreeVectorSpace (u,v,w) where + (u,v,w)^*^(u',v',w') = (u^*^u', v^*^v', w^*^w') + vmap f (u,v,w) = (vmap f u, vmap f v, vmap f w) +instance ( FreeVectorSpace u, FreeVectorSpace v, FreeVectorSpace w, FreeVectorSpace x + , Scalar x ~ Scalar v, Scalar v ~ Scalar u, Scalar v ~ Scalar w ) + => FreeVectorSpace (u,v,w,x) where + (u,v,w,x)^*^(u',v',w',x') = (u^*^u', v^*^v', w^*^w', x^*^x') + vmap f (u,v,w,x) = (vmap f u, vmap f v, vmap f w, vmap f x) + diff --git a/Data/VectorSpace/Free/FiniteSupportedSequence.hs b/Data/VectorSpace/Free/FiniteSupportedSequence.hs --- a/Data/VectorSpace/Free/FiniteSupportedSequence.hs +++ b/Data/VectorSpace/Free/FiniteSupportedSequence.hs @@ -12,22 +12,31 @@ {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE CPP #-} {-# LANGUAGE ConstrainedClassMethods #-} +{-# LANGUAGE MultiWayIf #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE UnicodeSyntax #-} module Data.VectorSpace.Free.FiniteSupportedSequence ( - FinSuppSeq (..) + FinSuppSeq (..) + , SparseSuppSeq (..) + , SemisparseSuppSeq (..) ) where import Data.AffineSpace import Data.VectorSpace +import Data.VectorSpace.Free.Class import Data.Basis import qualified Data.Foldable as Foldable +import qualified Data.Vector.Generic as Arr import qualified Data.Vector.Unboxed as UArr import qualified Data.Vector.Generic.Mutable as MArr import GHC.Exts (IsList(..)) +import Control.Arrow (first, second) +import Control.Monad (forM_) @@ -89,6 +98,10 @@ type Scalar (FinSuppSeq n) = n μ*^FinSuppSeq v = FinSuppSeq $ UArr.map (μ*) v +instance (Num n, AdditiveGroup n, UArr.Unbox n) => FreeVectorSpace (FinSuppSeq n) where + FinSuppSeq v^*^FinSuppSeq w = FinSuppSeq (UArr.zipWith (*) v w) + vmap f (FinSuppSeq v) = FinSuppSeq $ UArr.map f v + instance (Num n, AdditiveGroup n, UArr.Unbox n) => InnerSpace (FinSuppSeq n) where FinSuppSeq v<.>FinSuppSeq w = UArr.sum (UArr.zipWith (*) v w) @@ -106,3 +119,155 @@ instance (UArr.Unbox n, Show n) => Show (FinSuppSeq n) where show = ("fromList "++) . show . toList + + + + + +-- | Sparsely supported sequences (what other languages would call /sparse vectors/) +-- are sequences consisting of lots of zeroes, with finitely many +-- nonzeroes scattered around. Only these nonzero elements are stored. +data SparseSuppSeq n = SparseSuppSeq { + sparseNonzeroes :: UArr.Vector (Int,n) + } + +instance (Num n, UArr.Unbox n) => AffineSpace (SparseSuppSeq n) where + type Diff (SparseSuppSeq n) = SparseSuppSeq n + (.-.) = (^-^) + (.+^) = (^+^) + +instance (Num n, UArr.Unbox n) => AdditiveGroup (SparseSuppSeq n) where + zeroV = SparseSuppSeq $ UArr.empty + SparseSuppSeq u ^+^ SparseSuppSeq v = SparseSuppSeq w + where w = Arr.unfoldrN (Arr.length u + Arr.length v) seekws (0,0) + seekws (pu,pv) = case (u Arr.!? pu, v Arr.!? pv) of + (Just (ju,uj), Just (jv,vj)) + -> if | ju>jv -> Just ((jv, vj), (pu, pv+1)) + | ju Just ((ju, uj), (pu+1, pv)) + | otherwise -> Just ((ju, uj+vj), (pu+1, pv+1)) + (Just (ju,uj), Nothing) + -> Just ((ju, uj), (pu+1, pv)) + (Nothing, Just (jv,vj)) + -> Just ((jv, vj), (pu, pv+1)) + (Nothing, Nothing) -> Nothing + negateV (SparseSuppSeq v) = SparseSuppSeq $ Arr.map (second negate) v + +instance (Num n, UArr.Unbox n) => VectorSpace (SparseSuppSeq n) where + type Scalar (SparseSuppSeq n) = n + μ *^ SparseSuppSeq v = SparseSuppSeq $ Arr.map (second (*μ)) v + +instance (Num n, UArr.Unbox n) => FreeVectorSpace (SparseSuppSeq n) where + SparseSuppSeq u ^*^ SparseSuppSeq v = SparseSuppSeq w + where w = Arr.unfoldrN (Arr.length u `min` Arr.length v) seekws (0,0) + seekws (pu,pv) = case (u Arr.!? pu, v Arr.!? pv) of + (Just (ju,uj), Just (jv,vj)) + -> if | ju>jv -> seekws (pu, pv+1) + | ju seekws (pu+1, pv) + | otherwise -> Just ((ju, uj*vj), (pu+1, pv+1)) + _ -> Nothing + vmap f (SparseSuppSeq v) = SparseSuppSeq $ Arr.map (second f) v + +instance (Num n, AdditiveGroup n, UArr.Unbox n) => InnerSpace (SparseSuppSeq n) where + v <.> w = case v ^*^ w of SparseSuppSeq vw -> Arr.foldl' (\acc (_,q) -> acc+q) 0 vw + +instance (Num n, UArr.Unbox n) => HasBasis (SparseSuppSeq n) where + type Basis (SparseSuppSeq n) = Int + basisValue i = SparseSuppSeq $ UArr.singleton (i,1) + decompose (SparseSuppSeq v) = UArr.toList v + decompose' (SparseSuppSeq v) i = goBisect 0 (Arr.length v) + where goBisect jb jt + | jb==jt = 0 + | otherwise = case first (`compare`i) $ v Arr.! jm of + (LT,_) -> goBisect (jm+1) jt + (EQ,q) -> q + (GT,_) -> goBisect jb jm + where jm = (jb+jt)`div`2 + +instance (UArr.Unbox n, Eq n, Num n) => IsList (SparseSuppSeq n) where + type Item (SparseSuppSeq n) = n + fromListN n xs = SparseSuppSeq $ Arr.unfoldrN n go (0,xs) + where go (_,[]) = Nothing + go (j,0:xs) = go (j+1,xs) + go (j,x:xs) = Just ((j,x), (j+1,xs)) + fromList l = fromListN (length l) l + toList (SparseSuppSeq xs) = go 0 0 + where go i j = case xs Arr.!? j of + Just (i',x) | i==i' -> x : go (i+1) (j+1) + Nothing -> [] + _ -> 0 : go (i+1) j + + + + +-- | Like 'SparseSuppSeq', this type of number-sequence ignores zeroes and only stores +-- nonzero elements with positional information, but it does this not for every single +-- entry separately: only the first position of each contiguous /chunk/ of nonzeroes +-- is tracked. It is thus more suited for vectors that are in some places dense +-- but still have lots of zeroes. +-- +-- The drawback is that random access (i.e. 'decompose'') has complexity 𝓞(𝑛) +-- – instead of 𝓞(1) for 'FinSuppSeq', or 𝓞(log 𝑛) for 'SparseSuppSeq' – +-- so this type should only be used for “abstract vector operations”. +data SemisparseSuppSeq n = SemisparseSuppSeq { + chunkSparseNonzeroes :: UArr.Vector n + , sparseNonzeroLocation :: UArr.Vector (Int, Int) + -- ^ Start index of block, + -- size of block of consecutive nonzeroes + } + + +asSemisparse :: UArr.Unbox n => SparseSuppSeq n -> SemisparseSuppSeq n +asSemisparse (SparseSuppSeq v) = SemisparseSuppSeq (Arr.map snd v) + $ Arr.unfoldrN (Arr.length v) mkIndex 0 + where mkIndex :: Int -> Maybe ((Int, Int), Int) + mkIndex i + | Just (j,_) <- v Arr.!? i = case mkIndex $ i+1 of + Just ((j',l),n) | j'==j+1 -> Just ((j,l+1), n) + _ -> Just ((j,1), i+1) + | otherwise = Nothing + +fromSemisparse :: ∀ n . UArr.Unbox n => SemisparseSuppSeq n -> SparseSuppSeq n +fromSemisparse (SemisparseSuppSeq v ssIx) = SparseSuppSeq . (`Arr.zip`v) $ Arr.create (do + ix <- MArr.new $ Arr.length v + Arr.foldM_ (\i (j,l) -> do + forM_ [0..l-1] -- TODO: faster loop and unsafeWrite + $ \k -> MArr.write ix (i+k) (j+k) + return $ i+l + ) 0 ssIx + return ix + ) + +instance (Num n, UArr.Unbox n) => AffineSpace (SemisparseSuppSeq n) where + type Diff (SemisparseSuppSeq n) = SemisparseSuppSeq n + (.-.) = (^-^) + (.+^) = (^+^) + +instance (Num n, UArr.Unbox n) => AdditiveGroup (SemisparseSuppSeq n) where + zeroV = SemisparseSuppSeq UArr.empty UArr.empty + u ^+^ v = -- TODO: faster, direct implementation + asSemisparse $ fromSemisparse u ^+^ fromSemisparse v + negateV (SemisparseSuppSeq v vis) = SemisparseSuppSeq (Arr.map negate v) vis + +instance (Num n, UArr.Unbox n) => VectorSpace (SemisparseSuppSeq n) where + type Scalar (SemisparseSuppSeq n) = n + μ *^ SemisparseSuppSeq v ix = SemisparseSuppSeq (Arr.map (μ*) v) ix + +instance (Num n, UArr.Unbox n) => FreeVectorSpace (SemisparseSuppSeq n) where + u ^*^ v = -- TODO: faster, direct implementation + asSemisparse $ fromSemisparse u ^*^ fromSemisparse v + vmap f (SemisparseSuppSeq v ix) = SemisparseSuppSeq (Arr.map f v) ix + +instance (Num n, AdditiveGroup n, UArr.Unbox n) => InnerSpace (SemisparseSuppSeq n) where + v <.> w = fromSemisparse v <.> fromSemisparse w + +instance (Num n, UArr.Unbox n) => HasBasis (SemisparseSuppSeq n) where + type Basis (SemisparseSuppSeq n) = Int + basisValue i = SemisparseSuppSeq (UArr.singleton 1) (UArr.singleton (i,1)) + decompose = decompose . fromSemisparse + decompose' v = decompose' $ fromSemisparse v + +instance (UArr.Unbox n, Eq n, Num n) => IsList (SemisparseSuppSeq n) where + type Item (SemisparseSuppSeq n) = n + fromListN n = asSemisparse . fromListN n + fromList = asSemisparse . fromList + toList = toList . fromSemisparse diff --git a/free-vector-spaces.cabal b/free-vector-spaces.cabal --- a/free-vector-spaces.cabal +++ b/free-vector-spaces.cabal @@ -2,7 +2,7 @@ -- documentation, see http://haskell.org/cabal/users-guide/ name: free-vector-spaces -version: 0.1.2.0 +version: 0.1.4.0 synopsis: Instantiate the classes from the vector-space package with types from linear description: The package offers efficient vector types — where vector means /element of a free vector space/, i.e. @@ -50,10 +50,10 @@ exposed-modules: Data.VectorSpace.Free , Data.VectorSpace.Free.FiniteSupportedSequence , Data.VectorSpace.Free.Sequence - -- other-modules: + other-modules: Data.VectorSpace.Free.Class -- other-extensions: build-depends: base >=4.6 && <5.1, - vector-space >=0.8 && <0.13, + vector-space >=0.11 && <0.13, MemoTrie, linear >=1.18 && <1.23, lens >= 4 && < 5,