np-linear 0.1.1.1 → 0.2.0.0
raw patch · 4 files changed
+54/−47 lines, 4 filesdep ~basedep ~reflection
Dependency ranges changed: base, reflection
Files
- np-linear.cabal +3/−3
- src/Algebra/Linear.hs +15/−28
- src/Algebra/Linear/Integral.hs +2/−0
- src/Algebra/Module/Free.hs +34/−16
np-linear.cabal view
@@ -1,5 +1,5 @@ name: np-linear-version: 0.1.1.1+version: 0.2.0.0 synopsis: Linear algebra for the numeric-prelude framework -- description: license: BSD3@@ -19,12 +19,12 @@ other-modules: Auxiliary build-depends:- base >= 4.5 && < 4.7,+ base >= 4.5 && < 4.8, binary >= 0.6.3 && < 0.8, containers >= 0.5 && < 0.6, numeric-prelude >= 0.3 && < 0.5, -- bifunctors >= 4.1 && < 0.5,- reflection >= 1.3 && < 1.5,+ reflection >= 1.3 && < 1.6, tagged == 0.7.* hs-source-dirs: src default-language: Haskell2010
src/Algebra/Linear.hs view
@@ -83,23 +83,6 @@ inverseImage a = solveUpperTriangular u . matrixVector b where (u,b,_) = reduce a --- Gives the row-reduced matrix, together with the determinant of the applied row operations.--- rowReduce :: forall k. (Algebra.Field.C k,Eq k) => Matrix k -> (Matrix k,k)--- rowReduce [] = ([],1)--- rowReduce xs@([] : _) = (xs,1)--- rowReduce vs = case findAmong ((/= 0) . headE "rr1" . snd) ivs of--- Nothing -> recurse (map (tail . snd) ivs)--- Just ((i,v₀),rest) -> first (v₀' :) . second (* s) $ recurse vs' where--- λ = headE "rr2" v₀--- v₀' = map (/ λ) v₀--- vs' = map (tail . f . snd) $ rest--- f v = zipWith (\ x y -> x - c * y) v v₀' where--- c = headE "rr3" v--- s = (if i == 0 then id else negate) λ--- where--- ivs = zip [0 ..] vs :: [(Int,Vector k)]--- recurse = first (map (0 :)) . rowReduce- invert :: (Algebra.Field.C k,Eq k) => Matrix k -> Maybe (Matrix k) invert m = fmap strip . process . (\ (r,_,_) -> r) . reduce . adorn $ m where n = length m@@ -119,13 +102,7 @@ strip = map (drop n) determinant :: (Algebra.Field.C k,Eq k) => Matrix k -> k-determinant [] = one-determinant a- | m == n = product (diagonal rr) * σ- where- (rr,_,σ) = reduce a- m = length a- n = length (head a)+determinant x = case reduce x of (_,_,σ) -> σ adjoint :: (Algebra.Lattice.C k,Algebra.Field.C k,Eq k) => Matrix k -> Maybe (Matrix k) adjoint a = (abs (determinant a) *>) <$> invert a@@ -172,24 +149,26 @@ reduce [] = ([],[],1) reduce xs@([] : _) = (xs,identity (length xs),1) reduce vs = case nonZero of- [] -> (\ (x,u,σ) -> (map (0 :) x,u,σ)) $ reduce (map tail vs)+ [] -> (\ (x,u,_σ) -> (map (0 :) x,u,0)) $ reduce (map tail vs) (v@(v₀ : _),i) : [] -> let (h,u,σ) = reduce (map (tail . fst) startZero)+ sign = if i == 1 then 1 else -1 in ( (map (/ v₀) v :) . map (0 :) $ h , normalisation v₀ `matrixProduct` rowSwap n (1,i) `matrixProduct` shift u- ,v₀ * σ+ ,sign * v₀ * σ ) (v@(v₀ : _),i) : rest -> let (reduced,translates) = unzip . flip map rest $ \ (x@(x₀ : _),j) -> let c = x₀ / v₀ in ((zipWith (\ vᵢ xᵢ -> xᵢ - c * vᵢ) v x,j),(j,c)) (h,u,σ) = reduce $ v : map fst reduced ++ map fst startZero+ permutation = i : map snd reduced ++ map snd startZero in ( h , u- `matrixProduct` permute (i : map snd reduced ++ map snd startZero)+ `matrixProduct` permute permutation `matrixProduct` affine n i (map (second negate) translates)- , σ+ , permutationSign permutation * σ ) where (startZero,nonZero) = partition ((==) 0 . head . fst) . flip zip [1 ..] $ vs@@ -217,6 +196,14 @@ permute :: (Algebra.Ring.C k) => [Int] -> Matrix k permute is = matrixFromFunction n n (\ i j -> δ (is !! (i - 1)) j) where n = length is++permutationSign :: (Algebra.Ring.C k) => [Int] -> k+permutationSign [] = one+permutationSign (x : xs) = ε * permutationSign xs where+ inversions = length $ filter (x >) xs+ ε+ | even inversions = one+ | otherwise = negate one -- Example x :: Matrix Rational
src/Algebra/Linear/Integral.hs view
@@ -33,6 +33,8 @@ -- Compute the Hermite normal form of the matrix, -- together with the unimodular basis transformation matrix.+-- In particular, if (h,u) = hermite l, then u l = h, u is unimodular,+-- and h is upper triangular. hermite :: Matrix Integer -> (Matrix Integer,Matrix Integer) hermite [] = ([],[]) hermite xs@([] : _) = (xs,identity (length xs))
src/Algebra/Module/Free.hs view
@@ -6,7 +6,7 @@ {-# LANGUAGE DeriveGeneric #-} module Algebra.Module.Free (- Free+ T(Free) , first , second @@ -18,13 +18,17 @@ , matrix , vector , fromVector+ , linear+ , filter+ , filterZeroes+ , toList_ , Enumerable , enumerate , coefficient , basisVector , fromList- , FreeBasis(FreeBasis)+ , Basis(Basis) ) where @@ -34,21 +38,23 @@ import qualified Algebra.Module import qualified Algebra.ModuleBasis import qualified Algebra.Ring-import NumericPrelude+import NumericPrelude hiding (filter) -- import Data.Bifunctor (Bifunctor,first,second) import Data.Binary (Binary)-import Data.List (intercalate)+import Data.List (intercalate,foldl') import qualified Data.Map.Strict as Map import Data.Proxy (Proxy(Proxy)) import GHC.Generics (Generic) -newtype Free k t+newtype T k t = Free (Map.Map t k) deriving (Generic) --- The natural Bifunctor Free instance is not possible, because of the constraints.+type Free k t = T k t++-- The natural Bifunctor Free instance is not possible, because of the constraints on 'second'. first :: (k₁ -> k₂) -> Free k₁ t -> Free k₂ t first f (Free m) = Free (fmap f m) second :: (Algebra.Additive.C k,Ord t₂) => (t₁ -> t₂) -> Free k t₁ -> Free k t₂@@ -75,17 +81,20 @@ basisVector :: b -> BasisElement b m -> m fromList :: (ModuleBasis b m,Algebra.Module.C (Scalar b m) m) => b -> [(BasisElement b m,Scalar b m)] -> m-fromList b = sum . map (\ (x,c) -> c *> basisVector b x)+fromList b = foldl' (+) zero . map (\ (x,c) -> c *> basisVector b x) -data FreeBasis- = FreeBasis+toList_ :: Free k t -> [(t,k)]+toList_ (Free m) = Map.toList m -instance (Algebra.Ring.C k,Ord t) => ModuleBasis FreeBasis (Free k t) where- type Scalar FreeBasis (Free k t) = k- type BasisElement FreeBasis (Free k t) = t- coefficient FreeBasis x (Free m) = Map.findWithDefault zero x m- basisVector FreeBasis x = Free (Map.singleton x one)+data Basis+ = Basis +instance (Algebra.Ring.C k,Ord t) => ModuleBasis Basis (Free k t) where+ type Scalar Basis (Free k t) = k+ type BasisElement Basis (Free k t) = t+ coefficient Basis x (Free m) = Map.findWithDefault zero x m+ basisVector Basis x = Free (Map.singleton x one)+ class (Ord t) => Enumerable p t where enumerate :: p -> [t] @@ -112,13 +121,22 @@ m matrix :: forall k t₁ t₂ p₁ p₂. (Enumerable p₁ t₁,Enumerable p₂ t₂,Algebra.Ring.C k) => p₁ -> p₂ -> Free k (t₁,t₂) -> Matrix k-matrix p₁ p₂ m = [[coefficient FreeBasis (x₁,x₂) m | x₁ <- enumerate p₁] | x₂ <- enumerate p₂]+matrix p₁ p₂ m = [[coefficient Basis (x₁,x₂) m | x₁ <- enumerate p₁] | x₂ <- enumerate p₂] vector :: forall k t p. (Enumerable p t,Algebra.Ring.C k) => p -> Free k t -> Vector k-vector p v = [coefficient FreeBasis x v | x <- enumerate p]+vector p v = [coefficient Basis x v | x <- enumerate p] fromVector :: forall k t p. (Enumerable p t,Algebra.Ring.C k) => p -> Vector k -> Free k t fromVector p = Free . Map.fromList . zip (enumerate p)++linear :: (Algebra.Module.C k a) => (t -> a) -> Free k t -> a+linear f (Free m) = sum . map (\ (t,c) -> c *> f t) . Map.toList $ m++filterZeroes :: (Algebra.Additive.C k,Eq k,Ord t) => Free k t -> Free k t+filterZeroes = filter $ const $ (/=) zero++filter :: (Algebra.Additive.C k,Ord t) => (t -> k -> Bool) -> Free k t -> Free k t+filter f (Free m) = Free $ Map.filterWithKey f m -- data X2 = S1 | S2 deriving (Eq,Ord) -- instance Enumerable X2 where