packages feed

sparse-linear-algebra 0.1.0.2 → 0.1.0.3

raw patch · 2 files changed

+129/−158 lines, 2 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

- Math.Linear.Sparse: decode :: (Int, Int) -> Int -> (Rows, Cols)
- Math.Linear.Sparse: encode :: (Int, Int) -> (Rows, Cols) -> Int
- Math.Linear.Sparse: findWithDefault0IM :: Num a => Key -> IntMap a -> a
- Math.Linear.Sparse: matMatU :: Num a => SpMatrix a -> SpMatrix a -> SpMatrix a
- Math.Linear.Sparse: nzRowU :: SpMatrix a -> Key -> Int
- Math.Linear.Sparse: ones :: Num a => Int -> [a]
- Math.Linear.Sparse: (@@) :: Num a => SpMatrix a -> (Key, Key) -> a
+ Math.Linear.Sparse: (@@) :: Num a => SpMatrix a -> (IxRow, IxCol) -> a
- Math.Linear.Sparse: candidateRows :: IntMap (IntMap a) -> Key -> Key -> Maybe [Key]
+ Math.Linear.Sparse: candidateRows :: IntMap (IntMap a) -> IxRow -> IxCol -> Maybe [Key]
- Math.Linear.Sparse: extractCol :: SpMatrix a -> Int -> SpVector a
+ Math.Linear.Sparse: extractCol :: SpMatrix a -> IxCol -> SpVector a
- Math.Linear.Sparse: extractColSM :: SpMatrix a -> Int -> SpMatrix a
+ Math.Linear.Sparse: extractColSM :: SpMatrix a -> IxCol -> SpMatrix a
- Math.Linear.Sparse: extractRow :: SpMatrix a -> Int -> SpVector a
+ Math.Linear.Sparse: extractRow :: SpMatrix a -> IxRow -> SpVector a
- Math.Linear.Sparse: extractRowSM :: SpMatrix a -> Int -> SpMatrix a
+ Math.Linear.Sparse: extractRowSM :: SpMatrix a -> IxRow -> SpMatrix a
- Math.Linear.Sparse: extractSubmatrixSM :: SpMatrix a -> (Int, Int) -> (Int, Int) -> SpMatrix a
+ Math.Linear.Sparse: extractSubmatrixSM :: SpMatrix a -> (IxRow, IxCol) -> (IxRow, IxCol) -> SpMatrix a
- Math.Linear.Sparse: firstNonZeroColumn :: IntMap a -> Key -> Bool
+ Math.Linear.Sparse: firstNonZeroColumn :: IntMap a -> IxRow -> Bool
- Math.Linear.Sparse: givens :: SpMatrix Double -> Int -> Int -> SpMatrix Double
+ Math.Linear.Sparse: givens :: SpMatrix Double -> IxRow -> IxCol -> SpMatrix Double
- Math.Linear.Sparse: lookupSM :: SpMatrix a -> Key -> Key -> Maybe a
+ Math.Linear.Sparse: lookupSM :: SpMatrix a -> IxRow -> IxCol -> Maybe a
- Math.Linear.Sparse: lookupWD_IM :: Num a => IntMap (IntMap a) -> (Key, Key) -> a
+ Math.Linear.Sparse: lookupWD_IM :: Num a => IntMap (IntMap a) -> (IxRow, IxCol) -> a
- Math.Linear.Sparse: lookupWD_SM :: Num a => SpMatrix a -> (Key, Key) -> a
+ Math.Linear.Sparse: lookupWD_SM :: Num a => SpMatrix a -> (IxRow, IxCol) -> a
- Math.Linear.Sparse: ncols :: SpMatrix a -> Int
+ Math.Linear.Sparse: ncols :: SpMatrix a -> Cols
- Math.Linear.Sparse: nrows :: SpMatrix a -> Int
+ Math.Linear.Sparse: nrows :: SpMatrix a -> Rows

Files

sparse-linear-algebra.cabal view
@@ -1,5 +1,5 @@ name:                sparse-linear-algebra-version:             0.1.0.2+version:             0.1.0.3 synopsis:            Numerical computation in native Haskell  description:         Currently it provides iterative linear solvers, matrix decompositions, eigenvalue computations and related utilities. Please see README.md for details homepage:            https://github.com/ocramz/sparse-linear-algebra
src/Math/Linear/Sparse.hs view
@@ -213,7 +213,6 @@   --- * IntMap implementation  instance Set IM.IntMap where   liftU2 = IM.unionWith@@ -256,8 +255,6 @@   ---- ** instances for SpVector instance Functor SpVector where   fmap f (SV n x) = SV n (fmap f x) @@ -376,19 +373,17 @@   show (SV d x) = "SV (" ++ show d ++ ") "++ show (IM.toList x)  --- | lookup+-- | lookup an index in a SpVector (returns 0 if lookup fails)  lookupDenseSV :: Num a => IM.Key -> SpVector a -> a lookupDenseSV i (SV _ im) = IM.findWithDefault 0 i im  -findWithDefault0IM :: Num a => IM.Key -> IM.IntMap a -> a-findWithDefault0IM = IM.findWithDefault 0    --- | SV manipulation + -- | Tail elements tailSV :: SpVector a -> SpVector a tailSV (SV n sv) = SV (n-1) ta where@@ -400,9 +395,7 @@   --- | concatenate SpVector--+-- | concatenate two sparse vectors concatSV :: SpVector a -> SpVector a -> SpVector a concatSV (SV n1 s1) (SV n2 s2) = SV (n1+n2) (IM.union s1 s2') where   s2' = IM.mapKeys (+ n1) s2@@ -417,7 +410,6 @@   -- | promote a SV to SM- svToSM :: SpVector a -> SpMatrix a svToSM (SV n d) = SM (n, 1) $ IM.singleton 0 d @@ -449,8 +441,6 @@                       smData :: IM.IntMap (IM.IntMap a)} deriving Eq  ---- ** Instances for SpMatrix instance Show a => Show (SpMatrix a) where   show sm@(SM _ x) = "SM " ++ sizeStr sm ++ " "++ show (IM.toList x) @@ -504,7 +494,7 @@   --- ** MATRIX METADATA+-- ** Matrix metadata  -- type synonyms type Rows = Int@@ -513,7 +503,7 @@ type IxRow = Int type IxCol = Int --- *** predicates+-- *** Predicates -- |Are the supplied indices within matrix bounds? validIxSM :: SpMatrix a -> (Int, Int) -> Bool validIxSM mm = inBounds02 (dim mm)@@ -536,24 +526,29 @@   --- -- internal projection functions, do not export:+-- | Internal data (do not export) immSM :: SpMatrix t -> IM.IntMap (IM.IntMap t) immSM (SM _ imm) = imm +-- | (Number of rows, Number of columns) dimSM :: SpMatrix t -> (Rows, Cols) dimSM (SM d _) = d +-- | Number of rows times number of columns nelSM :: SpMatrix t -> Int nelSM (SM (nr,nc) _) = nr*nc --- | nrows, ncols : size accessors-nrows, ncols :: SpMatrix a -> Int+-- | Number of rows+nrows :: SpMatrix a -> Rows nrows = fst . dim-ncols = snd . dim +-- | Number of columns+ncols :: SpMatrix a -> Cols+ncols = snd . dim   +-- *** SpMatrix information data SMInfo = SMInfo { smNz :: Int,                        smSpy :: Double} deriving (Eq, Show) @@ -569,17 +564,16 @@  -- *** Non-zero elements in a row -nzRowU :: SpMatrix a -> IM.Key -> Int-nzRowU s i = maybe 0 IM.size (IM.lookup i $ immSM s)- nzRow :: SpMatrix a -> IM.Key -> Int nzRow s i | inBounds0 (nrows s) i = nzRowU s i-          | otherwise = error "nzRow : index out of bounds"+          | otherwise = error "nzRow : index out of bounds" where+              nzRowU :: SpMatrix a -> IM.Key -> Int+              nzRowU s i = maybe 0 IM.size (IM.lookup i $ immSM s)    --- *** bandwidth bounds (min, max)+-- *** Mandwidth bounds (min, max)  bwMinSM :: SpMatrix a -> Int bwMinSM = fst . bwBoundsSM@@ -605,12 +599,13 @@   --- ** SPARSE MATRIX BUILDERS+-- ** Sparse matrix builders +-- | Zero SpMatrix of size (m, n) zeroSM :: Int -> Int -> SpMatrix a zeroSM m n = SM (m,n) IM.empty  -+-- | Insert an element in a preexisting Spmatrix at the specified indices insertSpMatrix :: IxRow -> IxCol -> a -> SpMatrix a -> SpMatrix a insertSpMatrix i j x s   | inBounds02 d (i,j) = SM d $ insertIM2 i j x smd @@ -636,9 +631,8 @@     --- | to List --- |Populate missing entries with 0+-- |Convert SpMatrix to list and populate missing entries with 0 toDenseListSM :: Num t => SpMatrix t -> [(IxRow, IxCol, t)] toDenseListSM m =   [(i, j, m @@ (i, j)) | i <- [0 .. nrows m - 1], j <- [0 .. ncols m- 1]]@@ -653,10 +647,8 @@  -- *** Identity matrix eye :: Num a => Int -> SpMatrix a-eye n = mkDiagonal n (ones n)+eye n = mkDiagonal n (replicate n 1) -ones :: Num a => Int -> [a]-ones n = replicate n 1     @@ -683,11 +675,11 @@   -encode :: (Int, Int) -> (Rows, Cols) -> Int-encode (nr,_) (i,j) = i + (j * nr)+-- encode :: (Int, Int) -> (Rows, Cols) -> Int+-- encode (nr,_) (i,j) = i + (j * nr) -decode :: (Int, Int) -> Int -> (Rows, Cols)-decode (nr, _) ci = (r, c) where (c,r ) = quotRem ci nr+-- decode :: (Int, Int) -> Int -> (Rows, Cols)+-- decode (nr, _) ci = (r, c) where (c,r ) = quotRem ci nr   @@ -697,10 +689,10 @@   --- ** SUB-MATRICES-+-- ** Sub-matrices -extractSubmatrixSM :: SpMatrix a -> (Int, Int) -> (Int, Int) -> SpMatrix a+-- | Extract a submatrix given the specified index bounds+extractSubmatrixSM :: SpMatrix a -> (IxRow, IxCol) -> (IxRow, IxCol) -> SpMatrix a extractSubmatrixSM (SM (r, c) im) (i1, i2) (j1, j2)   | q = SM (m', n') imm'   | otherwise = error $ "extractSubmatrixSM : invalid indexing " ++ show (i1, i2) ++ ", " ++ show (j1, j2) where@@ -727,21 +719,21 @@     | otherwise = error $ "toSV : incompatible dimensions " ++ show (m,n)  --- *** Extract jth column-extractColSM :: SpMatrix a -> Int -> SpMatrix a+-- *** Extract j'th column+extractColSM :: SpMatrix a -> IxCol -> SpMatrix a extractColSM sm j = extractSubmatrixSM sm (0, nrows sm - 1) (j, j)  -- |", and place into SpVector-extractCol :: SpMatrix a -> Int -> SpVector a+extractCol :: SpMatrix a -> IxCol -> SpVector a extractCol m j = toSV $ extractColSM m j  --- *** Extract ith row-extractRowSM :: SpMatrix a -> Int -> SpMatrix a+-- *** Extract i'th row+extractRowSM :: SpMatrix a -> IxRow -> SpMatrix a extractRowSM sm i = extractSubmatrixSM sm (i, i) (0, ncols sm - 1)  -- |", and place into SpVector-extractRow :: SpMatrix a -> Int -> SpVector a+extractRow :: SpMatrix a -> IxRow -> SpVector a extractRow m i = toSV $ extractRowSM m i  @@ -750,7 +742,7 @@   --- ** MATRIX STACKING+-- ** Matrix stacking  -- | Vertical stacking vertStackSM, (-=-) :: SpMatrix a -> SpMatrix a -> SpMatrix a@@ -778,18 +770,18 @@   --- ** MATRIX ELEMENT LOOKUP+-- ** Matrix element lookup -lookupSM :: SpMatrix a -> IM.Key -> IM.Key -> Maybe a+lookupSM :: SpMatrix a -> IxRow -> IxCol -> Maybe a lookupSM (SM _ im) i j = IM.lookup i im >>= IM.lookup j  -- | Looks up an element in the matrix with a default (if the element is not found, zero is returned) -lookupWD_SM, (@@) :: Num a => SpMatrix a -> (IM.Key, IM.Key) -> a+lookupWD_SM, (@@) :: Num a => SpMatrix a -> (IxRow, IxCol) -> a lookupWD_SM sm (i,j) =   fromMaybe 0 (lookupSM sm i j) -lookupWD_IM :: Num a => IM.IntMap (IM.IntMap a) -> (IM.Key, IM.Key) -> a+lookupWD_IM :: Num a => IM.IntMap (IM.IntMap a) -> (IxRow, IxCol) -> a lookupWD_IM im (i,j) = fromMaybe 0 (IM.lookup i im >>= IM.lookup j)  -- | Zero-default lookup, infix form@@ -806,11 +798,13 @@   --- *** MISC SpMatrix OPERATIONS+-- *** Misc. SpMatrix operations +-- | Left fold over SpMatrix foldlSM :: (a -> b -> b) -> b -> SpMatrix a -> b foldlSM f n (SM _ m)= foldlIM2 f n m +-- | Indexed left fold over SpMatrix ifoldlSM :: (IM.Key -> IM.Key -> a -> b -> b) -> b -> SpMatrix a -> b ifoldlSM f n (SM _ m) = ifoldlIM2' f n m @@ -819,25 +813,16 @@   --- | mapping  ------- | folding---- count sub-diagonal nonzeros+-- |Count sub-diagonal nonzeros countSubdiagonalNZSM :: SpMatrix a -> Int countSubdiagonalNZSM (SM _ im) = countSubdiagonalNZ im  --- | filtering- -- extractDiagonalSM :: (Num a, Eq a) => SpMatrix a -> SpVector a -- extractDiagonalSM (SM (m,n) im) = mkSpVectorD m $ extractDiagonalIM2 im --- | extract with default 0+-- | Extract the diagonal as a SpVector (with default 0) extractDiagonalDSM :: Num a => SpMatrix a -> SpVector a extractDiagonalDSM mm = fromListDenseSV n $ foldr ins [] ll  where   ll = [0 .. n - 1]@@ -850,8 +835,7 @@     ---  filtering the index subset that lies below the diagonal-+-- |Filter the index subset that lies below the diagonal (used in the QR decomposition, for example) subdiagIndicesSM :: SpMatrix a -> [(IM.Key, IM.Key)] subdiagIndicesSM (SM _ im) = subdiagIndices im @@ -859,8 +843,7 @@   --- ** sparsify : remove 0s (!!!)-+-- ** sparsify : remove almost-0 elements (i.e. if |x| < eps) sparsifyIM2 :: IM.IntMap (IM.IntMap Double) -> IM.IntMap (IM.IntMap Double) sparsifyIM2 = ifilterIM2 (\_ _ x -> abs x >= eps) @@ -869,8 +852,8 @@   --- ** ROUNDING operations (!!!)-                              +-- ** Rounding operations (!!!)+-- | Round almost-0 and almost-1 to 0 and 1 respectively roundZeroOneSM :: SpMatrix Double -> SpMatrix Double roundZeroOneSM (SM d im) = sparsifySM $ SM d $ mapIM2 roundZeroOne im @@ -884,12 +867,9 @@   --- * ALGEBRAIC PRIMITIVE OPERATIONS----- ** Matrix transpose-+-- * Primitive algebra operations +-- | transposeSM, (#^) : Matrix transpose transposeSM, (#^) :: SpMatrix a -> SpMatrix a transposeSM (SM (m, n) im) = SM (n, m) (transposeIM2 im) @@ -897,21 +877,13 @@   --- | A^T B-(#^#) :: SpMatrix Double -> SpMatrix Double -> SpMatrix Double-a #^# b = transposeSM a #~# b  --- | A B^T-(##^) :: SpMatrix Double -> SpMatrix Double -> SpMatrix Double-a ##^ b = a #~# transposeSM b    ----- ** matrix action on a vector+-- ** Matrix action on a vector  {-  FIXME : matVec is more general than SpVector's :@@ -922,7 +894,7 @@   --- |matrix on vector+-- |Matrix-on-vector matVec, (#>) :: Num a => SpMatrix a -> SpVector a -> SpVector a matVec (SM (nr, nc) mdata) (SV n sv)   | nc == n = SV nr $ fmap (`dot` sv) mdata@@ -930,7 +902,7 @@  (#>) = matVec --- |vector on matrix (FIXME : transposes matrix: more costly than `matVec`, I think)+-- |Vector-on-matrix (FIXME : transposes matrix: more costly than `matVec`, I think) vecMat, (<#) :: Num a => SpVector a -> SpMatrix a -> SpVector a   vecMat (SV n sv) (SM (nr, nc) mdata)   | n == nr = SV nc $ fmap (`dot` sv) (transposeIM2 mdata)@@ -954,19 +926,16 @@  -- ** Matrix-matrix product --- unsafe matMat-matMatU :: Num a => SpMatrix a -> SpMatrix a -> SpMatrix a-matMatU m1 m2 =-  SM (nrows m1, ncols m2) im where-    im = fmap (\vm1 -> (`dot` vm1) <$> transposeIM2 (immSM m2)) (immSM m1)-- matMat, (##) :: Num a => SpMatrix a -> SpMatrix a -> SpMatrix a matMat m1 m2   | c1 == r2 = matMatU m1 m2   | otherwise = error $ "matMat : incompatible matrix sizes" ++ show (d1, d2) where       d1@(r1, c1) = dim m1       d2@(r2, c2) = dim m2+      matMatU :: Num a => SpMatrix a -> SpMatrix a -> SpMatrix a+      matMatU m1 m2 =+        SM (nrows m1, ncols m2) im where+          im = fmap (\vm1 -> (`dot` vm1) <$> transposeIM2 (immSM m2)) (immSM m1)       (##) = matMat@@ -982,7 +951,8 @@   --- ** Matrix-matrix product, sparsified (prunes all elements `x` for which `abs x <= eps`)+-- ** Matrix-matrix product, sparsified+-- | Removes all elements `x` for which `| x | <= eps`) matMatSparsified, (#~#)  :: SpMatrix Double -> SpMatrix Double -> SpMatrix Double matMatSparsified m1 m2 = sparsifySM $ matMat m1 m2 @@ -991,8 +961,22 @@   +-- ** Sparsified matrix products +-- | A^T B+(#^#) :: SpMatrix Double -> SpMatrix Double -> SpMatrix Double+a #^# b = transposeSM a #~# b ++-- | A B^T+(##^) :: SpMatrix Double -> SpMatrix Double -> SpMatrix Double+a ##^ b = a #~# transposeSM b++++++ -- * Predicates  -- |is the matrix orthogonal? i.e. Q^t ## Q == I@@ -1083,7 +1067,7 @@ non-zero but A has zeros in row k for all columns less than j. -} -givens :: SpMatrix Double -> Int -> Int -> SpMatrix Double+givens :: SpMatrix Double -> IxRow -> IxCol -> SpMatrix Double givens mm i j    | validIxSM mm (i,j) && isSquareSM mm =        sparsifySM $ fromListSM' [(i,i,c),(j,j,c),(j,i,-s),(i,j,s)] (eye (nrows mm))@@ -1095,12 +1079,12 @@     b = mm @@ (i, j)   -- element to zero out  -- |Is the `k`th the first nonzero column in the row?-firstNonZeroColumn :: IM.IntMap a -> IM.Key -> Bool+firstNonZeroColumn :: IM.IntMap a -> IxRow -> Bool firstNonZeroColumn mm k = isJust (IM.lookup k mm) &&                           isNothing (IM.lookupLT k mm)  -- |Returns a set of rows {k} that satisfy QR.C1-candidateRows :: IM.IntMap (IM.IntMap a) -> IM.Key -> IM.Key -> Maybe [IM.Key]+candidateRows :: IM.IntMap (IM.IntMap a) -> IxRow -> IxCol -> Maybe [IM.Key] candidateRows mm i j | IM.null u = Nothing                      | otherwise = Just (IM.keys u) where   u = IM.filterWithKey (\irow row -> irow /= i &&@@ -1356,7 +1340,7 @@  data LinSolveMethod = CGS_ | BICGSTAB_ deriving (Eq, Show)  --- | linear solve with _random_ starting vector+-- | Linear solve with _random_ starting vector linSolveM ::   PrimMonad m =>     LinSolveMethod -> SpMatrix Double -> SpVector Double -> m (SpVector Double)@@ -1368,7 +1352,7 @@     case method of CGS_ -> return $ _xBicgstab (bicgstab aa b x0 x0)                    BICGSTAB_ -> return $ _x (cgs aa b x0 x0) --- | linear solve with _deterministic_ starting vector (every component at 0.1) +-- | Linear solve with _deterministic_ starting vector (every component at 0.1)  linSolve ::   LinSolveMethod -> SpMatrix Double -> SpVector Double -> SpVector Double linSolve method aa b@@ -1383,8 +1367,7 @@       (m, n) = dim aa       nb     = dim b --- <\> : sets default solver method -+-- | <\> : linSolve using BiCGSTAB method and by default (<\>) :: SpMatrix Double -> SpVector Double -> SpVector Double       (<\>) = linSolve BICGSTAB_    @@ -1416,7 +1399,7 @@   --- * PRETTY PRINTING+-- ** Pretty printing of SpVector and SpMatrix   sizeStr :: SpMatrix a -> String@@ -1495,39 +1478,9 @@   ------- | rounding to 0 or 1 within some predefined numerical precision---almostZero, almostOne :: Double -> Bool-almostZero x = abs x <= eps-almostOne x = x >= (1-eps) && x < (1+eps)--withDefault :: (t -> Bool) -> t -> t -> t-withDefault q d x | q x = d-                  | otherwise = x--roundZero, roundOne :: Double -> Double-roundZero = withDefault almostZero 0-roundOne = withDefault almostOne 1--with2Defaults :: (t -> Bool) -> (t -> Bool) -> t -> t -> t -> t-with2Defaults q1 q2 d1 d2 x | q1 x = d1-                            | q2 x = d2-                            | otherwise = x--roundZeroOne :: Double -> Double-roundZeroOne = with2Defaults almostZero almostOne 0 1----+-- ** Control primitives for bounded iteration with convergence check  -- | transform state until a condition is met- modifyUntil :: MonadState s m => (s -> Bool) -> (s -> s) -> m s modifyUntil q f = do   x <- get@@ -1535,7 +1488,8 @@   put y   if q y then return y          else modifyUntil q f     -  ++-- | Keep a moving window buffer (length 2) of state `x` to assess convergence, stop when either a condition on that list is satisfied or when max # of iterations is reached   loopUntilAcc :: Int -> ([t] -> Bool) -> (t -> t)  -> t -> t loopUntilAcc nitermax q f x = go 0 [] x where   go i ll xx | length ll < 2 = go (i + 1) (y : ll) y @@ -1544,7 +1498,7 @@                            else go (i + 1) (take 2 $ y:ll) y                 where y = f xx --- | keep state `x` in a moving window of length 2 to assess convergence, stop when either a condition on that list is satisfied or when max # of iterations is reached+-- | Keep a moving window buffer (length 2) of state `x` to assess convergence, stop when either a condition on that list is satisfied or when max # of iterations is reached (runs in State monad) modifyInspectN ::   MonadState s m =>     Int ->           -- iteration budget@@ -1583,13 +1537,11 @@   ---+-- | iterate until convergence is verified or we run out of a fixed iteration budget untilConverged :: MonadState a m => (a -> SpVector Double) -> (a -> a) -> m a untilConverged fproj = modifyInspectN 100 (normDiffConverged fproj) --- convergence check (FIXME)+-- | convergence check (FIXME) normDiffConverged :: (Foldable t, Functor t) =>      (a -> SpVector Double) -> t a -> Bool normDiffConverged fp xx = normSq (foldrMap fp (^-^) (zeroSV 0) xx) <= eps@@ -1599,8 +1551,7 @@   --- run `niter` iterations and append the state `x` to a list `xs`, stop when either the `xs` satisfies a predicate `q` or when the counter reaches 0-+-- | run `niter` iterations and append the state `x` to a list `xs`, stop when either the `xs` satisfies a predicate `q` or when the counter reaches 0 runAppendN :: ([t] -> Bool) -> (t -> t) -> Int -> t -> [t] runAppendN qq ff niter x0 | niter<0 = error "runAppendN : niter must be > 0"                           | otherwise = go qq ff niter x0 [] where@@ -1609,7 +1560,7 @@     if n <= 0 || q xs then xs                       else go q f (n-1) x (x : xs) --- ", NO convergence check +-- | ", NO convergence check  runAppendN' :: (t -> t) -> Int -> t -> [t] runAppendN' ff niter x0 | niter<0 = error "runAppendN : niter must be > 0"                         | otherwise = go ff niter x0 [] where@@ -1625,13 +1576,39 @@   +-- ** Rounding operations +-- | Rounding rule+almostZero, almostOne :: Double -> Bool+almostZero x = abs x <= eps+almostOne x = x >= (1-eps) && x < (1+eps) +withDefault :: (t -> Bool) -> t -> t -> t+withDefault q d x | q x = d+                  | otherwise = x --- | random matrices and vectors+roundZero, roundOne :: Double -> Double+roundZero = withDefault almostZero 0+roundOne = withDefault almostOne 1 --- dense+with2Defaults :: (t -> Bool) -> (t -> Bool) -> t -> t -> t -> t+with2Defaults q1 q2 d1 d2 x | q1 x = d1+                            | q2 x = d2+                            | otherwise = x +-- | Round to respectively 0 or 1 within some predefined numerical precision eps+roundZeroOne :: Double -> Double+roundZeroOne = with2Defaults almostZero almostOne 0 1++++++++-- *** Random matrices and vectors++-- |Dense SpMatrix randMat :: PrimMonad m => Int -> m (SpMatrix Double) randMat n = do   g <- MWC.create@@ -1639,7 +1616,8 @@   let ii_ = [0 .. n-1]       (ix_,iy_) = unzip $ concatMap (zip ii_ . replicate n) ii_   return $ fromListSM (n,n) $ zip3 ix_ iy_ aav-  ++-- | Dense SpVector   randVec :: PrimMonad m => Int -> m (SpVector Double) randVec n = do   g <- MWC.create@@ -1649,8 +1627,7 @@   --- sparse-+-- | Sparse SpMatrix randSpMat :: Int -> Int -> IO (SpMatrix Double) randSpMat n nsp | nsp > n^2 = error "randSpMat : nsp must be < n^2 "                 | otherwise = do@@ -1660,7 +1637,7 @@   jj <- replicateM nsp (MWC.uniformR (0, n-1) g :: IO Int)   return $ fromListSM (n,n) $ zip3 ii jj aav -+-- | Sparse SpVector randSpVec :: Int -> Int -> IO (SpVector Double) randSpVec n nsp | nsp > n = error "randSpVec : nsp must be < n"                 | otherwise = do@@ -1674,15 +1651,13 @@   --- | misc utils--+-- *** Misc. utilities  -- | integer-indexed ziplist denseIxArray :: [b] -> [(Int, b)] denseIxArray xs = zip [0..length xs-1] xs  --- ", 2d arrays+-- | ", 2d arrays denseIxArray2 :: Int -> [c] -> [(Int, Int, c)] denseIxArray2 m xs = zip3 (concat $ replicate n ii_) jj_ xs where   ii_ = [0 .. m-1]@@ -1693,16 +1668,18 @@  -- folds +-- | foldr over the results of a fmap foldrMap :: (Foldable t, Functor t) => (a -> b) -> (b -> c -> c) -> c -> t a -> c foldrMap ff gg x0 = foldr gg x0 . fmap ff +-- | strict left fold foldlStrict :: (a -> b -> a) -> a -> [b] -> a foldlStrict f = go   where     go z []     = z     go z (x:xs) = let z' = f z x in z' `seq` go z' xs --- indexed fold?+-- | indexed right fold ifoldr :: Num i =>      (a -> b -> b) -> b -> (i -> c -> d -> a) -> c -> [d] -> b   ifoldr mjoin mneutral f  = go 0 where@@ -1710,8 +1687,7 @@   go _ _ [] = mneutral  --- bounds checking-+-- *** Bounds checking type LB = Int type UB = Int @@ -1733,11 +1709,6 @@   -------