linearmap-category 0.1.0.0 → 0.1.0.1

raw patch · 3 files changed

+124/−68 lines, 3 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• Math/LinearMap/Category.hs +0/−21
• Math/VectorSpace/Docile.hs +123/−46
• linearmap-category.cabal +1/−1

@@ -196,27 +196,6 @@   --- | For real matrices, this boils down to 'transpose'.---   For free complex spaces it also incurs complex conjugation.---   --- The signature can also be understood as------ @--- adjoint :: (v +> w) -> (DualVector w +> DualVector v)--- @--- --- Or------ @--- adjoint :: (DualVector v +> DualVector w) -> (w +> v)--- @--- --- But /not/ @(v+>w) -> (w+>v)@, in general (though in a Hilbert space, this too is--- equivalent, via 'riesz' isomorphism).-adjoint :: (LSpace v, LSpace w, Scalar v ~ Scalar w)-               => (v +> DualVector w) -+> (w +> DualVector v)-adjoint = arr fromTensor . transposeTensor . arr asTensor-   

@@ -34,6 +34,7 @@ import Data.Set (Set) import Data.Ord (comparing) import Data.List (maximumBy, unfoldr)+import qualified Data.Vector as Arr import Data.Foldable (toList) import Data.Semigroup @@ -91,15 +92,15 @@                         --   the functional list.      -> ([(Int,v)] -> Forest (Int, DualVector v))                         -- ^ Suitable definition of 'dualBasisCandidates'.-cartesianDualBasisCandidates dvs abss vcas = go 0 sorted+cartesianDualBasisCandidates dvs abss vcas = go 0 0 sorted  where sorted = sortBy (comparing $ negate . snd . snd)                        [ (i, (av, maximum av)) | (i,v)<-vcas, let av = abss v ]-       go k ((i,(av,_)):scs)-          | k<n   = Node (i, dv) (go (k+1) [(i',(zeroAt j av',m)) | (i',(av',m))<-scs])-                                : go k scs+       go k nDelay scs@((i,(av,_)):scs')+          | k<n   = Node (i, dv) (go (k+1) 0 [(i',(zeroAt j av',m)) | (i',(av',m))<-scs'])+                                : go k (nDelay+1) (bringToFront (nDelay+1) scs)         where (j,_) = maximumBy (comparing snd) $ zip jfus av               dv = dvs !! j-       go _ _ = []+       go _ _ _ = []                jfus = [0 .. n-1]        n = length dvs@@ -108,6 +109,11 @@        zeroAt _ [] = []        zeroAt 0 (_:l) = (-1/0):l        zeroAt j (e:l) = e : zeroAt (j-1) l+       +       bringToFront :: Int -> [a] -> [a]+       bringToFront i l = case splitAt i l of+           (_,[]) -> []+           (f,s:l') -> s : f++l'  instance (Fractional'' s, SemiInner s) => SemiInner (ZeroDim s) where   dualBasisCandidates _ = []@@ -117,26 +123,36 @@ (<.>^) :: LSpace v => DualVector v -> v -> Scalar v f<.>^v = (applyDualVector$f)$v -orthonormaliseDuals :: (SemiInner v, LSpace v, Fractional'' (Scalar v))-                          => [(v, DualVector v)] -> [(v,DualVector v)]-orthonormaliseDuals [] = []-orthonormaliseDuals ((v,v'₀):ws)-          = (v,v') : [(w, w' ^-^ (w'<.>^v)*^v') | (w,w')<-wssys]- where wssys = orthonormaliseDuals ws-       v'₁ = foldl' (\v'i (w,w') -> v'i ^-^ (v'i<.>^w)*^w') v'₀ wssys-       v' = v'₁ ^/ (v'₁<.>^v)+orthonormaliseDuals :: (SemiInner v, LSpace v, RealFrac' (Scalar v))+                          => Scalar v -> [(v, DualVector v)] -> [(v,DualVector v)]+orthonormaliseDuals _ [] = []+orthonormaliseDuals ε ((v,v'₀):ws)+        | abs ovl > ε  = (v,v') : [(w, w' ^-^ (w'<.>^v)*^v') | (w,w')<-wssys]+        | otherwise    = (v,zeroV) : wssys+ where wssys = orthonormaliseDuals ε ws+       v'₁ = foldl' (\v'i (w,w') -> v'i ^-^ (v'i<.>^w)*^w') (v'₀ ^/ (v'₀<.>^v)) wssys+       v' = v'₁ ^/ ovl+       ovl = v'₁<.>^v -dualBasis :: (SemiInner v, LSpace v, Fractional'' (Scalar v)) => [v] -> [DualVector v]-dualBasis vs = snd <$> orthonormaliseDuals (zip' vsIxed candidates)+dualBasis :: (SemiInner v, LSpace v, RealFrac' (Scalar v)) => [v] -> [DualVector v]+dualBasis vs = snd <$> orthonormaliseDuals epsilon (zip' vsIxed candidates)  where zip' ((i,v):vs) ((j,v'):ds)         | i<j   = zip' vs ((j,v'):ds)         | i==j  = (v,v') : zip' vs ds        zip' _ _ = []-       candidates = sortBy (comparing fst) . findBest-                             $ dualBasisCandidates vsIxed-        where findBest [] = []-              findBest (Node iv' bv' : _) = iv' : findBest bv'+       candidates+         | Just bestCandidates <- findBest n $ dualBasisCandidates vsIxed+             = sortBy (comparing fst) bestCandidates+        where findBest 0 _ = Just []+              findBest _ [] = Nothing+              findBest n (Node (i,v') bv' : alts)+               | v'<.>^(lookupArr Arr.! i) /= 0+               , Just best' <- findBest (n-1) bv'+                            = Just $ (i,v') : best'+               | otherwise  = findBest n alts        vsIxed = zip [0..] vs+       lookupArr = Arr.fromList vs+       n = Arr.length lookupArr  instance SemiInner ℝ where   dualBasisCandidates = fmap ((`Node`[]) . second recip)@@ -181,16 +197,13 @@          combineBaseis _ forbidden ([], bv) = combineBaseis True forbidden ([],bv)  -instance ∀ s u v . ( LSpace u, FiniteDimensional (DualVector u), SemiInner (DualVector u)-                   , SemiInner v, FiniteDimensional v-                   , Scalar u ~ s, Scalar v ~ s, RealFrac' s )+instance ∀ s u v . ( SimpleSpace u, SimpleSpace v, Scalar u ~ s, Scalar v ~ s )            => SemiInner (Tensor s u v) where   dualBasisCandidates = map (fmap (second $ arr transposeTensor . arr asTensor))                       . dualBasisCandidates                       . map (second $ arr asLinearMap) -instance ∀ s u v . ( SemiInner u, FiniteDimensional u, Scalar u ~ s-                   , SemiInner v, FiniteDimensional v, Scalar v ~ s, RealFrac' s )+instance ∀ s u v . ( SimpleSpace u, SimpleSpace v, Scalar u ~ s, Scalar v ~ s )            => SemiInner (LinearMap s u v) where   dualBasisCandidates = sequenceForest                       . map (second pseudoInverse) -- this is not efficient@@ -254,7 +267,6 @@   --   library).   uncanonicallyFromDual :: DualVector v -+> v   uncanonicallyToDual :: v -+> DualVector v-     instance (Num''' s) => FiniteDimensional (ZeroDim s) where@@ -307,7 +319,7 @@ #define FreeFiniteDimensional(V, VB, dimens, take, give)        \ instance (Num''' s, LSpace s)                            \             => FiniteDimensional (V s) where {            \-  data SubBasis (V s) = VB;                             \+  data SubBasis (V s) = VB deriving (Show);             \   entireBasis = VB;                                      \   enumerateSubBasis VB = toList $ Mat.identity;      \   subbasisDimension VB = dimens;                       \@@ -392,7 +404,7 @@        = [ u⊗v | u <- enumerateSubBasis bu, v <- enumerateSubBasis bv ]   subbasisDimension (TensorBasis bu bv) = subbasisDimension bu * subbasisDimension bv   decomposeLinMap muvw = case decomposeLinMap $ curryLinearMap $ muvw of-         (bu, mvwsg) -> first (TensorBasis bu) . go id $ mvwsg []+         (bu, mvwsg) -> first (TensorBasis bu) . go $ mvwsg []    where (go, _) = tensorLinmapDecompositionhelpers   decomposeLinMapWithin (TensorBasis bu bv) muvw                = case decomposeLinMapWithin bu $ curryLinearMap $ muvw of@@ -418,13 +430,13 @@  tensorLinmapDecompositionhelpers       :: ( FiniteDimensional v, LSpace w , Scalar v~s, Scalar w~s )-      => ( DList w -> [v+>w] -> (SubBasis v, DList w)+      => ( [v+>w] -> (SubBasis v, DList w)          , SubBasis v -> DList w -> [v+>w] -> DList (v+>w)                         -> (Bool, (SubBasis v, DList w)) ) tensorLinmapDecompositionhelpers = (go, goWith)-   where go _ [] = decomposeLinMap zeroV-         go prevdc (mvw:mvws) = case decomposeLinMap mvw of-              (bv, cfs) -> snd (goWith bv prevdc mvws (mvw:))+   where go [] = decomposeLinMap zeroV+         go (mvw:mvws) = case decomposeLinMap mvw of+              (bv, cfs) -> snd (goWith bv cfs mvws (mvw:))          goWith bv prevdc [] prevs = (False, (bv, prevdc))          goWith bv prevdc (mvw:mvws) prevs = case decomposeLinMapWithin bv mvw of               Right cfs -> goWith bv (prevdc . cfs) mvws (prevs . (mvw:))@@ -439,6 +451,9 @@              \\nif it cannot decompose in the given basis, do so in a proper\              \\nsuperbasis of the given one (so that any vector that could be\              \\ndecomposed in the old basis can also be decomposed in the new one)."+  +deriving instance (Show (SubBasis u), Show (SubBasis v))+             => Show (SubBasis (Tensor s u v))   instance ∀ s u v .@@ -475,7 +490,10 @@   uncanonicallyFromDual = fmap uncanonicallyFromDual >>> arr asTensor              >>> transposeTensor >>> arr fromTensor >>> fmap uncanonicallyFromDual   +deriving instance (Show (SubBasis (DualVector u)), Show (SubBasis v))+             => Show (SubBasis (LinearMap s u v)) + infixr 0 \$  -- | Inverse function application, aka solving of a linear system:@@ -486,8 +504,10 @@ -- f '$' f '\$' u  ≡  u -- @ -- --- If @f@ does not have full rank, the behaviour is undefined (but we expect--- it to be reasonably well-behaved or even give a least-squares solution).+-- If @f@ does not have full rank, the behaviour is undefined. However, it+-- does not need to be a proper isomorphism: the+-- first of the above equations is still fulfilled if only @f@ is /injective/+-- (overdetermined system) and the second if it is /surjective/. --  -- If you want to solve for multiple RHS vectors, be sure to partially -- apply this operator to the linear map, like@@ -502,18 +522,53 @@ -- @ -- [f '\$' v₁, f '\$' v₂, ...] -- @-(\$) :: ( FiniteDimensional u, FiniteDimensional v, SemiInner v-        , Scalar u ~ Scalar v, Fractional' (Scalar v) )+(\$) :: ∀ u v . ( SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v )           => (u+>v) -> v -> u-(\$) m = fst . \v -> recomposeSB mbas [v'<.>^v | v' <- v's]- where v's = dualBasis $ mdecomp []-       (mbas, mdecomp) = decomposeLinMap m+(\$) m+  | du > dv    = (unsafeRightInverse m $)+  | du < dv    = (unsafeLeftInverse m $)+  | otherwise  = let v's = dualBasis $ mdecomp []+                     (mbas, mdecomp) = decomposeLinMap m+                 in fst . \v -> recomposeSB mbas [v'<.>^v | v' <- v's]+ where du = subbasisDimension (entireBasis :: SubBasis u)+       dv = subbasisDimension (entireBasis :: SubBasis v)      -pseudoInverse :: ( FiniteDimensional u, FiniteDimensional v, SemiInner v-                 , Scalar u ~ Scalar v, Fractional' (Scalar v) )+pseudoInverse :: ∀ u v . ( SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v )           => (u+>v) -> v+>u-pseudoInverse m = recomposeContraLinMap (fst . recomposeSB mbas) v's+pseudoInverse m+  | du > dv    = unsafeRightInverse m+  | du < dv    = unsafeLeftInverse m+  | otherwise  = unsafeInverse m+ where du = subbasisDimension (entireBasis :: SubBasis u)+       dv = subbasisDimension (entireBasis :: SubBasis v)++-- | If @f@ is injective, then+-- +-- @+-- unsafeLeftInverse f . f  ≡  id+-- @+unsafeLeftInverse :: ∀ u v . ( SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v )+                     => (u+>v) -> v+>u+unsafeLeftInverse m = unsafeInverse (m' . (fmap uncanonicallyToDual $ m))+                         . m' . arr uncanonicallyToDual+ where m' = adjoint $ m :: DualVector v +> DualVector u++-- | If @f@ is surjective, then+-- +-- @+-- f . unsafeRightInverse f  ≡  id+-- @+unsafeRightInverse :: ∀ u v . ( SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v )+                     => (u+>v) -> v+>u+unsafeRightInverse m = (fmap uncanonicallyToDual $ m')+                          . unsafeInverse (m . (fmap uncanonicallyToDual $ m'))+ where m' = adjoint $ m :: DualVector v +> DualVector u++-- | Invert an isomorphism. For other linear maps, the result is undefined.+unsafeInverse :: ( SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v )+          => (u+>v) -> v+>u+unsafeInverse m = recomposeContraLinMap (fst . recomposeSB mbas) v's  where v's = dualBasis $ mdecomp []        (mbas, mdecomp) = decomposeLinMap m @@ -525,7 +580,7 @@        let (bas, compos) = decomposeLinMap $ sampleLinearFunction $ applyDualVector $ dv        in fst . recomposeSB bas $ compos [] -sRiesz :: (FiniteDimensional v, InnerSpace v) => DualSpace v -+> v+sRiesz :: FiniteDimensional v => DualSpace v -+> v sRiesz = LinearFunction $ \dv ->        let (bas, compos) = decomposeLinMap $ dv        in fst . recomposeSB bas $ compos []@@ -543,6 +598,7 @@             . showsPrec 7 (sRiesz$dv)  instance Show (LinearMap ℝ (V0 ℝ) ℝ) where showsPrec = showsPrecAsRiesz+instance Show (LinearMap ℝ ℝ ℝ) where showsPrec = showsPrecAsRiesz instance Show (LinearMap ℝ (V1 ℝ) ℝ) where showsPrec = showsPrecAsRiesz instance Show (LinearMap ℝ (V2 ℝ) ℝ) where showsPrec = showsPrecAsRiesz instance Show (LinearMap ℝ (V3 ℝ) ℝ) where showsPrec = showsPrecAsRiesz@@ -561,22 +617,22 @@  instance Show (LinearMap s v (V0 s)) where   show _ = "zeroV"-instance (FiniteDimensional v, InnerSpace v, Scalar v ~ ℝ, Show v)+instance (FiniteDimensional v, v ~ DualVector v, Scalar v ~ ℝ, Show v)               => Show (LinearMap ℝ v (V1 ℝ)) where   showsPrec p m = showParen (p>6) $ ("ex .< "++)                        . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._x)) $ m)-instance (FiniteDimensional v, InnerSpace v, Scalar v ~ ℝ, Show v)+instance (FiniteDimensional v, v ~ DualVector v, Scalar v ~ ℝ, Show v)               => Show (LinearMap ℝ v (V2 ℝ)) where   showsPrec p m = showParen (p>6)               $ ("ex.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._x)) $ m)          . (" ^+^ ey.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._y)) $ m)-instance (FiniteDimensional v, InnerSpace v, Scalar v ~ ℝ, Show v)+instance (FiniteDimensional v, v ~ DualVector v, Scalar v ~ ℝ, Show v)               => Show (LinearMap ℝ v (V3 ℝ)) where   showsPrec p m = showParen (p>6)               $ ("ex.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._x)) $ m)          . (" ^+^ ey.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._y)) $ m)          . (" ^+^ ez.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._z)) $ m)-instance (FiniteDimensional v, InnerSpace v, Scalar v ~ ℝ, Show v)+instance (FiniteDimensional v, v ~ DualVector v, Scalar v ~ ℝ, Show v)               => Show (LinearMap ℝ v (V4 ℝ)) where   showsPrec p m = showParen (p>6)               $ ("ex.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._x)) $ m)@@ -635,3 +691,24 @@   unsafeFromFullUnboxVect arrv = arr (unsafeFromFullUnboxVect arrv :: LinearMap s u v)                                         ++-- | For real matrices, this boils down to 'transpose'.+--   For free complex spaces it also incurs complex conjugation.+--   +-- The signature can also be understood as+--+-- @+-- adjoint :: (v +> w) -> (DualVector w +> DualVector v)+-- @+-- +-- Or+--+-- @+-- adjoint :: (DualVector v +> DualVector w) -> (w +> v)+-- @+-- +-- But /not/ @(v+>w) -> (w+>v)@, in general (though in a Hilbert space, this too is+-- equivalent, via 'riesz' isomorphism).+adjoint :: (LSpace v, LSpace w, Scalar v ~ Scalar w)+               => (v +> DualVector w) -+> (w +> DualVector v)+adjoint = arr fromTensor . transposeTensor . arr asTensor

@@ -2,7 +2,7 @@ -- documentation, see http://haskell.org/cabal/users-guide/  name:                linearmap-category-version:             0.1.0.0+version:             0.1.0.1 synopsis:            Native, complete, matrix-free linear algebra. description:         The term /numerical linear algebra/ is often used almost                      synonymous with /matrix modifications/. However, what's interesting