diff --git a/src/Data/LinearMap.hs b/src/Data/LinearMap.hs --- a/src/Data/LinearMap.hs +++ b/src/Data/LinearMap.hs @@ -1,30 +1,30 @@ -{-# LANGUAGE TypeOperators, FlexibleContexts, TypeFamilies #-} -{-# OPTIONS_GHC -Wall -fno-warn-orphans #-} --- {-# OPTIONS_GHC -funbox-strict-fields #-} --- {-# OPTIONS_GHC -ddump-simpl-stats -ddump-simpl #-} +{-# LANGUAGE TypeOperators, FlexibleContexts, TypeFamilies, GeneralizedNewtypeDeriving, StandaloneDeriving #-} +-- {-# OPTIONS_GHC -Wall -fno-warn-orphans #-} ---------------------------------------------------------------------- -- | -- Module : Data.LinearMap -- Copyright : (c) Conal Elliott 2008 -- License : BSD3 --- +-- -- Maintainer : conal@conal.net -- Stability : experimental --- +-- -- Linear maps ---------------------------------------------------------------------- module Data.LinearMap - ( (:-*) , linear, lapply, atBasis, idL, (*.*) - , liftMS, liftMS2, liftMS3 - , liftL, liftL2, liftL3 - ) where + ( (:-*) , linear, lapply, atBasis, idL, (*.*) + , inLMap, inLMap2, inLMap3 + , liftMS, liftMS2, liftMS3 + , liftL, liftL2, liftL3 + ) + where -import Control.Applicative ((<$>),Applicative,liftA2,liftA3) +import Control.Applicative (Applicative,liftA2,liftA3) import Control.Arrow (first) import Data.MemoTrie ((:->:)(..)) -import Data.AdditiveGroup (Sum(..),inSum2, AdditiveGroup(..)) +import Data.AdditiveGroup (Sum(..), AdditiveGroup(..)) import Data.VectorSpace (VectorSpace(..)) import Data.Basis (HasBasis(..), linearCombo) @@ -36,35 +36,33 @@ -- | An optional additive value type MSum a = Maybe (Sum a) --- nsum :: MSum a --- nsum = Nothing - jsum :: a -> MSum a jsum = Just . Sum +type LMap' u v = MSum (Basis u :->: v) + -- | Linear map, represented as an optional memo-trie from basis to -- values, where 'Nothing' means the zero map (an optimization). -type u :-* v = MSum (Basis u :->: v) +newtype u :-* v = LMap { unLMap :: LMap' u v } +deriving instance (HasTrie (Basis u), AdditiveGroup v) => AdditiveGroup (u :-* v) + +-- Before version 0.7, u :-* v was a type synonym, resulting in a subtle +-- ambiguity: u:-*v == u':-*v' does not imply that u==u', since Basis +-- might map different types to the same basis (e.g., Float & Double). +-- See . +-- See also . + -- TODO: Try a partial trie instead, excluding (known) zero elements. -- Then 'lapply' could be much faster for sparse situations. Make sure to -- correctly sum them. It'd be more like Jason Foutz's formulation -- -- which uses in @IntMap@. - --- PROBLEM: u :-* v is a type synonym, and Basis is an associated type synonym, resulting in a subtle --- ambiguity: u:-*v == u':-*v' does not imply that u==u', since Basis --- might map different types to the same basis (e.g., Float & Double). --- See --- --- Work in progress. See NewLinearMap.hs - - -- | Function (assumed linear) as linear map. linear :: (HasBasis u, HasTrie (Basis u)) => (u -> v) -> (u :-* v) -linear f = jsum (trie (f . basisValue)) +linear f = LMap (jsum (trie (f . basisValue))) atZ :: AdditiveGroup b => (a -> b) -> (MSum a -> b) atZ f = maybe zeroV (f . getSum) @@ -72,29 +70,36 @@ -- atZ :: AdditiveGroup b => (a -> b) -> (a -> b) -- atZ = id --- | Evaluate a linear map on a basis element. I've loosened the type to --- work around a typing problem in 'derivAtBasis'. --- atBasis :: (AdditiveGroup v, HasTrie (Basis u)) => --- (u :-* v) -> Basis u -> v -atBasis :: (HasTrie a, AdditiveGroup b) => MSum (a :->: b) -> a -> b -m `atBasis` b = atZ (`untrie` b) m +inLMap :: (LMap' r s -> LMap' t u) -> ((r :-* s) -> (t :-* u)) +inLMap = unLMap ~> LMap +inLMap2 :: (LMap' r s -> LMap' t u -> LMap' v w) + -> ((r :-* s) -> (t :-* u) -> (v :-* w)) +inLMap2 = unLMap ~> inLMap + +inLMap3 :: (LMap' r s -> LMap' t u -> LMap' v w -> LMap' x y) + -> ((r :-* s) -> (t :-* u) -> (v :-* w) -> (x :-* y)) +inLMap3 = unLMap ~> inLMap2 + -- | Apply a linear map to a vector. lapply :: ( VectorSpace v, Scalar u ~ Scalar v , HasBasis u, HasTrie (Basis u) ) => (u :-* v) -> (u -> v) -lapply = atZ lapply' +lapply = atZ lapply' . unLMap --- Handy for 'lapply' and '(*.*)'. +-- | Evaluate a linear map on a basis element. +atBasis :: (AdditiveGroup v, HasTrie (Basis u)) => + (u :-* v) -> Basis u -> v +LMap m `atBasis` b = atZ (`untrie` b) m + +-- | Handy for 'lapply' and '(*.*)'. lapply' :: ( VectorSpace v, Scalar u ~ Scalar v , HasBasis u, HasTrie (Basis u) ) => (Basis u :->: v) -> (u -> v) lapply' tr = linearCombo . fmap (first (untrie tr)) . decompose - - --- Identity linear map -idL :: (HasBasis u, HasTrie (Basis u)) => +-- | Identity linear map +idL :: (HasBasis u, HasTrie (Basis u)) => u :-* u idL = linear id @@ -108,35 +113,40 @@ (v :-* w) -> (u :-* v) -> (u :-* w) -- Simple definition, but only optimizes out uv == zero --- --- (*.*) vw = (fmap.fmap) (lapply vw) +-- vw *.* uv = LMap ((fmap.fmap.fmap) (lapply vw) (unLMap uv)) + +(*.*) vw = (inLMap.fmap.fmap.fmap) (lapply vw) + +-- Eep: +-- (*.*) = inLMap.fmap.fmap.fmap.lapply + + -- Instead, use Nothing/zero if /either/ map is zeroV (exploiting linearity -- when uv == zeroV.) --- Nothing *.* _ = Nothing --- _ *.* Nothing = Nothing --- Just (Sum vw) *.* Just (Sum uv) = Just (Sum (lapply' vw <$> uv)) +-- LMap Nothing *.* _ = LMap Nothing +-- _ *.* LMap Nothing = LMap Nothing +-- LMap (Just (Sum vw)) *.* LMap (Just (Sum uv)) = LMap (Just (Sum (lapply' vw <$> uv))) --- (*.*) = liftA2 (\ (Sum vw) (Sum uv) -> Sum (lapply' vw <$> uv)) +-- (*.*) = liftA2 (\ (LMap (Sum vw)) (LMap (Sum uv)) -> LMap (Sum (lapply' vw <$> uv))) --- (*.*) = (liftA2.inSum2) (\ vw uv -> lapply' vw <$> uv) -(*.*) = (liftA2.inSum2) (\ vw uv -> lapply' vw <$> uv) +-- (*.*) = (liftA2.inSum2.inLMap2) (\ vw uv -> lapply' vw <$> uv) --- (*.*) = (liftA2.inSum2) (\ vw -> fmap (lapply' vw)) +-- (*.*) = (liftA2.inSum2.inLMap2) (\ vw -> fmap (lapply' vw)) --- (*.*) = (liftA2.inSum2) (fmap . lapply') +-- (*.*) = (liftA2.inSum2.inLMap2) (fmap . lapply') -- It may be helpful that @lapply vw@ is evaluated just once and not -- once per uv. 'untrie' can strip off all of its trie constructors. -- Less efficient definition: --- +-- -- vw `compL` uv = linear (lapply vw . lapply uv) --- +-- -- i.e., compL = inL2 (.) --- +-- -- The problem with these definitions is that basis elements get converted -- to values and then decomposed, followed by recombination of the -- results. @@ -186,3 +196,43 @@ (a -> b -> c -> d) -> (MSum (f a) -> MSum (f b) -> MSum (f c) -> MSum (f d)) liftL3 = liftMS3 . liftA3 + +{- + + +infixr 9 *.* +-- | Compose linear maps +(*.*) :: ( HasBasis u, HasTrie (Basis u) + , HasBasis v, HasTrie (Basis v) + , VectorSpace w + , Scalar v ~ Scalar w ) => + (v :-* w) -> (u :-* v) -> (u :-* w) + +-- Simple definition, but only optimizes out uv == zero +-- +-- (*.*) vw = (fmap.fmap) (lapply vw) + +-- Instead, use Nothing/zero if /either/ map is zeroV (exploiting linearity +-- when uv == zeroV.) + +-- Nothing *.* _ = Nothing +-- _ *.* Nothing = Nothing +-- Just (Sum vw) *.* Just (Sum uv) = Just (Sum (lapply' vw <$> uv)) + +-- (*.*) = liftA2 (\ (Sum vw) (Sum uv) -> Sum (lapply' vw <$> uv)) + +-- (*.*) = (liftA2.inSum2) (\ vw uv -> lapply' vw <$> uv) +(*.*) = (liftA2.inSum2) (\ vw uv -> lapply' vw <$> uv) + +-- (*.*) = (liftA2.inSum2) (\ vw -> fmap (lapply' vw)) + +-- (*.*) = (liftA2.inSum2) (fmap . lapply') + + +-} + + +----- + +(~>) :: (a' -> a) -> (b -> b') -> ((a -> b) -> (a' -> b')) +(f ~> h) g = h . g . f diff --git a/src/Data/Maclaurin.hs b/src/Data/Maclaurin.hs --- a/src/Data/Maclaurin.hs +++ b/src/Data/Maclaurin.hs @@ -75,7 +75,7 @@ (b -> c) -> (a :> b) -> (a :> c) fmapD f = lf where - lf (D b0 b') = D (f b0) (liftL lf b') + lf (D b0 b') = D (f b0) ((inLMap.liftL) lf b') (<$>>) = fmapD @@ -84,7 +84,7 @@ (b -> c -> d) -> (a :> b) -> (a :> c) -> (a :> d) liftD2 f = lf where - lf (D b0 b') (D c0 c') = D (f b0 c0) (liftL2 lf b' c') + lf (D b0 b') (D c0 c') = D (f b0 c0) ((inLMap2.liftL2) lf b' c') -- | Apply a /linear/ ternary function over derivative towers. @@ -95,7 +95,7 @@ liftD3 f = lf where lf (D b0 b') (D c0 c') (D d0 d') = - D (f b0 c0 d0) (liftL3 lf b' c' d') + D (f b0 c0 d0) ((inLMap3.liftL3) lf b' c' d') -- TODO: Can liftD2 and liftD3 be defined in terms of a (<*>>) similar to @@ -173,8 +173,8 @@ distrib op = (#) where u@(D u0 u') # v@(D v0 v') = - D (u0 `op` v0) ( liftMS (inTrie ((# v) .)) u' ^+^ - liftMS (inTrie ((u #) .)) v' ) + D (u0 `op` v0) ( (inLMap.liftMS) (inTrie ((# v) .)) u' ^+^ + (inLMap.liftMS) (inTrie ((u #) .)) v' ) -- TODO: I think this distrib is exponential in increasing degree. Switch @@ -240,7 +240,7 @@ , AdditiveGroup (Scalar u)) => (u -> u) -> ((a :> u) -> (a :> Scalar u)) -> (a :> u) -> (a :> u) -f >-< f' = \ u@(D u0 u') -> D (f u0) (liftMS (f' u *^) u') +f >-< f' = \ u@(D u0 u') -> D (f u0) ((inLMap.liftMS) (f' u *^) u') -- TODO: express '(>-<)' in terms of '(@.)'. If I can't, then understand why not. diff --git a/src/Data/OldLinearMap.hs b/src/Data/OldLinearMap.hs new file mode 100644 --- /dev/null +++ b/src/Data/OldLinearMap.hs @@ -0,0 +1,188 @@ +{-# LANGUAGE TypeOperators, FlexibleContexts, TypeFamilies #-} +{-# OPTIONS_GHC -Wall -fno-warn-orphans #-} +-- {-# OPTIONS_GHC -funbox-strict-fields #-} +-- {-# OPTIONS_GHC -ddump-simpl-stats -ddump-simpl #-} +---------------------------------------------------------------------- +-- | +-- Module : Data.LinearMap +-- Copyright : (c) Conal Elliott 2008 +-- License : BSD3 +-- +-- Maintainer : conal@conal.net +-- Stability : experimental +-- +-- Linear maps +---------------------------------------------------------------------- + +module Data.LinearMap + ( (:-*) , linear, lapply, atBasis, idL, (*.*) + , liftMS, liftMS2, liftMS3 + , liftL, liftL2, liftL3 + ) where + +import Control.Applicative ((<$>),Applicative,liftA2,liftA3) +import Control.Arrow (first) + +import Data.MemoTrie ((:->:)(..)) +import Data.AdditiveGroup (Sum(..),inSum2, AdditiveGroup(..)) +import Data.VectorSpace (VectorSpace(..)) +import Data.Basis (HasBasis(..), linearCombo) + + +-- Linear maps are almost but not quite a Control.Category. The type +-- class constraints interfere. They're almost an Arrow also, but for the +-- constraints and the generality of arr. + +-- | An optional additive value +type MSum a = Maybe (Sum a) + +-- nsum :: MSum a +-- nsum = Nothing + +jsum :: a -> MSum a +jsum = Just . Sum + +-- | Linear map, represented as an optional memo-trie from basis to +-- values, where 'Nothing' means the zero map (an optimization). +type u :-* v = MSum (Basis u :->: v) + +-- TODO: Try a partial trie instead, excluding (known) zero elements. +-- Then 'lapply' could be much faster for sparse situations. Make sure to +-- correctly sum them. It'd be more like Jason Foutz's formulation +-- +-- which uses in @IntMap@. + + +-- PROBLEM: u :-* v is a type synonym, and Basis is an associated type synonym, resulting in a subtle +-- ambiguity: u:-*v == u':-*v' does not imply that u==u', since Basis +-- might map different types to the same basis (e.g., Float & Double). +-- See +-- +-- Work in progress. See NewLinearMap.hs + + +-- | Function (assumed linear) as linear map. +linear :: (HasBasis u, HasTrie (Basis u)) => + (u -> v) -> (u :-* v) +linear f = jsum (trie (f . basisValue)) + +atZ :: AdditiveGroup b => (a -> b) -> (MSum a -> b) +atZ f = maybe zeroV (f . getSum) + +-- atZ :: AdditiveGroup b => (a -> b) -> (a -> b) +-- atZ = id + +-- | Evaluate a linear map on a basis element. I've loosened the type to +-- work around a typing problem in 'derivAtBasis'. +-- atBasis :: (AdditiveGroup v, HasTrie (Basis u)) => +-- (u :-* v) -> Basis u -> v +atBasis :: (HasTrie a, AdditiveGroup b) => MSum (a :->: b) -> a -> b +m `atBasis` b = atZ (`untrie` b) m + +-- | Apply a linear map to a vector. +lapply :: ( VectorSpace v, Scalar u ~ Scalar v + , HasBasis u, HasTrie (Basis u) ) => + (u :-* v) -> (u -> v) +lapply = atZ lapply' + +-- Handy for 'lapply' and '(*.*)'. +lapply' :: ( VectorSpace v, Scalar u ~ Scalar v + , HasBasis u, HasTrie (Basis u) ) => + (Basis u :->: v) -> (u -> v) +lapply' tr = linearCombo . fmap (first (untrie tr)) . decompose + + + +-- Identity linear map +idL :: (HasBasis u, HasTrie (Basis u)) => + u :-* u +idL = linear id + + +infixr 9 *.* +-- | Compose linear maps +(*.*) :: ( HasBasis u, HasTrie (Basis u) + , HasBasis v, HasTrie (Basis v) + , VectorSpace w + , Scalar v ~ Scalar w ) => + (v :-* w) -> (u :-* v) -> (u :-* w) + +-- Simple definition, but only optimizes out uv == zero +-- +-- (*.*) vw = (fmap.fmap) (lapply vw) + +-- Instead, use Nothing/zero if /either/ map is zeroV (exploiting linearity +-- when uv == zeroV.) + +-- Nothing *.* _ = Nothing +-- _ *.* Nothing = Nothing +-- Just (Sum vw) *.* Just (Sum uv) = Just (Sum (lapply' vw <$> uv)) + +-- (*.*) = liftA2 (\ (Sum vw) (Sum uv) -> Sum (lapply' vw <$> uv)) + +-- (*.*) = (liftA2.inSum2) (\ vw uv -> lapply' vw <$> uv) +(*.*) = (liftA2.inSum2) (\ vw uv -> lapply' vw <$> uv) + +-- (*.*) = (liftA2.inSum2) (\ vw -> fmap (lapply' vw)) + +-- (*.*) = (liftA2.inSum2) (fmap . lapply') + + +-- It may be helpful that @lapply vw@ is evaluated just once and not +-- once per uv. 'untrie' can strip off all of its trie constructors. + +-- Less efficient definition: +-- +-- vw `compL` uv = linear (lapply vw . lapply uv) +-- +-- i.e., compL = inL2 (.) +-- +-- The problem with these definitions is that basis elements get converted +-- to values and then decomposed, followed by recombination of the +-- results. + +liftMS :: (AdditiveGroup a) => + (a -> b) + -> (MSum a -> MSum b) +-- liftMS _ Nothing = Nothing +-- liftMS h ma = Just (Sum (h (z ma))) + +liftMS = fmap.fmap + +liftMS2 :: (AdditiveGroup a, AdditiveGroup b) => + (a -> b -> c) -> + (MSum a -> MSum b -> MSum c) +liftMS2 _ Nothing Nothing = Nothing +liftMS2 h ma mb = Just (Sum (h (fromMS ma) (fromMS mb))) + +liftMS3 :: (AdditiveGroup a, AdditiveGroup b, AdditiveGroup c) => + (a -> b -> c -> d) -> + (MSum a -> MSum b -> MSum c -> MSum d) +liftMS3 _ Nothing Nothing Nothing = Nothing +liftMS3 h ma mb mc = Just (Sum (h (fromMS ma) (fromMS mb) (fromMS mc))) + +fromMS :: AdditiveGroup u => MSum u -> u +fromMS Nothing = zeroV +fromMS (Just (Sum u)) = u + + +-- | Apply a linear function to each element of a linear map. +-- @liftL f l == linear f *.* l@, but works more efficiently. +liftL :: (Functor f, AdditiveGroup (f a)) => + (a -> b) -> MSum (f a) -> MSum (f b) +liftL = liftMS . fmap + +-- | Apply a linear binary function (not to be confused with a bilinear +-- function) to each element of a linear map. +liftL2 :: (Applicative f, AdditiveGroup (f a), AdditiveGroup (f b)) => + (a -> b -> c) + -> (MSum (f a) -> MSum (f b) -> MSum (f c)) +liftL2 = liftMS2 . liftA2 + +-- | Apply a linear ternary function (not to be confused with a trilinear +-- function) to each element of a linear map. +liftL3 :: ( Applicative f + , AdditiveGroup (f a), AdditiveGroup (f b), AdditiveGroup (f c)) => + (a -> b -> c -> d) + -> (MSum (f a) -> MSum (f b) -> MSum (f c) -> MSum (f d)) +liftL3 = liftMS3 . liftA3 diff --git a/vector-space.cabal b/vector-space.cabal --- a/vector-space.cabal +++ b/vector-space.cabal @@ -1,5 +1,5 @@ Name: vector-space -Version: 0.6.2 +Version: 0.7.1 Cabal-Version: >= 1.2 Synopsis: Vector & affine spaces, linear maps, and derivatives (requires ghc 6.9 or better) Category: math